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

Carlos "Guâno" Grohmann carlos.grohmann at gmail.com
Thu Sep 28 14:52:59 EDT 2006


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




More information about the grass-stats mailing list