[R-sig-Geo] Predict and plot spatial data with ggplot

[Roger Bivand]

On Thu, 24 Jun 2010, Felipe Carrillo wrote:

> Thanks Roger for letting me know about this list.

Please re-formulate your question using data that ships with packages on 
CRAN. Trying things on data other than your own very often solves the 
problem - you get to look at the problem from a different perspective.

Please do not use ggplot until you have sorted out data handling. I for 
example have never used ggplot, so cannot answer questions involving its 
use (my bad of course, but life is tooo short...). Use sp plot and image 
methods until you are sure that you know that the data handling issues are 
resolved, simplifying as much as necessary does no harm.

Roger

>
> 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
> 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
> # Based on data from WOLVES_99 want to predict data for WOLVES_01.
> #Not sure how to proceed from here to fit the 'mod1' model to all the 5 states and show the prediction on map.
> # I am checking if predict()  would do it but so far no success.
> # I wanted to use ggplot2 but spplot would be fine. Thanks
>
>  
> Felipe D. Carrillo
> Supervisory Fishery Biologist
> Department of the Interior
> US Fish & Wildlife Service
> California, USA
>
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no