Métodos de interpolación utilizando datos de precipitación CHIRPS, para el departamento de Norte de Santander

Ana María Montaño Hernández

Junio 24 de 2020

Geomática Básica

1. Introducción

El análisis espacial se centra en el estudio de los componentes del espacio, definiendo sus elementos y la forma como éstos se comportan bajo determinadas condiciones, valiéndose de herramientas técnicas que ayuden a comprender la dinámica del espacio (Madrid A. & Ortiz L., 2018). Su objetivo es derivar nueva información a partir de datos existentes que le permitan tomar decisiones, para ello aplica operaciones geográficas, estadísticas y matemáticas a los datos que se representan en el mapa (Plataforma ArcGis., s.f)

Además, el análisis espacial es aplicable a una gran cantidad de disciplinas que poseen la necesidad de recolectar datos en puntos específicos. El conjunto de datos discretos debe ser convertido en una representación continua de la variable, y para pasar de mediciones puntuales a superficies continuas que expresen una propiedad continua de la naturaleza, es necesario aplicar un proceso llamado interpolación (Planet botany UWC., s.f). La interpolación es un proceso que permite obtener un ráster (superficie continua) y predecir valores para las celdas de este a partir de un número limitado de puntos de datos de muestra. Para generar las superficies ráster existen dos tipos de modelos: los modelos deterministas que emplean una función matemática para predecir los valores desconocidos; como el IDW, y los modelos geo estadísticos que asumen que las superficies no obedecen una fórmula fija y emplean las leyes de la probabilidad; como kriging. (Planet botany UWC., s.f)

Por otra parte, la adquisición de datos puntuales se puede realizar con distintas variables, en este caso, se emplea la variable discreta de la precipitación, que se mide en estaciones pluviométricas o con imágenes de satélite. La interpolación permite solucionar la problemática de que no existen estaciones en todas partes, estimando valores en otros sitios de interés por medio de la generación de superficies continuas. Esto permite comprender cómo es la distribución de la lluvia en el departamento, para su posterior análisis y toma de decisiones.

Por lo tanto, este trabajo pretende por medio del análisis espacial y la interpolación, utilizar estadística espacial para determinar lo que puede contar un conjunto de puntos de muestra sobre la distribución de la precipitación en el departamento de Norte de Santander, utilizando datos CHIRPS.

2. Descripción de la zona de estudio

Zona de estudio: Departamento de Norte de Santander.

Norte de Santander es uno de los 32 departamentos de Colombia, ubicado en la zona nororiental del país. Geográficamente, se localiza entre los 06º56’42’ y 09º18’01’’ de latitud norte y los 72º01’13’’ y 73º38’25’’ de longitud oeste.(Gobernación de Norte de Santander, 2018).

Tiene una extensión de 22.130 km2, que equivalen al 1.91% del millón ciento cincuenta y nueve mil ochocientos setenta y un kilómetros cuadrados (1.159.871,41 km2) del territorio nacional. Limita por el norte y el oriente con la República Bolivariana de Venezuela, por el sur con los departamentos de Arauca y Boyacá y por el occidente con los departamentos de Santander y Cesar.(Gerencia Departamental Norte de Santander, 2015) También, cabe resaltar que hace parte de la Región Andina de Colombia y de la región centro oriente, según el Departamento Nacional de Planeación.

Según la Contraloría General de la República (2015), Norte de Santander se divide administrativa y políticamente en 40 municipios y 6 subregiones; posee 96 corregimientos, 1719 veredas, 38 caseríos y 24 asentamientos indígenas. La división político administrativa del departamento se puede observar en la ilustración 1.

Ilustración 1. Mapa geográfico del departamento de Santander

Departamento del Norte de Santander 2003. Fuente: IGAC - Instituto Geográfico Agustín Codazzi

Según el Censo Nacional de Población y Vivienda realizado por el DANE en el año 2018, Norte de Santander posee una población de 1.346.806 habitantes, lo cual representaría aproximadamente un 2,84% del total de la población del país. Además, se calcula que el 62,39% de la población departamental se concentra en el área metropolitana, donde habitan alrededor de 824.045 personas y tan sólo en la ciudad de Cúcuta se concentra el 47,77% de la población.

Climáticamente, presenta un clima tropical modificado por la altitud y las variaciones del relieve, lo que ocasiona diferencias entre los elementos climáticos. La variación altitudinal va desde alturas menores a 100 m.s.n.m en las cuencas hidrográficas de los ríos Catatumbo y Magdalena; hasta alturas mayores a los 4.000 m.s.n.m en los páramos (TodaColombia, 2019). Por otra parte, según el IDEAM las mayores temperaturas medias anuales son de 28°C y se registran en los municipios de Tbú y Cúcuta, ubicados en la zona fronteriza con Venezuela. Las menores temperaturas medias pueden ser inferiores a los 14°C y se registran en sectores altos de municipios como Mutiscua, Silos y Cacota. Haciendo referencia a los pisos térmicos, Norte de Santander presenta desde piso térmico cálido a partir de los 50 m.s.n.mm hasta el páramo alto a más de 3.600 m.s.n.m. Esta variedad climática es consecuencia del dominio de relieves de montaña (Gobernación de Norte de Santander, 2018).

En cuanto a la precipitación, el departamento se caracteriza por tener un régimen bimodal, lo que indica la ocurrencia de dos períodos de lluvia y dos períodos secos cada año. La primera época de lluvias se extiende por lo general, de marzo a junio y la segunda época va de septiembre a noviembre.Sin embargo, esto puede variar un poco en cada municipio; por ejemplo en Pamplona, Cúcuta, Sardinara y Tibú la primera temporada de lluvias está centrada en el mes de abril, mientras que en San Calixto se centra en el mes de julio. Los municipios con mayor precipitación son Tibú con 2.554 mm/año y Toledo con más de 3.000 mm/año, ubicados en el noreste en la cuenca del Catatumbo y al sureste en la cuenca del Margua respectivamente. Por otra parte, las zonas con menores precipitaciones se encuentran en la cuenca baja del río Pamplonita y Zulia, como los municipios de Cúcuta, Villa del Rosario y Los Patio; cada uno con promedios anuales de 1.000 mm. (Gobernación de Norte de Santander, 2018).

En la figura 1, se muestra el comportamiento de la precipitación promedio en el departamento, donde se observa su régimen bimodal característico, teniendo los dos picos de mayores lluvias en los meses de abril-mayo y octubre. Esta figura se elaboró tomando promedios climatológicos (1981-2010) de una base de datos de las estaciones meteorológicas en Colombia.

Figura 1. Precipitación de Norte de Santander

Fuente: propia

3. Descripción de datos y métodos

3.1 Datos:

En el presente informe se van a emplear datos de precipitación CHIRPS. La precipitación es una variable discreta que se mide en estaciones meteorológicas y pluviométricas, es decir, se mide en puntos determinados que tienen asociada una latitud y longitud geográfica (Arias et al., s.f). Para aplicar la interpolación, es posible asumir que entre dos estaciones meteorológicas la precipitación varía en forma continua.

Por su parte, CHIRPS: Climate Hazards Group InfraRed Precipitation with Station data. (Grupo de Riesgos Climáticos de precipitación infrarroja con datos de estación); es un conjunto de datos de precipitación global que se emplea hace más de 35 años. Desde 1981, incorpora imágenes de satélite con una alta resolución de 0.05°, equivalente a 5,5 km aproximadamente, y datos de estaciones meteorológicas, con el fin de crear series de tiempo de lluvia para posteriores análisis de tendencias y monitoreo estacional de sequías (UC Santa Bárbara, s.f)

Estos datos CHIRPS son elaborados por The Climate Hazards Center (El Centro de Peligros Climáticos), de la Universidad de California, Santa Barbara. El CHC es una alianza de científicos multidisciplinarios que utilizan modelos climáticos y de cultivos, observaciones satelitales de la tierra, y conjuntos de datos; para predecir y monitorear sequías (UC Santa Bárbara, s.f). Además, estos datos se crearon en colaboración con científicos del Centro de Observación y Ciencia de los Recursos de la Tierra (EROS) del USGS (UC Santa Bárbara, s.f)

Por otra parte, los científicos de la Universidad de California presentan convenios con todos los países para recibir diariamente los datos de las estaciones meteorológicas, y poder integrar los datos a los de las imágenes satelitales para posteriormente crear un ráster de precipitación. La estructura ráster es útil para realizar monitoreo en agricultura y caracterizar sequías.

Adicionalmente, CHIRPS produce datos de precipitación de forma diaria, pentadal (cada cinco días) y mensual. En el presente informe se toman datos CHIRPS de forma pentadal, que corresponden al período del 26 al 30 de abril del 2020.

Cabe mencionar, que en el informe se emplea un shapefile del área de interés, en este caso, el departamento de Norte de Santander. Es un objeto de la clase SpatialPolygonsDataFrame, está en coordenadas WGS 84, tiene 40 municipios y 9 variables. Tanto los datos CHIRPS como el shapefile del departamento se encuentran con el mismo sistema de coordenadas.

En cuanto a la obtención y procesamiento de los datos; primero, se descargan los datos correspondientes al último pentadal del mes de abril. Se les realiza un pre-procesamiento ya que son datos de tipo ráster y se necesitan convertir a vector en geometría de puntos, para poder aplicar los distintos métodos de interpolación. Se descargan datos en formato ráster .tif y el formato de compresión es .gz.

Por otra parte, debido a que los datos CHIRPS son datos globales, es necesario recortarlos de tal manera que se obtengan solo los de la zona de interés. Para ello, se usa la función crop de la librería raster. Se aplica además la función mask para colocar valores nulos a todo aquello que quede fuera de los límites del departamento.

Después de realizar esto se obtiene un Ráster con las siguientes características: 49 filas, 32 columnas, 1568 celdas, una resolución de 0.05° y coordenadas geográficas en WGS84. Además, los valores de precipitación en el departamento en los últimos cinco días el mes de abril varían entre 2.060457 y 38.14587 mm.

Mapa 1. Fuente:propia

En el mapa 1 se observan las celdas que tienen un tamaño de 5,5km. Además, se puede identificar que las mayores precipitaciones ocurren en los municipios de La Esperanza, Cáchira y Sálazar, con valores entre 25 – 35 mm; mientras que, en los demás municipios la precipitación es baja, con valores entre 5 y 15 mm aproximadamente.

Mapa 1 utilizando librería leaflet

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

Otra posibilidad que ofrece R para representar los datos CHIRPS (tomando el 70% de los datos) de manera más interactiva, se ilustra a continuación en el mapa 2.

Mapa 2, Precipitación en Norte de Santander

Luego, se utiliza la función rasterToPoints de la libería Raster para obtener como resultado un objeto espacial y se pasan los datos Ráster a puntos, para poder aplicar la interpolación. Se obtiene un objeto de la clase: SpatialPointsDataFrame con 721 puntos, que corresponden al centro de cada celda. El atributo corresponde a la lluvia y tiene coordenadas WGS84. En el mapa 3 se observan los 721 puntos que se van a emplear en la interpolación.

Mapa 3. Fuente:propia

3.2 Métodos

ArcGis (2016) afirma que la interpolación es un proceso que permite predecir valores para celdas que no fueron muestreadas en un ráster a partir de un número limitado de puntos de datos de muestra. Se puede usar para predecir valores desconocidos para cualquier punto geográfico, como elevación, lluvia, concentraciones químicas y niveles de ruido. Además, se basa en el supuesto de que los objetos distribuidos espacialmente están correlacionados espacialmente; en otras palabras, las cosas que están juntas tienden a tener características similares. (Planet botany UWC., s.f) El producto de la interpolación es un modelo de superficie continua que permite analizar el fenómeno de interés en toda la superficie. Además, la interpolación se basa en la creación de una estructura regular a partir de datos irregularmente distribuidos.

Estos modelos de superficie permiten almacenar información en un SIG, y proporcionan un número infinito de valores continuos en el área de interés que sirven para medir variables que cambian continuamente en el espacio, como es el caso de la precipitación. Estos modelos pueden ser de tipo ráster o TIN.En el presente informe se emplea el modelo ráster.

Por otra parte, existen dos modelos para generar superficies ráster: Los modelos deterministas y los estadísticos.

3.2.1 Modelos deterministas

Las técnicas de interpolación determinísticas producen una superficie continua a partir de los puntos de muestra basándose por ejemplo, en el grado de similitud que existe entre ellos (ArcGis, 2019) Estos modelos se caracterizan por usar interpoladores determinísticos es decir, se basan en que el fenómeno de interés se comporta según una fórmula conocida y además, no emplean elementos probabilísticos. Entre estos métodos se encuentran: IDW (Distancia inversa ponderada), Natural Neighboor, de spline, entre otros.

3.2.1.1 Método de interpolación IDW (Distancia inversa ponderada)

Es un modelo determinístico que infiere que cada uno de los puntos empleados tiene una influencia que va disminuyendo con la distancia desde su ubicación muestreada, por lo tanto, la ponderación se realiza principalmente con los puntos más cercanos a la celda a la cual se le está calculando la precipitación. En este método, es posible utilizar la totalidad de los puntos o una cantidad específica que entre en el radio de búsqueda (Planet botany UWC., s.f). IDW generalmente se aplica a datos muy variables y la superficie interpolada que se genera, por lo general va a presentar valores menores al dato máximo de la muestra, y valores mayores al dato menor.

Según ArcGis, 2019 IDW se considera un interpolador exacto, es decir, esta técnica de interpolación predice un valor que es idéntico al valor medido en una ubicación que fue muestreada.

Además, IDW calcula el valor para cada nodo (punto de intersección o conexión de varios puntos de información) o coordenada dada dentro de la cuadrícula mediante una media ponderada de los puntos de influencia seleccionados (Olaya, 2014).

A continuación, se observa la fórmula que se emplea en IDW:

Fórmula de IDW. Fuente: Gis geography

Donde,

Zi, corresponde a los valores que se conocen.

n, corresponde a la cantidad de valores que entraron en el radio de búsqueda establecido.

〖di〗^p, corresponde a la distancia entre cada punto y el sitio de interés. P, corresponde a la potencia.

Asimismo, de la fórmula se puede afirmar que si la distancia es mayor, existirá un número más pequeño en el denominador y el efecto sobre el valor que se desea interpolar será menor, que en distancias más cortas. Por otro lado, Según Planet botany UWC, otros parámetros que influyen en el resultado de IDW son:

• (Power) o potencia: Esta variable define la influencia de los puntos vecinos usando una tasa exponencial, donde va disminuyendo la importancia de los puntos, a medida en que se encuentren más lejos. Por lo tanto, una alta potencia aumenta la influencia relativa de los puntos más cercanos, mientras que una baja potencia pondrá énfasis en los datos más alejados. Su valor oscila entre 1 y 10, pero el más utilizado en idp= 2.0

• (Search Radius) o radio de búsqueda: Se define como el tamaño del círculo, que está alrededor de cada uno de los nodos, dentro del cual va a buscar los puntos del conjunto de datos original y estos se van a ponderar, de acuerdo a su distancia. En la figura 2, se observa que se genera un radio de búsqueda alrededor de cada nodo de la cuadrícula y en este se van a seleccionar los puntos que van a influir en el cálculo. Este radio de búsqueda puede ser fijo o variable, limitando el número de puntos que se van a utilizar para calcular cada celda interpolada.

• (Barrier) o barrera: Se puede emplear una línea o un polígono como barrera, limitando la búsqueda de puntos que entran en la ponderación. Los polígonos y líneas empleadas pueden representar cualquier interrupción en el paisaje.

Figura 2. Radio de búsqueda

Radio de búsqueda. Fuente: Planet botany

3.2.2 Modelos estadísticos

Las técnicas estadísticas o también llamadas técnicas geoestadísticas (porque se está aplicando a datos geoespaciales), se basan en modelos que incluyen un autocorrelación, es decir, estudian la relación estadística que existe entre los puntos medidos (Planet botany UWC., s.f). Además, estos modelos no sólo producen una superficie de predicción, o superficie continua, sino que proporcionan una medida de la precisión o certeza de esas predicciones (Planet botany UWC., s.f).

Los modelos estadísticos incluyen a el método de kriging ordinario y kriging universal.

3.2.2.1 Método de kriging

Kriging es un método estadístico que genera una superficie continua o estimada, a partir de un conjunto de puntos dispersos en la zona de interés, cada uno de ellos con un valor z. Este modelo supone que la distancia entre los puntos muestreados refleja una correlación espacial entre ellos, que puede emplearse para explicar la variación de los valores en la superficie (ArcGis, 2016).A su vez, el grado de correlación espacial se rige por la Primera Ley de la Geografía o Ley de Tobler, que establece que las cosas que están cerca en el espacio se van a parecer más entre sí, que aquellas que se encuentran a mayor distancia.

Adicionalmente, Kriging es una técnica o modelo estadístico que utiliza un semi-variograma para desarrollar un modelo de superficie continua (Planet botany UWC., s.f).

Además, se puede afirmar que usa interpoladores estocásticos o probabilísticos que asumen que las superficies no obedecen una fórmula fija, que se rigen por las leyes de la probabilidad con una distribución aleatoria de los datos, y estudian tanto el efecto que tiene la cercanía de los puntos como la aleatoriedad entre ellos.

Por otra parte, según Olaya 2014, para obtener una interpolación óptima se deben cumplir los siguientes requisitos:

• El error de predicción debe ser mínimo.

• Los puntos cercanos deben tener pesos mayores que los lejanos.

• La presencia de un punto cercano en una dirección dada debe restar influencia a puntos en la misma dirección pero más lejanos.

• Puntos muy cercanos con valores muy similares deben agruparse, de tal forma que no aparezca sesgo por sobremuestreo.

• La estimación del error debe hacerse en función de la estructura de los puntos, no de los valores.

La fórmula general de kriging, similar a la de IDW, también realiza una suma ponderada de los datos, y viene dada por:

Fórmula de kriging. Fuente: Columbia Public Health

Donde,

Z(Xi), corresponde al valor medido en la ubicación i.

λi , corresponde a una ponderación desconocida para el valor medido en la ubicación i.

Xo, corresponde a la ubicación de la predicción.

N, corresponde a la cantidad de valores medidos.

A diferencia de IDW, la ponderación de los datos en kriging no solo está basada en la distancia entre los puntos medidos y la ubicación de la predicción, sino en la disposición espacial de dichos puntos (ArcGis, 2016).

Por otro lado, según ArcGis 2016, para desarrollar el método de kriging hay que seguir los siguientes pasos:

• Análisis estadístico exploratorio de los datos
• Modelado de variogramas
• Creación de la superficie
• Análisis de la varianza

Variografía

Se conoce como variografía al modelado espacial que se le aplica a un conjunto de datos. Por su parte, un variograma se define como una herramienta estadística que permite definir la correlación espacial de los datos, es decir, los valores de dependencia estadística de estos, y se fundamentan en el concepto de semivarianza. La semivarianza expresa la medida de la auto correlación espacial de una variable determinada, entre dos puntos (Olaya, 2014)

También, puede definirse un semivariograma como la función básica que describe la variabilidad espacial de un fenómeno de interés y permite responder a la pregunta ¿Qué tan parecidos son los puntos en el espacio a medida que estos se encuentran más alejados? . Lo anterior, relacionando las distancias entre los puntos y sus semivarianzas (Gallardo, 2006)

Se debe comenzar con un gráfico de semivariograma empírico. Este gráfico va a expresar cuanto varía la propiedad z, o el valor z, según la distancia que existe entre los puntos, por medio de una fórmula que calcula la diferencia cuadrada entre estas dos distancias (ArcGis, 2016).

Debido a que van a existir gran número de pares de puntos (ubicaciones) será difícil graficar estos resultados y, por tanto, estos pares se agrupan en intervalos. El semivariograma empírico que se obtiene va a representar de forma gráfica los valores promediados de semivarianza en el eje Y, y la distancia en el eje X, y se obtiene una nube de puntos como la mostrada en la figura 3.

Figura 3. Semivariograma empírico

Ejemplo de un semivariograma empírico.Fuente: ArcGis, 2016

Ajuste del modelo

Luego de obtener un semivariograma empírico, es necesario seleccionar y ajustarle un modelo, ya sea representado por una función o por una curva continua. Lo anterior, para asegurar que las predicciones de kriging tengan varianzas positivas.

Entre los modelos se encuentran:

• Circular

• Esférico

• Exponencial

• Gaussiano

• Lineal

Igualmente, el modelo que se seleccione va a tener influencia en la predicción de los valores a interpolar, y cada uno de estos está diseñado para ajustarse de forma precisa a diferentes fenómenos y variables (ArcGis, 2016).

Por otra parte, según Gallardo 2006, el ajuste a un modelo permite extraer una serie de parámetros que son los que van a ser usados para la interpolación geostadística de kriging, que definen el grado y escala de variación espacial y se observan en la figura 4.

Figura 4. Parámetros de un semivariograma

Componentes del semivariograma.Fuente: ArcGis, 2016

- Rango: Después de una cierta distancia, el modelo comienza a nivelarse. Se conoce como rango a la distancia a la cual el modelo comienza a aplanarse. También, puede definirse como la distancia hasta la cual hay correlación espacial y valor en el que se alcanza la máxima varianza.

- Nugget: Teóricamente, conforme la distancia de separación tiende a cero, el valor de la semivarianza también tiende a este valor. Sin embargo, a distancias muy pequeñas, el semivariograma presenta un efecto nugget que se evidencia en valores de semivarianza mayores que cero. Por lo tanto, representa una variabilidad que no puede explicarse mediante la estructura espacial y corresponde al valor en que el modelo intercepta el eje Y.

- Sill: Es el valor en el eje Y, que corresponde al valor en el cual el modelo alcanza el rango.

- Partial sill: Corresponde a la diferencia entre el sill y el nugget.

- Nugget: Conforme la distancia tiende a cero, el valor de la semivarianza tiende a este valor. Representa una variabilidad que no puede explicarse mediante la estructura espacial

Finalmente, cuando ya se tiene el análisis de dependencia de los datos y el modelo ajustado a estos, se procede a realizar la predicción de los valores desconocidos. Como ya se mencionó, kriging realiza ponderaciones a partir de los valores medidos y aquellos que estén más cerca tienen mayor influencia.

3.2.3 Método de Polígonos de Thiessen

Según Planet botany UWC, este método se puede relacionar con el análisis espacial de datos mediante el análisis de vecindad, donde se genera una red de polígonos de Thiessen a partir de las ubicaciones de los puntos, creando un diagrama de Voronoi.Esto debido a que determina zonas de influencia partiendo de la cercanía de distancias que existe entre los puntos de muestreo.

Los polígonos de Thiessen son un método de interpolación, que realiza construcciones geométricas y se basa en la distancia euclidiana. Se puede definir distancia euclidiana, como un número positivo que indica la separación que tienen dos puntos en el espacio, siguiendo los teoremas de la geometría de Euclides. Matemáticamente, se define la distancia euclidiana d(A,B) entre los puntos A y B, ubicados sobre una recta, como la raíz cuadrada del cuadrado de las diferencias de sus coordenadas x, garantizando que la distancia entre dos puntos sea siempre una cantidad positiva. Y que la distancia entre A y B sea igual a la distancia entre B y A (Lifeder, 2019).

Asimismo, estos polígonos se crean al unir los puntos entre sí y trazando mediatrices. Las intersecciones entre las mediatrices generan una serie de polígonos alrededor de los puntos de muestreo, de manera que el perímetro de cada polígono sea equidistante a los puntos vecinos y se le designa un área de influencia a cada uno de ellos (Andrade, 2011). Además, este modelo permite establecer relaciones matemáticas entre los puntos generando zonas de influencia, y marcando parámetros como las densidades y las equidistancias. Los lados de cada polígono se encuentran a la misma distancia de un punto que del otro.

El método de los polígonos de Thiessen sirve par determinar la precipitación media de un área de interés, y se aplica cuando se conocen las variaciones en los datos de precipitación en cada estación, teniendo además el supuesto de que existe una topografía plana en la zona. (Acasi, 2016).

Analizando la variable de precipitación, este método asigna a cada estación pluviométrica un peso proporcional a su área de influencia (Ingeniería civil, s.f). La influencia se define con los siguientes pasos:

• Las estaciones que se encuentren contiguas se van a conectar mediante líneas rectas, conformando una red irregular de triángulos (TIN) y sin existir líneas interceptadas.

• A cada una de las líneas anteriores se le traza una mediatriz, que se prolonga hasta que se corte con las demás mediatrices. Cabe mencionar, que una mediatriz es una línea recta perpendicular a un segmento de recta y que parte de su punto medio (Acasi, 2016).

• Los puntos de intersección entre estas mediatrices representan los puntos del polígono de Thiessen, cuya superficie constituye el área de influencia de la estación pluviométrica relacionada con dicho polígono.

• Finalmente, el área de cada uno de los polígonos se debe calcular para poder relacionarla y obtener el valor de la precipitación, con la siguiente fórmula:

Fórmula para el cálculo de la precipitación con polígonos de Thiessen. Fuente: Acasi, 2016

Donde,

P, corresponde a la precipitación media del área de interés.

Pi, a la precipitación de cada estación i.

Ai, el área del polígono correspondiente a cada estación i.

4. Presentación de resultados

4.1 Método de interpolación IDW (Distancia inversa ponderada)

Primero, se aplicó el método de interpolación determinístico IDW, tomando la totalidad de los puntos de muestreo y se obtuvo una superficie continua que sigue la tendencia de la distribución de la precipitación en el departamento de Norte de Santander y se puede observar en el mapa 4.

Mapa 4.IDW. Fuente: propia

Para el IDW se uso una potencia de 2.0, es decir, se toma un valor bajo de potencia donde no está disminuyendo la influencia de los puntos más distantes, al momento de hacer el ponderado. Además, una potencia de 2 aumenta la influencia general que tienen los valores conocidos. En este caso, los valores están localizados y no se promedian tanto como una potencia de 1. (GIS geography, 2020)

A continuacuón, se observa nuevamente el mapa 4, pero empleando la librería leaflet, permitiendo una mejor visualización e interpretación de los resultados obtenidos.

Mapa 4 utilizando libreria leaflet

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

Por otro lado, en la figura 5, se registran los resultados de la relación que existe entre los valores predichos y los valores observados o medidos. Además,el análisis de la nube de puntos permite determinar que tan ajustada y adecuada quedó la predicción con el método de IDW, sabiendo que en una interpolación perfecta, la nube de puntos debe seguir la línea negra.

Figura 5. Relación entre valor predicho y valor observado. Fuente: propia

El error cuadrático medio (RMSE) que se obtuvo fue de 2.393775, indicando la relación existente entre los valores observados, obtenidos por los datos CHIRPS, y los valores predichos con el método de interpolación IDW.

Validación cruzada

Se realizó la validación cruzada para el método IDW, tomando un intervalo de confianza del 95 %. Además, se realizó dicha validación con solamente el 20% de los datos, que el programa tomó como 100 puntos, esto con el objetivo de evitar que se vuelva un proceso largo y tedioso.

El resultado obtenido se observa en el mapa 5.

Mapa 5. Fuente:propia

4.2 Método Kriging

Luego, se procede a aplicar el método de interpolación estadístico Kriging.

Semivariograma

Primero, es necesario crear un semivariograma empírico que relacione la semivarianza con la distancia existente entre cada par de puntos. Después, se le ajustó el modelo de semivariograma que mejor explicaba la correlación espacial existente en estos datos.

El semivariograma obtenido se observa en la figura 6.

Figura 6. Semivariograma. Fuente: propia

Los parámetros que se le asignaron al semivariograma fueron:

• Rango: 150000 m

• Nugget: 0.2

• Sill: 3

Además, el semivariograma se ajustó a un modelo lineal (mat).

Posteriormente, al obtener un modelo ajustado adecuadamente a los datos se procede a aplicar la fórmula de kriging, explicada en la sección de Datos y Métodos, y se obtiene una superficie continua que explica de forma apropiada el comportamiento de la variable precipitación en el departamento. La interpolación obtenida se observa en el Mapa 6.

Mapa 6. Kriging.Fuente:propia

El método kriging hace posible generar un mapa de varianza, que ayude a explicar posteriormente el error y la varianza existente entre los datos, para poder conocer el grado de ajuste con que el modelo representa los datos. En el mapa 7 se observan los resultados obtenidos, mostrando intervalos que van desde 1-2 hasta varianzas de 9-10 mm2. Para el mapa de varianza se analizaron solamente el 20% de los datos (100 puntos).

Mapa 7. Mapa de varianza

Mapa 7. Fuente:propia

Asimismo, se genera un mapa sobre el intervalo de confianza, que en este caso, corresponde a 95%, para explicar la incertidumbre de los resultados obtenidos. Para la elaboración del mapa 8, se tomó el 20% de los datos, correspondiente a 100 puntos de muestreo.

Mapa 8. Intervalo de confianza del 95%

Mapa 8. Fuente: propia

4.3 Polígonos de Thiessen

Finalmente, se aplicó el método de Polígonos de Thiessen.

Para realizar la interpolación y obtener polígonos con formas distintas a un cuadradro, se tomaron el 70% de los puntos, que corresponde a 504 datos de precipitación. Lo anterior permite realizar un análisis más adecuado de los datos.

La superficie generada se observa en el mapa 9.

Mapa 9. Fuente: propia

5. Análisis de resultados

5.1 Método IDW

El método de IDW solo tiene en cuenta la distancia entre los puntos, pero no la posición de estos. Es decir, un punto situado a determinada distancia hacia el Norte va a tener la misma influencia que uno situado la misma distancia, pero hacia el Sur. Esto puede causar error en la obtención de los valores, ya que variables climáticas como la precipitación se ven afectadas por factores como la latitud, altitud y temperatura, los cuales varían con la posición (Olaya, 2014). Igualmente, este método basado en distancia no genera valores que se encuentren fuera del rango de valores de los datos que se utilizaron para hacer la interpolación. Eso causa efectos indeseados en caso de que el muestreo de datos no recoja los puntos característicos de la superficie interpolada, y, por tanto, se puede incurrir en sub o sobrestimaciones de la precipitación en determinadas partes del departamento.

RMSE:

EL RMSE (Error cuadrático medio) es una medida de error absoluta que calcula el cuadrado de las desviaciones y mide la cantidad de error que existe entre dos conjuntos de datos, comparando un valor predicho y un valor observado o conocido. (Oracle, s.f) Por lo tanto, el RMSE cuantifica cuán diferente es un conjunto de valores. RMSE va a tomar la diferencia para cada valor de IDW y el valor observado. Se calcula mediante la siguientes fórmula:

Fórmula de RMSE

Donde,

n, corresponde al número de valores observados.

Pi, cada valor predicho i.

Oi, cada valor observado i.

En la gráfica de la figura 5 se observa la relación entre los valores observados y los estimados mediante la interpolación IDW.

Valor predicho: Valor estimado con el método IDW.

Valor observado: Datos de precipitación de CHIRPS.

El valor estimado está representado por la línea negra, mientras que los puntos negros que se observan son los valores observados. El error cuadrático medio estará dado por la distancia que existe entre estos puntos y la línea modelo, y la línea roja permite observar de forma general que tan alejados están estos valores. Siempre se obtendrán valores positivos.Se puede analizar que los datos que se encuentran entre valores de precipitación de 0 y 15 mm aproximadamente, tienen menor error. A medida que la precipitación aumenta, los puntos se encuentran más alejados de la línea recta, aumentando el error en las zonas donde cae mayor precipitación. El promedio de estos errores es de 2,4 mm. Teniendo mayor influencia los errores de los altos valores de precipitación.

Por lo tanto, el valor de RMSE obtenido de 2.393775 indica que existe un error medio de 2.4 mm entre ambos conjuntos de datos, lo cual relacionado con los valores de precipitación, que tienen un rango de 2mm a 38 mm, es un error bajo y asegura que la técnica de IDW está proporcionando una interpolación que se ajusta adecuadamente al conjunto de datos. Además, cuanto más pequeño es el valor de RMSE, más cercanos son los valores predichos y los observados, existe mayor relación y precisión del método IDW (GIS Geography, 2020)

En cuanto a la validación cruzada, representada en el mapa 5 , esta consiste en tomar todos los datos originales (721 puntos) y crear a partir de ellos dos conjuntos separados: un primer conjunto de entrenamiento (y prueba), y un segundo conjunto de validación. A cada uno de estos puntos se le debe realizar la interpolación, por tanto, puede resultar un proceso largo y tedioso. En este caso, se toma el 20% de los datos, que el programa tomó como 100 puntos, indicando que para realizar la validación cruzada se realizaron aproximadamente 100 interpolaciones.

Una vez finalizadas las interpolaciones, se calcula la precisión y el error para el modelo producido (Delgado, 2018). La validación cruzada permite saber que tanto son los errores, y cuál es el intervalo de confianza. En este caso, se toma un intervalo de confianza del 95% y se observa que, en la mayoría del departamento, el error varía entre 0.0 a 0.2 mm, pero existen zonas donde el error es de 0.4 a 0.6 mm o incluso de 0.8 a 1.0 mm. Es posible deducir que el error que produce el método de IDW es muy bajo y se ajusta a la tendencia que tienen los datos de precipitación en Norte de Santander.

5.2 Método de kriging

El modelo seleccionado para el ajuste del variograma fue el lineal (mat), y de este se puede analizar que la variabilidad espacial aumenta linealmente con la distancia. Además, al ser lineal, el usuario debe seleccionar un valor de sill y rango, de forma arbitraria, ya que no se representan directamente en el gráfico (Gis Geography, 2019).El rango elegido de 150000 m no indica con certeza en este modelo, que hasta ese valor existe correlación espacial.

Además, en la figura 6 se puede observar que los datos quedaron bien ajustados al modelo seleccionado ya que la línea recta que se produce, pasa por la mayoría de los puntos, o cerca a ellos, con diferencias muy bajas.

La representación del semivariograma para la precipitación en Norte de Santander, se observa en la figura 5. Como se puede apreciar, esta figura no representa todos los intervalos de distancia posibles, ya que solo llega a 110000 m, siendo el departamento de 145000 x 145000 m aproximadamente (área de 21000 km2). Este aspecto es importante, ya que como regla general sólo se representa en el semivariograma distancias aproximadas a la mitad de la dimensión del área de interés. La razón se puede encontrar en que el número de pares a las distancias mayores desciende y la semivarianza puede llegar a mostrar un comportamiento erróneo (Gallardo, 2006).

Asimismo, el variograma nos da bastante información del comportamiento espacial de la precipitación, y permite deducir que los valores de la muestra separados por intervalos de distancia entre 0 y 20000 m son más parecidos que aquellos separados por 60000 m o 80000 m. A su vez, la semivarianza aumenta a medida que aumenta la distancia que separa las muestras hasta que se llega a un intervalo de distancia aproximada de 100000 m, donde la semivarianza alcanza un valor máximo. Lo anterior, nos permite confirmar la Primera Ley de la geografía en que se basa kriging, la cual afirma que las cosas que están cerca en el espacio, se parecen mucho.

En cuanto al valor de nugget, teóricamente el valor de la función debe ser cero en el origen, pero en este caso se le asigna un valor de 0.2, que corresponde a la intersección en el eje Y. Este error nugget o también llamada, discontinuidad en el origen, puede relacionarse con errores de medición o con fuentes espaciales de variación en distancias que son menores a los valores del muestreo. Por su parte, los errores de medición son inherentes a los dispositivos de medición, que, en este caso, corresponden a las estaciones pluviométricas y errores relacionados con el análisis de las imágenes de satélite. Asimismo, variables naturales como la precipitación, varían en el espacio en un rango de valores determinado, pero la variación a escalas muy pequeñas, que tiendan a cero, ocasionan la aparición de errores en el modelo.

Por otra parte, modelos geoestadísticos como kriging, no sólo tienen la capacidad de producir una superficie de predicción sino que proporcionan alguna medida de la certeza o exactitud de las predicciones (ArcGis, 2016). Esto se observa en el mapa 7 y 8, que muestran la varianza y el intervalo de confianza del 95%, respectivamente.

Por su parte, se puede definir la varianza como una medida de dispersión que representa la variabilidad de una serie de datos respecto a su media. La varianza corresponde a la desviación estándar elevada al cuadrado. En el mapa de varianza (mapa 7) se puede analizar que los valores de varianza obtenidos no son muy altos (1 a 2 mm2 ) excepto en los bordes, donde se observan valores mayores que varían de (3 a 4 mm2 ) hasta 8 a 9 mm2 en el extremo sur-oriental. Se podría afirmar que en ese extremo hay mayor error en los datos porque se encuentran más dispersos con respecto a la media.

Por otra parte, si se saca la raíz cuadrada, por ejemplo en los datos que tienen varianza de 2 mm2 ,se obtiene la desviación estándar que corresponde a 1,4 mm, y esto quiere decir que, en media, la diferencia entre los valores interpolados y los reales será de aproximadamente 1,4 mm. Así, valores bajos de varianza representan un buen indicativo de que el modelo quedó ajustado adecuadamente a los datos.

En cuanto al intervalo de confianza del 95%; cuando se realiza estimación por intervalos de confianza se está analizando tanto la etimación puntual de la variable como un intervalo que mide el margen de error de esta estimación. De cada dato muestreado es posible obtenerse un intervalo de confianza diferente. Un intervalo de confianza del 95% nos indica que dentro del rango dado se encuentra el valor real de la variable con 95% de certeza (Quintela, 2019). También, el intervalo de confianza describe la variabilidad entre la medida obtenida en la interpolación y la medida real de la precipitación (el valor real) (Candia y Caiozzi, 2005). En el mapa 8, se observa que los valores oscilan entre 2.0 y 3.0 mm aproximadamente en todo el departamento, sin embargo, existen zonas en los extremos donde las variaciones son más altas alcanzando valores de 5.0 mm. Se determina que estos intervalos de confianza muestran una variabilidad baja y dan certeza de los datos, afirmando que el error producido por Kriging es muy bajo.

5.3 Polígonos de Thiessen

En un primer intento, se habían obtenido polígonos de Thiessen con forma cuadrada, sin embargo, este resultado puede que no represente adecuadamente la forma como se comporta la precipitación en el departamento, ya que al estar las estaciones distribuidas regularmente, los polígonos que forman las intersecciones de las mediatrices, son todos iguales. Por tanto, se realizó un segundo intento donde se toman el 70% de los datos, que corresponde al mapa 9, donde se puede observar que se tienen polígonos con diversas formas, variando desde cuadrados, rectángulos hasta trapecios, hexágonos y formas irregulares.

Los polígonos irregulares que se forman, se debe a que en este método se está delimitando la región correspondiente a cada estación pluviométrica y se están uniendo las estaciones cercanas con segmentos de recta. Además, en este mapa no se emplea la totalidad de los puntos para poder obtener polígonos distintos y observar mejor el comportamiento de la variable.

En este caso, se puede afirmar que cada polígono obtenido está asociado a cada dato del muestreo, ya sea de estación o de imagen de satélite, y los valores que se predicen se relacionan espacialmente con cada estación cercana. Los polígonos de mayor área van a tener mayor influencia en la predicción de los valores desconocidos.

Por otra parte, el método de Thiessen se puede emplear para estimar o calcular la precipitación de un área de interés, principalmente cuando existe una distribución no uniforme de las estaciones pluviométricas y no considera las influencias orográficas. En cuanto a la precipitación en Norte de Santander, se puede observar que los puntos (que corresponden tanto a datos obtenidos en estaciones como con imágenes de satélite) están distribuidos de forma regular, y por tanto, este método no es el más indicado para analizar el comportamiento espacial de la precipitación en el departamento. Además, es necesario tener en cuenta que el relieve en todo el departamento no es uniforme, y por tanto, no se cumple con el condicionante de que la topografía del área de estudio sea plana, y no se están reflejando los efectos altitudinales.

Sin embargo, en el mapa 9 se observa que la distribución de la lluvia es muy similar con respecto a las superficies obtenidas con el método de IDW y kriging.

También, es importante considerar que en los puntos más externos, aquellos que delimitan el departamento, se podrán encontrar errores en los valores interpolados, ya que los polígonos no cuentan con algunos puntos cercanos que influirán en la forma del polígono obtenido.

5.4 Comparación de métodos.

La tarea más dificil al momento de crear una superficie continua por medio de la interpolación, es elegir la técnica más adecuada y que se ajuste mejor a la distribución de los puntos del conjunto de datos que se tienen. Además, es de gran importancia realizar una investigación del comportamiento espacial de la variable, antes de seleccionar el mejor método de interrpolación. En el mapa 3, se puede observar que los datos de precipitación que se tomaron para la interpolación, se distribuyen de forma regular sobre el departamento.

Al comparar los tres métodos de interpolación aplicados, se puede observar que de forma general los tres producen superficies continuas similares que parecen capturar la variación que tienen los datos de muestra usados. Los 3 se comportan razonablemente bien de manera visual si se compara con la distribución de los datos de precipitación en el mapa 1. Se consideran datos con sentido ya que siguen la tendencia natural de la precipitación medida en el departamento.

Sin embargo, se puede afirmar que kriging es más preciso que los demás, debido a que emplea un procedimiento geoestadístico avanzado donde, su algoritmo de cálculo no solo relaciona la distancia entre los puntos y las ubicaciones circundantes, como en el caso de IDW, sino que realiza un trabajo más arduo debido a que emplea más parámetros y tiene en cuenta la naturaleza y correlación espacial de los datos, obteniendo a su vez, resultados con mejor resolución. Se puede observar en el mapa 6, que al usar kriging, varía un poco las zonas que tienen precipitación entre 0 y 5 mm, mientras que las demás regiones son más parecidas.

Adicionalmente, opino que Kriging puede resultar en el método más adecuado para emplear en proyectos posteriores, ya que si se compara la varianza obtenida en kriging, que es aproximadamente 2mm2 (desviación estándar de 1.4 mm), con el error medio cuadrático obtenido en IDW, que corresponde a 2,4 mm; el error o dispersión de los datos obtenidos en Kriging es menor, existiendo mayor relación y cercanía entre los valores interpolados y la tendencia natural de la precipitación.

Finalmente, de los tres métodos se puede observar que el municipio con mayor precipitación es La Esperanza con 35 a 40 mm aproximadamente, y la menor precipitación se registra en Toledo, Tibú y Puerto Santander, con valores entre 0 y 5 mm. Asimismo, comparando los valores obtenidos, con los de la figura 1, se observan datos relativamente bajos de precipitación, teniendo en cuenta que el mes de abril se encuentra entre el primer período de mayores lluvias del departamento, con valores cercanos a los 200 mm. Esto se puede deber a que los últimos 5 días del mes tuvieron bajas precipitaciones, y además, los datos que se tomaron para realizar el climograma de la figura 1, corresponden a promedios climatológicos de los años 1980-2010 y actualmente, diez años después se está viviendo una problemática de calentamiento global y cambio climático, que puede influenciar en estos valores.

6. Conclusiones

• La interpolación es una técnica viable para predecir el valor desconocido de un punto geográfico, como la precipitación, ya que muestrear todas las ubicaciones del área de estudio para medir esta variable, puede resultar difícil y costoso.

• El método de polígonos de Thiessen tiene su principal aplicación en las áreas de la hidrología y la climatología, ya que permite calcular la precipitación de un área de interés, analizando la distribución espacial de las estaciones pluviométricas.

• El uso de variogramas en el método de kriging, es una herramienta muy importante para estudiar la relación espacial de conjuntos de datos relacionados con variables agronómicas y ambientales.

• IDW tienen un error cuadrático medio de 2,4. Un valor que indica que existe adecuada relación entre los valores medidos en las estaciones y las imágenes de satélite, y los valores predichos por el método IDW.

• El método de kriging se puede considerar más preciso que el IDW, ya que las ponderaciones que se realizan en el cálculo provienen de un semivariograma que se desarrolló observando la naturaleza espacial de los datos, y no solo de la distancia entre ellos. Además, tiene una varianza de los datos, menor que el error de IDW.

7. Referencias

Acasi, D. 2016. Calcular la precipitación media con polígonos de Thiessen. Recuperado de: https://www.youtube.com/watch?v=Q9OA6HDsdsM

Andrade, E. 2011.Polígonos de Thiessen. Recuperado de: https://civilgeeks.com/2011/09/24/poligonos-de-thiessen/

ArcGis. 2016. Comparar métodos de interpolación. Recuperado de: https://desktop.arcgis.com/es/arcmap/10.3/tools/spatial-analyst-toolbox/comparing-interpolation-methods.htm

ArcGis. 2019. Deterministic methods for spatial interpolation. Recuperado de: https://desktop.arcgis.com/es/arcmap/latest/extensions/geostatistical-analyst/deterministic-methods-for-spatial-interpolation.htm

Arias, D; Cerne, B; Gulizia, C; Zazulie, N. s.f. Primer encuentro. Recuperado de: http://difusion.df.uba.ar/ConectarIgualdad/presentacionCS%20de%20la%20ATmosfera.pdf

Candia, R y Caiozzi, G. 2005. Intervalos de confianza. Recuperado de: https://scielo.conicyt.cl/scielo.php?script=sci_arttext&pid=S0034-98872005000900017

Delgado, R. 2018.Introducción a la Validación Cruzada (k-fold Cross Validation) en R. Recuperado de: http://rstudio-pubs-static.s3.amazonaws.com/405322_6d94d05e54b24ba99438f49a6f8662a9.html

Gallardo, A. 2006. Geoestadística. Ecosistemas. 2006/3. Recuperado de: https://fjferrer.webs.ull.es/Bibliog/Biblio/Geoestadistica.pdf

Ingeniería civil. s.f. Métodos para el Cálculo de la Precipitación Media en una Cuenca. Recuperado de:http://ingenieriacivil.tutorialesaldia.com/metodos-para-el-calculo-de-la-precipitacion-media-en-una-cuenca/

Lifeder, 2019. Distancia euclidiana: concepto, fórmula, cálculo, ejemplo. Recuperado de: https://www.lifeder.com/distancia-euclidiana/

Madrid, A. & Ortiz, L.2018. Análisis espacial. Recuperado de: http://www.bdigital.unal.edu.co/1239/3/02CAPI01.pdf

Plataforma ArcGis. -. Resolver problemas con análisis espacial. Recuperado de: https://www.sigsa.info/productos/esri/plataforma-arcgis/analisis-espacial

Planet botany UWC. Introducción al análisis espacial. Recuperado de: http://planet.botany.uwc.ac.za/nisl/GIS/spatial/chap_1_01.htmv

UC Santa Bárbara.-. Climate Hazards Center. UC Santa Bárbara. Recuperado de: https://www.chc.ucsb.edu/

UC Santa Bárbara.-. CHIRPS: Rainfall Estimates from Rain Gauge and Satellite Observations. Recuperado de: https://chc.ucsb.edu/data/chirps

Olaya, V. 2014. Sistemas de información geográfica. Recuperado de: https://www.icog.es/TyT/files/Libro_SIG.pdf

Gis Geography, 2019. Semi-Variogram: Nugget, Range and Sill. Recuperado de: https://gisgeography.com/semi-variogram-nugget-range-sill/

GIS Geography, 2020. How to Calculate Root Mean Square Error (RMSE) in Excel. Recuperado de: https://gisgeography.com/root-mean-square-error-rmse-gis/

DANE. 2018. Resultados del Censo Nacional de Población y Vivienda 2018. Recuperado de: https://www.dane.gov.co/files/censo2018/informacion-tecnica/presentaciones-territorio/190806-CNPV-presentacion-Norte-de-Santander.pdf

Departamento del Norte de Santander. 2003. Norte de Santander. Recuperado de: https://repository.urosario.edu.co/bitstream/handle/10336/8681/1020763251-13.pdf?sequence=13

Contraloría General de la República.2015. PERFIL DE LA GESTION FISCAL DEL DEPARTAMENTO DE NORTE DE SANTANDER. Recuperado de: https://www.contraloria.gov.co/documents/487635/520928/Perfil+gesti%C3%B3n+NORTE+DE+SANTANDER.pdf/135607d2-4743-4488-9b83-8b2b5c1053a2?version=1.0

Gobernación del departamento de Norte de Santander. 2018. Información General Norte de Santander. Recuperado de: http://www.nortedesantander.gov.co/Gobernaci%C3%B3n/Nuestro-Departamento/Informaci%C3%B3n-General-Norte-de-Santander#:~:text=Norte%20de%20Santander%20es%20uno,sobre%20la%20frontera%20con%20Venezuela.&text=Limita%20al%20norte%20y%20al,oeste%20con%20Santander%20y%20Cesar.

Quintela, A. 2019. 7.4 Intervalos de confianza. Recuperado de: https://bookdown.org/aquintela/EBE/intervalos-de-confianza.html

ANEXOS

Datos CHIRPS

ANEXO 1. Métodos de interpolación

ANEXO 2. Interpolación espacial