[R] Plotting Data on a Map

Felipe Carrillo mazatlanmexico at yahoo.com
Wed Jun 23 22:57:35 CEST 2010


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?



----- 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]]






More information about the R-help mailing list