Re: [R] Plotting Data on a Map

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

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 <tomhopper@gmail.com>
>To: Felipe Carrillo <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 <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 <tomhopper_at_gmail.com>
>>>To: 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 <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 <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 <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 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 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(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 <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 ggplot2_at_googlegroups.com
>>>>> To unsubscribe: email ggplot2+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 ggplot2_at_googlegroups.com
>>>To unsubscribe: email ggplot2+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. Received on Wed 23 Jun 2010 - 20:55:22 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 Wed 23 Jun 2010 - 22:10: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