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

[Roger Bivand]

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

> Just another doubt? what the surface with degree=0 mean? it look like
> a membrane that passes through all the points..

A flat surface at the mean?

Roger

> 
> carlos
> 
> 
> 
> 
> On 9/28/06, Carlos Guâno Grohmann <carlos.grohmann at gmail.com> wrote:
> > Thanks Roger, that's what I needed (BTW, sorry for the other mail, I
> > should sent it to the list, not just for you).
> >
> > I guess I will user the surfaces from gstat and anova with surf.ls.
> >
> > thanks again
> >
> > Carlos
> >
> >
> >
> > On 9/28/06, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > > 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
> > >
> > >
> >
> >
> > --
> > +-----------------------------------------------------------+
> >               Carlos Henrique Grohmann - Guano
> >   Geologist M.Sc  - Doctorate Student at IGc-USP - Brazil
> > Linux User #89721  - carlos dot grohmann at gmail dot com
> > +-----------------------------------------------------------+
> > _________________
> > "Good morning, doctors. I have taken the liberty of removing Windows
> > 95 from my hard drive."
> > --The winning entry in a "What were HAL's first words" contest judged
> > by 2001: A SPACE ODYSSEY creator Arthur C. Clarke
> >
> 
> 
> 

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