Interpolación espacial con datos de precipitación CHIRPS en el amazonas

1. Introducción

El análisis espacial es una técnica ampliamente usada en topografía y geografía la cual consiste en la recopilación de información en un punto de muestreo específico, esto permitiendo la modelación o interpolación de datos puntuales para volverlos datos continuos, adicionalmente, para este proceso se suele usar el Sistema de Información Geográfica SIG el cual proporciona herramientas de análisis espacial para calcular estadísticas en el geoprocesamiento de datos por medio de celdas en un raster. La interpolación espacial se refiere principalmente al uso de puntos con valores cuantitativos dados, con los cuales se estima el valor de otros datos que son desconocidos. Es decir, es el paso de mediciones puntuales a representaciones continuas de variables que expresan propiedades continuas de la naturaleza. Usualmente, es usada para mediciones de la precipitación, el pH, humedad de suelo, concentración química en el suelo, etc. Este método también abarca varios tipos de interpolación como son:

• Distancia Inversa Ponderada (IDW)
• Red Irregular Triangulada (TIN)
• Kriging
• Polígonos Thiessen
• Natural Neighbourhood Interpolation (NN)

La interpolación espacial también tiene como propósito ser una opción viable para los objetos que están distribuidos en el espacio y que tienen relación mutua.¹ Es decir, que dependiendo de los puntos de los datos muestra se puede determinar el valor de puntos cercanos a este. En otras palabras, lo que se realiza es la predicción de valores para celdas en un raster a partir de valores de puntos muestra. Aunque este método no ha sido considerado 100% confiable, ya que, la información de los puntos muestra, es decir, los datos originales, pueden tender a degradarse al brindar datos a otras celdas por el nivel de probabilidad que este puede llegar a alzanzar, de igual forma, esto también depende de la distancia a la que se realiza esa predicción de los valores para saber una aproximación con una alta confianza del dato.² Una de las principales razones por las cuales el método de interpolación espacial ha sido tan usado y práctico es por los niveles de dificultad que se pueden presentar al realizar estudios de variables como el pH del suelo, ya que, se necesitaría de mucho dinero y también de personas para realizar un estudio detallado de, por ejemplo, toda la superficie del suelo de un país. En cambio, con el uso de la interpolación, no es necesario tener tantos estudios en una zona para hacer una determinación aproximada de los valores que pueden tener puntos sin presencia de datos.Por último, pero no menos importante, para la generación de superficies raster existen dos métodos generales: el uso de técnicas deterministas y técnicas estadísticas (geoestadísticas). Las técnicas deterministas se basan en el uso de una función matemática para poder conocer valores desconocidos y adicionalmente dar el valor de estas características. Por otro lado, las técnicas estadísticas o geoestadísticas, producen límites de confianza para la precisión de una predicción, aunque son un poco más difíciles de ejecutar es necesario para poder obtener un mayor numero de parámetros.³

En este informe, se muestran tres métodos de interpolación espacial, usados en el departamento del Amazonas en el país de Colombia con el uso de datos de precipitación CHIRPS con el objetivo de brindar la información necesaria para la comprensión de los diferentes conceptos e información que abarca el análisis espacial y uno de sus usos. Fueron elegidos los datos de precipitación ya que gráficamente proporcionan información útil a diferencia de la humedad y elevación del suelo, ya que, el amazonas es un departamento con una alta y uniforme humedad y la elevación del terreno es en general muy baja y sin muchas variaciones, a diferencia de la precipitación ya que presenta variaciones en diferentes zonas del departamento y en distintas épocas.

2. Descripción de la zona de estudio

El departamento del amazonas fue fundado en el año 1991 y cuenta con 2 municipios, Puerto Nariño y Leticia (su capital), y además, cuenta con 9 corregimientos. Este departamento se ha caracterizado por su multivariedad de culturas tanto en la época de la conquista como en las que se mantienen hoy en día. De igual modo, con el paso del tiempo el clima fue transformándose de sabana a uno mucho más humedo formando poco a poco la gran selva que hoy que existe actualmente. Por otro lado, las practicas que empleaban los humanos en el transcurso del tiempo fue evolucionando hasta la pesca y la agricultura desarrollando la cultura de la selva tropical. El amazonas es ampliamente reconocido por ser la zona forestal más grande del mundo que igualmente es compartida con los países de Venezuela, Brasil, Ecuador, Perú, Guyana, Suriname y Bolivia por esta razón es conocido por el nombre alterno del pulmón del mundo. El departamento del amazonas en Colombia es considerada la zona más forestal y, de igual modo, este terreno comprende las selvas tropicales de la llanura amazónica de relieve plano y la saliente del Vaupés, que posee colinas de poca altura y muy antiguas. El clima del departamento usualmente es cálido y húmedo con una temperatura media constante de 25.4°C y la precipitación varía entre los 2660 mm/año. El amazonas también es conocido ya que posee parte de los principales ríos de Colombia como Caquetá, Putumayo, Guaviare, Amazonas, Vichada, Apaporis, Inírida y Vaupés. Por la gran humedad que se presenta en este lugar no se permite la construcción de vías terrestres e igualmente por el temor al impacto ambiental. El departamento del amazonas es el más grande de Colombia ocupando 109665km2 de extensión territorial del país, además, su población cuenta con más de 76000 habitantes por estudios tomados en el 2018, y cuenta con países vecinos como Brasil y Perú brindando una gran riqueza cultural. Adicionalmente, también limita con los departamentos de Caquetá, y Vaupés. El amazonas está situado a 00° 07’ 08” de latitud norte, a 04° 13’ 19” latitud sur y a 74° 23’ 21” de longitud oeste y comprende una superficie plana desde los 80 msnm en cercanías de Leticia, hasta los 320 msnm hacia la zona norte cerca al Caquetá. ⁴ Según el IDEAM, los datos registrados de las estaciones meteorológicas, la Amazonía colombiana presenta precipitaciones en un régimen monomodal indicando los valores más bajos en los meses de diciembre y enero y los máximos en mayo y junio, mostrando una diferencia en los datos dependiendo de la zona sur o norte de departamento. ⁵

3. Descripción de datos y métodos

3.1. Datos

En este informe fueron empleados los datos Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) los cuales se basan en datos de precipitación infrarojo combinados con datos de estaciones, elaborados por un grupo llamado Amenazas climáticas. Este grupo de la Universidad de California, Santa Bárbara tienen un grupo de datos de precipitacion de todo el mundo de más de 35 años, es decir, se pueden adquirir datos desde el año 1981 hasta el año en curso 2020. La gran mayoría de los datos son adquiridos a una resolución de 0.05°, sabiendo que por cada grado hay una distancia de 111.11 km aproximadamente, entonces, 0.05° corresponderian a 5.55 km de distancia. Adicionalmente, a estas imágenes satelitales también se incorporan datos de estaciones climatológicas de todo el mundo gracias a convenios con la Universidad. También, debido al uso de un algoritmo y la union de estos datos se crea un raster de precipitacion, con el motivo de analizar sequías por medio de tendencias y monitoreos estacionales. Estos datos pueden ser adquiridos en diferentes periodos, como anuales, mensuales, cada cinco dias y diarios. Al descargar los datos hay tres opciones de formato raster: .bils, .netcdf y .tif. Los datos descargados en este informe fueron en formato .tif gracias a que se tiene más familiaridad con este tipo de formatos. ⁶ En este informe fueron usados los datos de precipitación de los últimos cinco días del mes de abril de este año, que fueron publicados el 18 de mayo. Al cargar los datos se obtuvo información del archivo como la siguiente:

class      : RasterLayer 
dimensions : 500, 1800, 9e+05  (nrow, ncol, ncell)
resolution : 0.2, 0.2  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : chirps.v2.0.2020.04.6 
values     : -9999, 401.1604  (min, max)

Esto indica que la clase del archivo cargado es raster, y que tiene 9e+05 celdas, es decir que el número de celdas es de 900000, de igual forma, las filas son en total 500 y las columnas son 1800 y los datos se encuentran en coordenadas geográficas. En este caso, las celdas tienen una resolución de 0.2°, un poco menos de lo mencionado anteriormente (0.05°) debido a que fue realizado un proceso adicional en el ploteo de los mapas, ya que, el número de celdas que habia era muy alto, entonces, fue usada la función aggregate del paquete raster para realizar el cambio de tamaño. A continuación, se muestra el mapa del amazonas con las celdas de los datos de precipitación graficados en un mapa interactivo, pero con la resolución original de 0.05°:

En la obtención de los puntos para la interpolación en este informe, se pasaron los datos raster a puntos, lo cual permitió que se obtuvieran los puntos del centro de cada una de las celdas que corresponden a la precipitación del departamento del amazonas.En el siguiente mapa se puede observar la representación gráfica de los puntos en cada celda representados en color coral,

De igual forma, otro representación más de la precipitación en esta zona puede ser el siguiente mapa ya que se pueden ver directamente los datos de precipitación en mm en cada punto, asimismo, entre más azul es el punto representado, es mayor el nivel de precipitación:

3.2. Métodos

En este informe fueron usados tres métodos de interpolación, los cuales fueron IDW, polígonos de Thiessen y por último el método de Kriging.

• IDW

Inverse Distance Weighted Interpolation o según sus siglas en inglés IDW, se refiere a la ponderación basada en el inverso de la distancia. Este método deterministico define una región cerca a la zona de interés donde se quieren tener valores interpolados por medio de un radio de busqueda, si lo datos están mucho más alejados de ese radio de busqueda se considera que la influencia de estos es menor, es decir, que la confiabilidad disminuye en relación con el aumento de la distancia. Para hallar la estimación de lo valores a interpolar se realizan promedios de los valores muestra en la las celdas vecinas a las celdas de procesamiento, esto, teniendo en cuenta la distancia del punto al centro de la celda. Para este tipo de interpolación es usado un algoritmo que es aplicado a datos que tienen una alta variabilidad, aunque, la superficie interpolada debe ser menor que el valor máximo local de la zona y mayor que el mínimo local. En otras palabras, este método se basa en el inverso de la distancia, y esta distancia es elevada a una potencia matemática; esto permitiendo fundamentarse en la distancia desde el punto de salida. Entre mayor sea el valor de la potencia, se tienen mayor énfasis en los puntos que estan a una distancia menor. Es decir, los datos más cercanos tienen más influencia y la superficie tendrá más detalles, aunque, el valor más común usado en la potencia es de 2. ^{7 8} La formula base que es empleada para el cálculo del valor del punto de interés es la siguiente:

\(z_{p}\): Valor del punto de interés

\(n\): Número de puntos que quedan dentro del radio de busqueda

\(z_i\): Valor de cada punto que es conocido y que están dentro del radio de busqueda

\(d_{i}^{p}\): Distancia entre los puntos y \(^{p}\) significa la potencia

A continuación, se muestra una ilustración del método de IDW tomado de la página esri ⁹

Donde, el circulo amarillo representa el radio de busqueda, el punto en el centro es el punto de interés, los puntos rojos representan los puntos dentro del radio de busqueda y los que son usados para determinar el valor del punto de interés, por último, los puntos negros indican los puntos que debido a su distancia desde el centro del radio de busqueda, no son incluidos en el cálculo.

• Polígonos Thiessen

También llamados Diagramas de Voronoi, se basan en organizar en polígonos las zonas que contienen una mayor densidad de puntos que están más cerca que otros, es decir, en cada polígono hay un conjunto de puntos con menor distancia entre cada uno, a diferencia de otros puntos más lejanos. No todas las zonas con mayor densidad de puntos están en simetría con la superficie de estudio formando figuras rectangulares, debido a que, no todas las superficies son planas o no todos los datos son graficados de forma igual; con respecto a la altura se pueden presentar variaciones en la elevación como montañas o colinas mostrando los datos en polígonos de más de 4 lados y con diferentes longitudes. Después de tener estos polígonos se calcula el área de estos asignando un punto central a cada uno mediante interpolación. ¹⁰

• Kriging

Es un método geoestadístico que asume los fenomenos que ocurren en la superficie terrestre y que tienen un componente aleatorio. En este método, se tiene en cuenta la relación directa que existe entre los puntos que son medidos, gracias a estos datos se pueden obtener superficies de predicción y también una mejor presición de las predicciones de los valores.En este método es usada la distancia o dirección de los puntos, para representar una correlación espacial que es usada para explicar la variación de la superficie según la variable empleada. El Kriging emplea varios pasos, como son, análisis estadístico exploratorio de los datos, modelado de variogramas, creación de la superficie y como alternativa opcional la exploración de la superficie de varianza. Este método es principalmente utilizado cuando se conoce que hay una influencia de la distancial correlacionada espacialmente de los datos. ^{11 12}

El kriging tiene cierta similitud con el método de interpolación IDW, ya que, pondera los valores medidos circundantes para calcular una predicción de una zona sin mediciones. La formula empleada en kriging es:

\[\hat{Z}(s_0) = \sum_{i=1}^{n}\lambda_{i}Z(s_i)\]

\(Z(s_i)\): Valor medido en la ubicación \(i\)

\(\lambda_{i}\): Ponderación desconocida para el valor medido en la ubicación \(i\)

\(s_0\): Ubicación del punto de interés

\(n\): Cantidad de valores medidos

En idw \(\lambda_{i}\) depende de la distancia a la ubicación del punto de interés. Con este método de kriging, las ponderaciones se basan en la distancia entre otros puntos, en la ubicación del punto de estudio y también de la disposición en el espacio de los puntos medidos.

Para el modelado de variogramas se empieza con un gráfico del semivariograma empírico el cual representa la distancia entre pares de los puntos en función de la semivarianza, la cual, representa la medida de la diferencia que existe entre esas distancias. Los puntos que se representan en este gráfico de semivariograma se deben ajustar con el modelo de variograma el cual es esencial para realizar la interpolación kriging. Este modelo de variograma tiene diversas formas como exponencial, esférica, Gaussiana o Matern.^13]

4. Presentación de resultados

A continuación, se muestran los resultados obtenidos al programar los mapas, gráficas y formulas de los tres métodos de interpolación mencionados anteriormente, usando los datos CHIRPS de precipitación en la zona del departamento del Amazonas.

4.1. Interpolación mediante el método IDW

En la siguiente imagen se puede observar el mapa del amazonas con el método de interpolación de IDW donde se observa en colores marrones las zonas interpoladas con menor precipitación en el departamento y en colores verdes oscuro las mayores precipitaciones de la zona. Cabe mencionar que para el siguiente mapa fue usada como potencia 2.

Adicionalmente, también se muestra el siguiente mapa interactivo indicando de igual forma la interpolación por el método IDW. A diferencia del anterior mapa, este muestra el número original de celdas raster que fueron obtenidas al descargar los datos de precipitación de CHIRPS, esto, con el fin de tener una mejor interpretación de los datos. Las zonas que presentan un color verde más claro indican menores precipitaciones y los colores verdes azulados un poco más oscuros indican mayor precipitación.

En el siguiente gráfico se quiso explicar de una manera distinta los resultados obtenidos con esta interpolación, dando como resultado el valor de precipitación observado (eje \(x\)) y los valores estimados mediante interpolación (eje \(y\)).

Conjunto con el gráfico, se obtuvo el siguiente valor de error medio cuadrático (RMSE) que corresponde a una medida de error absoluta que calcula el cuadrado de las desviaciones y mide la cantidad de error que existe entre dos variables con distintos datos, comparando un valor predicho y un valor observado o conocido. Este dato se calcula con la siguiente fórmula:

\(N\): Número de valores observados

\(y_i\): Valor estimado

\(\hat{y_i}\): Valor observado

Finalmente, en el siguiente mapa se puede observar el nivel de confianza que presentan los datos interpolados con el método de IDW, donde, los valores con color amarillo claro indican un menor intervalo de error, es decir, una mayor confianza de los datos, y los valores con colores más rojizos indican una menor confianza de los datos.

4.2. Interpolación mediante el método de polígonos Thiessen

A continuación, se muestra el mapa correspondiente a la interpretación de los datos de precipitación por medio de la interpolación con el método de los polígonos thiessen. Como se puede observar, los polígonos pueden ser visibles, adicionalmente, fue usado el número total de celdas raster del archivo original descargado de los datos CHIRPS de precipitación, para que la visualización de los polígonos fuera mucho más clara y más concluyente. De igual forma, como en el mapa anterior, los polígonos con colores marrones indican zonas con la menor precipitación y los polígonos con colores verde oscuro indican mayor precipitación en esa zona.

Este mapa adicional muestra de una forma interactiva la disposición de los polígonos y el número de celdas que se pueden encontrar en cada uno de ellos dependiendo de la cercanía a la que se quiera observar el mapa. En este caso, los colores de las celdas raster del mapa indican que los más rojizos presentan menores precipitaciones y los colores azul oscuro indica la zonas con mayor precipitación según la interpolación realizada. De igual forma, entre más cerca se esté de la zona el número de puntos automaticamente van acomodandose a la zona, y si también se pueden filtrar los datos por “Entrenamiento” y “Validación”.

4.3 Interpolación por el método de kriging

La siguiente gráfica indica el variograma, donde fue ajustado usando el modelo de variograma exponencial. Como se puede ver, los datos de distancia varían entre 0 y 600000 metros y la semivarianza va desde 0 hasta 120. El modelo exponencial fue el más adecuado para graficar los datos, ya que, con el modelo gaussiano solo era permitido un rango máximo de 200 datos.

Después de realizar el variograma, fue posible graficar el mapa del amazonas, que se muestra a continuación, con los datos de precipitacion por interpolación con el método de kriging. Al igual que algunis mapas anteriores, el color marrón indica menores precipitaciones y los colores verde oscuro indican mayor precipitación en esa zona. Para este mapa también fueron usadas la totalidad de celdas raster descargadas en CHIRPS con el propósito de que la visualización del gráfico fue mas ‘suave’ y más entendible.

En este mapa, se puede ver más claramente y de forma interactiva los datos obtenidos empleando este método de interpolación.

Una forma de visualización de los datos interpolados por el método de kriging según su varianza e incertidumbre se muestra en el siguiente mapa, donde se aprecia el nivel de variación de los datos de precipitación según la zona del departamento.

Por último, en el siguiente mapa se muestran los intervalos de confianza de los datos de precipitación obtenidos después de realizar kriging. Este nivel de confianza que se muestra se valora en un 95% de confianza total.

5. Análisis de resultados

Según los dos primeros mapas obtenidos para la interpolación IDW, se puede observar que los valores indicados de precipitación son congruentes con los valores que se indican en los datos originales de CHIRPS descargados, mostrando simetría con los mapas mostrados en la descripción de los datos. En los datos representados en la gráfica, se puede ver que los puntos tienen una alta variación con respecto a la línea guía que se muestra, indicando que los valores obtenidos por la interpolación con el método IDW fue efectiva en algunas zonas, pero se presentaron otras donde la variación a los datos esperados fue alta. Si el 100 % de los datos coincidieran, los puntos graficados estarían alineados con la línea punteada roja que se muestra. Adicionalmente, el valor del error medio cuadrático indica que el rango de error es muy alto, es decir que, la comparación entre los datos esperados y los datos obtenidos con la interpolación IDW tuvieron un rango de variación alto. Finalmente, con el gráfico realizado del nivel de confianza, se puede observar que la mayoría de los datos están en un intervalo de confianza muy alto, es decir, que el error es muy bajo, indicando que, aunque se presentaron algunos errores, en general los resultados de la interpolación fueron muy satisfactorios. Aunque en los dos primeros mapas parecieran no tener muchas variaciones con respecto al original, al ver los datos puestos en una tabla y comparados segun su nivel de error se puede comprobar que los datos varían bastante con respecto a los datos esperados.

En la interpolación mediante el método de polígonos Thiessen, se puede observar que se empiezan a presentar variaciones en los datos con respecto a la interpolación de IDW. Aunque, al comparar visualmente los datos que tiene cada uno de los polígonos del primer mapa y los datos de cada celda de los datos esperados, se puede ver cierta similitud, sin embargo, hay algunas pequeñas zonas que no coinciden del todo con los datos esperados. Se puede decir que es un método muy acertado pero, considerar algun otro método de comparación de los datos sería mucho más concluyente. No obstante, los datos del segundo gráfico son mucho más acertados si se observan los datos por cada celda, adicionalmente, los mapas interactivos al tener la opción de poder ver los datos más de cerca, se puede realizar una comparación más específica dependiendo de la zona del departamento. Sin embargo, al usar el método de interpolación con polígonos thiessen no es tan simple comparar los resultados gráficos con solo los mapas, sería muy gratificante poder tener otro método de comparación de los datos como gráficas para tener una conclusión certera de los datos.

Por último, en el variograma al usar el método exponencial, para la interpolación de kriging, se indicó al inicio de la gráfica que la función era creciente y simétrica con la línea exponencial, pero, después de una distancia de aproximadamente 150000 metros los puntos cambian su semivarianza constantemente aumentando y disminuyendo sin obtener datos concluyentes al método exponencial usado. Sin embargo, si se observa la continuidad de los puntos, se puede ver cierta homogeneidad de la gráfica donde se tiende a ver los puntos con ciertas formas parabólicas, indicando la relación de los datos segun la semivarianza. De igual forma, con respecto a los mapas se puede observar que los datos graficados fueron concluyentes y muy exactos en comparación con los datos de precipitación esperados, asimismo, el mapa de varianza mostró un intervalo bajo de variabilidad de los datos, es decir que efectivamente, los datos de la interpolación fueron congruentes. Por otro lado, el último mapa el cual indica el nivel de confianza de los datos obtenidos no mostró datos tan satisfactorios con un nivel de confianza en un rango medio en gran parte del mapa, y en los bordes del mapa se puede ver que los niveles de variación en los datos van de 6 a 8 mm. Dado esto, se puede decir que aunque los datos se acercaron a los esperados, este método no es totalmente certero en la aproximación de los datos.

6. Conclusiones

• El amazonas es una zona con variación climática media de la precipitación, debido a sus niveles de humedad y a su gran zona boscosa, presentando alteraciones desde 20 mm hasta 70 mm de precipitación en solo 5 días.

• En comparación con los tres métodos estudiados, el que tuvo más datos concluyentes, exactos y precisos fue el método de kriging. Aunque la variación de algunos datos fue de 6 a 8 mm los datos fueron más acertados y, adicionalmente, el variograma y los diferentes mapas realizados daban una mayor idea de qué tan similares fueron los datos en comparación con los originales.

• Según lo anterior se concluye que las técnicas de interpolación geoestadísticas son mucho más precisas y certeras en comparación con las técnicas deterministas.

• El método que brindó menor información en los datos fue el método de los polígonos de Thiessen debido a que el hecho de comparar la efectividad es muy subjetiva si no se tiene un gráfico adicional que brinde la información precisa acerca de cuáles son las variaciones en comparación con los datos reales.

• Aunque en el método de IDW se observó que hubo alta variación en los datos al observar el gráfico y el resultado del error medio cuadrático fue alto, se puede decir que el nivel de confianza de los datos obtenidos fue bastante satisfactorio en gran parte del departamento.

• La interpolación espacial es una técnica muy útil que brinda resultados en general muy precisos, evitando trabajos en campo de zonas muy grandes que cuestan mucho dinero y tiempo.

• Con los resultados obtenidos al emplear los tres métodos, se pudo observar que en general los datos fueron muy precisos en especial en la zona central del departamento del amazonas presentando los menores niveles de variabilidad.

7. Referencias

¹ Esri. ArcGIS Desktop. 2016. Understanding interpolation analysis.

² Planet Botany UWC. What is a spatial interpolation?.

³ Planet Botany UWC. Why interpolate?.

⁴ Gobernación del amazonas. Nuestro departamento.

⁵ Instituto SINCHI. Sistema de Información Ambiental Territorial de la Amazonía Colombiana. Clima.

⁶ Climate Hazard Center. CHIRPS: Rainfall Estimates from Rain Gauge and Satellite Observations. University of California, Santa Bárbara.

⁷ Planet Botany UWC. Inverse Distance Weighting (IDW).

⁸ Esri. ArcGIS Desktop. 2016. IDW

⁹ Esri. ArcGIS Desktop. 2016. How IDW works

¹⁰ Planet Botany UWC.Spatial analysis of categorical data using Neigbourhood Analysis (e.g. generation of soil maps).

¹¹ Planet Botany UWC.Statistical techniques using a semi-variogram for developing continuous surface models (Kriging).

¹² Planet Botany UWC. How kriging works

¹³ R for Spatial Statistics. Colorado State University, Oregon State University, and Humboldt State University. R Variograms & Kriging.

Olaya, V,. 2014. Sistemas de información Geografica. Creative commons.Pp 311-336.