Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/

Daniela Ballari (PhD)
daniela.ballari@ucuenca.edu.ec
http://www.researchgate.net/profile/Daniela_Ballari/publications

Curso completo en: http://rpubs.com/daniballari/analisis_espacial_indice



Objetivo de la práctica. Introducir los conceptos básicos de geoestadística 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. February 27, 2015. http://www.css.cornell.edu/faculty/dgr2/teach/R/gs_short_ex.pdf

Una breve introducción a la geoestadística se encuentra aqui.

Temario

A- Cargar librerias
B- Cargar datos Meuse
C- Exploración de atributos
D- Estructura espacial local: autocorrelación
E- Interpolación con kriging
F- Otros métodos de interpolación no geoestadística

A-Cargar Librerías 

library(sp)
library(gstat)

B-Cargar datos Meuse

Utilizaremos el 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()

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

C-Exploración de atributos - 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]  375
sort(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] 1839
hist(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.265
head(meuse$zinc)
## [1] 1022 1141  640  257  269  281
head(meuse$logZn)
## [1] 3.009451 3.057286 2.806180 2.409933 2.429752 2.448706
#Luego para transformar a los valores originales se utiliza antilog. antilog (x)=10^x
head(10^(meuse$logZn)) 
## [1] 1022 1141  640  257  269  281

D-Estructura espacial local: autocorrelación 

#Crear Spatial Data Frame
coordinates(meuse) <- c("x", "y")
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(meuse) #Slot @proj4string no tiene asignado el sistema de referencia.
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  155 obs. of  13 variables:
##   .. ..$ cadmium: num [1:155] 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##   .. ..$ copper : num [1:155] 85 81 68 81 48 61 31 29 37 24 ...
##   .. ..$ lead   : num [1:155] 299 277 199 116 117 137 132 150 133 80 ...
##   .. ..$ zinc   : num [1:155] 1022 1141 640 257 269 ...
##   .. ..$ elev   : num [1:155] 7.91 6.98 7.8 7.66 7.48 ...
##   .. ..$ dist   : num [1:155] 0.00136 0.01222 0.10303 0.19009 0.27709 ...
##   .. ..$ om     : num [1:155] 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 [1:155] 50 30 150 270 380 470 240 120 240 420 ...
##   .. ..$ logZn  : num [1:155] 3.01 3.06 2.81 2.41 2.43 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ coords     : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 178605 329714 181390 333611
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
#Visualizar datos. Los puntos no están regularmente distribuidos, 
#sino que son mas densos en la cercanía del río.
plot(meuse, asp = 1, pch = 1) #asp=1 las dos escalas son iguales
#cargar otro conjunto de datos llamado meuse.riv que contiene lineas de las márgenes del río.
data(meuse.riv)
lines(meuse.riv)

#Otra visualización con símbolos proporcionales
plot(meuse, asp = 1, cex = 4 * meuse$zinc/max(meuse$zinc),pch = 1) 
#cex tamaño de circulos proporcionales al valor
lines(meuse.riv)

La dependencia espacial local significa que cuanto mas cerca estén los puntos en el espacio geográfico, mas cerca también lo están en el espacio de atributos. Esto se llama autocorrelación.

Los valores de un atributo pueden estar correlacionados con si mismos, y la fuerza de esta correlación depende de la distancia de separación entre puntos.

Cada par de puntos, estará separado por una distancia en el espacio geográfico y una semivarianza en el espacio de atributos.

#Calcule cuantos pares de puntos hay en el dataset meuse.
n <- length(meuse$logZn)
n * (n - 1)/2
## [1] 11935
#calcule la distancia y semivarianza entre los dos puntos primeros puntos del dataset
coordinates(meuse)[1, ]#punto 1
##      x      y 
## 181072 333611
coordinates(meuse)[2, ]#punto 2
##      x      y 
## 181025 333558
sep <- dist(coordinates(meuse)[1:2, ])
sep #distancia
##          1
## 2 70.83784
gamma <- 0.5 * (meuse$logZn[1] - meuse$logZn[2])^2 # semivarianza, unidades log(mg kg-1)^2
gamma
## [1] 0.001144082
#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.

#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) 
# 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. 
ve 
##     np       dist      gamma dir.hor dir.ver   id
## 1   41   72.24836 0.02649954       0       0 var1
## 2  212  142.88031 0.03242411       0       0 var1
## 3  320  227.32202 0.04818895       0       0 var1
## 4  371  315.85549 0.06543093       0       0 var1
## 5  423  406.44801 0.08025949       0       0 var1
## 6  458  496.09401 0.09509850       0       0 var1
## 7  455  586.78634 0.10656591       0       0 var1
## 8  466  677.39566 0.10333481       0       0 var1
## 9  503  764.55712 0.11461332       0       0 var1
## 10 480  856.69422 0.12924402       0       0 var1
## 11 468  944.02864 0.12290106       0       0 var1
## 12 460 1033.62277 0.12820318       0       0 var1
## 13 422 1125.63214 0.13206510       0       0 var1
## 14 408 1212.62350 0.11591294       0       0 var1
## 15 173 1280.65364 0.11719960       0       0 var1
plot(ve, plot.numbers = T, asp=1)

#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

# 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) 
vt
##   model psill range
## 1   Nug  0.01     0
## 2   Sph  0.12   850
plot(ve, pl = T, model = vt)

#Ajuste automático
#fit.variogram: ajuta el modelo de variograma a un variograma empírico.
va <- fit.variogram(ve, vt) 

va
##   model      psill    range
## 1   Nug 0.01004123   0.0000
## 2   Sph 0.11525698 967.2634
plot(ve, pl = T, model = va)

E-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.
data(meuse.grid) #malla de 40m x 40m, disponible con el dataset meuse.
coordinates(meuse.grid) <- c("x", "y")
gridded(meuse.grid) <- T #indica que el conjunto de datos es un raster

#Predicción
ok <- krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = va) 
## [using ordinary kriging]
#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.
ok$pred <- 10^(ok$var1.pred)#volver a valores originales

str(ok)
## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
##   ..@ data       :'data.frame':  3103 obs. of  3 variables:
##   .. ..$ var1.pred: num [1:3103] 2.83 2.88 2.83 2.78 2.94 ...
##   .. ..$ var1.var : num [1:3103] 0.0596 0.0471 0.0509 0.055 0.0335 ...
##   .. ..$ pred     : num [1:3103] 679 763 680 605 871 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 178460 329620
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cellsize         : Named num [1:2] 40 40
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cells.dim        : Named int [1:2] 78 104
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   ..@ grid.index : int [1:3103] 69 146 147 148 223 224 225 226 300 301 ...
##   ..@ coords     : num [1:3103, 1:2] 181180 181140 181180 181220 181100 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3103] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 178460 329620 181540 333740
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
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(20)),
             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)

F-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(log10(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(log10(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(log10(zinc) ~ 1, meuse, meuse.grid,nmax=6)
## [inverse distance weighted interpolation]
#Cambiando pesos
idw$idp05 = idw(log10(zinc) ~ 1, meuse, meuse.grid, idp = 0.5)$var1.pred
## [inverse distance weighted interpolation]
idw$idp5 = idw(log10(zinc) ~ 1, meuse, meuse.grid, idp = 5)$var1.pred
## [inverse distance weighted interpolation]
idw$idp10 = idw(log10(zinc) ~ 1, meuse, meuse.grid, idp = 10)$var1.pred
## [inverse distance weighted interpolation]
spplot(idw, c("var1.pred","idp05", "idp5", "idp10"),  col.regions=rev(heat.colors(50))
       ,main="IDW")

Guardar y recuperar archivos en memoria y espacio de trabajo

save.image(file = "geostat_b.RData")
load(file = "geostat_b.RData")

Ha llegado al final de este material!