[GRASS-stats] Re: [GRASS-user] grass/stats mapping/prediction question

Wed Dec 19 19:19:35 EST 2007

On Wednesday 19 December 2007, andrew haywood wrote:
> Thanks to all that replied.
>
> everything works well now. THank you all for such an elegant way to
> model!!! Much more fun to play with GRASS and a stats package at the same
> time.
>
> I have just received my copy of Open Source GIS A GRASS GIS Approach", 3rd
> edition and at a first glance much of the issues are covered in Chapter 10.
>
> Once again thanks.
>
> Andy

Glad to hear that you are making progress. Here is one other example 
illustrating some common GRASS-R-GRASS operations with the Spearfish data.

Cheers,

Dylan

> On 12/20/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > On Wed, 19 Dec 2007, Dylan Beaudette wrote:
> > > On Wednesday 19 December 2007, Roger Bivand wrote:
> > >> On Wed, 19 Dec 2007, Daniel McInerney wrote:
> > >>
> > >> The original questioner should have written to the grass-stats list in
> >
> > the
> >
> > >> first place - thanks for CC-ing. See below for inline comments:
> > >>> Hi Andy,
> > >>>
> > >>>>  I am unsure how to then move the 'outmap' back across to grass.
> > >>>>  How can i convert the df to a spatial grid object.
> > >>>
> > >>> AFAIK, 'predict' won't create a dataframe object.
> > >>> R should return FALSE for is.data.frame(outmap) and
> > >>> TRUE for is.numeric(outmap)
> > >>>
> > >>> You can slot the model output to the AttributeList of
> > >>
> > >> It hasn't been an AttributeList for a long time - the data slot *is* a
> > >> data frame. We *always* need the output of sessionInfo() to help -
> > >> update.packages() is most often very helpful too.
> > >>
> > >>> one of the SpatialGridDataFrames that you created when
> > >>> you read in a GRASS raster and then use writeRAST6 to
> > >>> write it back to GRASS.
> > >>>
> > >>> e.g. using 'anmax' from your example
> > >>>
> > >>> anmax$anmax <- outmap
> > >>> #if that doesn't work, you might try
> > >>> anmax$anmax <- as.numeric(outmap)
> > >>> writeRAST6(anmax, "NameOfNewGRASSRaster", "anmax")
> > >>>
> > >>> Regards,
> > >>> Daniel.
> > >>>
> > >>> cc: grass-stats at lists.osgeo.org
> > >>>
> > >>> andrew haywood wrote:
> > >>>>  Dear List,
> > >>>>
> > >>>>  I am having some problems analysing some ecoligical models in grass
> > >>>> using the spgrass package through R.
> > >>>>
> > >>>>  I have 130 plot locations where i have observed presence/absence of
> >
> > a
> >
> > >>>>  species. I have followed a similar framework to the BUGSITE
> >
> > modelling
> >
> > >>>>  example from Markus's 2003 grass gis handouts (Grass 5) I have no
> > >>>>  problems constructing the model based on the 130 plots and the
> > >>>>  environmental layers from grass.
> > >>
> > >> See the OSGeo tutorial September 2006:
> >
> > http://www.foss4g2006.org/contributionDisplay.py?contribId=46&sessionId=5
> >9&
> >
> > >> confId=1
> > >>
> > >> and the OSGeo Journal note:
> > >>
> > >> http://www.osgeo.org/files/journal/final_pdfs/OSGeo_vol1_GRASS-R.pdf
> > >>
> > >> for more up-to-date information.
> > >>
> > >>>>  However, I am having problems bringing all the maps through into R
> >
> > so I
> >
> > >>>>  can make a prediction map.
> > >>>>  The region isnt too large 1600 by 800 cells at 10m resolution
> > >>>>  I can bring all the environmental layers through to R using
> >
> > readRAST6()
> >
> > >>>>  which doesnt take too much time at all.
> > >>>>
> > >>>>  However i assume I must convert the spatial grid objects into
> > >>>> dataframes to apply the predicted model function.
> > >>
> > >> No, usually not at all, since the objects have a data.frame in the
> > >> data slot, and have the standard access methods.
> > >>
> > >>>>  So I then coerce them into dataframes using as.dataframe (this
> > >>>> takes ages) I then merge all the dataframes into a single dataframe.
> > >>>> (this takes ages)
> > >>>>
> > >>>>  I then apply the model predict to the new data frame.
> > >>>>
> > >>>>  I am unsure how to then move the 'outmap' back across to grass.
> > >>>>  How can i convert the df to a spatial grid object.
> > >>>>
> > >>>>  Im thinking i must be doing something wrong. As it quite quick to
> >
> > pull
> >
> > >>>>  through the layers . But seems to take quite a lot of processing to
> >
> > get
> >
> > >>>>  the layers into a datframe appropppriate for applying the
> >
> > predictions.
> >
> > >>>>  Any help would be greatly appreciated.
> > >>>>
> > >>>>  Andy
> > >>>>
> > >>>>
> > >>>> # pull through environmental layers
> > >>>> # FAST
> > >>>>  anmax <- readRAST6("anmax", ignore.stderr=TRUE)
> > >>>>  anmin <- readRAST6("anmin", ignore.stderr=TRUE)
> > >>>>  aspect <- readRAST6("aspect", ignore.stderr=TRUE)
> > >>>>  dem10_lidar <- readRAST6("dem10_lidar", ignore.stderr=TRUE)
> > >>
> > >> Wrong, do:
> > >>
> > >> mydata <- readRAST6(c("anmax", "anmin", "aspect", "dem10_lidar"),
> > >>    ignore.stderr=TRUE)
> > >>
> > >> then the data slot of the object is a data frame. Look at
> > >>
> > >> summary(mydata)
> > >>
> > >> for a sanity check.
> > >>
> > >>>> # coerce to dataframe
> > >>>> # SLOW
> > >>>>  mypred_anmaxDF<-as.data.frame(anmax)
> > >>>>  mypred_anminDF<-as.data.frame(anmin)
> > >>>>  mypred_aspectDF<-as.data.frame(aspect)
> > >>>>  mypred_dem10_lidarDF<-as.data.frame(dem10_lidar)
> > >>>>
> > >>>> # merge into single dataframe
> > >>>> # VERY SLOW
> > >>>>  merge_tmp<-merge(mypred_anmaxDF,mypred_anminDF)
> > >>>>  rm(mypred_anmaxDF,mypred_anminDF)
> > >>>>  merge_tmp1<-merge(merge_tmp,mypred_aspectDF)
> > >>>>  rm(merge_tmp,mypred_aspectDF)
> > >>>>  mypredDF<-merge(merge_tmp1,mypred_dem10_lidarDF)
> > >>>>
> > >>>>  #apply model
> > >>
> > >> What is tree? You may need to do extra steps depending on what
> >
> > class(tree)
> >
> > >> says - if you have used rpart() or some such, you may find that
> > >>
> > >> outmap   <- predict(tree,newdata=mydata, type="class")
> > >>
> > >> works,
> > >>
> > >> or
> > >>
> > >> outmap   <- predict(tree,newdata=as(mydata, "data.frame"),
> >
> > type="class")
> >
> > >> Maybe just assign into mydata straight away:
> > >>
> > >> mydata$outmap   <- predict(tree,newdata=as(mydata, "data.frame"),
> > >>    type="class")
> > >>
> > >> given ?predict.rpart saying:
> > >>
> > >> "If 'type="class"': (for a classification tree) a factor of
> > >> classifications based on the responses."
> > >>
> > >> which looks like a vector for a vector response in the formula to
> >
> > rpart().
> >
> > >> But do check what happens if there are NA in the newdata, because the
> > >> default predict() behaviour may be to drop those observations. Look at
> > >> summary(mydata).
> > >>
> > >> Some formula-using model fitting functions just work, like lm() and
> > >> the predict() method for lm objects.
> > >>
> > >>> From there, as Daniel wrote:
> > >>
> > >> writeRAST6(mydata, "rpartpred", "outmap")
> > >>
> > >> Hope this helps,
> > >>
> > >> Roger
> > >>
> > >>>>  outmap   <- predict(tree,newdata=mypredDF, type="class")
> > >
> > > An article on this in the OSGeo newsletter might be a nice way to
> >
> > document
> >
> > > simple modeling examples with GRASS and R
> >
> > Perhaps, and then there is the section in chapter 10 in "Open Source GIS
> > A GRASS GIS Approach", 3rd edition (my copy is still on its way, but from
> > the ToC, it looks as though pages 353-363 should be very helpful).
> >
> > In fact, your site is a convenient collection of resources, I ought to
> > have mentioned it in my reply!
> >
> > Roger
> >
> > > D
> >
> > --
> > 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

-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341