Re: [R] Plotting Data on a Map

From: Felipe Carrillo <mazatlanmexico_at_yahoo.com>
Date: Wed, 23 Jun 2010 13:57:35 -0700 (PDT)

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?

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@yahoo.com">mazatlanmexico@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@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@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@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@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@googlegroups.com">unsubscribe@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@googlegroups.com">ggplot2@googlegroups.com
>>>To 

> unsubscribe: email ggplot2+> href="mailto:unsubscribe@googlegroups.com">unsubscribe@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. Received on Wed 23 Jun 2010 - 20:59:36 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 24 Jun 2010 - 07:40:34 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive