[R] rgdal-maptools

Felipe Carrillo mazatlanmexico at yahoo.com
Thu Jun 24 03:42:04 CEST 2010


The shapefile data can be downloaded from the link below:
download all the six files and save them on your working directory
and make sure the dsn path is set to where the files are saved.
My shapefiles are saved on C:/Data.


https://secure.filesanywhere.com/fs/v.aspx?v=897263875a6472a99baa

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 to 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:/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 <- 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.

 
Felipe D. Carrillo
Supervisory Fishery Biologist
Department of the Interior
US Fish & Wildlife Service
California, USA






More information about the R-help mailing list