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