[R] Plotting Data on a Map

Wed Jun 23 23:47:21 CEST 2010

On Jun 23, 2010, at 4:57 PM, Felipe Carrillo wrote:

> For some reason the shapefile can't get attached.
> The shapefile is too large to put it in dput..Is there
> another way to do this?

If you dput or dump it and then label the output as <sth>.txt it will  
generally pass the server's scrutiny.

-- 
David.
>
>
>
> ----- Original Message ----
>> From: Felipe Carrillo <mazatlanmexico at yahoo.com>
>> To: Tom Hopper <tomhopper at gmail.com>
>> Cc: r-help at stat.math.ethz.ch; ggplot2 at googlegroups.com
>> Sent: Wed, June 23, 2010 1:52:29 PM
>> Subject: Re: [R] Plotting Data on a Map
>>
>> Hi:
> I am practicing with the attached shapefile and was wondering
> if I can
>> get some help. Haven't used 'rgdal' and 'maptools' much
> but it appears to be
>> a great way bring map data into R.
> Please take a look at the comments and let
>> me know if I need to
> explain better what I am trying to
>> accomplish.
>
> library(rgdal)
> library(maptools)
> library(ggplot2)
> dsn="C:/Documents
>> and Settings/fcarrillo/Desktop/Software/R Scripts
> and Datasets/NCTC/Data
>> Analisys II/Data 4 Data Analysis II March 2009-
> March2010/Data"
> wolves.map
>> <- readOGR(dsn=dsn,
>> layer="PNW_wolf_habitat_grid")
> class(wolves.map)
> dim(wolves.map)
> head(wolves.map,1)
> names(wolves.map)
> gpclibPermit()
> wolves.ds
>> <- fortify(wolves.map)
> head(wolves.ds);tail(wolves.ds)
> # Shapefile of 5
>> states
> wolves.plot <- ggplot(wolves.ds,aes(long,lat,group=group))
>> +
> geom_polygon(colour='white',fill='lightblue')
> wolves.plot
> # How to
>> show the state borders of Washington, Oregon, Idaho, Montana
> and Wyoming on
>> map?
>
> # Subset from wolves to create a logistic regression model based
>> on
> WOLVES_99 and then apply to all the 5 states
> # Remove the WOLVES_99
>> 2(two value) and use "one" for presence and
> "zero" for absence to end up with
>> 153 records.
> #wolfsub<-subset(wolves,subset=!(WOLVES_99==2))
> #wolfsub
>> <- subset(wolves.map,WOLVES_99!=2)
> wolfsub <-
>> wolves.map[!wolves.map$WOLVES_99 %in% 2,];wolfsub
> dim(wolfsub)
> # 42 =
>> Forest, 51 = Shrub, > 81 =
>> Agriculture
> wolfsub$Forest<-ifelse(wolfsub$MAJOR_LC==42,1,0)
> wolfsub$Shrub<-ifelse(wolfsub$MAJOR_LC==51,1,0)
> wolfsub$Agriculture<-ifelse(wolfsub$MAJOR_LC>81,1,0)
> names(wolfsub);dim(wolfsub)
> #
>> create the
>> model
> mod1<-glm(WOLVES_99~RD_DENSITY+Forest+Shrub
> +Agriculture,family=binomial,data=wolfsub)
> summary(mod1)
> wolfsub$pred99<-fitted(mod1)
> names(wolfsub)
> #fitted(mod1)
> wolfsub$pred99
>
> #
>> Add the wolfsub data to the map to see the map
> wolfsub <-
>> fortify(wolfsub);names(wolfsub)
> area_mod <- wolves.plot +
>> geom_polygon(data=wolfsub,aes(fill=????))
> #Want to fill by WOLVES_99 but is
>> gone #after fortify
> area_mod
> #Not sure how to proceed from here to fit the
>> 'mod1' model to all
> #the 5 states.
>
>
>
>
>>
>> From: Tom
>> Hopper <> href="mailto:tomhopper at gmail.com">tomhopper at gmail.com>
>> To: Felipe
>> Carrillo <> href="mailto:mazatlanmexico at yahoo.com">mazatlanmexico at yahoo.com 
>> >
>> Sent:
>> Tue, June 22, 2010 9:40:19 PM
>> Subject: Re: Plotting Data on a
>> Map
>>
>> Felipe,
>>
>>
>> I am just learning these
>> tools, too, so it may be a good learning opportunity for both of  
>> us. Please send
>> me the files and I will try to put it together and plot
>> it.
>>
>>
>> Regards,
>>
>>
>> Tom
>>
>>
>> On
>> Mon, Jun 21, 2010 at 11:57 PM, Felipe Carrillo <> ymailto="mailto:mazatlanmexico at yahoo.com 
>> "
>> href="mailto:mazatlanmexico at yahoo.com">mazatlanmexico at yahoo.com>
>> wrote:
>>
>> Hi Tom:
>>> I am just starting to use rgdal and
>> maptools but I have a long way to go. I went to a training
>>> a couple
>> of weeks ago and the instructor showed us a csv file and a  
>> shapefile with wolf
>> data from
>>> a national park in the midwest. I am trying to put all of
>> the csv data and some predicted data
>>> on a map using ggplot2 but I am
>> stuck with it. I am just trying to do this example because I want
>> to
>>> see if I can apply this example to fish. Let me know if
>> interested.
>>>
>>>
>>> Felipe D.
>> Carrillo
>>> Supervisory Fishery Biologist
>>> Department of the
>> Interior
>>> US Fish & Wildlife Service
>>> California,
>> USA
>>>
>>>
>>>
>>>>
>>>> From: Tom
>> Hopper <> href="mailto:tomhopper at gmail.com">tomhopper at gmail.com>
>>>> To:
>>> href="mailto:ggplot2 at googlegroups.com">ggplot2 at googlegroups.com
>>>> Sent:
>> Mon, June 21, 2010 2:27:14 PM
>>>> Subject: Re: Plotting Data on a
>> Map
>>>>
>>>>
>>>> Hadley,
>>
>>>>
>>>>
>>>> Thank
>> you!
>>>>
>>>>
>>>> I doing this, I've
>> encountered an encountered and unexpected difference between the   
>> graphs on my
>> Mac and my Windows machines. It appears that there is a default  
>> setting
>> different between the two
>> programs.
>>>>
>>>>
>>>> Using the same commands
>> on both Mac and Windows, I found that the polygon borders are  
>> plotted on the
>> Mac, but not on Windows, so on the Mac I see the country borders,  
>> but on Windows
>> I do not. On the Mac, it looks like the default for geom_polygon is  
>> something
>> like size = 0.01, colour = "gray50". On Windows, it's more like  
>> colour =
>> alpha("gray50", 0), though in fact the visibility of polygon  
>> borders seems to be
>> independent of both the size and colour
>> settings.
>>>>
>>>>
>>>> Regards,
>>>>
>>>>
>>>> Tom
>>>>
>>>>
>>>> On
>> Mon, Jun 21, 2010 at 3:00 AM, Hadley Wickham <> ymailto="mailto:hadley at rice.edu 
>> "
>> href="mailto:hadley at rice.edu">hadley at rice.edu>
>> wrote:
>>>>
>>>> Hi
>> Tom,
>>>>>
>>>>> Thanks for contribution! Ploting map
>> data is never easy and its always
>>>>> nice to see a success
>> story.
>>>>>
>>>>> Hadley
>>>>>
>>>>>
>>>>> On
>> Saturday, June 19, 2010, Tom Hopper <> href="mailto:tomhopper at gmail.com 
>> ">tomhopper at gmail.com>
>> wrote:
>>>>>> Well, it turns out that, in my haste to thank
>> Fernando, I posted code with some typos. I have also had a chance  
>> to test it on
>> my Mac, and found that fortify() was throwing an error ("Error in  
>> nchar(ID) :
>> invalid multibyte string 1"). I have added a call, currently  
>> commented out, to
>> Sys.setlocale() to fix this. I have tested the code below under R  
>> 2.11.1 on both
>> Windows XP and Mac OS X 10.5.8, with the latest packages installed.  
>> In fact, the
>> plot looks better in the Mac Quartz
>> window.
>>>>>>
>>>>>>
>>>>>>
>> Sorry for the previous
>> errors.
>>>>>>
>>>>>>
>>>>>>
>> # Download electricity generation data from
>> http://tonto.eia.doe.gov/cfapps/ipdbproject/iedindex3.cfm?tid=2&pid=2&aid=12&cid=&syid=2004&eyid=2008&unit=BKWH
>>>>>>
>>>>>>
>>>>>>
>> # Download new map data from
>> http://thematicmapping.org/downloads/world_borders.php
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # Bring the thematicmapping data into R
>>>>>>
>> library("rgdal")
>>>>>> world.map <- readOGR(dsn="C:/",
>> layer="TM_WORLD_BORDERS-0.3")
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # Save the map data as an R object
>>>>>> save(world.map,
>> "worldmap.rdata")
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # As needed, reload the data
>>>>>>
>> library(maptools)
>>>>>> gpclibPermit()
>>>>>>
>> load("worldmap.rdata")
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # Prepare the world.map data for ggplot2
>>>>>>
>> library(ggplot2)
>>>>>> # On some setups, fortify throws "Error
>> in nchar(ID)"
>>>>>> # in that case, use
>> setlocale
>>>>>> # Sys.setlocale("LC_ALL", locale =
>> "C")
>>>>>> world.ggmap <- fortify(world.map, region =
>> "NAME")
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # Load the electricity generation data and clean it up to match with
>> world.ggmape
>>>>>> elect.gen.tot <-
>> read 
>> .csv 
>> ("Total_Electricity_Net_Generation_(Billion_Kilowatthours).csv",  
>> sep =
>> ",", dec =
>> ".")
>>>>>>
>>>>>>
>>>>>>
>> names(elect.gen.tot) <- c("id", "y2004", "y2005", "y2006", "y2007",
>> "y2008")
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # make sure the id column is in the same case for both
>> sets
>>>>>> elect.gen.tot$id <-
>> tolower(elect.gen.tot$id)
>>>>>>
>>>>>>
>>>>>>
>> world.ggmap$id <-
>> tolower(world.ggmap$id)
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # merge the two data sets
>>>>>> world.ggmape <-
>> merge(world.ggmap, elect.gen.tot, by = "id", all =
>> TRUE)
>>>>>>
>>>>>>
>>>>>>
>> world.ggmape <- world.ggmape[order(world.ggmape$order),
>> ]
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> # NOTE: at this point, the electricity data country names do not  
>> match
>> 100%
>>>>>> # with the thematicmapping country names (column
>> "id").
>>>>>> # print out the country names
>> using
>>>>>> # table(world.ggmape$id)
>>>>>> #
>> and do a search for values of 1. Then find the corresponding
>> country
>>>>>> # name with values > 1 and rename the country
>> names in the electricity
>>>>>> # generation data to match
>> (e.g. "Tanzania" becomes "united republic of
>>>>>> #
>> tanzania").
>>>>>> # Once this data cleaning operation is
>> complete, repeat the above steps
>>>>>> # starting with reading
>> data into
>> elect.gen.tot.
>>>>>>
>>>>>>
>>>>>>
>> # Plot the data using ggplot2
>>>>>> world.plot <-
>> ggplot(data = world.ggmape, aes(x = long, y = lat, group =
>> group))
>>>>>>
>>>>>>
>>>>>>
>> world.plot + geom_polygon(aes(fill = y2007)) + labs(x =  
>> "Longitude", y =
>> "Latitude", fill = "TWh") + opts(title = "Net Electricity  
>> Generation in TWh for
>> 2007\nData from EIA International Energy Statistics,
>> 2008")
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> On Fri, Jun 18, 2010 at 11:39 AM, Tom Hopper <> ymailto="mailto:tomhopper at gmail.com 
>> "
>> href="mailto:tomhopper at gmail.com">tomhopper at gmail.com>
>> wrote:
>>>>>>
>> Fernando,
>>>>>>
>>>>>> That worked perfectly;
>> thank you!
>>>>>>
>>>>>> There were some
>> mismatches in the country names, as you noted, but after cleaning  
>> up the
>> electricity generation data everything looks great. I've updated  
>> the electricity
>> generation data with the cleaned version and posted a graph to >  
>> target="_blank" href="http://box.net">box.net.
>>>>>> the
>> graph:
>> http://www.box.net/tomhopper/1/22918943/452739712
>>>>>>
>>>>>>
>> Below, for reference, is the complete R
>> code.
>>>>>>
>>>>>> Thank you, and best
>> regards,
>>>>>>
>>>>>>
>> Tom
>>>>>>
>>>>>> # Download electricity
>> generation data from > href="http://tonto.eia.doe.gov/cfapps/ipdbproject/iedindex3.cfm?tid=2&pid=2&aid=12&cid=&syid=2004&eyid=2008&unit=BKWH 
>> "
>> target=_blank
>>> http://tonto.eia.doe.gov/cfapps/ipdbproject/iedindex3.cfm?tid=2&pid=2&aid=12&cid=&syid=2004&eyid=2008&unit=BKWH
>>>>>>
>> # Download new map data from > href="http://thematicmapping.org/downloads/world_borders.php 
>> " target=_blank
>>> http://thematicmapping.org/downloads/world_borders.php
>>>>>>
>>>>>>
>> # Bring the thematicmapping data into R
>>>>>>
>> library(maptools)
>>>>>> gpclibPermit()
>>>>>>
>> library("rgdal")
>>>>>> world.map <- readOGR(dsn="C:/",
>> layer="TM_WORLD_BORDERS-0.3")
>>>>>>
>>>>>> #
>> Save the map data as an R object for later use
>>>>>>
>> save(world.map,
>> "worldmap.rdata")
>>>>>>
>>>>>> # As needed,
>> reload the data
>>>>>> #
>> load("worldmap.rdata")
>>>>>>
>>>>>> # Prepare
>> the world.map data for ggplot2
>>>>>>
>> library(ggplot2)
>>>>>> world.ggmap <- fortify(world.map,
>> region = "NAME")
>>>>>>
>>>>>> # Load the
>> electricity generation data and clean it up to match with
>> world.ggmape
>>>>>> elect.gen.tot <-
>> read 
>> .csv("Total_Electricity_Net_Generation_(Billion_Kilowatthours).csv",
>>>>>>
>>      sep = ",", dec = ".")
>>>>>> names(elect) <- c("id",
>> "y2004", "y2005", "y2006", "y2007",
>> "y2008")
>>>>>>
>>>>>> elect.gen.tot$id <-
>> tolower(elect.gen.tot$id)
>>>>>>
>>>>>> #
>> merge the two data sets
>>>>>> world.ggmape <-
>> merge(world.ggmap, elect.gen.tot, by = "id", all = TRUE)
>>>>>>
>> world.ggmape <- world.ggmape[order(world.ggmape$order),
>> ]
>>>>>>
>>>>>> # NOTE: at this point, the
>> electricity data country names may not match 100%
>>>>>> # with
>> the thematicmapping country names (column "id").
>>>>>> # print
>> out the country names using
>>>>>> #
>> table(world.ggmape$id)
>>>>>> # and do a search for values of
>> 1. Then find the corresponding country
>>>>>> # name with
>> values > 1 and rename the country names in the
>> electricity
>>>>>> # generation data to match (e.g. "Tanzania"
>> becomes "United Republic of
>>>>>> #
>> Tanzania").
>>>>>> # Once this data cleaning operation is
>> complete, repeat the above steps
>>>>>> # starting with reading
>> data into elect.gen.tot.
>>>>>>
>>>>>> # Plot
>> the data using ggplot2
>>>>>> world.plot <- ggplot(data =
>> world, aes(x = long, y = lat, group = group))
>>>>>> world.plot
>> + geom_polygon(aes(fill = y2007)) +
>>>>>>      labs(x =
>> "Longitude", y = "Latitude", fill = "TWh") +
>>>>>>
>> opts(title = "Net Electricity Generation in TWh for 2007\nData from  
>> EIA
>> International Energy Statistics,
>> 2008")
>>>>>>
>>>>>>
>>>>>>
>> On Fri, Jun 18, 2010 at 2:10 AM, fernando <> ymailto="mailto:fgtaboada at gmail.com 
>> "
>> href="mailto:fgtaboada at gmail.com">fgtaboada at gmail.com>
>> wrote:
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> Hi
>> Tom,
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> I think fortify can help you in translating the sp object to
>> a
>>>>>> data.frame. Then you can use merge to join the two
>> tables. The code bellow
>>>>>> illustrates this, although I
>> think there are some problems in matching country
>>>>>> names.
>> Another issue is that if you add coord_map(), something strange
>> happens
>>>>>> with Antarctica (maybe a problem in shapefile
>> order?).
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>> --
>>>>>
>>>>>> You received this message because
>> you are subscribed to the ggplot2 mailing list.
>>>>>> Please
>> provide a reproducible example:
>> http://gist.github.com/270442
>>>>>>
>>>>>> To
>> post: email > href="mailto:ggplot2 at googlegroups.com">ggplot2 at googlegroups.com
>>>>>>
>> To unsubscribe: email ggplot2+> href="mailto:unsubscribe at googlegroups.com 
>> ">unsubscribe at googlegroups.com
>>>>>>
>> More options:
>> http://groups.google.com/group/ggplot2
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Assistant
>> Professor / Dobelman Family Junior Chair
>>>>> Department of
>> Statistics / Rice
>> University
>>>>> http://had.co.nz/
>>>>>
>>>> -- 
>>
>>>>
>>>> You received this message because you are
>> subscribed to the ggplot2 mailing list.
>>>> Please provide a
>> reproducible example: > >http://gist.github.com/270442
>>>>
>>>> To post:
>> email > href="mailto:ggplot2 at googlegroups.com">ggplot2 at googlegroups.com
>>>> To
>> unsubscribe: email ggplot2+> href="mailto:unsubscribe at googlegroups.com 
>> ">unsubscribe at googlegroups.com
>>>> More
>> options: > >http://groups.google.com/group/ggplot2
>>>>
>>>
>>
>
>
>
>>
>     [[alternative HTML version
>> deleted]]
>
>
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.