Introducción

El presente taller de Geoestadística tiene por objetivo determinar por medios determinísticos y geoestadísticos, utilizando modelos de regresión lineal, inverso a la distancia y modelos kriging, la predicción de las temperaturas máximas de la región de la Araucanía. Para aquello, se extrajeron los datos de las estaciones meteorológicas, ubicadas en la región, distribuidas en temporadas correspondiente al año 2020, y como variable predictiva, se ocuparon imágenes Modis de resolución espacial de 1km en las mismas épocas de los datos de las estaciones y como covariables, debido a la influencia en las temperaturas, un ráster de distancia a la costa y un ráster de elevaciones propio de la región.

Zona de estudio

La zona de estudio corresponde a la región de la Araucanía, Chile, comprendida entre las regiones del Biobío y Los Ríos, limitando al oriente con Argentina y al poniente con el Océano Pacífico, y está centrada en las coordenadas 38° 43’ 50’’ S y 72° 17’ 53’’ W (Figura 1). Su extensión total es de 31.842 km² y su altitud media es de 580 m.s.n.m. En la región es posible encontrarse con varios tipos de clima en la transición de norte a sur, tales como mediterráneo y oceánico lluvioso, así como también templado cálido con estación seca corta, templado cálido lluvioso con influencia mediterránea, templado frío lluvioso con influencia mediterránea, y clima de hielo de altura. En el contexto hidrográfico, la región cuenta con la presencia de tres grandes ríos: el Imperial, el Toltén y el Biobío. También cuenta con una serie de lagos entre los que se mencionan el Villarrica, Budi y el lago Caburga. En esta región además se encuentran los volcanes Llaima y el Villarrica los cuales son los más activos de Sudamérica.

Figura 1: Zona de estudio. Estaciones (en azul) y límite de la región de la Araucanía (en rojo).

Materiales

Se utilizaron imágenes Modis del año 2020 de invierno, otoño, primavera y verano, obtenidas y procesadas con Google Earth Engine (https://code.earthengine.google.com/) correspondiente a la colección MODIS/006/MOD11A1, LST_Day_1km. En cuanto a datos de terreno se obtuvieron desde estaciones temperaturas máximas obtenidas desde página https://agrometeorologia.cl/. La capa vectorial que representa el límite regional fue extraída de cartografía proveniente de http://www.ide.cl/.

Metodología y Resultados

Exploración de datos

Para la exploración de los datos, estos se reordenaron estableciendo sólo una columna de temperaturas y otra de periodos, de este modo se utiliza el paquete tidyr. Para realizar lo anterior primero se debe extraerr la parte gráfica de los puntos que representan las estaciones (columna geometry). Luego se aplica la función pivot. Lo siguiente es explorar los datos de las imágenes de MODIS, para ello se cargan desde una carpeta en común y se genera un Stack al que se le asignará un sistema de referencia en coordenadas geográficas WGS84 al igual que los datos vectoriales. Se debe extraer los datos a partir del Stack creado antes solamente en el lugar en donde se ubica cada estación. A los resultados de la extracción se les debe aplicar un factor de escala y luego convertirlos a grados Celsius. Luego de la conversión se acoplan los datos de dicha extracción en un solo SF al que se denominará “data_unida”. De estos datos se toman sólo los valores de temperatura, las estaciones y los periodos. Para esto nuevamente se omite la componente gráfica, luego se aplica la función pivot para agrupar tanto los períodos dados por las estaciones en terreno como los valores obtenidos a partir de la extracción desde las imágenes MODIS, el resultado se aprecia en la Tabla 1.

aws verano otono invierno primavera modis1 modis2 modis3 modis4
CARILLANCA 24.9 14.6 10.7 19.2 30.94667 13.715714 8.483333 22.70455
CABALLERIA 24.4 14.6 10.6 19.4 29.68368 13.270000 10.715000 22.40400
CUARTA_FAJA 24.7 14.6 10.7 19.5 28.58444 12.905556 7.890000 23.91333
DOMINGUEZ 18.9 14.5 10.9 14.8 23.97455 11.735714 7.456667 21.51500
EL_MEMBRILLO 23.6 14.8 9.5 19.0 26.62400 11.206667 5.980000 22.60333
EL_VERGEL 27.7 16.4 12.0 23.5 30.66545 13.980000 8.050000 26.24667
FAJA_MAISAN 22.7 15.1 10.9 18.0 20.91636 7.815714 3.180000 17.40750
GABY_RANQUILCO 27.7 16.1 10.9 22.1 35.10905 13.576667 8.650000 28.99133
HUISCAPI 23.7 13.9 9.7 18.3 27.75632 11.486364 7.010000 22.71571
LA_ISLA 27.3 16.6 10.2 21.6 33.05600 15.142500 10.230000 27.47353
LA_PROVIDENCIA 25.4 14.5 10.1 20.2 32.21211 13.333333 7.043333 23.32833
LAS_PALMAS 25.7 14.1 9.6 19.8 29.73200 14.550000 6.972857 21.65000
LOS_ARRAYANES 20.4 15.2 11.5 16.3 25.22000 12.388333 8.278571 20.78467
MANZANARES 28.2 16.7 12.0 23.4 29.34583 15.124286 7.525556 27.89526
MAQUEHUE 24.1 14.8 10.8 18.5 30.04529 13.666667 5.650000 23.49000
MARIMENUCO 24.2 13.4 4.3 20.7 25.68043 7.150000 -5.762000 17.88143
PUALA 24.9 14.9 8.7 20.6 25.18000 8.798889 3.420000 21.68882
QUIRIPIO 18.1 13.1 9.1 13.3 25.47353 11.776667 6.644000 19.65923
RADAL 24.8 14.5 10.6 19.2 27.96778 11.478000 6.663333 21.98333
RUCAMANQUE 23.4 13.2 9.2 17.9 24.86467 10.215454 3.455714 19.03200
SAN_ENRIQUE 23.9 15.1 10.1 19.4 22.73316 9.408182 4.234000 19.10412
SAN_LUIS 23.7 14.1 9.6 18.4 24.07211 10.818333 5.556000 20.96882
SAN_RAFAEL 27.0 15.5 10.7 22.0 33.28048 13.906000 5.036667 26.91526
SAN_SEBASTIAN 24.5 13.6 9.2 18.9 30.36400 13.864000 7.223333 23.18500
SANTA_ADELA 25.2 15.7 11.2 19.3 30.46294 12.507778 8.498000 24.12400
SANTA_INES 24.8 14.3 10.2 19.3 28.56000 11.900000 7.184286 22.77364
SURCO_Y_SEMILLA 26.2 14.5 10.2 20.8 29.00000 11.737692 7.078889 21.44000
TAPLON 24.4 14.5 10.9 19.3 29.37316 13.195000 6.660000 23.71222
TRANAPUENTE 19.4 15.6 11.8 15.5 25.79727 12.257143 7.931667 22.29727

Tabla 1: temperaturas registradas por estacion e imágenes MODIS.; modis1 (verano), modis2 (otoño), modis3 (invierno), modis4(primavera)

En la Figura 2 se aprecian los histogramas para cada periodo en el que podemos ver como se agrupan los valores de temperatura. En ellos, por ejemplo, prima una temperatura cuya moda es de 10°C aproximadamente en el invierno, subiendo a 15°C en otoño, en primavera hay un achatamiento en el histograma debido a que las temperaturas registradas rondan los 20°C y finalmente en verano en el cual hay nuevamente un alza en la temperatura alcanzando un promedio en los 25°C aproximadamente. De los resultados anteriores se generan gráficos para visualizar la distribución y relación de los datos (Figura 2).

Figura 2: histogramas de temperaturas en función de estaciones del año

En la Figura 3, se puede observar que en invierno se ve una buena concentración y baja variabilidad de temperaturas tanto para las registradas por las estaciones en tierra como para las registradas en MODIS. En cambio en verano, si bien hay una concentración de temperaturas en las estaciones de tierra de alrededor de los 24°C, los datos captados por MODIS son muy variables. Esto se corrobora en la Figura 4 en donde la curva que representa la distribución de los datos de verano es una curva más achatada.

Figura 3: relación entre las temperaturas registradas por MODIS y estaciones en tierra

Figura 4: histograma suavizado para los datos de temperatura de las estaciones.

Si comparamos los histogramas suavizados de la figura 4 (correspondientes a estaciones) con los de la figura 5 (correspondientes a los obtenidos por MODIS) podemos ver una discrepancia en las curvas para invierno, otoño y verano en lo que respecta a los valores de densidad, aunque visualmente son parecidos, los valores de MODIS son más bajos.

Figura 5: histograma suavizado para los datos de temperatura de MODIS.

Figura 6: datos MODIS para cada periodo.Arriba izq, verano; Arriba der, otoño; Abajo izq, invierno; Abajo der, primavera; temperaturas en grados Kelvin

Figura 7: mapas de temperaturas captadas por estaciones para los 4 periodos.

De igual manera se aplica un código similar pero para los datos obtenidos con MODIS, de esta manera se obtiene la vista de la Figura 8:

Figura 8: mapas de temperaturas captadas por MODIS para los 4 periodos.

Creacion de covariables raster de elevacion y distancias

Figura 9: interactivo de raster de elevación y distancias.

Análisis mediante Modelo de Regresión Lineal

Para realizar esto se utilizó un data frame (data_unida_est_modis) con los valores de Temperatura Máxima (Tmax) en cada una de las estaciones y los valores de Temperatura Superficial de Suelo (LST) extraído de los archivos modis en las coordenadas de las estaciones.

Archivo rasterstack de modis con las LST en los meses de enero, mayo, julio y noviembre.

Luego se trabajó de manera individual los datos de cada temporada. En cada periodo los pasos realizados fueron los siguientes:

  • Extraer del data frame los datos correspondientes al mes de interés.

  • Analizar los histogramas y gráficos de densidad para ver si la distribución de los datos presenta algún tipo de sesgo. Se debe trabajar con valores que presentan datos normales.

Adicionalmente se analizó la normalidad de los datos con las funciones qqnorm y shapiro test. En el caso de qqnorm debe mostrar una línea recta perfecta para que la distribución sea normal, mientras que con el shapiro test se analiza el p-value. Si el p-value es menor que 0.05 se rechaza la hipótesis de normalidad, por lo que los datos no presentarían una distribución normal y se deben ajustar los valores. En el caso que el valor de p-value sea mayor que 0.05 se acepta la hipótesis de normalidad, por lo que los datos presentan una distribución normal. En este caso para los meses de verano invierno se rechazó la hipótesis de normalidad, poseen un pequeño sesgo negativo, por lo que los datos se ajustaron la ecuación [1]. En las estaciones de otoño y primavera los valores de p-value son mayores a 0.05 por lo que los datos presentan una distribución normal.

                                √(max⁡(Tmax+1)-Tmax)                                               [1]

Donde, Tmax corresponde al valor de temperatura máxima del periodo.

http://www.datanovia.com/en/lessons/transform-data-to-normal-distribution-in-r/#non-normal-distributions

Luego se realizó el modelo de regresión para cada estación. En el caso de verano e invierno se utilizaron los datos ajustados. Los valores obtenidos para cada temporada se especifican en la Tabla 2.

tabla2

tabla2

Se procedió a calcular los valores de ajuste del modelo y los valores de los residuos. Una vez determinados estos valores de plotearon los gráficos de ajuste del modelo de regresión para cada periodo (Figura 10) y los gráficos de la calidad del modelo (Figura 11) en los cuales se observa como varían los residuos.

Es importante destacar, que en los gráficos de ajuste del modelo (Figura 10) se observa que los P-values son menores a 0.05, lo que implica que los coeficientes están bien estimados. En cambio, los modelos no superan los R2 a 50%, al ser bajos no explican la variabilidad de los datos, solo verano e invierno, están cercanos a esta cifra y se pueden ver la gráfica como los puntos están cercanos a la zona gris (grado de confiabilidad).

Figura 10: Gráfica de modelos de regresión lineal por periodos

Para analizar la calidad del modelo se gráficaron los residuos vs la temperatura superficial de suelo (LST). En estos gráficos se presenta una distribución heterogénea de los datos, no hay una dispersión sistemática. En invierno, se observa una agrupación de estaciones con temperaturas similares con residuales menores a 1, en cambio en otoño y primavera, la varianza es más amplia que el resto.

Figura 11. Gráficos de residuos vs temperatura superficial de suelo (LST) para cada periodo

Predicciones espaciales por Modelo de Regresion Lineal

Luego se realizó la predicción espacial utilizando el modelo de regresión calculado sobre el archivo raster de cada mes. En el caso de verano e invierno, una vez generada la predicción, se transformaron los valores para que las temperaturas se encontraran en el mismo rango de temperatura que las demás estaciones. Los mapas de predicción espacial se pueden observar en la Figura 12.

Figura 12. interactivo de predicción, modelo regresión lineal

Análisis de autocorrelación espacial

Índice Moran - Montecarlo

El índice de moran es una medida estadística que permite analizar la autocorrelación espacial de los datos, que por lo general se mueve entre los valores, -1 y 1. Las zonas con valores mayores que 0 indican una autocorrelación positiva, es decir que los valores similares están cercanos entre sí. Valores bajo 0 indican una autocorrelación negativo, lo que implica que las zonas tienen distintos valores que sus vecinos. Por ultimo, las zonas con valores 0 indican que no hay autocorrelación espacial. Es importante destacar que este índice va acompañado de una prueba de hipótesis (test de monte carlo), bajo el supuesto de normalidad. La hipótesis nula establece que no hay autocorrelación espacial. Esto se traduce en que si los valores obtenidos en el p-value son menores a 0.05 se rechaza la hipótesis de normalidad y por ende hay autocorrelación espacial de los datos. En este caso los valores obtenidos por el índice de moran (Figura 13) son positivos, pero muy bajos, lo que indica que presentan una autocorrelación pero no de manera significativa. Por otra parte, el test de monte carlo entrega como resultado que los p-value se encuentran todos bajo 0.05 lo que muestra que los datos presentan autocorrelación espacial.

Figura 13. Índice Moran - Montecarlo Arriba izq, verano; Arriba der, otoño; Abajo izq, invierno; Abajo der, primavera

Variogramas

El Variograma es una herramienta que permite analizar el comportamiento espacial de una variable sobre un dominio definido. Lo primero es determinar el variograma experimental para luego calcular el variograma ajustado utilizando diversos modelos de variogramas: esférico, circular, exponencial, gaussiano, etc. Con el objeto de analizar la autocorrelación espacial de las variables se realizaron variogramas experimentales y ajustados, lo que permitirá determinar como el parámetro de estudio va variando a diferencias distancias. En este caso en la Figura 14, se muestran variogramas con el modelo ajustado. Es importante destacar que es posible ver una autocorrelación espacial de los datos y a pesar de que hay algunos puntos que se alejan de la curva de ajuste, se lográn ajustar de forma general infiriendo la autocorrelación. Se sugiere seguir iterando el variograma para encontrar un mejor modelo de ajuste.

Figura 14. Variogramas con los modelos ajustados para los distintos periodos.

Variogramas Anisotrópicos

Los variogramas anisotrópicos permiten analizar la autocorrelación espacial en diferentes direcciones. En el gráfico que se presente la menor semi-varianza, los valores se parecerán más y por ende, habrá mayor autocorrelación a menor distancia. En este caso se realizaron variográmas anisotrópicos con dirección Norte (0º) y dirección E (90º) a los cuales no se le realizron ajuste por la baja cantidad de muestras. En el caso, en la Figura 15, se puede observar que en verano, otoño e invierno hay una leve diferencia entre los variogramas generados, siendo la de dirección norte (0º) la que presenta una menor semi-varianza. En el caso de primavera, la menor semi-varianza se observa en dirección E (90º). Es importante remarcar, que los apreciaciones realizadas son por medio de inspección visual, lo que puede generar un error de interpretación. Adicionalmente, se destaca que entre los cuatro periodos analizados, la de menor semi-varianza corresponde a invierno, lo que permite interpretar que es en este periodo en el que se presenta una mayor autocorrelación a una menor distancia.

No se consideraron encontrar modelos de ajuste por la escasez de muestras en la región,

Figura 15. Variogramas anisotrópicos de los periodos

Modelo Universal de Variación

El modelo universal de variación se puede dividir en una parte determinística y una parte de variación espacial. La parte determinística puede determinarse a partir de los modelos de regresión. Adicionalmente, los métodos de predicción espacial geoestadísticos, tales como Ordinary Kriging (OK) y Regression Kriging (RK) son derivados del modelo universal de variación. A continuación, se procederá a generar mapas de predicción espacial utilizando MRM, OK y RK.

Modelos de Regresión Multiple (MRM)

La regresión lineal múltiple trata de ajustar modelos lineales entre una variable dependiente y más de una variables independientes. En la Figura 16, se observa los mapas de predicción espacial realizado mediante Modelo de Regersión Limeal Multiple (MRM)

Figura 16. Mapas de predicción, modelos de regresión multiple por periodo

Modelos de Ordinary Kriging e Inversos a la distancia

Al efectuar la interpolación de los ordinario kriging, las validaciones cruzadas, nos dan r2 negativos,lo que se explica por la distribución de las muestras, que dentro de los dominios de estos, la varianza están cercanas a 0, pero hacia la periferia de la región alcanzan hasta el valor 4.

En relación al modelo de inverso a la distancia, es donde matematicamente se asume que los valores mas cercanos están mas relacionados que otro y determinados por su función.

A continuación se muestran los mapas obtenidos para Ordinary Kriging e Inverso a la distancia para cada mes

Ordinary Kriging - Inverso a la Distancia, Verano

Figura 17. Mapas de predicción ordinary kriging, inverso a la distancia en verano

Prueba corte de área con Ordinary Kriging

Para observar como se comportaba el modelo ordinary kriging, dentro de los dominios de las muestras, se efectuó un recorte de la región , generando un nuevo polígono de corte. Si bien se redujo la estimación de las temperaturas, los r2 aun se mantienen con valores negativos, ya que las muestras no cubren la zona por su distribución y distancia, originando extrapolaciones hacia la periferia de la región. Figura 19

Figura 18. Mapa de corte de dominio de las estaciones, ordinary kriging, verano.

Ordinary Kriging - Inverso a la Distancia, Otoño

Figura 19. Mapas de predicción ordinary kriging, inverso a la distancia en Otoño

Ordinary Kriging - Inverso a la Distancia, Invierno

Figura 20. Mapas de predicción ordinary kriging, inverso a la distancia en invierno

Ordinary Kriging - Inverso a la Distancia, Primavera

Figura 21. Mapas de predicción ordinary kriging, inverso a la distancia en primavera

Variograma de Ajuste de Residuos

Esta Herramienta permite ver como varían los residuos a medida que aumenta la distancia. Estos gráficos son necesarios para determinar si es posible aplicar el método geo-estadistico de Regression Kriging.

En este caso en particular, realizando una inspección visual se pueden observar en los variogramas de ajustes de residuos (Figura22, Figura23, Figura 24 y Figura 25), que en el mes de verano los puntos se encuentran mas cercanos a la curva de ajuste, mientras que en el mes de otoño los puntos se encuentran mas distantes a la curva de ajuste, lo que permite inferir que el modelo de verano es mejor que el de otoño. Esto se encuentra alineado con los R2 obtenidos mediante validación cruzada para el metodo de Regressión Kriging, los cuales pueden ser observados en la Figura. Estos R2 muestran que para verano se obtiene un valor de 0.81 y para el mes de otoño se obtuvo un R2 de 0,25 lo que permite corroborar que el modelo de verano es mejor que el de otoño.

Figura 22. Variograma de ajuste de los residuos correspondiente al mes de Verano

Figura 23. Variograma de ajuste de los residuos correspondiente al mes de Otoño

Figura 24. Variograma de ajuste de los residuos correspondiente al mes de Invierno

Figura 25. Variograma de ajuste de los residuos correspondiente al mes de Primavera

Mapas comparativos de predicciones

En esta sección muestran las distintas predicciones espaciales obtenidas para cada periodo. Los modelos utilizados son: MRM, RK y IDW.

En los mapas se puede observar que los rangos de temperatura son similares en los tres modelos. Partiucularmente, se puede apreciar en el modelo de inverso a la distancia los rangos de valores tienden a ser mas acotados.

Los mapas predictores de regresión kriging, que consideran los modelos ajustados de los residuos y la regresión múltiple, mejoran bastante explicando la variabilidad similarmente a los modelos de regresión múltiples, aunque estos mantienen en unos porcentajes mas altos que los de regresión kriging.

Se puede observar en las comparaciones de los modelos; en verano las temperaturas entre MRM y RK, son idénticas en grados de temperatura y forma de la interpolación, en cambio en el modelo IDW la distribución de las muestras va influir bastante en las formas de la interpolación que serán certeras mientras estén cercanas a la muestra. Aunque en las temperaturas más altas de MRM y RK como se aprecian en la figura 26 se escapan (entre 30°y 35°), y están cercanas a la cordillera, donde no hay toma de muestras, falseando la información. No así en IDW.

En otoño, no se logra explicar la variabilidad, considerando los r2 bajos de los modelos MRM y RK, entre ellos la diferencia es incluso distinta en el centro de la región, incluso hacia la cordillera la forma influye la covariable de elevación, por la forma que tiene la interpolación. En este periodo, el modelo IDW, se podría afirmar que posee más peso.

En invierno, el r2 del modelo RK es menor a MRM, que explica mejor la variabilidad en un 72 % e inclusive IDW tiene rangos de temperaturas similares a los dos modelos.

En primavera, el r2 de 0.84 del modelo de MRM, explica mejor el Modelo de RK, que origina un 72 % de explicación, aunque como ocurre en los distintos periodos genera también una extrapolación hacia la cordillera, y esto es por la falta de puntos de muestras en el sector. En cambio, donde hay mayor concentración de estaciones los valores de temperaturas son similares incluso con el modelo de IDW.

MRM - Regression Krigging - Inverso a la Distancia

Figura 26. Mapas comparativos de MRM, RK e Inverso a la distancia

Mapa Interactivo

Modelo de Regresión Multiple (MRM) - Regresion Kriging

Figura 27. Interactivo de MRM y Regresion kriging

Conclusiones

Los modelos de regresión simple en base al predictor LST poseen bajos r2, exceptuando verano e invierno cercanos a 50 %, y que mejoran considerablemente, cuando se le suman las covariables de distancia a la costa y elevación. Incluso en otoño que de un 18 % logra explicar la variabilidad en un 27 %, aun bajo si se considera ideal sobre un 50 %.

Los índices de Morán nos demuestran a nivel global una autocorrelación espacial, que nos lleva a generar modelos geoestadísticos, como Ordinario y regresión kriging, pero previa búsqueda de modelos de ajuste por medio de variogramas.

Por la cantidad y distribución de las estaciones en la región, se hace la tarea difícil encontrar modelos apropiados. Es como se observa en el periodo que corresponde a otoño, una distribución bien dispersa de la agrupación de los datos que conlleva a unos mapas predictivos bajos en r2.

En síntesis, la falta de más muestras, sobre todo en los extremos de la región, origina un sesgo en la predicción, que se debe considerar en un análisis posterior. Aunque con los puntos de muestras que se obtuvieron en algunos casos los modelos de regresión kriging se comportaron bien. En el caso del estudio, con los modelos de regresión múltiples se lograron mejores resultados, no así en otoño, donde se podría afirmar que la variable predictora LST haya influido en ese periodo.