[STATSGRASS] Trend surface analysis issue - sending results back
to GRASS
Carlos "Guâno" Grohmann
carlos.grohmann at gmail.com
Thu Sep 28 15:08:00 EDT 2006
Just another doubt? what the surface with degree=0 mean? it look like
a membrane that passes through all the points..
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
>
--
+-----------------------------------------------------------+
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