[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