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 <tomhop...@gmail.com>
>To: Felipe Carrillo <mazatlanmex...@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 <mazatlanmex...@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 <tomhop...@gmail.com>
>>>To: ggpl...@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 <had...@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 <tomhop...@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 <tomhop...@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 <fgtabo...@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 ggpl...@googlegroups.com
>>>>> To unsubscribe: email ggplot2+unsubscr...@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 ggpl...@googlegroups.com
>>>To unsubscribe: email ggplot2+unsubscr...@googlegroups.com
>>>More options: http://groups.google.com/group/ggplot2
>>>
>>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
