[GRASSLIST:10738] Re: krigging result looks weird

[Roger Bivand]

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

Some comments below:

On Sat, 4 Mar 2006, Carlos "Guâno" Grohmann wrote:

> 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 attach)... 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);
> 



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