Aplicación del Concepto de Interpolación Espacial al Departamento de Huila
Sergio García
10.12.2020
Introducción
Este cuaderno muestra la aplicación de tres técnicas de interpolación espacial para el departamento de Huila. Los conceptos aplicados en este cuaderno hacen parte del curso geomética básica de la Universidad Nacional de Colombia. La relización de este informe, estuvo apoyada en la guía dispuesta por el profesor Iván Lizarazo.
La recopilación de datos de interés agronómico en una zona particular es fundamental para el desarrollo de investigaciones o trabajos que busquen extraer información de forma rápida y con aproximaciones a la realidad. Por eso es tan importante aplicar un análisis espacial cuando tenemos pocos datos y queremos extrapolar a partir de esa información, otra que sea cercana y modele información desconocida. La forma de estimar esos datos desconocidos será utilizando cualquiera de los métodos de interpolación.
Estas interpolaciones utilizan un GIS (celdas) basado en ráster y hay una variedad de técnicas para derivar una predicción para cada ubicación. Cada una de estas técnicas se basa en diferentes supuestos (por ejemplo, si los datos se ajustan a una distribución normal) (Introducción al análisis espacial).
Si se toman puntos de información que se alejan unos de otros el paso a seguir sería ralizar una interpolación y no una extrapolación.
La interpolación espacial es utilizada para estimar o predecir valores desconocidos en una determinada zona a través de valores recopilados conocidos en un lugar específico. Estos valores sirven como recurso a la hora de obtener información y organizarla de forma clara de feñómenos naturales o características de una región especifica: elevación, lluvia, características físicas y químicas del suelo, entre otros.
Descripción de la zona de estudio
El departamento del Huila identificado con el código Dane 041, está localizado al suroccidente del país entre los 3º55’12” y 1º30’04” de latitud norte (entre el nacimiento del Rio Riachón, municipio de Colombia y el pico de la Fragua, municipio de Acevedo), y los 74º25’24” y 76º35’16” de longitud al oeste del meridiano de Greenwich ( entre el Alto de Las Oseras, municipio de Colombia y el páramo de Las Papas, municipio de San Agustín.) (Gobernación del Huila, 1969).
Según datos tomados del mapa físico-político de Colombia elaborado por el instituto Geográfico Agustín Codazzi, la superficie del Departamento es de 19.900 Km2 que representa tan solo un 1.8% de la superficie total del país. Comparada con la superficie de los demás departamentos, ocupa el 170 lugar, superando a Caldas, Atlántico, Quindío, Risaralda y Sucre (Gobernación del Huila, 1969).
Al norte limita con los departamentos de Cundinamarca y el Tolima al sur con los de Cauca y Caquetá, al oriente con los departamentos de Meta y Caquetá, y hacia el Occidente con los de Cauca y Tolima (Gobernación del Huila, 1969).
mapaSegún el Instituto Geográfico Agustín Codazzi, las tierras arables aptas para la agricultura y ganadería solo suman el 15,9%; es decir que se les ha ido la mano y han afectado la sostenibilidad de los recursos naturales.
Ante este panorama, el 83% de los suelos huilenses, catalogados por el IGAC como terrenos para el uso forestal y para la conservación ambiental, se ha reducido a más de la mitad: un alarmante 32,7%.
Los municipios en donde debería estar focalizada la producción del Huila son Palermo, Rivera, Campo Alegre, Hobo, Paicol, Gigante, Garzón, Altamira, Tarqui, La Plata, Pilalito e Isnos; pero la realidad es que sus 37 territorios están invadidos por este tipo de desarrollo (IGAC, 2016).
La economía del departamento del Huila se basa principalmente en la producción agrícola y ganadera, la explotación petrolera y el comercio. La agricultura se ha desarrollado y tecnificado en los últimos años y sus principales cultivos son café, algodón, arroz riego, fríjol, maíz tecnificado, maíz tradicional, sorgo, cacao, caña panelera, plátano, yuca, iraca y tabaco. Los campos de petróleo se encuentran en el norte del departamento y para la distribución de gas está conectado por el gasoducto Vasconia – Neiva en donde las reservas representan el 1.2% del total nacional (Gobernación del Huila, 2017).
En la cordillera Central se extrae plata y oro, este último se explota en 13 municipios. Otros minerales no preciosos que se explotan son cuarzo, calcita, marmol y azufre. La producción artesanal es muy laboriosa, especialmente la de cerámica y sombreros. La industria fabril está poco desarrollada, no obstante, en Neiva se han instalado fábricas de productos alimenticios, bebidas, jabones, cigarros y licores. Los centros de gran actividad comercial son Neiva, Garzón y en menor escala Gigante, La Plata y Campoalegre (Gobernación del Huila, 2017).
La fertilidad de las tierras del valle alto del Magdalena, es propicia para el cultivo intensivo de arroz (Gobernación del Huila, 2017).
Clima
La variación de temperatura en el departamento del Huila oscila entre 28º C experimentados en su parte más baja, especialmente en el valle de Neiva, extendiéndose hasta el extremo norte, atravesando toda la región semiárida de Yararaca (Tatacoa). En esta región se encuentran todos los climas y una gran variedad de suelos que facilitan la diversidad y extensión de la producción agrícola y ganadera; predomina el clima templado, con una temperatura media de 24 ºC; como puntos extremos están las cumbres montañosas del Nevado del Huila, que forma parte del Parque Nacional Natural que lleva su mismo nombre, donde la temperatura permanece bajo 0 ºC y las regiones cálidas de los valles de Neiva, Aipe y Villavieja, donde se encuentra el imponente desierto de La Tatacoa, con 35º C (Gobernación del Huila, 2017).
Descripción de datos y métodos
Para realizar el análisis de datos, es neceario descargar los datos de precipitación. Descargamos los datos de CHIRPS, cuyas siglas significan: Grupo de Riegos Climáticos de Precipitación Infrarroja.
Para este informe se utilizó un grupo de tres métodos, selecccionados dentro de varios métodos que existen para llevar a cobo el proceso de interpolación. Los métodos elegidos fueron Poligonos de Thiessen, IDW y Kriging.
Polígonos de Thiessen
Este método de interpolación es uno de los más simples y es utilizado cuando el conjunto de datos son cualitativos. Este método se forma a partir de la unión de puntos por medio de segmentos. La unión de todos estos segmentos describen una serie de polígonos.
Ponderación de Distancia Inversa (IDW)
Por medio de un conjunto de puntos dispersos, IDW, fija valores a cada celda en específico. “Cuando se le dan valores conocidos, la interpolación estima valores desconocidos” (GIS Geography). Para calcular el valor de una celda en especial que es deconocida y que está cercana a varios puntos conocidos se utiliza la siguiente fórmula:
Kriging
“Es un método geoestadístico de interpolación que ha probado ser útil y popular en muchos campos (Burgess y Webster, 1980). En la actualidad FUNDECOR utiliza este método de interpolación de modelos de elevación digital (mapas de curvas de nivel) para la planificación del aprovechamiento forestal de los bosques naturales. Dicho método provee, a partir de una muestra de puntos, ya sean regular o irregularmente distribuidos, valores estimados de aquellos sitios donde no hay información, sin sesgo y con una varianza mínima conocida” (FAO).
Presentación de resultados
Conjunto de puntos - Departamento de HUila
library(rgdal)## Loading required package: sp## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/Sergio/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/Sergio/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.library(raster)
library(sf)library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.3## -- Attaching packages ----------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.3 v dplyr 1.0.2
## v tidyr 1.1.2 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts -------------------------------------------------------- tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()library(tmap)## Warning: package 'tmap' was built under R version 4.0.3library(gstat)## Warning: package 'gstat' was built under R version 4.0.3library(sp)precip <- raster("C:/Users/Sergio/Downloads/GB/Trabajo Final/chirps-v2.0.2020.09.6.tif")Es conveniente conocer la información contenida en los archivos que manejamos. Así, para precip se requirió información con la instrucción:
precip## class : RasterLayer
## dimensions : 2000, 7200, 14400000 (nrow, ncol, ncell)
## resolution : 0.05, 0.05 (x, y)
## extent : -180, 180, -50, 50 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : C:/Users/Sergio/Downloads/GB/Trabajo Final/chirps-v2.0.2020.09.6.tif
## names : chirps.v2.0.2020.09.6(hui <- shapefile("C:/Users/Sergio/Downloads/GB/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")) ## class : SpatialPolygonsDataFrame
## features : 37
## extent : -76.62466, -74.41303, 1.552125, 3.843208 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 9
## names : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR, MPIO_CRSLC, MPIO_NAREA, MPIO_NANO, DPTO_CNMBR, Shape_Leng, Shape_Area
## min values : 41, 41001, ÃQUIRA, 1538, 80.44174325, 2017, HUILA, 0.517464409731, 0.00653203068962
## max values : 41, 41885, YAGUARÃ, Ordenanza 6 del 26 de Noviembre de 1965, 1584.96991931, 2017, HUILA, 3.07318677018, 0.128978175097precip.crop <- raster::crop(precip, extent(hui))Se enmascaró la capa raster:
precip.mask <- mask(x = precip.crop, mask = hui)Se verificó la que la instrucción se haya realizado correctamente:
precip.mask## class : RasterLayer
## dimensions : 46, 44, 2024 (nrow, ncol, ncell)
## resolution : 0.05, 0.05 (x, y)
## extent : -76.6, -74.4, 1.549999, 3.849999 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : chirps.v2.0.2020.09.6
## values : 2.736892, 48.57009 (min, max)Se trazó la capa ráster enmascarada:
plot(precip.mask, main= "CHIRPS lluvia en Huila desde el 26.09. al 30.09 en 2020 [mm]")
plot(hui, add=TRUE)
Existen otros tipos de visualizaciones de mapas. El siguiente mapa se hizo con leaflet.
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precip.mask),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(precip.mask, colors = pal, opacity = 0.6) %>%
addLegend(pal = pal, values = values(precip.mask),
title = "CHIRPS lluvia Huila desde el 26.09 al 30.09 en 2020 [mm]")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definitionLa conversión de ráster en puntos se puede realizar utilizando la función rasterToPoints de la biblioteca de ráster:
precip.points <- rasterToPoints(precip.mask, spatial = TRUE)precip.points## class : SpatialPointsDataFrame
## features : 596
## extent : -76.575, -74.425, 1.574999, 3.824999 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 1
## names : chirps.v2.0.2020.09.6
## min values : 2.73689198493958
## max values : 48.5700912475586names(precip.points) <- "rain"precip.points## class : SpatialPointsDataFrame
## features : 596
## extent : -76.575, -74.425, 1.574999, 3.824999 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 1
## names : rain
## min values : 2.73689198493958
## max values : 48.5700912475586str(precip.points)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 596 obs. of 1 variable:
## .. ..$ rain: num [1:596] 9.51 12.08 13.33 13.6 13.23 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:596, 1:2] -74.5 -74.5 -74.5 -74.6 -74.5 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "x" "y"
## ..@ bbox : num [1:2, 1:2] -76.57 1.57 -74.42 3.82
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"Cuando se realizaron la conversión de puntos, se procedió a disponer en un mapa las configuraciones efectuadas anteriormente:
plot(precip.mask, main= "CHIRPS lluvia desde el 26 al 30.09.2020 [mm]")
plot(hui, add=TRUE)
points(precip.points$x, precip.points$y, col = "red", cex = .6)
Se utilizó la instrucción rgadal para guardar el archivo con el fin de evitar que algún error adelante arruine lo hecho hasta este punto.
geojsonio::geojson_write(precip.points, file = "C:/Users/Sergio/Downloads/GB/Trabajo Final/chirps/ppoints.geojson")## Success! File is at C:/Users/Sergio/Downloads/GB/Trabajo Final/chirps/ppoints.geojson## <geojson-file>
## Path: C:/Users/Sergio/Downloads/GB/Trabajo Final/chirps/ppoints.geojson
## From class: SpatialPointsDataFrameLuego, se leyeron los puntos de precipitación.
precip.points <- geojsonio::geojson_read("C:/Users/Sergio/Downloads/GB/Trabajo Final/chirps/ppoints.geojson", what="sp")precip.points## class : SpatialPointsDataFrame
## features : 596
## extent : -76.575, -74.425, 1.574999, 3.824999 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 1
## names : rain
## min values : 2.73689198493958
## max values : 48.5700912475586En caso de que se pierda la conexión, podemos volver a leer el shapefile con el área de interés:
(hui <- shapefile("C:/Users/Sergio/Downloads/GB/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))## class : SpatialPolygonsDataFrame
## features : 37
## extent : -76.62466, -74.41303, 1.552125, 3.843208 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 9
## names : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR, MPIO_CRSLC, MPIO_NAREA, MPIO_NANO, DPTO_CNMBR, Shape_Leng, Shape_Area
## min values : 41, 41001, ÃQUIRA, 1538, 80.44174325, 2017, HUILA, 0.517464409731, 0.00653203068962
## max values : 41, 41885, YAGUARÃ, Ordenanza 6 del 26 de Noviembre de 1965, 1584.96991931, 2017, HUILA, 3.07318677018, 0.128978175097Se convirtió en una característica espacial:
nei_sf <- sf::st_as_sf(hui)Se disolvieron los límites internos:
(border_sf <-
nei_sf %>%
summarise(area = sum(MPIO_NAREA)))Se convirtió la característica espacial en un marco de datos espacial:
(border <- as(border_sf, 'Spatial'))## class : SpatialPolygonsDataFrame
## features : 1
## extent : -76.62466, -74.41303, 1.552125, 3.843208 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 1
## names : area
## value : 18141.65581079Otra conversión:
(nei.sf <- st_as_sf(hui) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))Se ejecutó la instrucción st_intersection. Se utilizó la biblioteca sf:
p.sf <- st_as_sf(precip.points)(precip.sf = st_intersection(nei.sf, p.sf))## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesReproyección:
p.sf.magna <- st_transform(precip.sf, crs=3116)nei.sf.magna <- st_transform(nei.sf, crs=3116)Una conversión:
(precip2 <- as(p.sf.magna, 'Spatial'))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
## but +towgs84= values preserved## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum Marco_Geocentrico_Nacional_de_Referencia in CRS definition## class : SpatialPointsDataFrame
## features : 596
## extent : 722040, 961398.5, 666034.3, 914733.3 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 3
## names : MUNIC, CODIGO, rain
## min values : ÃQUIRA, 41001, 2.73689198493958
## max values : YAGUARÃ, 41885, 48.5700912475586Se escribió de nuevo el conjunto de datos intermedio(para salvarlo):
shapefile(precip2, filename='C:/Users/Sergio/Downloads/GB/Trabajo Final/precip2.shp', overwrite=TRUE)precip2$rainfall <- round(precip2$rain, 1)Se leyeron los datos configurados:
precip2## class : SpatialPointsDataFrame
## features : 596
## extent : 722040, 961398.5, 666034.3, 914733.3 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 4
## names : MUNIC, CODIGO, rain, rainfall
## min values : ÃQUIRA, 41001, 2.73689198493958, 2.7
## max values : YAGUARÃ, 41885, 48.5700912475586, 48.6(nei2 <- as(nei.sf.magna, 'Spatial'))## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
## but +towgs84= values preserved## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum Marco_Geocentrico_Nacional_de_Referencia in CRS definition## class : SpatialPolygonsDataFrame
## features : 37
## extent : 716525, 962727.8, 663502.8, 916747.4 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 2
## names : MUNIC, CODIGO
## min values : ÃQUIRA, 41001
## max values : YAGUARÃ, 41885Al igual que antes, se guardó lo anteriormente hecho:
shapefile(nei2, filename='C:/Users/Sergio/Downloads/GB/Trabajo Final/nei2.shp', overwrite=TRUE)Se aseguró de que las extensiones coincidieran:
precip2@bbox <- nei2@bboxLuego, se dió paso a realizar el mapa con tmap:
tm_shape(nei2) + tm_polygons() +
tm_shape(precip2) +
tm_dots(col="rainfall", palette = "YlGnBu", midpoint = 0.1,
title="Precipitación muestreada \ n (en mm)", size=0.5) +
tm_text("rainfall", just="center", xmod=.6, size = 0.5) +
tm_legend(legend.outside=TRUE)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
Polígonos de Thiessen
Este tipo de interpolación se puede llevar a cabo con la función Spatstat.
library(spatstat)## Warning: package 'spatstat' was built under R version 4.0.3## Loading required package: spatstat.data## Warning: package 'spatstat.data' was built under R version 4.0.3## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## The following object is masked from 'package:raster':
##
## getData## Loading required package: rpart## Registered S3 method overwritten by 'spatstat':
## method from
## print.boxx cli##
## spatstat 1.64-1 (nickname: 'Help you I can, yes!')
## For an introduction to spatstat, type 'beginner'##
## Note: spatstat version 1.64-1 is out of date by more than 7 months; we recommend upgrading to the latest version.##
## Attaching package: 'spatstat'## The following object is masked from 'package:gstat':
##
## idw## The following objects are masked from 'package:raster':
##
## area, rotate, shift# Crea una superficie teselada
th <- as(dirichlet(as.ppp(precip2)), "SpatialPolygons")
# La función dirichlet no transfiere información de proyección
# requiriendo que esta información se agregue manualmente
crs(th) <- crs(precip2)
crs(nei2) <- crs(precip2)crs(th)## CRS arguments:
## +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defscrs(precip2)## CRS arguments:
## +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs# Mapear los datos
tm_shape(th.clp) +
tm_polygons(col="rainfall", palette="RdBu", midpoint=3.0,
title="Polígonos de Thiessen \ nPrecipitación prevista \ n (en mm)") +
tm_legend(legend.outside=TRUE)## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
Ponderación de Distancia Inversa (IDW)
# Cree una cuadrícula vacía donde n es el número total de celdas
grd <- as.data.frame(spsample(precip2, "regular", n=100000))
# Necesita averiguar cuál es el tamaño esperado del grd de salida
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE # Crear objeto SpatialPixel
fullgrid(grd) <- TRUE # Crear objeto SpatialGrid
# Agregue la información de proyección de P a la cuadrícula vacía
proj4string(grd) <- proj4string(precip2)## Warning in proj4string(precip2): CRS object has comment, which is lost in output## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
## but +towgs84= values preserved# Interpolar las celdas de la cuadrícula usando un valor de potencia de 2 (idp = 2.0)
P.idw <- gstat::idw(rainfall ~ 1, precip2, newdata=grd, idp=2.0)## Warning in proj4string(d$data): CRS object has comment, which is lost in output## Warning in proj4string(newdata): CRS object has comment, which is lost in output## [inverse distance weighted interpolation]## Warning in proj4string(newdata): CRS object has comment, which is lost in output## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
## but +towgs84= values preserved# Convertir a objeto ráster y luego recortar a HUI
r <- raster(P.idw)r## class : RasterLayer
## dimensions : 321, 312, 100152 (nrow, ncol, ncell)
## resolution : 789.6172, 789.6172 (x, y)
## extent : 716320.5, 962681.1, 663489.2, 916956.3 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : var1.pred
## values : 2.800302, 48.10019 (min, max)nei2## class : SpatialPolygonsDataFrame
## features : 37
## extent : 716525, 962727.8, 663502.8, 916747.4 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 2
## names : MUNIC, CODIGO
## min values : ÃQUIRA, 41001
## max values : YAGUARÃ, 41885r.m <- raster::mask(r, nei2)r.m ## class : RasterLayer
## dimensions : 321, 312, 100152 (nrow, ncol, ncell)
## resolution : 789.6172, 789.6172 (x, y)
## extent : 716320.5, 962681.1, 663489.2, 916956.3 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : var1.pred
## values : 2.800302, 48.10019 (min, max)# Mapa
tm_shape(r.m) +
tm_raster(n=10,palette = "RdBu", auto.palette.mapping = FALSE,
title="Distancia inversa ponderada \ nPrecipitación prevista \ n (en mm)") +
tm_shape(precip2) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)## Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
## for numeric data and stretch.palette for categorical data to control the palette
## mapping.
Se realizó una visualización mejor de la interpolación hecha:
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precip.mask),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
addLegend(pal = pal, values = values(r.m),
title = "IDW precipitación interpolada en Huila \ n De 26.09 a 30.09 en 2020 [mm]")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definitionAjuste de la interpolación
Se hizo un ajuste de interpolación de acuerdo con lo expuesto por el Profsor Iván Lizarazo: “La elección de la función de potencia puede ser subjetiva. Para ajustar la elección del parámetro de potencia, puede realizar una rutina de validación de dejar uno fuera para medir el error en los valores interpolados”(Lizarazo, 2020).
precip2## class : SpatialPointsDataFrame
## features : 596
## extent : 716525, 962727.8, 663502.8, 916747.4 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 4
## names : MUNIC, CODIGO, rain, rainfall
## min values : ÃQUIRA, 41001, 2.73689198493958, 2.7
## max values : YAGUARÃ, 41885, 48.5700912475586, 48.6P <- precip2IDW.out <- vector(length = length(P))Se efectuó la siguiente instrucción y al ejecutar el código generó un serie de respuestas muy largas y para efectos de presentación y estilo del informe, se prefirió mantenerlas ocultas.
#for (i in 1:length(P)) {
# IDW.out[i] <- gstat::idw(rainfall ~ 1, P[-i,], P[i,], idp=2.0)$var1.pred
#}# Traza las diferencias
OP <- par(pty="s", mar=c(4,3,0,0))
plot(IDW.out ~ P$rainfall, asp=1, xlab="Observado", ylab="Predicho", pch=16,
col=rgb(0,0,0,0.5))
abline(lm(IDW.out ~ P$rainfall), col="red", lw=2,lty=2)
abline(0,1)
par(OP)El error cuadrático medio (RMSE) se puede calcular a partir de IDW.out de la siguiente manera:
# Compute RMSE
sqrt( sum((IDW.out - P$rainfall)^2) / length(P))## [1] 3.555055Validación cruzada
Según la información suministrada podemos realizar otro tipo de superficie interpolada: “Además de generar una superficie interpolada, puede crear un mapa de intervalo de confianza del 95% del modelo de interpolación. Aquí crearemos un mapa de IC del 95% a partir de una interpolación IDW que usa un parámetro de potencia de 2 (idp = 2.0)” (Lizarazo, 2020).
La línea de código que sigue genera unos comentarios resultantes extensos por eso permanecen ocultos y mejor se escribieron como comentarios para presentar la línea de código.
# Implementación de una técnica navaja para estimar
# un intervalo de confianza en cada punto no muestreado.
# Crea la superficie interpolada
#img <- gstat::idw(rainfall~1, P, newdata=grd, idp=2.0)
#n <- length(P)
#Zi <- matrix(nrow = length(img$var1.pred), ncol = n)
# Retire un punto y luego interpole (haga esto n veces para cada punto)
#st <- stack()
#for (i in 1:n){
# Z1 <- gstat::idw(rainfall~1, P[-i,], newdata=grd, idp=2.0)
# st <- addLayer(st,raster(Z1,layer=1))
# Pseudovalor calculado Z en j
# Zi[,i] <- n * img$var1.pred - (n-1) * Z1$var1.pred
#}
# Estimador Jackknife del parámetro Z en la ubicación j
#Zj <- as.matrix(apply(Zi, 1, sum, na.rm=T) / n )
# Calcular (Zi * - Zj) ^ 2
#c1 <- apply(Zi,2,'-',Zj) # Calcule la diferencia
#c1 <- apply(c1^2, 1, sum, na.rm=T ) # Suma el cuadrado de la diferencia
# Calcule el intervalo de confianza
#CI <- sqrt( 1/(n*(n-1)) * c1)
# Crear ráster (CI / valor interpolado)
#img.sig <- img
#img.sig$v <- CI /img$var1.pred # Ploteado
tm_shape(r.m) + tm_raster(n=7,title="IDW \ n intervalo de confianza del 95% \ n (en mm)") +
tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)
Ajuste polinomial de primer orden
Para ajustar un modelo polinomial de primer orden de la forma precip = intersección + aX + bY
# Definir la ecuación polinomial de primer orden
f.1 <- as.formula(rainfall ~ X + Y)
# Suma X e Y a P
P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]
# Ejecutar el modelo de regresión
lm.1 <- lm( f.1, data=P)
# Utilice la salida del modelo de regresión para interpolar la superficie
dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) # Recorte el ráster interpolado a Texas
r <- raster(dat.1st)
r.m <- raster::mask(r, nei2)
# Trazar el mapa
tm_shape(r.m) +
tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,
title="Primer polinomio de primer orden \ nPrecipitación prevista \ n (en mm)") +
tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)## Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
## for numeric data and stretch.palette for categorical data to control the palette
## mapping.
Ajuste polinomial de segundo orden
Para ajustar un modelo polinomial de segundo orden de la forma precip = intersección + aX + bY + dX2 + eY2 + fXY
# Definir la ecuación polinomial de segundo orden
f.2 <- as.formula(rainfall ~ X + Y + I(X*X)+I(Y*Y) + I(X*Y))
# Suma X e Y a P
P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]
# Ejecutar el modelo de regresión
lm.2 <- lm( f.2, data=P)
# Utilice la salida del modelo de regresión para interpolar la superficie
dat.2nd <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.2, newdata=grd))) # Recorte el ráster interpolado a Texas
r <- raster(dat.2nd)
r.m <- raster::mask(r, nei2)
# Trazar el mapa
tm_shape(r.m) +
tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,midpoint = NA,
title="Ajuste polinomial de segundo orden \ nPrecipitación prevista \ n (en mm)") +
tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)## Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
## for numeric data and stretch.palette for categorical data to control the palette
## mapping.
Kriging
Ajustar el modelo de variograma
Primero, necesitamos crear un modelo de variograma. Tenga en cuenta que el modelo de variograma se calcula sobre los datos sin tendencia. Esto se implementa en el siguiente fragmento de código pasando el modelo de tendencia de primer orden (definido en un fragmento de código anterior como el objeto fórmula f.1) a la función variograma.
# Definir la ecuación polinomial de primer orden
f.1 <- as.formula(rainfall ~ X + Y)
# Calcule el variograma de muestra; tenga en cuenta que el modelo de tendencia f.1 es uno de los
# parámetros pasados a variogram (). Esto le dice a la función que cree el
# variograma sobre los datos sin tendencia.
var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=100000, width=8990)
# Calcule el modelo de variograma pasando los valores de nugget, umbral y rango
# para ajustar.variogram () a través de la función vgm ().
dat.fit <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
vgm(psill=3, model="Mat", range=150000, nugget=0.0))## Warning in fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE, : No
## convergence after 200 iterations: try different initial values?## El modelo puede ser Exp, Gau, Sph, Mat
## Ver paquete gstat - función vgm - para obtener los nombres de estos modelos
## Repasa la teoría de kriging para comprender su significadodat.fit# La siguiente gráfica nos permite evaluar el ajuste
plot(var.smpl, dat.fit, xlim=c(0,130000))
Genera superficie Kriged
A continuación, utilice el modelo de variograma dat.fit para generar una superficie interpolada kriged. La función krige nos permite incluir el modelo de tendencia, lo que nos ahorra tener que eliminar la tendencia de los datos, krige los residuos y luego combinar los dos rásteres. En cambio, todo lo que tenemos que hacer es pasar a krige la fórmula de tendencia f.1.
# Definir el modelo de tendencia
f.1 <- as.formula(rainfall ~ X + Y)
# Realice la interpolación de krige (tenga en cuenta el uso del modelo de variograma
# creado en el paso anterior)
dat.krg <- krige( f.1, P, grd, dat.fit)## Warning in proj4string(d$data): CRS object has comment, which is lost in output## Warning in proj4string(newdata): CRS object has comment, which is lost in output## [using universal kriging]## Warning in proj4string(newdata): CRS object has comment, which is lost in output## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
## but +towgs84= values preserved# Convierta la superficie kriged en un objeto ráster para recortar
r <- raster(dat.krg)
r.m <- raster::mask(r, nei2)r.m## class : RasterLayer
## dimensions : 321, 312, 100152 (nrow, ncol, ncell)
## resolution : 789.6172, 789.6172 (x, y)
## extent : 716320.5, 962681.1, 663489.2, 916956.3 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : var1.pred
## values : 3.245264, 42.41171 (min, max)# Mapa
tm_shape(r.m) +
tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,
title="Kriging universal \ nPrecipitación prevista \ n (en mm)") +
tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)## Warning: The argument auto.palette.mapping is deprecated. Please use midpoint
## for numeric data and stretch.palette for categorical data to control the palette
## mapping.
Una mejor visualización de la superficie interpolada:
library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("red", "orange", "yellow", "blue", "darkblue"), values(precip.mask),
na.color = "transparent")
leaflet() %>% addTiles() %>%
addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
addLegend(pal = pal, values = values(r.m),
title = "Kriging interpolación precipitación en Huila \ n De 26.09 a 30.09 en 2020 [mm]")## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definitionGenere los mapas de intervalo de confianza y varianza
El objeto dat.krg almacena no solo los valores interpolados, sino también los valores de varianza. Estos se pueden pasar al objeto ráster para el mapeo de la siguiente manera:
r <- raster(dat.krg, layer="var1.var")
r.m <- raster::mask(r, nei2)
tm_shape(r.m) +
tm_raster(n=7, palette ="Reds",
title="Interpolación de Kriging \ nMapa de varianza \ n (en mm cuadrados)") +tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)
Un mapa más fácilmente interpretable es el mapa del intervalo de confianza del 95% que se puede generar a partir del objeto de varianza de la siguiente manera (los valores del mapa deben interpretarse como el número de mm por encima y por debajo de la cantidad de lluvia estimada).
r <- sqrt(raster(dat.krg, layer="var1.var")) * 1.96
r.m <- raster::mask(r, nei2)
tm_shape(r.m) +
tm_raster(n=7, palette ="Reds",
title="Interpolación de Kriging \ nMapa de IC del 95% \ n (en mm)") +tm_shape(P) + tm_dots(size=0.2) +
tm_legend(legend.outside=TRUE)
sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252 LC_CTYPE=Spanish_Colombia.1252
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] spatstat_1.64-1 rpart_4.1-15 nlme_3.1-148
## [4] spatstat.data_1.5-2 gstat_2.0-6 tmap_3.2
## [7] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
## [10] purrr_0.3.4 readr_1.3.1 tidyr_1.1.2
## [13] tibble_3.0.3 ggplot2_3.3.2 tidyverse_1.3.0
## [16] raster_3.3-13 rgdal_1.5-16 sp_1.4-2
## [19] rgl_0.100.54 leaflet_2.0.3 RColorBrewer_1.1-2
## [22] sf_0.9-6
##
## loaded via a namespace (and not attached):
## [1] leafem_0.1.3 colorspace_1.4-1 deldir_0.2-3
## [4] ellipsis_0.3.1 class_7.3-17 base64enc_0.1-3
## [7] fs_1.5.0 dichromat_2.0-0 httpcode_0.3.0
## [10] rstudioapi_0.11 farver_2.0.3 fansi_0.4.1
## [13] lubridate_1.7.9 xml2_1.3.2 splines_4.0.2
## [16] codetools_0.2-16 knitr_1.29 polyclip_1.10-0
## [19] jsonlite_1.7.1 tmaptools_3.1 broom_0.7.0
## [22] dbplyr_1.4.4 png_0.1-7 rgeos_0.5-5
## [25] shiny_1.5.0 compiler_4.0.2 httr_1.4.2
## [28] backports_1.1.9 Matrix_1.2-18 lazyeval_0.2.2
## [31] assertthat_0.2.1 fastmap_1.0.1 cli_2.0.2
## [34] later_1.1.0.1 leaflet.providers_1.9.0 htmltools_0.5.0
## [37] tools_4.0.2 gtable_0.3.0 glue_1.4.2
## [40] geojson_0.3.4 V8_3.4.0 Rcpp_1.0.5
## [43] cellranger_1.1.0 vctrs_0.3.4 crul_1.0.0
## [46] leafsync_0.1.0 crosstalk_1.1.0.1 lwgeom_0.2-5
## [49] xfun_0.16 rvest_0.3.6 mime_0.9
## [52] miniUI_0.1.1.1 lifecycle_0.2.0 goftest_1.2-2
## [55] XML_3.99-0.5 jqr_1.2.0 zoo_1.8-8
## [58] scales_1.1.1 spatstat.utils_1.17-0 hms_0.5.3
## [61] promises_1.1.1 parallel_4.0.2 yaml_2.2.1
## [64] curl_4.3 stringi_1.5.3 maptools_1.0-2
## [67] e1071_1.7-3 manipulateWidget_0.10.1 intervals_0.15.2
## [70] rlang_0.4.7 pkgconfig_2.0.3 evaluate_0.14
## [73] lattice_0.20-41 tensor_1.5 htmlwidgets_1.5.1
## [76] tidyselect_1.1.0 magrittr_1.5 R6_2.4.1
## [79] geojsonio_0.9.2 generics_0.0.2 DBI_1.1.0
## [82] mgcv_1.8-31 foreign_0.8-80 pillar_1.4.6
## [85] haven_2.3.1 withr_2.2.0 units_0.6-7
## [88] stars_0.4-3 xts_0.12.1 abind_1.4-5
## [91] spacetime_1.2-3 modelr_0.1.8 crayon_1.3.4
## [94] KernSmooth_2.23-17 rmarkdown_2.3 grid_4.0.2
## [97] readxl_1.3.1 blob_1.2.1 FNN_1.1.3
## [100] reprex_0.3.0 digest_0.6.25 classInt_0.4-3
## [103] webshot_0.5.2 xtable_1.8-4 httpuv_1.5.4
## [106] munsell_0.5.0 viridisLite_0.3.0Análisis de resultados
Polígonos de Thiessen
Son 37 los municipios que hacen parte del depertamneto de Huila, cada uno caracterizado por precipitaciones diferentes a lo largo de un mes o año. A través de este método se obtuvo información del nivel de precpitaciones en cada municipio. Guadalupe, Suaza, Acevedo, Palestina y la parte sur de Pitalito son los municipios que presentaron una mayor cantidad de lluvias en comparación con los demás municipios (por encima de 40 mm desde el 26 de septiembre al 30 del mismo mes). Por otro lado, existen algunos municipios y parte de algunos de ellos en los que las precipitaciones se presentaron por debajo de 10 mm: Hobo, Yaguará, Tesalia, Gigante, Algeciras, Campoalegre, Rivera, Palermo, Neiva, Tello, Baraya, Villavieja, Aipe y colombia. Se analizaron varios niveles de precipitación en el departamento y se compararón con información sobre el mismo: el mapa de precipitación se aproxima a dar graficamente las caractteristicas de la condición de cada municipio. Por ejemplo: El desierto de la Tatacoa ubicado en el municipio de Villavieja es una de las regiones que está representada por un color casi blanco que indica la escaza frecuencia de precipitacion.
Ponderación de Distancia Inversa (IDW)
Comparado con el método anterior, a simple vista las informacion que se representa en el mapa no varía. Se llega a un mismo análisis donde los zonas de precipitación altas, medias y bajas son las mismas. Pero, el aspecto en el que se diferencia este método del anterior, está en que IDW proporciona un recurso de error cuadrático medio y un intervalo de confianza, que precisamente permite observar que tan confiable es el mapa representado.
Kriging
El modelo de variograma se calcula sobre los datos sin tendencia. De acuerdo a esto, el variograma que se analizó con los datos de precipitación dieron como resultado una linea de tendencia que se inclina hasta el punto de parecer plana, pero un conjunto de puntos en forma creciente. Luego, al generar la superficie Kriged esta tendencia está incluida y, analizando el mapa se identificaron zonas con una mayor descripción del nivel de precipitaciónes. Lo que hace Kriging a diferencia de los métodos anteriores es delimitar unas zonas bien definidas. De este modo hay mayor información en toda la región. Mientras que en los métodos anteriores habían demarcaciones de un solo valor de cantidad de lluvias, en este método estas zonas contienen a otras zonas con distintos niveles de precipitación.
Conclusiones
La variabilidad del clima actual en el que vivímos requiere de métodos como estos para realizar un seguimiento a través de los años, esto con el propósito de recopilar información, datos que más adelante sirvan como herramienta, para tomar deciisones de cual debe ser la vocación agrícola de una región.
Las representaciones gráficas de los datos recopilados resultan prácticos y abiertos al público al realizar una composición de la recoleción de mucha información.
Referencias
Lizarazo, I. (2020, 22 de mayo). Interpolation of precipitation data. https://bit.ly/3qqk67Y
Introduction to Spatial Analysis. https://bit.ly/2VEoTEG
Gobernación de Huila. (1969, 31 de noviembre). Identificación del departamento. https://bit.ly/39SXxTs
Instituto Geográfico Agustín Codazzi. (2016, 12 de mayo). Cultivos y ganado solo deberían estar presentes en el 15% del Huila. https://bit.ly/33PqbRE
Gobernación de Huila. (2017, 16 de febrero). Economía. https://bit.ly/3ghv2ju
Gobernación de Huila. (1969, 31 de diciembre). Relieve y Clima. https://bit.ly/3mLpWyJ
CHIRPS.https://bit.ly/2IhuaPz
GIS Geography. (2020, 17 de octubre). Inverse Distance Weighting (IDW) Interpolation. https://bit.ly/3gfc5Ot
Organización de las Naciones Unidas para la Alimentación y la Agricultura. (s.f.). 7.1.1 INTERPOLACIÓN DE MODELOS DE ELEVACIÓN DIGITAL (MED) MEDIANTE EL MÉTODO KRIGING. https://bit.ly/3lO3WSp