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

[Roger Bivand]

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")
>>
>>
>>
>>  ------------------------------------------------------------------------
>>
>>  _______________________________________________
>>  grass-user mailing list
>>  grass-user at lists.osgeo.org
>>  http://lists.osgeo.org/mailman/listinfo/grass-user
> _______________________________________________
> grass-stats mailing list
> grass-stats at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-stats
>

-- 
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