Re: Plotting Data on a Map

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Re: Plotting Data on a Map

Felipe Carrillo
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 <[hidden email]>
>To: Felipe Carrillo <[hidden email]>
>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 <[hidden email]> 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 <[hidden email]>
>>>To: [hidden email]
>>>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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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 [hidden email]
>>>>> To unsubscribe: email [hidden email]
>>>>> 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 [hidden email]
>>>To unsubscribe: email [hidden email]
>>>More options: http://groups.google.com/group/ggplot2
>>>
>>
>

     
        [[alternative HTML version deleted]]


______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Plotting Data on a Map

Felipe Carrillo
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 <[hidden email]>
> To: Tom Hopper <[hidden email]>
> Cc: [hidden email]; [hidden email]
> 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:[hidden email]">[hidden email]>
>To: Felipe
> Carrillo <> href="mailto:[hidden email]">[hidden email]>
>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:[hidden email]"
> href="mailto:[hidden email]">[hidden email]>
> 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:[hidden email]">[hidden email]>
>>>To:
> > href="mailto:[hidden email]">[hidden email]
>>>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:[hidden email]"
> href="mailto:[hidden email]">[hidden email]>
> 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:[hidden email]">[hidden email]>
> 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:[hidden email]"
> href="mailto:[hidden email]">[hidden email]>
> 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:[hidden email]"
> href="mailto:[hidden email]">[hidden email]>
> 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:[hidden email]">[hidden email]
>>>>>
> To unsubscribe: email ggplot2+> href="mailto:[hidden email]">[hidden email]
>>>>>
> 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:[hidden email]">[hidden email]
>>>To
> unsubscribe: email ggplot2+> href="mailto:[hidden email]">[hidden email]
>>>More
> options: > >http://groups.google.com/group/ggplot2
>>>
>>
>


 
>    
    [[alternative HTML version
> deleted]]




______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Plotting Data on a Map

David Winsemius

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 <[hidden email]>
>> To: Tom Hopper <[hidden email]>
>> Cc: [hidden email]; [hidden email]
>> 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:[hidden email]">[hidden email]>
>> To: Felipe
>> Carrillo <> href="mailto:[hidden email]">[hidden email]
>> >
>> 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:[hidden email]
>> "
>> href="mailto:[hidden email]">[hidden email]>
>> 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:[hidden email]">[hidden email]>
>>>> To:
>>> href="mailto:[hidden email]">[hidden email]
>>>> 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:[hidden email]
>> "
>> href="mailto:[hidden email]">[hidden email]>
>> 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:[hidden email]
>> ">[hidden email]>
>> 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:[hidden email]
>> "
>> href="mailto:[hidden email]">[hidden email]>
>> 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:[hidden email]
>> "
>> href="mailto:[hidden email]">[hidden email]>
>> 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:[hidden email]">[hidden email]
>>>>>>
>> To unsubscribe: email ggplot2+> href="mailto:[hidden email]
>> ">[hidden email]
>>>>>>
>> 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:[hidden email]">[hidden email]
>>>> To
>> unsubscribe: email ggplot2+> href="mailto:[hidden email]
>> ">[hidden email]
>>>> More
>> options: > >http://groups.google.com/group/ggplot2
>>>>
>>>
>>
>
>
>
>>
>     [[alternative HTML version
>> deleted]]
>
>
>
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Plotting Data on a Map

Brandon Hurr
In reply to this post by Felipe Carrillo
Upload it to GoogleDocs <http://docs.google.com> or
Dropbox<http://dropbox.com>and give us a public link to the dataset.

HTH,

Brandon

On Wed, Jun 23, 2010 at 21:57, Felipe Carrillo <[hidden email]>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?
>
>
>
> ----- Original Message ----
> > From: Felipe Carrillo <[hidden email]>
> > To: Tom Hopper <[hidden email]>
> > Cc: [hidden email]; [hidden email]
> > 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:[hidden email]">[hidden email]>
> >To: Felipe
> > Carrillo <> href="mailto:[hidden email]">
> [hidden email]>
> >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:
> [hidden email]"
> > href="mailto:[hidden email]">[hidden email]>
> > 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:[hidden email]">[hidden email]>
> >>>To:
> > > href="mailto:[hidden email]">[hidden email]
> >>>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:
> [hidden email]"
> > href="mailto:[hidden email]">[hidden email]>
> > 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:[hidden email]
> ">[hidden email]>
> > 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:
> [hidden email]"
> > href="mailto:[hidden email]">[hidden email]>
> > 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:
> [hidden email]"
> > href="mailto:[hidden email]">[hidden email]>
> > 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:[hidden email]">
> [hidden email]
> >>>>>
> > To unsubscribe: email ggplot2+> href="mailto:
> [hidden email]">[hidden email]
> >>>>>
> > 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:[hidden email]">[hidden email]
> >>>To
> > unsubscribe: email ggplot2+> href="mailto:[hidden email]">
> [hidden email]
> >>>More
> > options: > >http://groups.google.com/group/ggplot2
> >>>
> >>
> >
>
>
>
> >
>     [[alternative HTML version
> > deleted]]
>
>
>
>
> --
> 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 [hidden email]
> To unsubscribe: email [hidden email]<ggplot2%[hidden email]>
> More options: http://groups.google.com/group/ggplot2
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] 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.