Geoestadística con R. Parte 2. Geoestadística multivariada
Daniela Ballari
Septiembre de 2018
Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ
Objetivo de la práctica. Introducir los métodos de geoestadística multivariada utilizando el conjunto de datos de “Meuse”.
Esta práctica se basa en el documento: “An introduction to geostatistics with R/gstat” de D G Rossiter, Cornell University, School of Integrative Plant Science, Crop and Soil Sciences Section. June, 2018. http://www.css.cornell.edu/faculty/dgr2/teach/R/gs_intro_07Jun2018.pdf
Temario
A- Cargar librerias y Cargar datos Meuse
B- Exploración de atributos - univariado
C- Estructura espacial local: autocorrelación
D- Interpolación con kriging
E- Otros métodos de interpolación no geoestadística
F- Validación cruzada
G- Anisotropía
H- Exploración de atributos - bivariado
I- Exploración de atributos - multivariado
J- Kriging with external drift
K- Universal kriging
L- Visualizacion conjunta
Tipos de Kriging
Media constante
-simple kriging = media conocida -ordinary kriging = media desconocida -indicator krigning = datos categoricos / above - under threshold
Tendencia en la media
-co-kriging = covariables, solo disponible en ciertos puntos -universal kriging = tendencia en función de las coordenadas, covariables en todos los puntos de predicción. deterministic part is modeled as a function of coordinates -kriging with external drift= tendencia en funcion de covariables, covariables en todos los puntos de predicción. both components are predicted simultaneously. -Regresion kriging = se maneja por separado la componente deterministica y estocástica. the deterministic (regression) and stochastic (kriging) predictions are done separately
A-Cargar Librerías y Datos Meuse
Conjunto de datos “Meuse” de polución del suelo con muestras de concentración de metales pesados (paquete sp). Meuse está previamente cargado en R, de modo que solo lo llamamos con data()
library(sp)## Warning: package 'sp' was built under R version 3.4.4library(gstat)## Warning: package 'gstat' was built under R version 3.4.4data(meuse)
str(meuse) #estructura de datos meuse. Algunas variables son continuas y otras clasificadas (factors)## 'data.frame': 155 obs. of 14 variables:
## $ x : num 181072 181025 181165 181298 181307 ...
## $ y : num 333611 333558 333537 333484 333330 ...
## $ cadmium: num 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## $ copper : num 85 81 68 81 48 61 31 29 37 24 ...
## $ lead : num 299 277 199 116 117 137 132 150 133 80 ...
## $ zinc : num 1022 1141 640 257 269 ...
## $ elev : num 7.91 6.98 7.8 7.66 7.48 ...
## $ dist : num 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## $ om : num 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## $ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## $ dist.m : num 50 30 150 270 380 470 240 120 240 420 ...B-Exploración de atributos - Univariado- Analisis EDA (Exploratory Data Analysis)
meuse$zinc## [1] 1022 1141 640 257 269 281 346 406 347 183 189 251 1096 504
## [15] 326 1032 606 711 735 1052 673 402 343 218 200 194 207 180
## [29] 240 180 208 198 250 192 213 321 569 833 906 1454 298 167
## [43] 176 258 746 746 464 365 282 375 222 812 1548 1839 1528 933
## [57] 432 550 1571 1190 907 761 659 643 801 784 1060 119 778 703
## [71] 676 793 685 593 549 680 539 560 1136 1383 1161 1672 765 279
## [85] 241 317 545 505 420 332 400 553 577 155 224 180 226 186
## [99] 198 187 199 157 203 143 136 117 113 130 192 240 221 140
## [113] 128 166 191 232 203 722 210 198 139 253 703 832 262 142
## [127] 119 152 415 474 126 210 220 133 141 158 129 206 451 296
## [141] 189 154 169 403 471 612 601 783 258 214 166 496 342 162
## [155] 375sort(meuse$zinc)## [1] 113 117 119 119 126 128 129 130 133 136 139 140 141 142
## [15] 143 152 154 155 157 158 162 166 166 167 169 176 180 180
## [29] 180 183 186 187 189 189 191 192 192 194 198 198 198 199
## [43] 200 203 203 206 207 208 210 210 213 214 218 220 221 222
## [57] 224 226 232 240 240 241 250 251 253 257 258 258 262 269
## [71] 279 281 282 296 298 317 321 326 332 342 343 346 347 365
## [85] 375 375 400 402 403 406 415 420 432 451 464 471 474 496
## [99] 504 505 539 545 549 550 553 560 569 577 593 601 606 612
## [113] 640 643 659 673 676 680 685 703 703 711 722 735 746 746
## [127] 761 765 778 783 784 793 801 812 832 833 906 907 933 1022
## [141] 1032 1052 1060 1096 1136 1141 1161 1190 1383 1454 1528 1548 1571 1672
## [155] 1839hist(meuse$zinc, breaks = 16) #Distribución no simétrica, sesgada hacia la derecha
summary(meuse$zinc) # media mayor a la mediana (sesgada hacia la derecha, distribución no-normal)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 198.0 326.0 469.7 674.5 1839.0#Al ser una distribución no simétrica, se aplica logaritmo para transformar los valores y obtener una distribución simetrica (normal). Esto, además, reduce los posibles outliers.
meuse$logZn <- log10(meuse$zinc)
hist(meuse$logZn, breaks = 16)
summary(meuse$logZn)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.053 2.297 2.513 2.556 2.829 3.265head(meuse$zinc)## [1] 1022 1141 640 257 269 281head(meuse$logZn)## [1] 3.009451 3.057286 2.806180 2.409933 2.429752 2.448706head(10^(meuse$logZn))#Luego para transformar a los valores originales se utiliza antilog. antilog (x)=10^x## [1] 1022 1141 640 257 269 281C-Estructura espacial local: autocorrelación
3 Pasos - variograma empirico - variograma teorico o inicial - ajuste de variograma
coordinates(meuse) <- c("x", "y")#Crear Spatial Data Frame
plot(meuse, asp = 1, cex = 4 * meuse$zinc/max(meuse$zinc),pch = 1); data(meuse.riv);lines(meuse.riv) 
C.1 VARIOGRAMA EMPÍRICO O EXPERIMENTAL
¿Cómo la semivarianza y distancia se relacionan en toda el área de estudio? Es a través del Variograma empírico o experimental. Se define como la semivarianza media dividida para un rango de separación determinado. Por el principio de estacionariedad, la estructura espacial es la misma en toda la zona de estudio.
Graficar el variograma experimental de las concentraciones de Zinc (log). Es decir la media de las semivarianzas, respecto de la distancia media (bin o lag). En este caso de 90m, hasta una distancia máxima (cutoff) de 1300m.
ve <- variogram(logZn ~ 1, meuse, cutoff = 1300,width = 90)
plot(ve, plot.numbers = T, asp=1)
variogram: genera la nube de variograma. logZn ~ 1: logZn es dependiente de si misma –> autocorrelación. Por defecto (si no se especifica cutoff y with) el cutoff es 1/3 de la máxima distancia(diagonal del bbox). Esto es dividido en 15 clases igualmente espaciadas.
¿Cuál es la evidencia de dependencia espacial local? np=número de pares de puntos para cada una de las 15 clases, dist=distancia media, gamma=semivarianza media. A medida que la distancia aumenta, también lo hace la semivarianza, pero hasta una distancia determinada donde la semivarianza se estabiliza.
C.2 MODELO DE VARIOGRAMA O VARIOGRAMA TEÓRICO
Función que ajusta el variograma teórico. Hay diferentes tipos de funciones que pueden ser útiles:
show.vgms()
Métodos para seleccionar el modelo que mejor ajusta al variograma teórico. 1- Conocer el proceso espacial que gobierna a la variable 2- Ajuste visual 3- Ajuste automático y comparar la bondad del ajuste
C.3 AJUSTE VISUAL
Range o rango: separación o distancia entre pares de puntos en la cual ya no hay dependencia espacial. aprox 850m Nugget o pepita: semivarianza a la separación de 0m. aprox 0.01 Total-sill o meseta: semivarianza a la distancia del rango. aprox 0.13 Partial-sill o meseta parcial: total sill - nugget. aprox 0.12
vgm genera el modelo de variograma
vt <- vgm(psill = 0.12, model = "Sph", range = 850,nugget = 0.01)
plot(ve, pl = T, model = vt)
C.4 AJUSTE AUTOMÁTICO
fit.variogram: ajuta el modelo de variograma a un variograma empírico.
va <- fit.variogram(ve, vt)
plot(ve, pl = T, model = va)
Variograma inicial. Regla general para determinar parámetros automaticamente -Nugget = 0 -partial-sill = varianza total de los datos -Range = 1/4 de la distancia maxima, diagonal del bounding box
vt.i <- vgm(nugget=0, model="Sph", range=sqrt(diff(meuse@bbox["x",])^2 +
diff(meuse@bbox["y",])^2)/4,
psill=var(meuse$logZn))Ajuste
va <- fit.variogram(ve, vt.i)
plot(ve, pl = T, model = va)
D-Interpolación
Kriging ordinario
Usualmente kriging se utiliza para predecir los píxeles (o nodos) de una malla regular que cubre la zona de estudio.Kriging ordinario ordinary" significa que (1) la variable es modelada a partir de si misma; (2) la media espacial no es conocida a priori, sino estimada de los datos y no existente tendencia espacial. prediction (var1.pred) prediction variance (var1.var)
data(meuse.grid) #malla de 40m x 40m, disponible con el dataset meuse.
coordinates(meuse.grid) <- c("x", "y")
gridded(meuse.grid) <- TRUE #indica que el conjunto de datos es un raster
ok <- krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = va) ## [using ordinary kriging]ok$pred <- 10^(ok$var1.pred)#volver a valores originales
par(mfrow=c(2,1))
pts.s <- list("sp.points", meuse, col="white",pch=1, cex=4*meuse$zinc/max(meuse$zinc))
print(spplot(ok, "pred", asp=1, col.regions=rev(heat.colors(64)),
main="Predicción OK, log-ppm Zn",sp.layout = list(pts.s)),
split=c(1,1,2,1), more=TRUE)
pts.s <- list("sp.points", meuse, col="black", pch=20)
print(spplot(ok, zcol="var1.var",col.regions=rev(gray(seq(0,1,.01))), asp=1,
main="Varianza OK, log-ppm Zn^2",sp.layout = list(pts.s)),
split=c(2,1,2,1), more=FALSE)
E-Otros métodos de interpolación no-geoestadísticos
Polígonos de Thiessen: 1. Cambios abruptos de bordes. 2. Solo utiliza un punto para cada predicción.
Inverso de la distancia: 1. No hay una manera objetiva de seleccionar el peso (inverso o inveso cuadrado…). 2. No hay una manera objetiva de seleccionar el radio de interpolación.
En todos los casos, la distribución no regular puede sobre-enfatizar áreas (o lo contrario). La varianza debe ser estimada a partir de un conjunto de datos de validación, especificamente reservado para ello.
Poligonos de Thiessen
thiessen = krige(zinc ~ 1, meuse, meuse.grid, nmax = 1)## [inverse distance weighted interpolation]pts.s <- list("sp.points", meuse, col="white",pch=20)
spplot(thiessen, "var1.pred", asp=1, col.regions=rev(heat.colors(50)),
sp.layout = list(pts.s),main="Thiessen")
Inverso de la distancia (idw)
idw = idw(zinc ~ 1, meuse, meuse.grid)## [inverse distance weighted interpolation]spplot(idw, "var1.pred", asp=1, col.regions=rev(heat.colors(50)),
sp.layout = list(pts.s),main="IDW")
#Cambiando el número de vecinos
# idw$vecino<- krige(zinc ~ 1, meuse, meuse.grid,nmax=6)
#Cambiando pesos
# idw$idp05 = idw(zinc ~ 1, meuse, meuse.grid, idp = 0.5)$var1.pred
# idw$idp5 = idw(zinc ~ 1, meuse, meuse.grid, idp = 5)$var1.predF-Validación cruzada
thiessen.cv <- krige.cv(logZn ~ 1, meuse, nmax = 1)
idw.cv<- krige.cv(logZn ~ 1, meuse)
ok.cv <- krige.cv(logZn ~ 1, meuse, model = va)
thiessen.cv[1:5,]## coordinates var1.pred var1.var observed residual zscore fold
## 1 (181072, 333611) 3.057286 NA 3.009451 -0.04783475 NA 1
## 2 (181025, 333558) 3.009451 NA 3.057286 0.04783475 NA 2
## 3 (181165, 333537) 3.009451 NA 2.806180 -0.20327092 NA 3
## 4 (181298, 333484) 2.806180 NA 2.409933 -0.39624685 NA 4
## 5 (181307, 333330) 2.448706 NA 2.429752 -0.01895404 NA 5idw.cv[1:5,]## coordinates var1.pred var1.var observed residual zscore fold
## 1 (181072, 333611) 2.830957 NA 3.009451 0.1784941 NA 1
## 2 (181025, 333558) 2.797944 NA 3.057286 0.2593420 NA 2
## 3 (181165, 333537) 2.687633 NA 2.806180 0.1185473 NA 3
## 4 (181298, 333484) 2.605734 NA 2.409933 -0.1958009 NA 4
## 5 (181307, 333330) 2.505657 NA 2.429752 -0.0759050 NA 5ok.cv[1:5,]## coordinates var1.pred var1.var observed residual zscore
## 1 (181072, 333611) 2.941046 0.03408724 3.009451 0.068404690 0.37050151
## 2 (181025, 333558) 2.938452 0.03300608 3.057286 0.118833216 0.65409504
## 3 (181165, 333537) 2.735836 0.03416564 2.806180 0.070343758 0.38056671
## 4 (181298, 333484) 2.627050 0.04264164 2.409933 -0.217116957 -1.05142137
## 5 (181307, 333330) 2.425835 0.03281227 2.429752 0.003917734 0.02162802
## fold
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5ok.cv$obs <- 10^(ok.cv$observed)F.1 Scatterplots
par(mfrow=c(1,3))
print(plot(var1.pred~observed,thiessen.cv, main="Thiessen"), split=c(1,1,3,1), more=TRUE)## NULLprint(plot(var1.pred~observed,idw.cv, main="IDW"), split=c(2,1,3,1), more=TRUE)## NULLprint(plot(var1.pred~obs,ok.cv, main="OK"), split=c(3,1,2,1), more=FALSE)
## NULLF.2 Correlación, idealmente 1
cor(thiessen.cv$var1.pred,thiessen.cv$observed)## [1] 0.685839cor(idw.cv$var1.pred,idw.cv$observed)## [1] 0.7640391cor(ok.cv$var1.pred,ok.cv$obs)## [1] 0.7677732F.3 Media de residuos, idealmente 0
mean(thiessen.cv$residual)## [1] 0.003457601mean(idw.cv$residual)## [1] -0.005565866mean(ok.cv$residual)## [1] -0.0001470178F.4 Desviaciones estandar del error de interpolación (residuos). Idealmente pequeño
sd(thiessen.cv$residual)## [1] 0.2463506sd(idw.cv$residual)## [1] 0.2238086sd(ok.cv$residual)## [1] 0.1731182F.5 MSPE (mean square predictor error), idealmente pequeño
mean(thiessen.cv$residual^2)## [1] 0.06030903mean(idw.cv$residual^2)## [1] 0.0497981mean(ok.cv$residual^2)## [1] 0.02977656F.6 Error medio cuadrático (RMSE) es una medida general. Idealmente pequeño
sqrt(sum(thiessen.cv$residual^2)/length(thiessen.cv$residual))## [1] 0.245579sqrt(sum(idw.cv$residual^2)/length(idw.cv$residual))## [1] 0.2231549sqrt(sum(ok.cv$residual^2)/length(ok.cv$residual))## [1] 0.1725589F.7 Varianza de residuos, idealmente cercanos a cero
var(thiessen.cv$residual, na.rm=T)## [1] 0.06068862var(idw.cv$residual, na.rm=T)## [1] 0.05009028var(ok.cv$residual, na.rm=T)## [1] 0.02996989F.8 Cantidad de variación explicada por el modelo en cada método, idealmente cercano a 100
(1-var(thiessen.cv$residual, na.rm=T)/var(meuse$zinc))*100## [1] 99.99995(1-var(idw.cv$residual, na.rm=T)/var(meuse$zinc))*100## [1] 99.99996(1-var(ok.cv$residual, na.rm=T)/var(meuse$zinc))*100## [1] 99.99998G-Anisotropia o variogramas direccionales
https://rpubs.com/nnnagle/vario. Patron en los datos (suroeste - noreste). Esto viola el principio de estacionalidad (stationarity). La estructura espacial puede ser diferente en distintas direcciones. Por ello se genera un variograma en 4 direcciones: 0, 45, 90, 135
Variograma con anisotropía
vgm.aniso <- variogram(logZn ~ 1, meuse, alpha = c(0, 45, 90, 135))
plot(vgm.aniso)
Autocorrelación mas fuerte a los 45 y menos fuerte a los 135.
Mapa de variograma
vgm.map = variogram(logZn~1, meuse, cutoff = 1500, width = 100,map = TRUE)
# vgm.map = variogram(logZn~sqrt(dist), meuse, cutoff = 1500, width = 100,map = TRUE)
plot(vgm.map, threshold = 5)
Ajuste de variograma
vt.i.a <- vgm(nugget=0, model="Sph", range=sqrt(diff(meuse@bbox["x",])^2 +
diff(meuse@bbox["y",])^2)/4, psill=var(meuse$logZn), anis = c(45, 0.4))
#In two dimensions, two parameters define an anisotropy ellipse, say anis = c(30, 0.5). The first parameter, 30, refers to the main axis direction: it is the angle for the principal direction of continuity (measured in degrees, clockwise from positive Y, i.e. North). The second parameter, 0.5, is the anisotropy ratio, the ratio of the minor range to the major range (a value between 0 and 1). So, in our example, if the range in the major direction (North-East) is 100, the range in the minor direction (South-East) is 0.5 x 100 = 50.
# In three dimensions, anis = c(p,q,r,s,t). $p$ is the angle for the principal direction of continuity (measured in degrees, clockwise from Y, in direction of X), $q$ is the dip angle for the principal direction of continuity (measured in positive degrees up from horizontal), $r$ is the third rotation angle to rotateanis = c(30, 10, 0, 0.5, 0.3))
va.a <- fit.variogram(vgm.aniso, vt.i.a)
va.a## model psill range ang1 anis1
## 1 Nug 0.0111140 0.00 0 1.0
## 2 Sph 0.1102725 1742.85 45 0.4plot(vgm.aniso, pl = T, model = va.a)
A partir de aqui se procede con kriging pero utilizando este variograma
Predicción
ok <- krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = va.a) ## [using ordinary kriging]ok$pred <- 10^(ok$var1.pred)#volver a valores originales
par(mfrow=c(2,1))
pts.s <- list("sp.points", meuse, col="white",pch=1, cex=4*meuse$zinc/max(meuse$zinc))
print(spplot(ok, "pred", asp=1, col.regions=rev(heat.colors(64)),
main="Predicción OK, log-ppm Zn",sp.layout = list(pts.s)),
split=c(1,1,2,1), more=TRUE)
pts.s <- list("sp.points", meuse, col="black", pch=20)
print(spplot(ok, zcol="var1.var",col.regions=rev(gray(seq(0,1,.01))), asp=1,
main="Varianza OK, log-ppm Zn^2",sp.layout = list(pts.s)),
split=c(2,1,2,1), more=FALSE)
H-Exploración de atributos - Bivariado - Analisis EDA (Exploratory Data Analysis)
Análisis de correlación (https://www0.gsb.columbia.edu/premba/analytical/s7/s7_5.cfm)
Distancia al rio
cor(meuse$logZn, meuse$dist)## [1] -0.7394276plot(meuse$logZn ~ meuse$dist)
Regresión lineal (https://www0.gsb.columbia.edu/premba/analytical/s7/s7_6.cfm)
m.lzn.dist <- lm(logZn ~ dist, data = meuse)
summary(m.lzn.dist)##
## Call:
## lm(formula = logZn ~ dist, data = meuse)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.4889 -0.1583 -0.0007 0.1387 0.7286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.83759 0.02680 105.87 <2e-16 ***
## dist -1.17256 0.08631 -13.59 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2118 on 153 degrees of freedom
## Multiple R-squared: 0.5468, Adjusted R-squared: 0.5438
## F-statistic: 184.6 on 1 and 153 DF, p-value: < 2.2e-16plot(meuse$logZn ~ meuse$dist)
abline(m.lzn.dist)
I-Exploración de atributos - multivariado - Analisis EDA (Exploratory Data Analysis)
Modelo aditivo, asume que las variables actuan lineal e independientemente
meuse$logCu <- log10(meuse$copper)
m.lzn.lcu.dist <- lm(logZn ~ logCu + dist, data = meuse)
summary(m.lzn.lcu.dist)##
## Call:
## lm(formula = logZn ~ logCu + dist, data = meuse)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54110 -0.07904 -0.00324 0.08194 0.29847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.06843 0.10500 10.175 < 2e-16 ***
## logCu 1.02805 0.06032 17.042 < 2e-16 ***
## dist -0.41777 0.06736 -6.202 5.04e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1245 on 152 degrees of freedom
## Multiple R-squared: 0.8443, Adjusted R-squared: 0.8422
## F-statistic: 412.1 on 2 and 152 DF, p-value: < 2.2e-16Modelo de interacción
m.lzn.lcu.dist.i <- lm(logZn ~ logCu*dist, data = meuse)
summary(m.lzn.lcu.dist.i)##
## Call:
## lm(formula = logZn ~ logCu * dist, data = meuse)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53466 -0.07806 -0.00211 0.08016 0.29964
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.04958 0.11446 9.170 3.27e-16 ***
## logCu 1.04049 0.06736 15.447 < 2e-16 ***
## dist -0.24777 0.41044 -0.604 0.547
## logCu:dist -0.12031 0.28653 -0.420 0.675
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1249 on 151 degrees of freedom
## Multiple R-squared: 0.8445, Adjusted R-squared: 0.8414
## F-statistic: 273.3 on 3 and 151 DF, p-value: < 2.2e-16J-Kriging with external drift, “drift” being the value of the covariable
Predicción de concentración de metales de acuerdo con la distancia al río.
vr <- variogram(logZn ~ dist, location=meuse,cutoff=1300, width=90)
plot(vr, plot.numbers=T, main="Residuals, dist co-variable")
vked <- fit.variogram(vr, vgm(nugget=0, model="Sph", range=sqrt(diff(meuse@bbox["x",])^2 + diff(meuse@bbox["y",])^2)/4, psill=var(meuse$logZn)))
plot(vr, plot.numbers=T, model=vked, main="Dist co-variable")
ked <- krige(logZn ~ dist, locations=meuse,newdata=meuse.grid, model=vked)## [using universal kriging]ked$pred <- 10^(ked$var1.pred)#volver a valores originales
par(mfrow=c(2,1))
pts.s <- list("sp.points", meuse, col="white",pch=1, cex=4*meuse$zinc/max(meuse$zinc))
print(spplot(ked, "pred", asp=1, col.regions=rev(heat.colors(64)),
main="Predicción KED, log-ppm Zn",sp.layout = list(pts.s)),
split=c(1,1,2,1), more=TRUE)
pts.s <- list("sp.points", meuse, col="black", pch=20)
print(spplot(ked, zcol="var1.var",col.regions=rev(gray(seq(0,1,.01))), asp=1,
main="Varianza KED, log-ppm Zn^2",sp.layout = list(pts.s)),
split=c(2,1,2,1), more=FALSE)
K-Universal Kriging
Predicción de tendencia espacial geográfica (x e y) además de estructura local.
vr <- variogram(logZn ~ dist + x + y, location=meuse,cutoff=1300, width=90)
plot(vr, plot.numbers=T, main="Residuals, dist co-variable")
vuk <- fit.variogram(vr, vgm(nugget=0, model="Sph",range=sqrt(diff(meuse@bbox["x",])^2 + diff(meuse@bbox["y",])^2)/4, psill=var(meuse$logZn)))
plot(vr, plot.numbers=T, model=vuk, main="Dist co-variable")
uk <- krige(logZn ~ dist + x + y, locations=meuse,newdata=meuse.grid, model=vuk)## [using universal kriging]uk$pred <- 10^(uk$var1.pred)#volver a valores originales
par(mfrow=c(2,1))
pts.s <- list("sp.points", meuse, col="white",pch=1, cex=4*meuse$zinc/max(meuse$zinc))
print(spplot(uk, "pred", asp=1, col.regions=rev(heat.colors(64)),
main="Predicción UK, log-ppm Zn",sp.layout = list(pts.s)),
split=c(1,1,2,1), more=TRUE)
pts.s <- list("sp.points", meuse, col="black", pch=20)
print(spplot(uk, zcol="var1.var",col.regions=rev(gray(seq(0,1,.01))), asp=1,
main="Varianza uK, log-ppm Zn^2",sp.layout = list(pts.s)),
split=c(2,1,2,1), more=FALSE)
K-Visualización cojunta, en la misma escala
zmax <- max(ok$pred,ked$pred, uk$pred)
zmin <- min(ok$pred,ked$pred, uk$pred)
zmax <- round(zmax, 1) + 0.1
zmin <- round(zmin, 1) - 0.1
ramp <- seq(from=zmin, to=zmax, by=50)
p1 <- spplot(ok, "pred", asp=1, main="OK prediction",at=ramp, col.regions=bpy.colors(64))
p2 <- spplot(ked, "pred", asp=1, main="KED-prediction",at=ramp, col.regions=bpy.colors(64))
p3 <- spplot(uk, "pred", asp=1, main="UK-prediction",at=ramp, col.regions=bpy.colors(64))
plot(p1, split = c(1, 1, 3, 1), more = T)
plot(p2, split = c(2, 1, 3, 1), more = T)
plot(p3, split = c(3, 1, 3, 1), more = F)
Mapas de varianzas
summary(ok$var1.var)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01733 0.02517 0.02944 0.03289 0.03647 0.08100summary(ked$var1.var)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01978 0.02389 0.02584 0.02777 0.02980 0.05236summary(uk$var1.var)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01993 0.02382 0.02568 0.02761 0.02948 0.05155zmax <- max(ok$var1.var,ked$var1.var, uk$var1.var)
zmin <- min(ok$var1.var,ked$var1.var, uk$var1.var)
zmax <- round(zmax, 1) + 0.1
zmin <- round(zmin, 1) - 0.1
ramp <- seq(from=zmin, to=zmax, by=0.005)
p1 <- spplot(ok, "var1.var", asp=1, main="OK varianza",at=ramp, col.regions=cm.colors(64))
p2 <- spplot(ked, "var1.var", asp=1, main="KED-ffreq varianza",at=ramp, col.regions=cm.colors(64))
p3 <- spplot(uk, "var1.var", asp=1, main="UK-ffreq varianza",at=ramp, col.regions=cm.colors(64))
plot(p1, split = c(1, 1, 3, 1), more = T)
plot(p2, split = c(2, 1, 3, 1), more = T)
plot(p3, split = c(3, 1, 3, 1), more = F)
Actividad Individual
Compruebe con validación cruzada el método de kriging que mejor predice la variable interpolada. Compruebe además incorporando anisotropía en el variograma.