[R-sig-Geo] Regression Kriging cross validation

Michele Fiori mfiori at arpa.sardegna.it
Tue Apr 29 13:33:26 CEST 2014


Hi everyone,
I am working on rainfall interpolation using regression kriging method and I
need suggestions on how I can carry out a cross validation (leave-one-out)
for this elaboration. At first I tried to apply directly Krige.cv, similarly
to UK method (example for october: PP10uk.cv <- krige.cv(reg, prec2,
PP10.vgm)), but unfortunately when I applied Universal Kriging on the same
data, I realized that UK map was a little different from RK map.
So my question is: How could I manage universal kriging in order to make it
equivalent to regression kriging and use the above cross-validation, or is
there 
another different method to apply cross validation (leave-one-out) on
Regression Kriging interpolation?  
Below my code:
Many thanks

Michele Fiori

ARPAS - Environmental Protection Agency of Sardinia
MeteoClimatic Department - Meteorological Service

Viale Porto Torres 119 - 07100 Sassari, Italy
Tel + 39 079 258617
Fax + 39 079 262681
www.sardegnaambiente.it/arpas

####	Creating SpatialPixelDataFrame ("dem" - 250x250 m grid)
	....
####	Loading Precipitation data
   	prec2 <- read.table("prec2.txt", sep="\t", header =TRUE)
  	coordinates(prec2) <- c("x", "y")
   	proj4string(prec2) <- CRS("+init=epsg:32632")
####	Linear regression Model
  	mod.gen <- lm(PP10 ~ QUOTA_MARE + UTM_EST + UTM_NORD + DIST_MARE,
prec2)
  	step1 <- stepAIC(mod.gen, direction="both")
  	reg <- formula(step1)
   	PP10.lm <- lm(reg, prec2)
  	summary(PP10.lm)
   	prec2$residuals <- residuals(PP10.lm)
    	dem$predlm <- predict(PP10.lm, dem)
####	Variogram of residuals
	PP10.vgm <- vgm(nugget=51.46, model="Sph", range=38038.89,
psill=86.44)
####	Ordinary Kriging of residuals
	PP10.okr <- krige(PP10.lm$residuals ~ 1, prec2, dem, PP10.vgm,
maxdist=Inf)
	dem$varokr <- PP10.okr$var1.pred
####	Regression Kriging (Linear Regression + Ordinary Kriging of
residuals)
	dem$vark <- dem$predlm + dem$varokr
####	Universal kriging
	PP10.uk <- krige(reg, prec2, dem, PP10.vgm, maxdist=Inf)
	dem$varuk <- PP10.uk$var1.pred
	dem$difference <- dem$vark - dem$varuk
	spplot(dem[c("difference")], col.regions=terrain.colors(25),
contour=FALSE, cuts = 15)



More information about the R-sig-Geo mailing list