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

[Roger Bivand]

On Tue, 7 Mar 2006, Carlos "Guâno" Grohmann wrote:

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

Before trying alternatives for moving the raster back to GRASS, maybe we 
should zoom in and increase the number of breaks from the default. You can 
use the xlim= and ylim= arguments to view a part of the image. To get more 
breaks, give it for example col=topo.colors(80).

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

OK, the input data are points from the 90m grid, and the region read from 
GRASS is for a 30m grid. 

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

This is showing the data passed back to GRASS when just using gstat, so 
the problem is still there, right? Am I wrong that the "off" pixels are 
(roughly) on the 90m spots? What happens if you overplot the d.rast with a 
d.vect toposrtm? do the "off" pixels match? Does the same happen if you 
overplot in R?

We can try other intermediate files later on, but I've a feeling that 
something is happening to the predictions at the fit points. For a 
subscene, we could use overlay() to pull out the predictions at the data 
points, but we need to see the same artefact in R first - GRASS is using 
many more colours than R, so the details may be there anyway.

Best wishes,

Roger



> 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
> 
> _______________________________________________
> statsgrass mailing list
> statsgrass at grass.itc.it
> http://grass.itc.it/mailman/listinfo/statsgrass
> 

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