[GRASSLIST:10737] krigging result looks weird

GRASS GRASS user
1

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

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

(attachments)

kigging.jpg

2

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@nhh.no