El agua subterránea de California proporciona aproximadamente del 30 al 46 por ciento del suministro total de agua del estado, dependiendo de los años húmedos y secos. [CAg] A medida que aumenta la demanda de agua para usos domésticos, agrícolas e industriales, la gestión del agua subterránea se vuelve cada vez más importante. Diversas medidas de manejo necesitan conocer el comportamiento espacial y temporal del agua subterránea. Se necesitan mapas precisos de la profundidad del agua subterránea para predecir la dirección del flujo de la red y para monitorear la recarga del agua subterránea. En una región dispersa de observación de aguas subterráneas, como California, se pueden usar métodos geoestadísticos para determinar los valores de los puntos donde no se realizan mediciones. Las mediciones se tomaron en el otoño de 2016, entre octubre y diciembre.De las 4.727 mediciones originales, se mantuvieron 4.028 mediciones con ubicaciones únicas. Además, 42 ubicaciones con una profundidad negativa para la medición del agua subterránea se excluyeron del kriging ya que indican que el nivel del agua subterránea está por encima de la superficie del suelo en esas ubicaciones.

Profundidad a las aguas subterraneas y elevación de las aguas subterraneas.

La profundidad de las aguas subterráneas es la distancia vertical medida al agua en un pozo desde un punto de referencia definido. El DBGS es la profundidad a las aguas subterráneas menos la distancia desde el punto de referencia hasta la superficie del suelo. A menudo el punto de referencia es el borde superior de la carcasa del pozo, que está comúnmente por encima de la superficie del suelo. La profundidad de las aguas subterráneas y DBGS se expresa comúnmente como un número positivo. El DBGS puede ser un número negativo si el nivel de agua en una carcasa de pozo está por encima de la superficie del suelo.

La elevación de las aguas subterráneas o WSEL es la elevación de la superficie del suelo menos profundidad debajo de la superficie del suelo o DBGS. Los valores negativos de WSEL indican que los niveles de las aguas subterráneas están por debajo del nivel medio del mar.
  • Analisis Exploratorio

Mapa Profundidad de agua subterránea de California

Se procede a realizar el grafico de las mediciones de DGBS, en el estado de California para observar como se distribuyen:

Histograma de DGBS

Luego se realiza un histograma de log(DGBS) debido a que este podria aproximarse a una distribucion normal

  • Variogramas o semivariogramas

En geoestadistica, la correlacion espacial se modela ya sea por el variograma o la funcion de covarianza. Conociendo que el variograma es una herramienta que nos permite analizar el comportamiento espacial de una variable sobre un area definida obteniendo como resultado un variograma experimental que refleja la distancia maxima y la forma en que un punto tiene influencia sobre otro punto a diferentes distancias.

Se pretende hacer uso de esta herramienta para obtener un variograma experimental en los datos de profundidad de aguas subterraneas (DGBS), seguido de esto encontrar un modelo que se ajuste a este variograma y asi poder estimar un punto en especifico de profundidad con ayuda del metodo de kriging, el cual usa el modelo ajustado encontrado para definir un ponderador que sera aplicado a los puntos de profundidad de las aguas subterraneas (DGBS) y poder predecir puntos donde no se generan ubicaciones, es decir, puntos de profundidad.*

Grafica de variograma de DGBS

A continuación graficamos el variograma con respecto a la profundidad de superficie de aguas subterraneas:
  • Modelos de Variogramas
A continuacion se realizan diferentes modelos de variogramas para ajustar el variograma experimental:

Grafico de los diferentes modelos de variogramas

Considerando los diferentes modelos de variogramas para identificar cual de ellos se ajusta de mejor manera al variograma experimental en los puntos de profundidad se realizaron las sumas cuadraticas ponderadas para cada uno, obteniendose:
  • Variograma Exponencial
## [1] 21836.71
  • Variograma Gausiano
## [1] 68781.57
  • Variograma Esferico
## [1] 24817.44
  • Variograma Matern
## [1] 21836.71
  • Variograma Stein’s
## [1] 21974.03
  • Variograma Bessel
## [1] 26073.68
De esta manera se puede observar que el mejor modelo de variograma con menor suma cuadratica es el variograma Matern. 
#Informacion del variograma Matern
vMat
##   model     psill    range kappa
## 1   Nug 0.6357458   0.0000   0.0
## 2   Mat 0.7497980 176.3888   0.5
Para fines de evalucación de medidas de prediccion se procedio a dividir la informacion en datos de prueba que corresponden al 25% y datos experimetales que corresponden al 75% de los datos. 
#Separacion de data para ajuste y prediccion
trainIndices <- sample (1: length ( h20Edit ) , length ( h20Edit ) / 4 ,
                        replace = FALSE )
test <- h20Edit [ trainIndices , ] #25% datos
train <- h20Edit [ - trainIndices , ]  #75% datos
  • Kriging Ordinario
Para predecir los puntos de profundidad donde no hubo medicion se usara este metodo. Kriging ordinario 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 <- krige ( id = " logDist " , formula = log ( DGBS ) ~ 1,
              train , newdata = test , model = vMat )
## [using ordinary kriging]
Error Cuadratico Medio de OK: 
## [1] 0.4876096
Convertimos OK en un dataframe y creamos una nueva columna con las predicciones: 
UK_df <- as( OK , "data.frame" )
UK_df$Prediccion = exp(UK_df$X.logDist..pred)
colnames(UK_df)=c("Longitude","Latitude","LogPred","LogVar","Prediccion")
pander::pandoc.table(
head(UK_df))
## 
## -----------------------------------------------------------------
##   &nbsp;    Longitude   Latitude   LogPred   LogVar   Prediccion 
## ---------- ----------- ---------- --------- -------- ------------
##  **3834**    -121.3      36.72      4.275    0.7098     71.87    
## 
##  **1210**    -117.3      34.17      5.34     0.6866     208.6    
## 
##  **2840**    -122.5      41.77      3.026    0.7109     20.62    
## 
##  **4208**    -119.3      34.25      4.327    0.6798     75.68    
## 
##  **3631**    -122.1      39.77      3.993    0.6653     54.23    
## 
##  **374**     -119.1       34.1      4.109    0.6876     60.91    
## -----------------------------------------------------------------

Intervalo de prediccion para un punto usando OK

Escogemos de la tabla UK_df la variable “LogPred” un valor para calcular el intervalo de medida de predicción de un punto de profuncidad subterranea a un nivel de confianza del 95%: 
prediccionPunto= UK_df$LogPred[1]
prediccionPunto
## [1] 4.274919
exp(prediccionPunto)+c(-1,1)*1.96*sqrt(MSE_OK)
## [1] 70.50566 73.24296

Grafico de medidas de prediccion de DGBS con OK

  • Kriging Universal
Para la estimacion o predicción de medidas de profundidad de agua subterranea con Kriging Universal haremos un modelo de regresión, para este modelo se usaran las variables predictoras latitud y elevacion de las aguas subterraneas(WSEL), tambien mediante un grafico de dispersión podremos observar si existe una correlacion espacial entre las variables: 
##Modelo de regresion para estimar kriging universal
a <- data.frame(train)
mod1 <- lm(log(DGBS) ~ 1 + Latitude + WSEL, a)
summary(mod1)
## 
## Call:
## lm(formula = log(DGBS) ~ 1 + Latitude + WSEL, data = a)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3926 -0.7087  0.0931  0.7771  2.3316 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.148e+01  3.347e-01  34.293  < 2e-16 ***
## Latitude    -1.963e-01  8.963e-03 -21.901  < 2e-16 ***
## WSEL        -8.102e-05  1.316e-05  -6.155  8.5e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.039 on 2987 degrees of freedom
## Multiple R-squared:  0.169,  Adjusted R-squared:  0.1684 
## F-statistic: 303.7 on 2 and 2987 DF,  p-value: < 2.2e-16
pairs(a,col= 4)

Como podemos observar en el grafico de dispersión no hay una tendencia por parte de las variables, esto tambien implica a que las variables Latitud y WSEL no tengan una correlacion con la variable de estudio, provocando que el modelo no sea lo suficientemente bueno , es decir, las variables no explican de manera eficiente la profundidad del agua subterranea (DGBS) con un coeficiente de determinación muy bajo del 15%.

A pesar de ello, para fines de explicación de como funciona el kriging universal se toma este modelo para las predicciones espaciales. 
train_g_U <- gstat(id = "log_dist", formula = log(DGBS) ~1 + Latitude + WSEL, data = train)
vg_U <- variogram(train_g_U)
#plot(vg_U)
hchart(vg_U, "scatter", hcaes(x = dist, y = gamma))

Anosotropia o Variogramas direccionales

Tambien se procedio a realizar variogramas direccionales para observar el comportamiento en diferentes direcciones: 0,45,90 y 135 grados.

Variogramas direccionales:

Se puede deducir que el variograma con angulo de 135 grados, es aquel con mayor continuidad, es decir,la asimetria de los puntos de profundidad estarian definidas por un angulo de 135 grados,ya que esta tiende a tener un comportamiento mas parecido al variograma experimental.
  • Modelos de Variogramas
A continuacion se realizan diferentes modelos de variogramas para ajustar el variograma experimental:

  • Sumas cuadraticas ponderadas de los modelos de Variograma

  • Variograma Esferico

## [1] 67430.32
  • Variograma Matern
## [1] 14085.08
  • Variograma Stein’s
## [1] 14256.24
Aplicación de Kriging universal(UK): 
UK <- krige(id = "logDist", formula = log(DGBS) ~ 1 + Latitude + WSEL ,
            train , newdata = test , model = vMat2)
## [using universal kriging]
Error cuadratico medio del UK: 
MSE_UK <- mean((log(test$DGBS) - UK$logDist.pred )^2)
MSE_UK
## [1] 0.4740497
Se procede a calcular las predicciones: 
## 
## -----------------------------------------------------------------
##   &nbsp;    Longitude   Latitude   LogPred   LogVar   Prediccion 
## ---------- ----------- ---------- --------- -------- ------------
##  **3834**    -121.3      36.72      4.176     0.71      65.09    
## 
##  **1210**    -117.3      34.17      5.445    0.687      231.5    
## 
##  **2840**    -122.5      41.77      3.13     0.7111     22.86    
## 
##  **4208**    -119.3      34.25      4.338    0.6801     76.55    
## 
##  **3631**    -122.1      39.77      3.989    0.6658     53.99    
## 
##  **374**     -119.1       34.1      4.132    0.6878     62.32    
## -----------------------------------------------------------------

Intervalo de predicción de un punto usando UK

#intervalo de prediccion UK
prediccionPuntos=UK_df2$LogPred[1]
prediccionPuntos
## [1] 4.175822
exp(prediccionPuntos)+c(-1,1)*1.96*sqrt(MSE_UK)
## [1] 63.74384 66.44281

Grafico de medidas de prediccion de DGBS con UK

ggplot(california)+geom_sf()+
  geom_point(data =UK_df2,mapping = aes(x = Longitude, y = Latitude , color =Prediccion))

Referencias

Randolph, B. C. (2017). Extending kriging methods to large datasets with applications to California groundwater data. UCLA. ProQuest #ID: Randolph_ucla_0031N_16109. Merritt ID: ark:/13030/m53f9k41. Retrieved from https://escholarship.org/uc/item/23d1w7vv Obtencion de la base de datos : https://sgma.water.ca.gov/webgis/?appid=SGMADataViewer#gwlevels

California, G. d. (2013). Appendix E. California’s Groundwater. Obtenido de California Department WATER RESOURCES: #https://water.ca.gov/-/media/DWR-Website/Web-Pages/Programs/Groundwater-Management/Data-and-Tools/Files/Statewide-Reports/California-G#roundwater-Update-2013/California-Groundwater-Update-2013---Appendix-E.pdf

Vera, F. (8 de Julio de 2020). Youtube. Obtenido de Geoestadistica con R: https://www.youtube.com/watch?v=pHmbZAqU55I