Interpolación de datos de radiación solar para el departamento de Arauca

Resumen

La radiación solar resulta de gran interés agronómico debido a que es el factor que más influye en la fotosíntesis de las plantas, para obtener altas tasas de fotosíntesis se encuentra necesario maximizar la intercepción de la luz encontrando lugares idóneos de radiación solar, su máxima intercepción por parte de las plantas se observa en una mayor cantidad de materia seca, sobre todo en cultivos extensos y un crecimiento optimo en todas las etapas de crecimiento de la planta. La presente evalúa 3 métodos de interpolación espacial en R usando datos de radiación solar suministrados por World Clim donde se destacan metodos como kriging ordinal, ponderación de la distancia inversa (IDW) y poligonos de thieseen aplicados al departamento de Arauca en Colombia; Finalmente encontrara una muestra de los resultados y la interpretación de ellos para denotar la necesidad de radiación solar para los cultivos mas destacados del departamento.

Palabras clave: Radiación, distribución espacial, interpolación, kriging, IDW, poligonos de Thiessen.

1. Introducción

La energía que determina la dinámica de los procesos atmosféricos y climáticos emitida por el sol es conocida como radiación solar, la energía producida por el sol es radiación electromagnética proporcionada por reacciones de hidrogeno en el núcleo del sol por fusión nuclear y luego emitidas por la superficie solar, la radiación solar es esencial para los procesos productivos de un cultivo, sin embargo la exposición excesiva a esta fuente de energía resulta en estrés térmico que afecta notoriamente su desarrollo, afectando la planta desde su periodo vegetativo y de crecimiento, hasta su estado de cosecha;El sol emite energía en forma de radiación de onda corta principalmente en la banda del ultravioleta visible y el infrarrojo cercano, con longitudes de onda entre 0,2 y 3,0 micrómetros, el espectro electromagnético de este no tiene definidos limites superiores ni inferiores, y la energía de una fracción de radiación (fotón) es inversamente proporcional a su longitud de onda, por lo que a menor longitud de onda mayor contenido energético .

Para determinar que datos de radiación solar tenemos en Arauca es menester realizar una interpolación debido a que tomaremos unos datos distribuidos espacialmente en formato raster; Los objetos con una distribución espacial tienden a tener características similares entre más cerca se encuentren, es decir, sí tenemos un valor de radiación solar en un punto espacial de Arauca podemos proveer con un alto nivel de confianza, que este valor es similar o sigue una tendencia en un punto cercano al punto de muestra.

La presente evalua 3 métodos de interpolación espacial en R, usando datos de radiación solar suministrados por World Clim donde se destacan metodos como kriging ordinal, ponderación de la distancia inversa (IDW) y poligonos de thieseen aplicados al departamento de Arauca en Colombia; Finalmente encontrara una interpretación de los resultados y la interpretación de ellos para denotar la necesidad de radiación solar para los cultivos mas destacados del departamento.

2. Datos y metodos

2.1 Zona de estudio

En el norte de los llanos orientales, entre la cordillera oriental y la zona limítrofe con Venezuela se encuentra ubicado el departamento de Arauca, se localiza entre los 06°, 02’ 40” y 07° 06’ 13” de latitud norte y los 69° 25’ 54” y 72° 22’ 23” de longitud oeste, en el extremo norte de la Orinoquía colombiana, territorio denotado por la sierra nevada del cocuy y atravesado por los ríos meta, Arauca y Casanare. El epicentro agrícola, ganadero y petrolero del oriente colombiano comprende sistemas de selva, sabana y afluentes; Comprende los municipios de Tame, Arauquita, Fortul, Puerto Rondón, Saravena, Cravo Norte y Arauca como capital.Esta tierra ubicada en el extremo nororiental de la Orinoquia colombiana, posee un relieve mixto donde encontramos zonas montañosas las cuales no superan los 5000 m.s.n.m, sabanas que se levantan sobre los 1000 m.s.n.m y una temperatura entre los 28° y 30°.

Arauca posee tres conjuntos morfológicos: la cordillera oriental que cubre una quinta parte del departamento y comprende elevaciones de 500 hasta 5380 m.s.n.m; pie de monte cubierta de sabana y áreas ampliamente vegetativas, bosque ecuatorial y conformada por conos, abanicos, aluviales, terrazas de relieve plano e inclinado y la llanura aluvial; La formación orográfica mas destacada es la sierra nevada del cocuy la cual cuenta con accidentes geográficos como la piedra, el diamante, los altos nievecitas, osos y las cuchillas de Altamira y el salitre. Debido a su ubicación latitudinal, Arauca presenta un clima ecuatorial lluvioso, con influencia de los vientos alisios y la cordillera oriental; Presentando dos periodos climáticos: Entre abril y noviembre el clima predominante es lluvioso y entre noviembre y marzo predomina un clima seco sin llegar a alcanzar los 35°C.

library(sf)
library(cartography)

opar <- par(mar = c(0,0,1.2,0))

par(bg="grey75")
# plot municipalities
plot(st_geometry(nbi_munic_2), col = "#e4e9de", border = "red", 
     bg = "grey75", lwd = 0.5)

# plot labels
labelLayer(
  x = nbi_munic_2, 
  txt = "MUNICIPIO", 
  col= "yellow", 
  cex = 0.4, 
  font = 4,
  halo = TRUE, 
  bg = "grey25", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)

# map layout
layoutLayer(
  title = "Departamento de Arauca con sus respectivos municipios", 
  author = "Daniel Cortés", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal"
) 

2.2 Datos

En la presente se tomaran datos de radiación solar suministrados por WorldClim, usaremos la versión lanzada en enero del 2020 donde se encontraran 19 variables bioclimáticas a 4 resoluciones espaciales, cada descarga contiene un archivo zip donde se encontraran 12 archivos GeoTiff, cada archivo corresponde a un mes del año; Para este informe se usaron datos de radiación solar medidos en kJ m -2 día -1 ,a una resolución espacial de 5 minutos y se tomara un archivo Geotiff que contiene el promedio de radiación solar global para el mes de diciembre.

La medición de la radiación solar es de gran valor para un amplio rango de aplicaciones, en las áreas de ingeniería, arquitectura, agricultura, ganadería, salud humana y meteorología dentro de las cuales se destaca fundamentalmente el uso como fuente de energía alternativa. La luz solar es esencial para los mecanismos de crecimiento de las plantas, pero la exposición excesiva a esta fuente de energía resulta en estrés térmico que afecta su desarrollo, crecimiento y periodo productivo de la planta, especialmente en periodos de siembra y en plantas maduras durante todo el periodo vegetativo, afectando negativamente la etapa de cosecha del producto.

Arauca al ser una región caracterizada por su potencial productivo de plátano y tubérculos, se encuentran influenciados estos sistemas de producción por factores bióticos y abióticos, las altas y bajas temperaturas y la radiación solar son factores que afectan considerablemente la producción, para el caso del plátano, este es sembrado en condiciones muy variadas de radiación solar, sin embargo la falta de luz no interrumpe la emisión y desarrollo de las hojas, pero los limbos quedan blanquecidos debido a la ausencia de síntesis de clorofila, como consecuencia las vainas foliares y pseudotallos se alargan demasiado, la radiación también influye directamente en el proceso de maduración y composición química de los frutos de plátano.

A continuacion trataremos los datos que usaremos para interpolar: Iniciaremos limpiando la memoria

#rm(list=ls())

Cargaremos las librerias que necesitaremos para tratar los datos y interpolar

library(knitr)
library(rgdal)
library(raster)
library(sf)
library(tidyverse)
library(tmap)
library(gstat)
library(sp)
library(spatstat.data)
library(nlme)
library(rpart)
library(spatstat)

Leemos el conjunto de datos de radiación solar en formato ráster

## Leyendo el conjunto de datos de radiación solar en formato raster
carborg <- raster("C:\\Users\\johnedisoncortes\\Documents\\Inf3GMB\\rad.tif")

Carborg sera nuestros datos de radiación solar

## Qué adjuntamos en la variable carborg?
carborg

class      : RasterLayer 
dimensions : 2160, 4320, 9331200  (nrow, ncol, ncell)
resolution : 0.08333333, 0.08333333  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : C:/Users/johnedisoncortes/Documents/Inf3GMB/rad.tif 
names      : rad 
values     : 0, 45397  (min, max)

La variable “aoi” contiene el shapefile de nuestra área de interés, para este caso el departamento de Arauca.

(aoi <- shapefile("C:\\Users\\johnedisoncortes\\Documents\\Inf2GMB\\MGN2017_81_ARAUCA\\81_ARAUCA\\ADMINISTRATIVO\\MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 7 
extent      : -72.36662, -69.42756, 6.036228, 7.104381  (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  :         81,      81001,     ARAUCA,                                 1959,  945.13540708,      2017,     ARAUCA, 1.29412603421, 0.0772141267598 
max values  :         81,      81794,       TAME, Decreto Nal 677 de Abril 13 de  1987, 5787.94213047,      2017,     ARAUCA, 4.66801606458,  0.471626116664 

Datos de radiación solar para nuestra area de interes

carborg.crop <- raster::crop(carborg, extent(aoi))

Enmascaramos nuestra capa raster para nuestra area de interes

carborg.mask <- mask(x = carborg.crop, mask = aoi)

carborg.mask 

class      : RasterLayer 
dimensions : 13, 35, 455  (nrow, ncol, ncell)
resolution : 0.08333333, 0.08333333  (x, y)
extent     : -72.33333, -69.41667, 6, 7.083333  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : rad 
values     : 14780, 16991  (min, max)

Trazamos la capa raster

plot(carborg.mask, main=  "wordclim radiación solar para el departamento de Arauca [kj m^-2 dia^-1]")
plot(aoi, add=TRUE)

Plotearemos un mapa para obtener una mejor visualización de nuestros datos:

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("darkblue", "blue", "yellow", "orange", "red"), values(carborg.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(carborg.mask, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(carborg.mask),
    title = "wordclim sobre radiación solar para el departamento de Arauca [kj m^-2 dia^-1]")

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 +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded 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 +type=crsDiscarded datum World Geodetic System 1984 in CRS definition

Con la funcion rasterToPoints de la libreria raster podemos convertir el raster a puntos

carborg.points <- rasterToPoints(carborg.mask, spatial = TRUE)

carborg.points

class       : SpatialPointsDataFrame 
features    : 278 
extent      : -72.29167, -69.45833, 6.041667, 7.041667  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :   rad 
min values  : 14780 
max values  : 16991 

str(carborg.points)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 278 obs. of  1 variable:
  .. ..$ rad: num [1:278] 14933 15421 16890 16854 16833 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:278, 1:2] -71.8 -71.6 -70.9 -70.8 -70.7 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -72.29 6.04 -69.46 7.04
  .. ..- 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"

Ploteamos nuestra conversión de raster a puntos

plot(carborg.mask, main= "wordclim radiación solar para el departamento de Arauca [kj m^-2 dia^-1]")
plot(aoi, add=TRUE)

points(carborg.points$x, carborg.points$y, col = "blue", cex = 0.6)

Escribamos el objeto espacial de puntos en un archivo de disco. Solo para ahorrar tiempo y trabajo en caso de que algo salga mal.

Usaremos la biblioteca rgdal como se ilustra

geojsonio::geojson_write(carborg.points, file = "./Inf3GMB")

Success! File is at ./Inf3GMB.geojson

<geojson-file>
  Path:       ./Inf3GMB
  From class: SpatialPointsDataFrame

Leamos los puntos de precipitación. Consulte la documentación de geojsonio para comprender qué parámetros deben pasarse a la función geojson_read

carborg.points <- geojsonio::geojson_read("./Inf3GMB.geojson", what="sp")

carborg.points

class       : SpatialPointsDataFrame 
features    : 278 
extent      : -72.29167, -69.45833, 6.041667, 7.041667  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :   rad 
min values  : 14780 
max values  : 16991 

Preparación adicional de nuestros datos del area de interes (aoi)

(aoi <- shapefile("C:\\Users\\johnedisoncortes\\Documents\\Inf2GMB\\MGN2017_81_ARAUCA\\81_ARAUCA\\ADMINISTRATIVO\\MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 7 
extent      : -72.36662, -69.42756, 6.036228, 7.104381  (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  :         81,      81001,     ARAUCA,                                 1959,  945.13540708,      2017,     ARAUCA, 1.29412603421, 0.0772141267598 
max values  :         81,      81794,       TAME, Decreto Nal 677 de Abril 13 de  1987, 5787.94213047,      2017,     ARAUCA, 4.66801606458,  0.471626116664 

Necesitamos convertir nuestros datos en una caracteristica espacial

arauca_sf <-  sf::st_as_sf(aoi)

Disolveremos los limites internos

(border_sf <-
  arauca_sf %>%
  summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -72.36662 ymin: 6.036228 xmax: -69.42756 ymax: 7.104381
geographic CRS: WGS 84
      area                       geometry
1 23851.41 POLYGON ((-69.76161 6.53034...

Convertimos nuestra caracteristica espacial en datos espaciales

(border <- as(border_sf, 'Spatial')) 

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -72.36662, -69.42756, 6.036228, 7.104381  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :           area 
value       : 23851.40798475 

(arauca.sf <- st_as_sf(aoi) %>% mutate(MUNIC = MPIO_CNMBR, CODIGO = MPIO_CCDGO) %>% select(MUNIC, CODIGO))

Simple feature collection with 7 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -72.36662 ymin: 6.036228 xmax: -69.42756 ymax: 7.104381
geographic CRS: WGS 84
           MUNIC CODIGO                       geometry
1         ARAUCA  81001 POLYGON ((-70.68038 7.09393...
2 PUERTO RONDÓN  81591 POLYGON ((-70.87879 6.62133...
3    CRAVO NORTE  81220 POLYGON ((-70.40276 6.64639...
4      ARAUQUITA  81065 POLYGON ((-71.58441 7.04298...
5         FORTUL  81300 POLYGON ((-71.54473 6.75924...
6       SARAVENA  81736 POLYGON ((-71.76224 7.06602...
7           TAME  81794 POLYGON ((-71.73216 6.70303...

Realizamos la conversión

# conversion
p.sf <- st_as_sf(carborg.points)

Realizaremos una intersección de poligonos con puntos, conservando la información de ambos en un mismo lugar.

(carborg.sf = st_intersection(arauca.sf, p.sf))

although coordinates are longitude/latitude, st_intersection assumes that they are planar
attribute variables are assumed to be spatially constant throughout all geometries

Simple feature collection with 278 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -72.29167 ymin: 6.041667 xmax: -69.45833 ymax: 7.041667
geographic CRS: WGS 84
First 10 features:
        MUNIC CODIGO   rad                   geometry
6    SARAVENA  81736 14933 POINT (-71.79167 7.041667)
4   ARAUQUITA  81065 15421   POINT (-71.625 7.041667)
1      ARAUCA  81001 16890   POINT (-70.875 7.041667)
1.1    ARAUCA  81001 16854 POINT (-70.79167 7.041667)
1.2    ARAUCA  81001 16833 POINT (-70.70833 7.041667)
1.3    ARAUCA  81001 16876   POINT (-70.625 7.041667)
1.4    ARAUCA  81001 16893 POINT (-70.54167 7.041667)
6.1  SARAVENA  81736 16486 POINT (-71.95833 6.958333)
6.2  SARAVENA  81736 15884   POINT (-71.875 6.958333)
6.3  SARAVENA  81736 14996 POINT (-71.79167 6.958333)

Realizaremos dos reproyecciones, una de nuestros datos y otra de nuestro sitio de interes.

p.sf.magna <- st_transform(carborg.sf, crs=3116)

arauca.sf.magna <- st_transform(arauca.sf, crs=3116)

(carborg2 <- as(p.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definitionDiscarded datum Marco Geocentrico Nacional de Referencia in CRS definition

class       : SpatialPointsDataFrame 
features    : 278 
extent      : 1197609, 1511834, 1161866, 1271914  (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 +units=m +no_defs 
variables   : 3
names       :  MUNIC, CODIGO,   rad 
min values  : ARAUCA,  81001, 14780 
max values  :   TAME,  81794, 16991 

shapefile(carborg2, filename='C:\\Users\\johnedisoncortes\\Documents\\Inf3GMB\\carborg2.shp', layer= "carborg2", overwrite=TRUE)

carborg2

class       : SpatialPointsDataFrame 
features    : 278 
extent      : 1197609, 1511834, 1161866, 1271914  (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 +units=m +no_defs 
variables   : 3
names       :  MUNIC, CODIGO,   rad 
min values  : ARAUCA,  81001, 14780 
max values  :   TAME,  81794, 16991 

(arauca2 <- as(arauca.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definitionDiscarded datum Marco Geocentrico Nacional de Referencia in CRS definition

class       : SpatialPolygonsDataFrame 
features    : 7 
extent      : 1189337, 1515269, 1161263, 1278742  (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 +units=m +no_defs 
variables   : 2
names       :  MUNIC, CODIGO 
min values  : ARAUCA,  81001 
max values  :   TAME,  81794 

shapefile(arauca2, filename='C:\\Users\\johnedisoncortes\\Documents\\Inf3GMB\\arauca2.shp', layer= "arauca2", overwrite=TRUE)

Reemplazamos la extensión del límite del punto con la de Arauca

carborg2@bbox <- arauca2@bbox

Graficamos nuestra area de interes con los datos espaciales de radiación solar.

tm_shape(arauca2) + tm_polygons() +
  tm_shape(carborg2) +
  tm_dots(col="rad", palette = "RdBu", midpoint = 3.0,
             title="Muestreo de radiación solar [kj m^-2 dia^-1]", size=0.2) +
  tm_text("rad", just="center", xmod=.6, size = 0.5) +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output

2.3 Metodos

IDW

La ponderación de distancia inversa (IDW) es un método determinístico de interpolado en donde se da por hecho que la variable representada cartográficamente tiende a disminuir su influencia a mayor distancia desde la ubicación del punto de muestra, esto se logra a través de una combinación ponderada linealmente de un conjunto de puntos de muestreo, que generalmente suelen aplicarse a datos muy variables.

La potencia es la variable con la que podemos controlar la importancia de los puntos conocidos sobre los valores que queremos interpolar en función de su distancia desde su punto de salida.Es la variable que determina la tasa exponencial de caída de la influencia de los puntos vecinos cuanto mas lejos se encuentran del nodo de la cuadricula. Los valores asignados oscilan entre 1 y 10, especificando un valor más bajo dará mas importancia a los puntos que están mas lejos, obteniendo superficies mas suaves, por otro lado si asignamos una potencia alta se pone mas énfasis en los puntos mas cercanos por lo que el resultado será mas detallado y menos suave.

Krigin

Es un procedimiento geoestadístico que trabaja a partir de un conjunto de puntos dispersos para generar una superficie estimada, su principal diferencia de los otros métodos de interpolación es que se encuentra necesaria una investigación interactiva del comportamiento espacial de la variable representada. Interpolar con kriging presupone que la distancia entre los puntos de muestra refleja correlación espacial que se utiliza para explicar la variación en la superficie, al aplicar kriging ajustamos una función matemática a una cantidad especifica de puntos dentro de un radio determinado para conocer el valor de salida para cada ubicación.

Kriging esta basado en modelos estadísticos de autocorrelación, es decir, las relaciones estadísticas entre los puntos medidos, gracias a esto, podemos proporcionar alguna medida de certeza o precisión a la hora de realizar una estimación de un valor. Aplicando este método generamos una serie de procesos de análisis estadístico exploratorio, modelos de variogramas, creaciones de superficie, entre otros. Este método es adecuado cuando se sabe que hay una influencia de la distancia correlacionada espacialmente con los datos.

Poligonos de Thiessen

Los polígonos de Thiessen permiten establecer relaciones matemáticas entre elementos generando zonas de influencia con unas premisas matemáticas específicas, la relación parte de una nube de puntos sobe los que se generan una serie de polígonos, estos puntos se unen entre si y se proyectan mediatrices entre los segmentos de unión, siendo dichas mediatrices, los lados de los polígonos resultantes. La principal regla que se establece es que los lados de los polígonos generados, son equidistantes a los puntos vecinos y tratan de encontrar la menor distancia posible. Los lados de cada polígono se encuentran a la misma distancia de un punto que otro

3.Resultados

IDW

La función IDW se encuentra disponible en Spatstat como en Gstat, es importante mencionar que la función Dirichlet usada en esta interpolación (como la mayoría de las funciones del paquete spatstat) requiere un formato ppp, de allí la sintaxis as.ppp

La salida IDW es un raster. Requerimos crear una cuadricula ráster vacía, para luego interpolar los valores de radiación solar en cada celda de la cuadricula sin muestrear.

Utilizaremos un valor de potencia IDW de 2 (idp=2)

# Crearemos una cuadricula vacía donde n será el numero total de celdas
grd              <- as.data.frame(spsample(carborg2, "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  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

# Agregamos la información de proyección de C a la cuadrícula vacía
proj4string(grd) <- proj4string(carborg2)

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition

# Interpolamos las celdas de la cuadrícula usando un valor de potencia de 2 (idp = 2.0)
P.idw <- gstat::idw(rad ~ 1, carborg2, newdata=grd, idp=2.0)

CRS object has comment, which is lost in outputCRS object has comment, which is lost in output

[inverse distance weighted interpolation]

CRS object has comment, which is lost in outputDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition

# Convertimos a raster y limitamos nuestra area de estudio
r       <- raster(P.idw)

r

class      : RasterLayer 
dimensions : 190, 527, 100130  (nrow, ncol, ncell)
resolution : 618.7907, 618.7907  (x, y)
extent     : 1189078, 1515181, 1161195, 1278765  (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 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 14783.96, 16987.71  (min, max)

arauca2

class       : SpatialPolygonsDataFrame 
features    : 7 
extent      : 1189337, 1515269, 1161263, 1278742  (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 +units=m +no_defs 
variables   : 2
names       :  MUNIC, CODIGO 
min values  : ARAUCA,  81001 
max values  :   TAME,  81794 

r.m   <- raster::mask(r, arauca2)

r.m

class      : RasterLayer 
dimensions : 190, 527, 100130  (nrow, ncol, ncell)
resolution : 618.7907, 618.7907  (x, y)
extent     : 1189078, 1515181, 1161195, 1278765  (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 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 14783.96, 16987.71  (min, max)

# Plot
tm_shape(r.m) + 
  tm_raster(n=10,palette = "OrRd", auto.palette.mapping = FALSE,
            title="IDW de radiación solar prevista en [kj m^-2 dia^-1] ") + 
  tm_shape(carborg2) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

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 nuestro campo de interpolado

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("darkblue", "blue", "yellow", "orange", "red"), values(carborg.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(r.m),
    title = "IDW: Radiacíon solar interpolada en arauca medida en [kj m^-2 dia^-1] ")

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 +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded 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 +type=crsDiscarded datum World Geodetic System 1984 in CRS definition

Ajustamos la interpolación

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.

carborg2

class       : SpatialPointsDataFrame 
features    : 278 
extent      : 1189337, 1515269, 1161263, 1278742  (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 +units=m +no_defs 
variables   : 3
names       :  MUNIC, CODIGO,   rad 
min values  : ARAUCA,  81001, 14780 
max values  :   TAME,  81794, 16991 

c <- carborg2

IDW.out <- vector(length = length(c))

for (i in 1:length(c)) {
  IDW.out[i] <- gstat::idw(rad ~ 1, c[-i,], c[i,], idp=2.0)$var1.pred
}

Graficamos la diferencia entre los datos estimados y los datos observados

# Graficamos la diferencia
OP <- par(pty="s", mar=c(4,3,0,0))
  plot(IDW.out ~ c$rad, asp=1, xlab="Observed", ylab="Predicted", pch=16,
       col=rgb(0,0,0,0.5))
  abline(lm(IDW.out ~ c$rad), 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:

# Computar RMSE
sqrt( sum((IDW.out - c$rad)^2) / length(c))

[1] 251.2121

Validación cruzada

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). Encontrará en el codigo aplicado los metodos usados para realizar el proceso de interpolado.

# Implementamos una tecnica de navaja
# Estimando un intervalo de confianza en cada punto no muestreado

# Creamos la superficie de interpolado
img <- gstat::idw(rad~1, c, newdata=grd, idp=2.0)
n   <- length(c)
Zi  <- matrix(nrow = length(img$var1.pred), ncol = n)

# Removemos un punto para luego interpolar 
st <- stack()
for (i in 1:n){
  Z1 <- gstat::idw(rad~1, c[-i,], newdata=grd, idp=2.0)
  st <- addLayer(st,raster(Z1,layer=1))
  # calculamos la pseudovariable Z at j
  Zi[,i] <- n * img$var1.pred - (n-1) * Z1$var1.pred
}

# Estimador tipo navaja para el parametro z en la locación j
Zj <- as.matrix(apply(Zi, 1, sum, na.rm=T) / n )

# Calculamos (Zi* - Zj)^2
c1 <- apply(Zi,2,'-',Zj)            # Compute the difference
c1 <- apply(c1^2, 1, sum, na.rm=T ) # Sum the square of the difference

# Calculamos el intervalo de confianza
CI <- sqrt( 1/(n*(n-1)) * c1)

# Creamos (CI / Interpolación de la variable) raster
img.sig   <- img
img.sig$v <- CI /img$var1.pred 

Limitamos el raster de confianza a nuestra area de estudio

# Limitamos el raster de confianza a nuestra area de estudio
r <- raster(img.sig, layer="v")
r.m <- raster::mask(r, arauca2)

# Ploteamos nuestro mapa
tm_shape(r.m) + tm_raster(n=7,title="IDW\n95% confidence interval \n(in mm)") +
  tm_shape(c) + tm_dots(size=0.2) +
  
  tm_legend(legend.outside=TRUE)

Polinomio de primer orden

Este metodo sera incluido debido a que lo necesitaremos para generar nuestro interpolado por el metodo kriging.

# Define the 1st order polynomial equation
f.1 <- as.formula(rad ~ X + Y) 
 
# Add X and Y to P
c$X <- coordinates(c)[,1]
c$Y <- coordinates(c)[,2]

# Run the regression model
lm.1 <- lm( f.1, data=c)

# Use the regression model output to interpolate the surface
dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 

# Clip the interpolated raster to Texas
r   <- raster(dat.1st)
r.m <- raster::mask(r, arauca2)

# Plot the map
tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, 
            title="1st order polinomial fir \nPredicted precipitation \n(in mm)") +
  tm_shape(c) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

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 radsun = intersección + aX + bY + dX2 + eY2 + fXY

# Define the 2nd order polynomial equation
f.2 <- as.formula(rad ~ X + Y + I(X*X)+I(Y*Y) + I(X*Y))

# Add X and Y to P
c$X <- coordinates(c)[,1]
c$Y <- coordinates(c)[,2]
 
# Run the regression model
lm.2 <- lm( f.2, data=c)

# Use the regression model output to interpolate the surface
dat.2nd <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.2, newdata=grd))) 

# Clip the interpolated raster to Texas
r   <- raster(dat.2nd)
r.m <- raster::mask(r, arauca2)

# Plot the map
tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,midpoint = NA,
            title="2nd order polinomial fit\nPredicted precipitation \n(in mm)") +
  tm_shape(c) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

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

Definiremos la ecuación polinomial de primer orden, para luego calcular el variograma de muestra y posterior a ello ajustar nuestro variograma a travez de la funcion vgm()

# Definimos la ecuación polinomial de primer orden
f.1 <- as.formula(rad ~ X + Y) 

# Calculamos el variograma de muestra
# Esto le dice a la función que cree el variograma sobre los datos sin tendencia.
var.smpl <- variogram(f.1, c, cloud = FALSE, cutoff=120000, width=10000)

# 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=80000, nugget=0.0))

No convergence after 200 iterations: try different initial values?

dat.fit

# La siguiente grafica nos permite determinar el ajuste
plot(var.smpl, dat.fit, xlim=c(0,130000))

Generamos una superficie Kriged

A continuación, utilizaremos 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.

# Definimos el modelo de tendencia
f.1 <- as.formula(rad ~ X + Y) 

# Realizamos la interpolación de krige 
dat.krg <- krige( f.1, c, grd, dat.fit)

Convertiremos la superficie kriged en un objeto raster para recortar

# Conviertimos la superficie kriged en un objeto ráster para recortar
r <- raster(dat.krg)
r.m <- raster::mask(r, arauca2)

r.m

class      : RasterLayer 
dimensions : 190, 527, 100130  (nrow, ncol, ncell)
resolution : 618.7907, 618.7907  (x, y)
extent     : 1189078, 1515181, 1161195, 1278765  (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 +units=m +no_defs 
source     : memory
names      : var1.pred 
values     : 15312.72, 16968.03  (min, max)

# Plot the map
tm_shape(r.m) + 
  tm_raster(n=10, palette="Reds", auto.palette.mapping=FALSE, 
            title="K.O universal para Radiación solar [kj m^-2 dia^-1]") +
  tm_shape(c) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

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 nuestra superficie de interpolado

library(leaflet)
library(RColorBrewer)
pal <- colorNumeric(c("darkblue", "blue", "yellow", "orange", "red"), values(carborg.mask),
  na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(r.m, colors = pal, opacity = 0.6) %>%
  addLegend(pal = pal, values = values(r.m),
    title = "Interpolación por K.O en Arauca medido en [kj m^-2 dia^-1] ")

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 +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded 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 +type=crsDiscarded datum World Geodetic System 1984 in CRS definition

Generamos 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, arauca2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Interpolación por K.O/Mapa de varianza en n[kj m^-2 dia^-1] )") +tm_shape(c) + 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 [kj m^-2 dia^-1] por encima y por debajo de la cantidad de radiación solar estimada).

r   <- sqrt(raster(dat.krg, layer="var1.var")) * 1.96
r.m <- raster::mask(r, arauca2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Kriging Interpolation\n95% CI map \n(en [kj m^-2 dia^-1] )") +tm_shape(c) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

Thiessen poligonos

Los poligonos de Thiessen pueden ser creados usando la funcion dirichlet de spatstat

Primero crearemos una superficie teselada

# Create a tessellated surface
th  <-  as(dirichlet(as.ppp(carborg2)), "SpatialPolygons")

# The dirichlet function does not carry over projection information
# requiring that this information be added manually
crs(th) <- crs(carborg2)
crs(arauca2) <- crs(carborg2)

crs(th)

CRS arguments:
 +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000
+ellps=GRS80 +units=m +no_defs 

crs(carborg2)

CRS arguments:
 +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000
+ellps=GRS80 +units=m +no_defs 

La superficie teselada no almacena información de atributos de la capa de datos de puntos. Usaremos la función over () (del paquete sp) para unir los atributos de punto a la superficie teselada mediante una unión espacial. La función over () crea un marco de datos que necesitará ser agregado al objeto th creando así un objeto

th.z     <- over(th, carborg2, fn=mean)
th.spdf  <-  SpatialPolygonsDataFrame(th, th.z)

#Recortaremos la superficie teselada
th.clp   <- raster::intersect(arauca2,th.spdf)

# Mapeamos los datos
tm_shape(th.clp) + 
  tm_polygons(col="rad", palette="Reds", midpoint=3.0,
              title="Thiessen Polygons \nEstimación radiación solar \n(en [kj m^-2 dia^-1])") +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output

4.Interpretación de resultados

Los procesos de interpolación trabajados en la presente a partir de R lograron aclarar y denotar la concepción de análisis de datos geoespaciales, Aunque los 3 métodos de interpolación trabajados lograron arrojar resultados símiles, podemos afirmar que ciertos métodos fueron más precisos que otros. Los datos interpolados de radiación solar lograron mostrar como existen unos valores de radiación mas altos en los municipios de Arauquita, puerto rondón, Arauca y Cravo norte, en donde destaca el municipio cravo norte, el cual posee los valores más altos de radiación solar y es de los municipios donde la producción de plátano y tubérculos es mínima; Por otro lado, los municipios donde se concentra la mayor producción de plátano en el departamento de Arauca tienen valores más bajos de radiación solar como lo es Tame y Fortul.

4.1 IDW

Al aplicar el método IDW evidenciamos un proceso vigoroso de interpolado en cuanto este método fue capas de arrojar ploteos con transiciones de valores entre celdas no tan abruptos, es decir, siguen una tendencia clara, pese a usar un valor de potencia de 2. Nuestra diferencia de los valores estimados sobre los observados indica un cierto grado de precisión, sin embargo, la desviación de estos genera un valor de incertidumbre considerable; Sí bien puede verse que entre los dos hay una tendencia línea creciente, no existe una exactitud total puesto que se esta haciendo un proceso de interpolado en un área muy grande para los pocos puntos de muestra tomados.

4.2 Kriging

El método kriging a partir del semivariograma y el modelo de ajuste, lograron mostrar con eficacia que las observaciones estimadas tienen una correlación alta con las observaciones propuestas en campo, lo que lo hace de los mejores métodos de interpolado usados en este informe, gracias al algoritmo propuesto por kriging se generaron representaciones a escalas de colores que permitían ver fácilmente que zonas del departamento tenían valores más altos de radiación solar. El variograma suministrado para nuestros datos de radiación solar presenta un comportamiento logarítmico, pese a que este método de interpolado fue eficaz, el mapa de nivel de incertidumbre arrojo valores más altos en comparación al método IDW con valores oscilantes entre 2 y 3.

4.3 Poligonos de Thiessen

Interpolando con el método de polígonos de Thiessen podemos ver que los polígonos se encuentran a una misma medida, debido a que los puntos de muestra se encuentran uniformemente distribuidos espacialmente, pese a que no tenemos un indicativo de incertidumbre, podemos concluir que fue el método menos exacto de los trabajados, no es un método con una complejidad alta y presenta transiciones entre polígonos abruptas y poco sutiles.

Conclusiones

Finalmente podemos concluir que lo interesante de interpolar se encuentra en que los datos con una distribución espacial, se encuentran correlacionados entre si, es decir, los valores que están cerca tienden a tener características similares. El provecho de interpolar en cuanto a su precisión depende no solo del método que usaremos, sino además, de su resolución espacial y la densidad de puntos de muestra que tengamos. La presente evaluó 3 métodos de interpolación donde se destaca IDW, Kriging y polígonos de Thiessen, usando datos de radiación solar para el departamento de Arauca; Donde se pudo observar que el método mas fiable para los datos tratados fue el método de la ponderación de distancia inversa (IDW) puesto que nos suministró información de la distribución de valores de radiación solar a lo largo del departamento con el valor más bajo de incertidumbre respecto de los demás métodos trabajados. En el proceso se evidencio que los municipios de Tame, Fortul y Saravena poseen valores mas bajos de Radiación solar, estos departamentos, son quienes predominan productivamente en el departamento, mientras que los departamentos de Puerto rondón, Arauquita, Cravo norte y Arauca poseen los valores mas altos de radiación solar y a su vez son los municipios con una menor tasa de productividad agraria.

Referencias

Cuencar, J. Arauca Colombia guía turística. Adaptado de: https://ceo.uniandes.edu.co/images/Documentos/Gu%C3%ADa%20tur%C3%ADstica%20Arauca.pdf

IDEAM. (2014). La importancia de la radiación solar. Adaptado de: http://www.ideam.gov.co/web/tiempo-y-clima/radiacion-solar-ultravioleta

Barrera, J. Cardona, C. Cayón, D. (2011). El cultivo de plátano (MUSA AAB SIMMONDS) ecofisiología y manejo cultural sostenible. Adaptado de: https://editorialzenu.com/images/1467833541.pdf

ArcGIS. (2016). Vista general del conjunto de herramientas Interpolación. Adaptado de: https://desktop.arcgis.com/es/arcmap/10.3/tools/spatial-analyst-toolbox/an-overview-of-the-interpolation-tools.htm