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

[Dylan Beaudette]

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=59&
>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

D


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