Interpolación de datos de precipitación en Boyacá
Marly Daniela Quiroga Buitrago
24/06/2020
1. Introducción
La descripción espacial y temporal de variables climáticas es de gran utilidad para comprender el funcionamiento de procesos bio-físicos. La precipitación es una de las principales variables climáticas requeridas para la estimación de balances hídricos, recargas de frontera en la modelación de flujos de agua subterránea, evaluar procesos de erosión, así como para definir las condiciones climáticas actuales. Sin embargo, la representación espacial fiable de la precipitación, es particularmente difícil en zonas montañosas con escasa disponibilidad de estaciones climáticas en donde el efecto orográfico es grande (Hevesi et al., 1992; Huade et al., 2005).
Diversos métodos han sido desarrollados para predecir la distribución espacial de variables climáticas que difieren en su concepto y formulación matemática. Métodos globales de interpolación que permiten el procesamiento de datos topográficos y geográficos han sido ampliamente utilizados en diversas partes del mundo para generar mapas de precipitación (Ninyerola et al., 2000).
En términos generales, un método de interpolación es una herramienta que permite el cálculo de valores en puntos no muestreados, a partir de los valores recogidos en otra serie de puntos. El número de métodos distintos es muy amplio, y es importante reseñar que la bondad o no de uno u otro va ligada no solo al método en sí, sino también a la variable interpolada y al uso que se dé posteriormente a la capa resultante (Olaya, 2014). En este sentido, existen principalmente dos grandes grupos de modelos de interpolación, los determinísticos (o mecánicos), como polígonos de Thiessen, IDW e Isoyetas; y los geoestadísticos como Kriging y sus múltiples variacones (Mair & Fares 2011).
La interpolación de datos ofrece la ventaja de proyectar mapas o superficies continuas a partir de datos discretos. Dependiendo del tipo de datos analizados, su costo y dificultad de obtención determinan que tan valioso es finalmente el uso de la interpolación. Otro aspecto a mencionar, es que la precisión en el mapa generado, depende en gran medida de la estructura espacial de los datos, donde entre más fuerte la correlación espacial, mejor la calidad del mapeo (Kravchenko 2003). De esta manera, el objetivo de este trabajo es realizar un estudio comparativo de tres métodos de interpolación “polígonos de Thiessen, IDW y Kriging” para modelizar la precipitación durante los últimos cinco días del mes de abril de 2020 para el departamento de Boyacá a partir de datos de precipitación de CHIRPS (Centro de riesgos climáticos de Precipitación infrarroja con datos de estación).
2. Descripción de la zona de estudio
La zona de estudio corresponde al departamento de Boyacá.
Boyacá es un departamento ubicado en la cordillera oriental de los Andes; localizado entre los 04°39′10″ y los 07°03′17″ de latitud norte y los 71°57′49″ y los 74°41’35″ de longitud oeste. Cuenta con una superficie de 23,189 km², lo que representa el 2,03 % del territorio nacional. Limita por el norte con los departamentos de Santander y Norte de Santander, por el este con los departamentos de Arauca, Casanare y con el país vecino de Venezuela, por el sur con Cundinamarca, y por el oeste con el rio Magdalena, Caldas y Antioquia (Aranza, Y. et al. 2016).
El departamento de Boyacá cuenta con una población de 1’276.407 Habitantes (Proyección DANE 2015), está dividido en 123 municipios (Fig 1.), 44 corregimientos, 185 inspecciones de policía, así como, numerosos caseríos y sitios poblados. Su capital es Tunja.
Figura 1. División politico administrativa de Boyacá. Fuente IGAC, 2002.
En este territorio se presenta una diversidad de accidentes geográficos que forman las regiones del valle del río Magdalena, la cordillera Oriental, el Altiplano Cundiboyacense y el piedemonte de los llanos orientales. Gracias a ello, en el departamento se presentan todos los pisos térmicos con temperaturas desde los 35 °C en Puerto Boyacá, hasta temperatura bajo cero grados, en la Sierra Nevada de Güicán. El departamento presenta altitudes que varían entre 145-3100 m.s.n.m a excepción de los nevados que están a 5.490 m.s.n.m (Gobernación Boyacá, 2012).
El sistema hidrográfico de Boyacá (Fig 2.) está constituido por numerosos ríos y quebradas que nacen en la cordillera Oriental, afluentes directos o de alguno de los tributarios de los ríos Magdalena, Meta y Arauca; entre los ríos que confluyen a la cuenca del Magdalena están el Ermitaño, Negro, Minero, Suárez, Sutamanchán, Sáchica, Chíquiza, Iguaque, Arcabuco o Pómeca, Ubazá, Riachuelo, Moniquirá, Chicamocha, Chorrera, Tuta, Pesca, Tota, Saguera, Sasa, Cambas, Loblanco, Rechiminiga, Chitano, Susacón. También cuenta con lagunas de gran belleza natural, entre ellas son notables las de Tota, Sochagota y Fúquene (Aranza, Y. et al. 2016).
Figura 2. Hidrografía de Boyacá. Fuente IGAC, 2002.
La agricultura en el departamento de Boyacá, es una de las principales actividades económicas junto con la producción pecuaria y la minería (Fig 3). Boyacá, rico en recursos naturales, con una infraestructura energética, vial y de servicios adecuada y una posición geográfica privilegiada ha centrado su actividad económica en la agricultura tradicional, con un alto porcentaje de la producción dedicada a abastecer el mercado interno, se caracteriza por ser uno de los departamentos con mayor producción papera en el país, aunque también se cultiva maíz, cebada, cebolla, caña, café, cacao, frutales, hortalizas y muchos más (Secretaría de Agricultura, 2019).
Figura 3. Economía de Boyacá. Fuente Secretaria de Fomento Agropecuario
En cuanto a clima, el departamento presenta una alta variabilidad territorial en la distribución de la precipitación (Fig 4.). Las áreas con menores lluvias, de 500 a 1000 mm anuales, se localizan a lo largo del altiplano cundiboyacense y a partir de esta zona las lluvias se incrementan tanto al oriente como al occidente del departamento. Al oriente ocurren las lluvias orográficas correspondientes a la vertiente oriental de la cordillera oriental, registrando valores en el rango de 2000 a más de 5000 mm. Al occidente, en dirección al valle del Magdalena, las lluvias oscilan entre 2000 y 3000 mm al año (Gobernación Boyacá, 2012)..
El régimen de lluvias es variado, al centro y occidente del departamento se observa un comportamiento bimodal, con una época seca principal a principios de año y una ligeramente menos marcada, a mediados de año. Las temporadas de lluvia en estas dos zonas tienen lugar en los meses de marzo-mayo y septiembre-noviembre. Al oriente del departamento, el régimen es el característico de la Orinoquia, con una única estación seca de diciembre a marzo y una temporada lluviosa que se prolonga durant el resto del año con máximos de lluvia en los meses de junio-julio. El número de días con lluvia oscila entre 100 y 150 sobre el altiplano cundiboyacense y aumenta hasta 250 en los municipios localizados sobre la vertiente de la Orinoquia. El resto del departamento registra entre 150 y 250 días lluviosos al año (Gobernación Boyacá, 2012)..
Figura 4. Precipitación total anual del departamento de Boyacá. Fuente IDEAM, 2019
3. Descripción de datos y métodos
3.1. Datos: Se utilizaron datos de precipitación de CHIRPS (Centro de riesgos climáticos de Precipitación infrarroja con datos de estación). Estos son un conjunto de datos de precipitaciones cuasi global de más de 30 años, con un alcance de 50 ° S-50 ° N (y todas las longitudes). Se basa en enfoques previos de técnicas de interpolación “inteligentes” y estimaciones de alta resolución con datos de la estación in situ para crear series de tiempo de lluvia cuadriculadas o estimaciones de precipitación diarias, pentadales y mensuales desde 1981 al presente, para el análisis de tendencias y el monitoreo estacional de la sequía (Funk, C. et al. 2015).
Chirps combina los datos de la estación para producir un producto de información preliminar con una latencia de 2 días y un producto final con una latencia promedio de 3 semanas. Casi todos los datos tienen una resolución espacial de 0.05 ° × 0.05 ° grados que equivale a 5.5km x 5.5km. Estos datos se proporcionan en formatos ráster como lo es NetCDF, GeoTiff y Esri BIL, sus unidades son mm por período de tiempo, por ejemplo, mm por día, mm por pentad, mm por mes. Otro aspesto importante es que maneja tres dominios espaciales principales que son Global (7200 × 2000 píxeles, 180 ° W a 180 ° E, 50 ° N a 50 ° S), África (1500 × 1600 píxeles, 20 ° W a 55 ° E, 40 ° N a 40 ° S) y Centroamérica-Caribe (720 × 350 píxeles, 93 ° W a 57 ° W, 23.5 ° N a 6 ° N)(Funk, C. et al. 2015).
Los datos de precipitación utilizados este trabajo se descargaron del sitio FTP de CHIRPS. Los datos pertenecen al dominio Global_pentad del directorio tifs, específicamente al pentad 6 del mes de abril de 2020, es decir, los últimos cinco días que son del 26 al 30 del respectivo mes. Se escogió el directorio tifs porque es un formato ráster con facilidad de manejar.
Teniendo en cuenta que los datos de precipitación son un conjunto de datos de cobertura global y su CRS es un sistema de coordenadas geográficas, se carga un archivo shape que represente el área de interés, en este caso Boyacá, que también esta en coordenadas geográficas, para recortar y trazar los datos de precipitación de cultivos por extensión de Boyacá (Fig 5).
Figura 5. Precipitación del departamento de Boyacá del 26 al 30 de abril, 2020..
Algo adicional que se realizó fue modificar el tamaño de las celdas del ráster antes de convertirlo a puntos para obtener un mapa de puntos mucho más legible, esto se hace mediante la función aggregate de la librería ráster, donde con fact = 3 dará como resultado un nuevo objeto Raster con 3 * 3 = 9 veces menos celdas. Es decir, si antes las celdas ráster tenían una resolución de 0.05° con un total de 14’400.000 celdas, ahora es de 0.15° con 1’600.800 celdas.
Posteriormente, como los datos que se tienen están en datum WGS84, los cuales están en coordenadas geográficas, no se puede hacer interpolación o no se recomienda. Para hacer la interpolación se hace la conversión a coordenadas planas (más recomendable). Adicional a la transformación de coordenadas o reproyección, se tienen que colocar los datos en el formato apropiado para que las librerías de R puedan trabajar, las librerías que hay trabajan bien con el tipo de estructura SpatialPointsDataFrame, en otras palabras, es la conversión de ráster a puntos utilizando la función * rasterToPoints* de la librería ráster, obteniendo finalmente un mapa de puntos de precipitación para el departamento de Boyacá, estos datos muestran que en ese periodo de tiempo la precipitación estuvo entre 4.7-93 mm y se puede observar en la Figura 6.
Figura 6. Mapa de puntos de la precipitación del departamento de Boyacá del 26 al 30 de abril, 2020.
3.2. Métodos: A continuación se realiza una breve descripción de los tres métodos de interpolación empleados en el desarrollo de este informe:
3.2.1. Polígonos de Thiessen
Es el método más básico y simple de interpolación, siendo especialmente apropiado cuando las variables son cualitativas. Está basado en las áreas de influencias de las estaciones, que se crean al unirlas entre sí, trazando las mediatrices de los segmentos de unión. Las intersecciones de estas mediatrices determinan una serie de polígonos en un espacio bidimensional alrededor de un conjunto de puntos de control o estaciones. A partir de los polígonos creados, simplemente se asignan pesos como segmentos del área de influencia sobre elk área total (Moreano, 2008).
El área definida por cada polígono es el área de influencia de la estación. La lluvia total estimada es el resultado de la sumatoria del producto de los datos de precipitación promedio de las lluvias registradas en cada estación por el área relativa de influencia de las estaciones (Arias, 2001). Este método tiene la ventaja de aprovechar la información de estaciones vecinas, sin embargo su principal falencia radica en que sólo considera el posicionamiento de las estaciones y una superficie plana de influencia para cada una, sin considerar las diferencias topográficas que se pueden presentar. Además, asume que la precipitación de la estación es la misma de la zona que representa geométricamente, lo cual no siempre es cierto(Pizarro, et.al, 2003). La ecuación general del método polígonos de thiessen es:
3.2.2. Interpolación ponderada de distancia inversa (IDW)
El método IDW está basado principalmente en la inversa de la distancia elevada a una potencia matemática para estimar el valor de una celda. En este método los puntos de muestreo se ponderan durante la interpolación de tal manera que la influencia de un punto en relación con otros disminuye con la distancia desde el punto desconocido que se desea crear, es decir, la influencia de la ponderación decae mientras la distancia hacia el punto nuevo se incrementa (Díaz, 2008; Olaya, 2014). La ecuación general es:
La precisión del modelo de interpolación se define por la medida en que el valor predicho de una ubicación coincide con el valor real en esa ubicación. Una manera de estimar el error o incertidumbre de los resultados es mediante validación cruzada.
La validación cruzada, conocida como el método “jacknife” remueve consecutivamente un valor de los datos y el valor ausente es interpolado con los datos remanentes para luego ser comparado con el dato removido u observado. Comparar un porcentaje de las observaciones reales con las observaciones interpoladas, entre más similares sean estos valores existe un aumento en la efectividad de predicción (Mateu y Morrell, 2003; Henríquez et al., 2013)
3.2.3. Kriging
El kriging es un método de interpolación estocástico, exacto, aplicable tanto de forma global como local. Se trata de un método complejo con una fuerte carga geo-estadística, del que existen además diversas variantes (Olaya. 2014). El Kriging es un conjunto de técnicas geoestadísticas de interpolación que predicen valores de la variable de interés en lugares no muestreado, proporciona un mapa de predicciones en todas las ubicaciones de interés y una varianza estimada del error de predicción. En otras palabras,en vez de decir solamente cuánta lluvia hay en lugares específicos, kriging también le dice la probabilidad de cuánta lluvia habrá en un lugar específico.
El método Kriging cuantifica la estructura espacial de los datos mediante el uso de variogramas llamados algunas veces semivariogramas debido a su similitud en el cálculo y los predice mediante la interpolación, usando estadística. Se basa en la suposición de que entre dos puntos cercanos debe existir algún tipo de correlación que no existe entre puntos lejanos, de forma que asigna peso, se asume que los datos más cercanos a un punto conocido tienen mayor peso o influencia sobre la interpolación, influencia que va disminuyendo conforme se aleja del punto de interés (Villatoro, 2008).
El semivariograma se trata de un gráfico en él que se representa la semivarianza en función de la distancia a la que se encuentran los datos, describe la variabilidad espacial del fenómeno de interés y responde a la siguiente pregunta ¿Qué tan parecidos son los puntos en el espacio a medida que estos se encuentran más alejados? y proporciona información del comportamiento espacial de la variable (Núñez, 2014).
La ecuación general del método krigin es:
4. Resultados
4.1. Interpolación Polígonos de Thiessen: Para esta interpolación se subió el mapa de puntos y el archivo shape que representa el área de interés “Boyacá”, se crea una superficie teselada recortada a los límites del departamento para obtener el resultado de la interpolación que se observa en la figura 7. Revisar el Anexo 2. Interpolación por el método Polígonos de Thiessen de datos de precipitación en Boyacá
Figura 7. Interpolación polígonos de Thiessen de precipitación del departamento de Boyacá del 26 al 30 de abril, 2020.
4.2. Interpolación IDW: Se siguió el código que está en el Anexo 3. Interpolación ponderada de distancia inversa (IDW) de datos de precipitación en Boyacá y el resultado de la interpolación se observa en la figura 8.
Figura 8. Interpolación IDW de precipitación del departamento de Boyacá del 26 al 30 de abril, 2020.
Validación cruzada
Rutina de validación de omisión para medir el error en los valores interpolados y se trazan las diferencias (fig 9). El error cuadrático medio (RMSE) es [1] 10.52289.
Figura 9 Validación cruzada de la interpolación IDW.
Intervalos de Confianza: En la figura 10 se observa el mapa de intervalo de confianza del 95% del modelo de interpolación utilizando un parámetro de potencia de 2 (idp = 2.0).
Figura 10. Intervalo de confianza del 95% de la interpolación IDW.
4.3. Interpolación Kriging: Se ajustó el modelo de variograma Gaussiano a los datos de precipitación como se observa en la figura 11:
Figura 11. Variograma de modelo Gaussiano de precipitación del departamento de Boyacá del 26 al 30 de abril, 2020.
A partir del variograma Gaussiano y definir un modelo de tendencia, se generó la superficie interpolada kriged como se evidencia a continuación:
Figura 12. Interpolación método Kriging de precipitación del departamento de Boyacá del 26 al 30 de abril, 2020.
Se generó el mapa de varianza y el mapa de intervalo de confianza de la interpolación kriging a partir del objeto dar.krg.
Figura 13. Mapa de varianza de la interpolación Kriging.
Figura 14. Intervalo de confianza del 95% de la interpolación Kriging.