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

[Roger Bivand]

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

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