Objetivo de la practica. Introducir los metodos de geoestadistica multivariada utilizando el conjunto de datos de “Meuse”.

Esta practica 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. Febr, 2022. http://www.css.cornell.edu/faculty/dgr2/_static/files/R_PDF/gs_intro_sf.pdf

Otros materiales del mismo autor: http://www.css.cornell.edu/faculty/dgr2/tutorials/index.html

Hengl, T. (2009). A practical guide to geostatistical mapping. http://spatial-analyst.net/book/system/files/Hengl_2009_GEOSTATe2c1w.pdf

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 funcion de las coordenadas, covariables en todos los puntos de prediccion. deterministic part is modeled as a function of coordinates
kriging with external drift= tendencia en funcion de covariables, covariables en todos los puntos de prediccion. both components are predicted simultaneously.
Regresion kriging = se maneja por separado la componente deterministica y estocastica. the deterministic (regression) and stochastic (kriging) predictions are done separately

A-Cargar Librerias y Datos Meuse

Conjunto de datos “Meuse” de polucion del suelo con muestras de concentracion de metales pesados (paquete sp). Meuse esta previamente cargado en R, de modo que solo lo llamamos con data()

library(sp)
library(gstat)
library(sf)

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_FileGDB.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_FileGDB.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_MSSQLSpatial.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_MSSQLSpatial.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_OCI.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_OCI.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_PG.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_PG.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_FileGDB.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_FileGDB.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_MSSQLSpatial.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_MSSQLSpatial.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_OCI.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_OCI.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_PG.dll
## 126: No se puede encontrar el módulo especificado.

## Warning in CPL_gdal_init(): GDAL Error 1: Can't load requested DLL: C:\Program Files\GeoDa Software\ogr_PG.dll
## 126: No se puede encontrar el módulo especificado.

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

data(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 ...

meuse.sf <- st_as_sf(meuse, coords = c("x","y"))
st_crs(meuse.sf) <- 28992

B-Exploracion de atributos

Univariado- Analisis EDA (Exploratory Data Analysis)

meuse.sf$zinc

##   [1] 1022 1141  640  257  269  281  346  406  347  183  189  251 1096  504  326
##  [16] 1032  606  711  735 1052  673  402  343  218  200  194  207  180  240  180
##  [31]  208  198  250  192  213  321  569  833  906 1454  298  167  176  258  746
##  [46]  746  464  365  282  375  222  812 1548 1839 1528  933  432  550 1571 1190
##  [61]  907  761  659  643  801  784 1060  119  778  703  676  793  685  593  549
##  [76]  680  539  560 1136 1383 1161 1672  765  279  241  317  545  505  420  332
##  [91]  400  553  577  155  224  180  226  186  198  187  199  157  203  143  136
## [106]  117  113  130  192  240  221  140  128  166  191  232  203  722  210  198
## [121]  139  253  703  832  262  142  119  152  415  474  126  210  220  133  141
## [136]  158  129  206  451  296  189  154  169  403  471  612  601  783  258  214
## [151]  166  496  342  162  375

hist(meuse.sf$zinc, breaks = 16) #Distribucion no simetrica, sesgada hacia la derecha

#Al ser una distribucion no simetrica, se aplica logaritmo para transformar los valores y obtener una distribucion simetrica (normal). Esto, ademas, reduce los posibles outliers.
meuse.sf$logZn <- log10(meuse.sf$zinc)
hist(meuse.sf$logZn, breaks = 16)

# summary(meuse.sf$logZn)
# head(meuse.sf$zinc)
# head(meuse.sf$logZn)
# head(10^(meuse.sf$logZn))#Luego para transformar a los valores originales se utiliza antilog. antilog (x)=10^x

C-Autocorrelacion

Estructura espacial local

3 Pasos
- variograma empirico
- variograma teorico o inicial
- ajuste de variograma

data(meuse.riv, package="sp")
meuse.riv.sf <- st_sfc(st_linestring(meuse.riv, dim="XY"), crs = st_crs(meuse.sf))

plot(meuse.sf["zinc"],
reset = FALSE,
nbreaks = 64, pch = 20,
cex=4*meuse.sf$zinc/max(meuse.sf$zinc),
main = "Zn concentration [ppm]")
plot(meuse.riv.sf, add = TRUE)

C.1 VARIOGRAMA EMPIRICO

Como la semivarianza y distancia se relacionan en toda el area de estudio? Es a traves del Variograma empirico o experimental. Se define como la semivarianza media dividida para un rango de separacion 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 maxima (cutoff) de 1300m.

ve <- variogram(logZn ~ 1, meuse.sf, 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 –> autocorrelacion. Por defecto (si no se especifica cutoff y with) el cutoff es 1/3 de la maxima distancia(diagonal del bbox). Esto es dividido en 15 clases igualmente espaciadas.

Cual es la evidencia de dependencia espacial local? np=numero de pares de puntos para cada una de las 15 clases, dist=distancia media, gamma=semivarianza media. A medida que la distancia aumenta, tambien lo hace la semivarianza, pero hasta una distancia determinada donde la semivarianza se estabiliza.

C.2 VARIOGRAMA TEORICO

Funcion que ajusta el variograma teorico. Hay diferentes tipos de funciones que pueden ser utiles:

show.vgms()

Metodos para seleccionar el modelo que mejor ajusta al variograma teorico. 1- Conocer el proceso espacial que gobierna a la variable 2- Ajuste visual 3- Ajuste automatico y comparar la bondad del ajuste

C.3 AJUSTE VISUAL

Range o rango: separacion o distancia entre pares de puntos en la cual ya no hay dependencia espacial. aprox 850m Nugget o pepita: semivarianza a la separacion 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 AUTOMATICO

fit.variogram: ajuta el modelo de variograma a un variograma empirico.

va <- fit.variogram(ve, vt) 
plot(ve, pl = T, model = va)

Variograma inicial. Regla general para determinar parametros 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(st_bbox(meuse.sf)[c("xmax", "xmin")])^2 +
                 diff(st_bbox(meuse.sf)[c("ymax", "ymin")])^2)/4, 
                 psill=var(meuse.sf$logZn))

Ajuste

va <- fit.variogram(ve, vt.i) 
plot(ve, pl = T, model = va)

D-Interpolacion

Kriging ordinario

Usualmente kriging se utiliza para predecir los pixeles (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.
meuse.grid.sf <- st_as_sf(meuse.grid, coords = c("x", "y"))
st_crs(meuse.grid.sf) <- st_crs(meuse.sf)

ok <- krige(logZn ~ 1, locations = meuse.sf, newdata = meuse.grid.sf, model = va)

## [using ordinary kriging]

ok$pred <- 10^(ok$var1.pred)#volver a valores originales

plot(ok["var1.pred"], pch=15,
nbreaks = 64,
main="OK prediction [log10(Zn ppm)]")

plot(ok["var1.var"], pch=15,
nbreaks = 64,
pal = cm.colors,
main="OK prediction variance [log10(Zn ppm)^2]")

E-Metodos no-geoestadisticos

Poligonos de Thiessen: 1. Cambios abruptos de bordes. 2. Solo utiliza un punto para cada prediccion.

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 interpolacion.

En todos los casos, la distribucion no regular puede sobre-enfatizar areas (o lo contrario). La varianza debe ser estimada a partir de un conjunto de datos de validacion, especificamente reservado para ello.

Poligonos de Thiessen

thiessen = krige(zinc ~ 1, meuse.sf, meuse.grid.sf, nmax = 1)

## [inverse distance weighted interpolation]

plot(thiessen["var1.pred"], pch=15,main="Thiessen")

Inverso de la distancia (idw)

idw = idw(zinc ~ 1, meuse.sf, meuse.grid.sf)

## [inverse distance weighted interpolation]

plot(idw["var1.pred"], pch=15,main="IDW")

#Cambiando el numero de vecinos
# idw$vecino<- krige(zinc ~ 1, meuse.sf, meuse.grid.sf,nmax=6)

#Cambiando pesos
# idw$idp05 = idw(zinc ~ 1, meuse.sf, meuse.grid.sf, idp = 0.5)$var1.pred
# idw$idp5 = idw(zinc ~ 1, meuse.sf, meuse.grid.sf, idp = 5)$var1.pred

F-Validacion cruzada

thiessen.cv <- krige.cv(logZn ~ 1, meuse.sf, nmax = 1)
idw.cv<- krige.cv(logZn ~ 1, meuse.sf)
ok.cv <- krige.cv(logZn ~ 1, meuse.sf, model = va)

thiessen.cv[1:5,]

## Simple feature collection with 5 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 181025 ymin: 333330 xmax: 181307 ymax: 333611
## Projected CRS: Amersfoort / RD New
##   var1.pred var1.var observed    residual zscore fold              geometry
## 1  3.057286       NA 3.009451 -0.04783475     NA    1 POINT (181072 333611)
## 2  3.009451       NA 3.057286  0.04783475     NA    2 POINT (181025 333558)
## 3  3.009451       NA 2.806180 -0.20327092     NA    3 POINT (181165 333537)
## 4  2.806180       NA 2.409933 -0.39624685     NA    4 POINT (181298 333484)
## 5  2.448706       NA 2.429752 -0.01895404     NA    5 POINT (181307 333330)

idw.cv[1:5,]

## Simple feature collection with 5 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 181025 ymin: 333330 xmax: 181307 ymax: 333611
## Projected CRS: Amersfoort / RD New
##   var1.pred var1.var observed   residual zscore fold              geometry
## 1  2.830957       NA 3.009451  0.1784941     NA    1 POINT (181072 333611)
## 2  2.797944       NA 3.057286  0.2593420     NA    2 POINT (181025 333558)
## 3  2.687633       NA 2.806180  0.1185473     NA    3 POINT (181165 333537)
## 4  2.605734       NA 2.409933 -0.1958009     NA    4 POINT (181298 333484)
## 5  2.505657       NA 2.429752 -0.0759050     NA    5 POINT (181307 333330)

ok.cv[1:5,]

## Simple feature collection with 5 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 181025 ymin: 333330 xmax: 181307 ymax: 333611
## Projected CRS: Amersfoort / RD New
##   var1.pred   var1.var observed     residual      zscore fold
## 1  2.941046 0.03408724 3.009451  0.068404690  0.37050151    1
## 2  2.938452 0.03300608 3.057286  0.118833216  0.65409504    2
## 3  2.735836 0.03416564 2.806180  0.070343758  0.38056671    3
## 4  2.627050 0.04264164 2.409933 -0.217116957 -1.05142137    4
## 5  2.425835 0.03281227 2.429752  0.003917734  0.02162802    5
##                geometry
## 1 POINT (181072 333611)
## 2 POINT (181025 333558)
## 3 POINT (181165 333537)
## 4 POINT (181298 333484)
## 5 POINT (181307 333330)

ok.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)

## NULL

print(plot(var1.pred~observed,idw.cv, main="IDW"), split=c(2,1,3,1), more=TRUE)

## NULL

print(plot(var1.pred~obs,ok.cv, main="OK"), split=c(3,1,2,1), more=FALSE)

## NULL

F.2 Correlacion

Idealmente 1

cor(thiessen.cv$var1.pred,thiessen.cv$observed)

## [1] 0.685839

cor(idw.cv$var1.pred,idw.cv$observed)

## [1] 0.7640391

cor(ok.cv$var1.pred,ok.cv$obs)

## [1] 0.7677732

F.3 Media de residuos

Idealmente 0

mean(thiessen.cv$residual)

## [1] 0.003457601

mean(idw.cv$residual)

## [1] -0.005565866

mean(ok.cv$residual)

## [1] -0.0001470178

F.4 Desviaciones estandar de residuos

Desviaciones estandar del error de interpolacion (residuos)

Idealmente pequeno

sd(thiessen.cv$residual)

## [1] 0.2463506

sd(idw.cv$residual)

## [1] 0.2238086

sd(ok.cv$residual)

## [1] 0.1731182

F.5 MSPE (mean square predictor error)

Idealmente pequeno

mean(thiessen.cv$residual^2)

## [1] 0.06030903

mean(idw.cv$residual^2)

## [1] 0.0497981

mean(ok.cv$residual^2)

## [1] 0.02977656

F.6 Error medio cuadratico (RMSE)

Idealmente pequeno

sqrt(sum(thiessen.cv$residual^2)/length(thiessen.cv$residual))

## [1] 0.245579

sqrt(sum(idw.cv$residual^2)/length(idw.cv$residual))

## [1] 0.2231549

sqrt(sum(ok.cv$residual^2)/length(ok.cv$residual))

## [1] 0.1725589

F.7 Varianza de residuos

Idealmente cercanos a cero

var(thiessen.cv$residual, na.rm=T)

## [1] 0.06068862

var(idw.cv$residual, na.rm=T)

## [1] 0.05009028

var(ok.cv$residual, na.rm=T)

## [1] 0.02996989

F.8 Variacion explicada

Cantidad de variacion explicada por el modelo en cada metodo

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.99998

G-Anisotropia

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 anisotropia

vgm.aniso <- variogram(logZn ~ 1, meuse.sf, alpha = c(0, 45, 90, 135))
plot(vgm.aniso)

Autocorrelacion mas fuerte a los 45 y menos fuerte a los 135.

Mapa de variograma

vgm.map = variogram(logZn~1, meuse.sf, 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(st_bbox(meuse.sf)[c("xmax", "xmin")])^2 +
                 diff(st_bbox(meuse.sf)[c("ymax", "ymin")])^2)/4, psill=var(meuse.sf$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.4

plot(vgm.aniso, pl = T, model = va.a)

A partir de aqui se procede con kriging pero utilizando este variograma

Prediccion

ok <- krige(logZn ~ 1, locations = meuse.sf, newdata = meuse.grid.sf, model = va.a)

## [using ordinary kriging]

ok$pred <- 10^(ok$var1.pred)#volver a valores originales

plot(ok["var1.pred"], pch=15,
nbreaks = 64,
main="OK prediction [log10(Zn ppm)]")

plot(ok["var1.var"], pch=15,
nbreaks = 64,
pal = cm.colors,
main="OK prediction variance [log10(Zn ppm)^2]")

H-Exploracion - Bivariado

Analisis EDA - Exploracion de atributos - Bivariado

Analisis de correlacion (https://www0.gsb.columbia.edu/premba/analytical/s7/s7_5.cfm)

Distancia al rio

cor(meuse.sf$logZn, meuse.sf$dist)

## [1] -0.7394276

plot(meuse.sf$logZn ~ meuse.sf$dist)

Regresion lineal (https://www0.gsb.columbia.edu/premba/analytical/s7/s7_6.cfm)

m.lzn.dist <- lm(logZn ~ dist, data = meuse.sf)
summary(m.lzn.dist)

## 
## Call:
## lm(formula = logZn ~ dist, data = meuse.sf)
## 
## 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-16

plot(meuse.sf$logZn ~ meuse.sf$dist)
abline(m.lzn.dist)

I-Exploracion - Multivariado

Analisis EDA (Exploratory Data Analysis) - Exploracion de atributos - multivariado

Modelo aditivo, asume que las variables actuan lineal e independientemente

meuse.sf$logCu <- log10(meuse.sf$copper)
m.lzn.lcu.dist <- lm(logZn ~ logCu + dist, data = meuse.sf)
summary(m.lzn.lcu.dist)

## 
## Call:
## lm(formula = logZn ~ logCu + dist, data = meuse.sf)
## 
## 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-16

Modelo de interaccion

m.lzn.lcu.dist.i <- lm(logZn ~ logCu*dist, data = meuse.sf)
summary(m.lzn.lcu.dist.i)

## 
## Call:
## lm(formula = logZn ~ logCu * dist, data = meuse.sf)
## 
## 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-16

J-Kriging with external drift

“drift” being the value of the covariable

Prediccion de concentracion de metales de acuerdo con la distancia al rio.

vr <- variogram(logZn ~ dist, location=meuse.sf,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(st_bbox(meuse.sf)[c("xmax", "xmin")])^2 + diff(st_bbox(meuse.sf)[c("ymax", "ymin")])^2)/4, psill=var(meuse.sf$logZn)))
plot(vr, plot.numbers=T, model=vked, main="Dist  co-variable")

ked <- krige(logZn ~ dist, locations=meuse.sf,newdata=meuse.grid.sf, model=vked)

## [using universal kriging]

ked$pred <- 10^(ked$var1.pred)#volver a valores originales

plot(ked["var1.pred"], pch=15,
nbreaks = 64,
main="KED prediction [log10(Zn ppm)]")

plot(ked["var1.var"], pch=15,
nbreaks = 64,
pal = cm.colors,
main="KED prediction variance [log10(Zn ppm)^2]")

K-Universal Kriging

Prediccion de tendencia espacial geografica (x e y) ademas de estructura local.

meuse.sf$x<- st_coordinates(meuse.sf)[,"X"]
meuse.sf$y<- st_coordinates(meuse.sf)[,"Y"]

meuse.grid.sf$x<- st_coordinates(meuse.grid.sf)[,"X"]
meuse.grid.sf$y<- st_coordinates(meuse.grid.sf)[,"Y"]

vr <- variogram(logZn ~ dist + x + y, location=meuse.sf,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(st_bbox(meuse.sf)[c("xmax", "xmin")])^2 + diff(st_bbox(meuse.sf)[c("ymax", "ymin")])^2)/4, psill=var(meuse.sf$logZn)))
plot(vr, plot.numbers=T, model=vuk, main="Dist  co-variable")

uk <- krige(logZn ~ dist + x + y, locations=meuse.sf,newdata=meuse.grid.sf, model=vuk)

## [using universal kriging]

uk$pred <- 10^(uk$var1.pred)#volver a valores originales

plot(uk["var1.pred"], pch=15,
nbreaks = 64,
main="KU prediction [log10(Zn ppm)]")

plot(uk["var1.var"], pch=15,
nbreaks = 64,
pal = cm.colors,
main="KU prediction variance [log10(Zn ppm)^2]")

Visualizacion 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(as(ok, "Spatial"), "pred", asp=1, main="OK prediction",at=ramp, col.regions=bpy.colors(64))
p2 <- spplot(as(ked, "Spatial"), "pred", asp=1, main="KED-prediction",at=ramp, col.regions=bpy.colors(64))
p3 <- spplot(as(uk, "Spatial"), "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.08100

summary(ked$var1.var)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01978 0.02389 0.02584 0.02777 0.02980 0.05236

summary(uk$var1.var)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01993 0.02382 0.02568 0.02761 0.02948 0.05155

zmax <- 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(as(ok, "Spatial"), "var1.var", asp=1, main="OK varianza",at=ramp, col.regions=cm.colors(64))
p2 <- spplot(as(ked, "Spatial"), "var1.var", asp=1, main="KED-ffreq varianza",at=ramp, col.regions=cm.colors(64))
p3 <- spplot(as(uk, "Spatial"), "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 validacion cruzada el metodo de kriging que mejor predice la variable interpolada. Compruebe ademas incorporando anisotropia en el variograma.

Geoestadistica con R. Parte 2. Geoestadistica multivariada

Daniela Ballari

Tipos de Kriging

A-Cargar Librerias y Datos Meuse

B-Exploracion de atributos

C-Autocorrelacion

C.1 VARIOGRAMA EMPIRICO

C.2 VARIOGRAMA TEORICO

C.3 AJUSTE VISUAL

C.4 AJUSTE AUTOMATICO

D-Interpolacion

Kriging ordinario

E-Metodos no-geoestadisticos

Poligonos de Thiessen

Inverso de la distancia (idw)

F-Validacion cruzada

F.1 Scatterplots

F.2 Correlacion

F.3 Media de residuos

F.4 Desviaciones estandar de residuos

F.5 MSPE (mean square predictor error)

F.6 Error medio cuadratico (RMSE)

F.7 Varianza de residuos

F.8 Variacion explicada

G-Anisotropia

Variograma con anisotropia

Mapa de variograma

Prediccion

H-Exploracion - Bivariado

I-Exploracion - Multivariado

J-Kriging with external drift

K-Universal Kriging

Visualizacion cojunta

Mapas de varianzas

Actividad Individual