[STATSGRASS] Re: [GRASSLIST:10737] krigging result looks weird

[Carlos "Guâno" Grohmann]

> Can we move this to STATGRASS, probably a more suitable list?

OK.

>
> > Hello again
> >
> > I was playing with R, I get a SRTM (90m) raster, converted it to
> > vector, put it into R, calculate variogram and krigged a new map, with
> > 30m resolution. The final maps looks "pixelate" (see it in www.igc.usp.br/pessoais/guano/tempo/krigging.jpg)... Maybe
> > the parameters I used?
> >
> > here's what I did:
> >
> > 1-r.to.vect (90m)
> > 2-g.region res=30
> >
> > R:
> >
> > library(spgrass6) ; library(spatial);
> > srtm <- getSites6sp("toposrtm");
> > coords<-coordinates(srtm);
> > class(srtm);
> > 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));
> > class(mask_SG);
> >
>
> Why do the trend surface outside gstat?
> Is any other resampling going on?
> Is the artefact in the trend surface residuals?
> Is the artefact visible in image(resid.img)?
> Is the artefact visible in image(OK_pred, "var1.pred")?
>
> > trend <- surf.ls(1, x=coords[,1], y=coords[,2], z=srtm$cat);
> > resid1<-residuals(trend);
> > resid2<-as.data.frame(resid1);
> > resid.trend<-SpatialPointsDataFrame(coords,resid2);
> > resid.img <- SpatialGridDataFrame(grd, resid.trend);
> >
> > library(gstat);
> >
> > variog1 <- variogram(resid.trend$resid1 ~ 1, locations=resid.trend,
> > width=90, cutoff=1200);
> > plot(variog1);
> >
> > vrg.eye<-(vgm(psill=7500, model="Gau", range=420, nugget=0));
> > plot(variog1, model=vrg.eye);
> >
> > vrg.fit<-fit.variogram(variog1,vrg.eye);
> > vrg.fit;
> >
> >   model     psill    range
> > 1   Nug  173.7465   0.0000
> > 2   Gau 7661.1776 442.5015
> >
> > plot(variog1, model=vrg.fit);
> > vrg.eye2<-(vgm(psill=7650, model="Gau", range=440, nugget=170));
> > plot(variog1, model=vrg.eye2);
> >
> > OK_pred <- krige(srtm$cat ~ 1, loc=srtm, newdata=mask_SG,
> > model=vrg.eye2, maxdist=270);
> > names(OK_pred);
> > writeRast6sp(OK_pred, "srtm.krig", zcol="var1.pred", NODATA=-9999);
> >
>
>


OK. I did it only in gstat now:

>variog1 <-variogram(srtm$cat~coords[,1]+coords[,2], loc=srtm, srtm);

I don't see the artifact in

>image(OK_pred, "var1.pred") (although I don't know how to zoom
in...the area is bigger than in the sample image I send)

The only resample was when I changed the GRASS region from 90m
resolution to 30m. The boundaries are the same.

image in: www.igc.usp.br/pessoais/guano/tempo/krigging.jpg

thanks

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