[STATSGRASS] Trend surface analysis issue - sending results back to GRASS

[Roger Bivand]

On Thu, 28 Sep 2006, Carlos "Guâno" Grohmann wrote:

> Roger, I'm still lost :)
> 
> let's see:
> 
> 
> >library(spgrass6) ; library(spatial); library(gstat); library(lattice);
> >G <- gmeta6();
> >grd <- GridTopology(cellcentre.offset=c(G$west+(G$ewres/2),
> G$south+(G$nsres/2)), cellsize=c(G$ewres, G$nsres),
> cells.dim=c(G$cols, G$rows));
> >mask_SG <- SpatialGridDataFrame(grd, data=list(k=rep(1,
> G$cols*G$rows)), proj4string=CRS(G$proj4));
> 
> >pts_trend<- getSites6sp("points_trend2");
> 
> > summary(pts_trend)
> Object of class SpatialPointsDataFrame
> Coordinates:
>                 min       max
> coords.x1  293529.5  297411.7
> coords.x2 8873992.2 8878313.7
> Is projected: TRUE
> proj4string : [+proj=utm +south +zone=24 +a=6378160 +rf=298.25
> +no_defs +towgs84=-60.0,-2.0,-41.0]
> Number of points: 31
> Data attributes:
>       cat             x                y                alt
>  Min.   : 1.0   Min.   :293530   Min.   :8873992   Min.   :514.4
>  1st Qu.: 8.5   1st Qu.:294619   1st Qu.:8874820   1st Qu.:522.8
>  Median :16.0   Median :295493   Median :8875555   Median :533.0
>  Mean   :16.0   Mean   :295542   Mean   :8875925   Mean   :536.8
>  3rd Qu.:23.5   3rd Qu.:296673   3rd Qu.:8877029   3rd Qu.:550.4
>  Max.   :31.0   Max.   :297412   Max.   :8878314   Max.   :571.3
> 
> ******************************************************************
> now gstat:
> 
> >tbv_trend2_gstat <- gstat(id="tr2", formula = alt ~ 1, data =
> pts_trend2, degree=2);
> 
> > names(tbv_trend2_gstat)
> [1] "data"  "model" "call"

?predict.gstat

OK, the Wiki uses the one-step krige() instead of gstat() then predict() 
of the output object, I guess in your case:

tr2_pred <- predict(tbv_trend2_gstat, newdata=mask_SG)
names(tr2_pred)

Instead of anova(), I suggest leave-one-out cross-validation:

pe <- gstat.cv(tbv_trend2_gstat, debug.level=0, random=FALSE)$residual
sqrt(mean(pe^2)) # RMS_LOO_PE

which may take more time for many points. Have a look at ?gstat.cv for 
some alternatives.

Please let me know if this helps,

Roger

> 
> 
> And I am here... I don't know how to display the results or how to get
> it back in GRASS, how to do analysis of variance among trend surfaces
> of different degrees... If I use surf.ls(), I can get the anova table
> but I don't know how to get the image or anything else..
> 
> 
> Carlos
> 
> 
> 
> 
> 
> 
> On 9/27/06, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > On Wed, 27 Sep 2006, Carlos "Guâno" Grohmann wrote:
> >
> > > Hello all
> > >
> > > Since the new sp and spGRASS6 packages, I haven't done any trand
> > > surface analysis yet, so I'm a little lost here..
> > >
> > > in the old day, I would do this:
> > >
> > > g1r<-surf.ls(1,gradiente$east,gradiente$north,gradiente$num1);
> > >  g1g<-trmat.G(G,g1r);
> > >  rast.put(G,lname="gradiente.trend1",layer=g1g);
> > >
> > > now that trmat.G doesen't exists anymore, how do I do to send my trend
> > > sufaces back to GRASS?
> >
> > Use the standard gstat() approach, with the degree= argument. The gstat
> > package and surf.ls() in MASS use the same coordinate shrinking approach
> > and give the same results.
> >
> > Roger
> >
> > >
> > > thanks
> > >
> > > Carlos
> > >
> > >
> >
> > --
> > 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
> >
> >
> 
> 
> 

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