Interpolacion de datos espaciales

INTRODUCCIÓN

La interpolación sirve para predecir valores para las celdas de un ráster a partir de una cantidad limitada de puntos de datos de muestra. Esta herramienta se puede utilizar para prever valores que no se conocen de cualquier dato de un punto geográfico como por ejemplo la elevación y precipitación.

Las herramientas de interpolación de superficie crean una superficie continua (o predicción) a partir de valores de punto de muestra. Lo anterior se hace viable cuando se hace la suposición de que los objetos distribuidos espacialmente están correlacionados espacialmente, o sea, las cosas que están cerca, tienden a parecerce. Así, los puntos cercanos a los puntos de muestra, tienen más posibilidades de ser más similares que los que están más alejados. Para esto, existen diferentes maneras de derivar una predicción para cada ubicación; cada método se denomina modelo. Con cada modelo, se realizan diferentes suposiciones acerca de los datos, y algunos modelos son más aplicables para determinados datos y cada uno produce predicciones usando diferentes cálculos.

En el presente trabajo, se realizarán 3 métodos de interpolación los cuales se describen a continuación:

Kriggin: Interpola una superficie de ráster a partir de puntos IDW: Interpola una superficie de ráster a partir de puntos utilizando una técnica de distancia inversa ponderada

IDW: Interpola una superficie de ráster a partir de puntos utilizando una técnica de distancia inversa ponderada

Polígonos de Thiessen: Se basa en considerar que los puntos de una zona coinciden con el valor medido más próximo y generan, como consecuencia, zonas de influencia siendo el resultado la unión de polígonos que generan distintas zonas temáticas

library(tidyverse)
library(sf)

DESCRIPCIÓN DE LA ZONA DE ESTUDIO

El departamento de Boyacá está situado en el centro del país, en la cordillera oriental de los Andes; localizado entre los 04°39′10″ y los 07°03′17″ de latitud norte y los 71°57′49″ y los 74°41’35″ de longitud oeste. Cuenta con una superficie de 23 189 km², lo que representa el 2,03 % del territorio nacional. Limita por el norte con los departamentos de Santander y norte de Santander, por el este con los departamentos de Arauca, Casanare y con el país vecino de Venezuela, por el sur con Cundinamarca, y por el oeste con Cundinamarca, Caldas y Antioquia.

a continuacion podra ver la ubicacion en el mapa nacional:

colmap<-  read_sf("C:/Users/David Perdomo/Desktop/Geomatica/COL_adm/COL_adm1.shp")
st_crs(colmap)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

ggplot() + geom_sf(data = colmap) + coord_sf(crs=st_crs(3978))

boyaca <-  colmap %>%   filter(NAME_1 == "Boyacá")
ggplot() + geom_sf(data = boyaca) 

Hidrografia de Boyacá

El sistema hidrográfico de Boyacá está constituido por numerosos ríos y quebradas que nacen en la cordillera Oriental, afluentes directos o de alguno de los tributarios de los ríos Magdalena, Meta y Arauca; entre los ríos que confluyen a la cuenca del Magdalena están el Ermitaño, Negro, Minero, Suárez, Sutamanchán, Sáchica, Chíquiza, Iguaque, Arcabuco o Pómeca, Ubazá, Riachuelo, Moniquirá, Chicamocha, Chorrera, Tuta, Pesca, Tota, Saguera, Sasa, Cambas, Loblanco, Rechiminiga, Chitano, Susacón.

Actividades Economicas de Boyacá:

La economía de Boyacá se basa principalmente en la producción agrícola y ganadera, la explotación de minerales, la industria siderúrgica, el comercio y el turismo. La agricultura se ha desarrollado y tecnificado en los últimos años; los principales cultivos son papa, maíz, cebolla, trigo, cebada, caña panelera, yuca. Existen en Boyacá alrededor de 6.800 explotaciones mineras, en su mayoría de carácter artesanal, principalmente de carbón coquizable, arcilla, roca fosfórica, arena, yeso, mármol y caliza de hierro; el carbón extraído técnicamente tienen reservas importantes en 19 municipios.

El petróleo se encuentra en el territorio Vásquez; en el occidente se explotan las esmeraldas consideradas las más famosas del mundo por su calidad y belleza. La producción artesanal es muy laboriosa especialmente en cerámica, tejidos de lana de oveja y fique, tagua, tapices, instrumentos musicales y cestería, entre otros

Clima de Boyacá:

El Clima del departamento de Boyacá está determinado principalmente por su relieve y posición geográfica, dando como resultado, que se presenten todos los pisos térmicos desde el nival, en la sierra nevada del Cocuy, hasta el cálido, en el Territorio Vásquez y la vertiente oriental de la cordillera Oriental.

El departamento presenta una alta variabilidad territorial en la distribución de la precipitación. Las áreas con menores lluvias, de 500 a 1000 mm anuales, se localizan a lo largo del altiplano cundiboyacense y a partir de esta zona las lluvias se incrementan tanto al oriente como al occidente del departamento. Al oriente ocurren las lluvias orográficas correspondientes a la vertiente oriental de la cordillera oriental, registrando valores en el rango de 2000 a más de 5000 mm. Al occidente, en dirección al valle del Magdalena, las lluvias oscilan entre 2000 y 3000 mm al año.

Divición Administrativa de Boyacá:

El departamento de Boyacá está dividido en 123 municipios, los cuales se agrupan en 12 provincias: Sugamuxi, Tundama, Valderrama, Norte, Gutierrez, Centro, Neira, Oriente, Lengupa, Ricaurte, Marquez y Occidente.

DATOS Y METODOS:

Tecnicas de Interpolacion a usar:

Kriging: Es una herramienta de interpolación que basicamente interpola una superficie de ráster a partir de puntos. Kriging presupone que la distancia o la dirección entre los puntos de muestra reflejan una correlación espacial que puede utilizarse para explicar la variación en la superficie. Además ajusta una función matemática a una cantidad especificada de puntos o a todos los puntos dentro de un radio específico para determinar el valor de salida para cada ubicación. Pondera los valores medidos circundantes para calcular una predicción de una ubicación sin mediciones, estas ponderaciones están basadas no solo en la distancia entre los puntos medidos y la ubicación de la predicción, sino también en la disposición espacial general de los puntos medidos. Para utilizar la disposición espacial en las ponderaciones, la correlación espacial debe estar cuantificada. Por lo tanto, en un kriging ordinario, la ponderación depende de un modelo ajustado a los puntos medidos, la distancia a la ubicación de la predicción y las relaciones espaciales entre los valores medidos alrededor de la ubicación de la predicción.

IDW: La interpolación mediante distancia inversa ponderada es útil para determiar los valores de celda mediante una combinación ponderada linealmente de un conjunto de puntos de muestra. La ponderación es una función de la distancia inversa. La superficie que se interpola debe ser la de una variable dependiente de la ubicación. Así mismo, este método presupone que la variable que se representa cartográficamente disminuye su influencia a mayor distancia desde su ubicación de muestra, es decir, asume que los valores más cercanos están más relacionados que otros con su función. IDW es un método de interpolación espacial muy flexible ya que se puede configurar la interpolación de diferentes maneras, pues se puede espicificar el radio de búsqueda y la interpolación sólo utilizará el número de puntos conocidos dentro del radio de búsqueda. También se pueden establecer barreras. Si hay crestas en un perfil de elevación, o barreras de ruido, entonces estos son ejemplos apropiados para usar una barrera. Esta barrera de polilínea impide buscar los puntos de entrada de la muestra.

Polígonos de Thiessen: La interpolación de áreas utiliza el concepto de polígonos de Thiessen para dibujar límites en torno a las áreas que tienen el mismo valor. Este tipo de interpolador se llama interpolador proximal. Esta técnica construye polígonos mediante la asignación por cada punto de un área de interés un valor del punto de muestra que sea más cercano. Crea superficies escalonadas, donde los valores pueden cambiar drásticamente en distancias cortas, es decir, un valor estimado está representado dentro del área de cada polígono. Una vez descargado el archivo CHIRPS, leerlo sin comprimir:

library(tmap)
library(gstat)
library(sp)
library(knitr)
library(rgdal)
library(raster)
library(sf)
library(tidyverse)

precip <- raster("C:/Users/David Perdomo/Desktop/Geomatica/tif/chirps-v2.0.2020.10.6.tif") 

Como ya habrá notado, este archivo representa la precipitación acumulada en los últimos días de octubre de 2020, es decir, la precipitación del 26 al 31 de octubre.

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/David Perdomo/Desktop/Geomatica/tif/chirps-v2.0.2020.10.6.tif 
names      : chirps.v2.0.2020.10.6 

Tenga en cuenta que este es un conjunto de datos con cobertura global. Su CRS es un sistema de coordenadas geográficas.

Carguemos un shapefile que represente nuestra área de interés

(aoi <- shapefile("C:/Users/David Perdomo/Desktop/Geomatica/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 123 
extent      : -74.66496, -71.94885, 4.655196, 7.055557  (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  :         15,      15001,    ALMEIDA,                               1537,   25.35271459,      2017,    BOYACÁ, 0.204497978002, 0.00206924119091 
max values  :         15,      15897, ZETAQUIRÁ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233,      2017,    BOYACÁ,  2.40391520946,   0.123609862522 

Tenga en cuenta que los sistemas de referencia de coordenadas para ambos conjuntos de datos son los mismos.

Ahora, recortemos los datos de precipitación:

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

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

precip.mask 

class      : RasterLayer 
dimensions : 48, 54, 2592  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -74.65, -71.95, 4.649999, 7.049999  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : chirps.v2.0.2020.10.6 
values     : 12.18011, 95.86112  (min, max)

Trace la capa ráster enmascarada:

plot(precip.mask, main= "CHIRPS lluvia en Boyacá 2020 [mm]")
plot(aoi, add=TRUE)

Se puede hacer un mejor mapa con Leaflet usando la función addRasterImage:

library(leaflet)

package 㤼㸱leaflet㤼㸲 was built under R version 4.0.3

library(RColorBrewer)

package 㤼㸱RColorBrewer㤼㸲 was built under R version 4.0.3

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 en Boyacá 2020 [mm]")

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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_1984 in CRS definitionSome values were outside the color scale and will be treated as NA

NA

La 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    : 755 
extent      : -74.625, -71.975, 4.674999, 7.024999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       : chirps.v2.0.2020.10.6 
min values  :      12.1801061630249 
max values  :      95.8611221313477 

names(precip.points) <- "rain"
precip.points

class       : SpatialPointsDataFrame 
features    : 755 
extent      : -74.625, -71.975, 4.674999, 7.024999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :             rain 
min values  : 12.1801061630249 
max values  : 95.8611221313477 

str(precip.points)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 755 obs. of  1 variable:
  .. ..$ rain: num [1:755] 63.1 78.5 84.6 95.9 53.7 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:755, 1:2] -72.2 -72.1 -72.1 -72 -72.3 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] -74.62 4.67 -71.97 7.02
  .. ..- 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"

Grafiquemos los puntos:

plot(precip.mask, main= "CHIRPS precipitaciones 2020 [mm]")
plot(aoi, add=TRUE)

points(precip.points$x, precip.points$y, col = "red", cex = .5)

geojsonio::geojson_write(precip.points, file = "./chirps/ppoints.geojson")

Success! File is at ./chirps/ppoints.geojson

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

Leamos los puntos de precipitación.

precip.points <- geojsonio::geojson_read("./chirps/ppoints.geojson", what="sp")
precip.points

class       : SpatialPointsDataFrame 
features    : 755 
extent      : -74.625, -71.975, 4.674999, 7.024999  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :             rain 
min values  : 12.1801061630249 
max values  : 95.8611221313477 

¿Cuál es el punto de convertir datos de precipitación ráster en datos de precipitación puntual?

Bueno, hay pocas estaciones meteorológicas de la Organización Meteorológica Mundial (OMM) en nuestro país Por lo tanto, CHIRPS puede ser una opción para obtener los datos de precipitación puntual necesarios para obtener una superficie de precipitación continua.

ADICIONAL:

En caso de que se pierda la conexión, podemos volver a leer el shapefile con el área de interés:

(aoi <- shapefile("C:/Users/David Perdomo/Desktop/Geomatica/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 123 
extent      : -74.66496, -71.94885, 4.655196, 7.055557  (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  :         15,      15001,    ALMEIDA,                               1537,   25.35271459,      2017,    BOYACÁ, 0.204497978002, 0.00206924119091 
max values  :         15,      15897, ZETAQUIRÁ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233,      2017,    BOYACÁ,  2.40391520946,   0.123609862522 

Necesitamos convertirlo en una característica espacial:

boya_sf <-  sf::st_as_sf(aoi)
(border_sf <-
  boya_sf %>%
  summarise(area = sum(MPIO_NAREA)))

Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
geographic CRS: WGS 84
      area                       geometry
1 23138.05 POLYGON ((-73.03051 4.98230...

Convierta la característica espacial en un marco de datos espacial:

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

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -74.66496, -71.94885, 4.655196, 7.055557  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 1
names       :           area 
value       : 23138.04812969 

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

Simple feature collection with 123 features and 2 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
geographic CRS: WGS 84
First 10 features:
        MUNIC CODIGO                       geometry
1       TUNJA  15001 POLYGON ((-73.34014 5.58308...
2     ALMEIDA  15022 POLYGON ((-73.36793 5.01349...
3   AQUITANIA  15047 POLYGON ((-72.76242 5.63856...
4    ARCABUCO  15051 POLYGON ((-73.50487 5.84347...
5      BELÉN  15087 POLYGON ((-72.91692 6.08612...
6      BERBEO  15090 POLYGON ((-73.0677 5.27048,...
7  BETÉITIVA  15092 POLYGON ((-72.81796 5.97422...
8     BOAVITA  15097 POLYGON ((-72.64907 6.43640...
9  BOYACÃ\u0081  15104 POLYGON ((-73.34806 5.47411...
10   BRICEÑO  15106 POLYGON ((-73.89118 5.73749...

Ahora, intentaremos una tarea st_intersection. Esto se hace usando la biblioteca sf:

p.sf <- st_as_sf(precip.points)
(precip.sf = st_intersection(boya.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 755 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -74.625 ymin: 4.674999 xmax: -71.975 ymax: 7.024999
geographic CRS: WGS 84
First 10 features:
        MUNIC CODIGO     rain                 geometry
122   CUBARÃ\u0081  15223 63.14175 POINT (-72.175 7.024999)
122.1 CUBARÃ\u0081  15223 78.48558 POINT (-72.125 7.024999)
122.2 CUBARÃ\u0081  15223 84.64105 POINT (-72.075 7.024999)
122.3 CUBARÃ\u0081  15223 95.86112 POINT (-72.025 7.024999)
122.4 CUBARÃ\u0081  15223 53.68439 POINT (-72.275 6.974999)
122.5 CUBARÃ\u0081  15223 60.33939 POINT (-72.225 6.974999)
122.6 CUBARÃ\u0081  15223 59.96722 POINT (-72.175 6.974999)
122.7 CUBARÃ\u0081  15223 62.04287 POINT (-72.125 6.974999)
122.8 CUBARÃ\u0081  15223 73.22034 POINT (-72.075 6.974999)
122.9 CUBARÃ\u0081  15223 84.24580 POINT (-72.025 6.974999)

Dos tareas de reproyección:

p.sf.magna <- st_transform(precip.sf, crs=3116)
boya.sf.magna <- st_transform(boya.sf, crs=3116)
(precip2 <- as(p.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS definition

class       : SpatialPointsDataFrame 
features    : 755 
extent      : 939360, 1232381, 1008778, 1269088  (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  :    ALMEIDA,  15001, 12.1801061630249 
max values  : ZETAQUIRÁ,  15897, 95.8611221313477 

Nuevamente, podemos escribir el conjunto de datos intermedio

shapefile(precip2, filename='./chirps/precip2.shp', overwrite=TRUE)
precip2$rainfall <- round(precip2$rain, 1)
precip2

class       : SpatialPointsDataFrame 
features    : 755 
extent      : 939360, 1232381, 1008778, 1269088  (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  :    ALMEIDA,  15001, 12.1801061630249,     12.2 
max values  : ZETAQUIRÁ,  15897, 95.8611221313477,     95.9 

(boya2 <- as(boya.sf.magna, 'Spatial'))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS definition

class       : SpatialPolygonsDataFrame 
features    : 123 
extent      : 934932.9, 1235253, 1006592, 1272398  (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  :    ALMEIDA,  15001 
max values  : ZETAQUIRÁ,  15897 

shapefile(boya2, filename='./chirps/boya2.shp', overwrite=TRUE)

Asegúrese de que las dos extensiones coincidan:

precip2@bbox <- boya2@bbox

Trazar los datos con tmap:

tm_shape(boya2) + tm_polygons() +
  tm_shape(precip2) +
  tm_dots(col="rainfall", palette = "RdBu", midpoint = 3.0,
             title="Precipitación muestreada \n(en mm)", size=0.2) +
  tm_text("rainfall", just="center", xmod=.6, size = 0.5) +
  tm_legend(legend.outside=TRUE)

CRS object has comment, which is lost in output

INTERPOLACIÓN DE DATOS DE PRECIPITACIÓN

1. polígonos de Thiessen

Los polígonos de Thiessen (o interpolación de proximidad) se pueden crear utilizando la función dirichlet de Spatstat.

library(spatstat.data)

package 㤼㸱spatstat.data㤼㸲 was built under R version 4.0.3

library(nlme)

package 㤼㸱nlme㤼㸲 was built under R version 4.0.3
Attaching package: 㤼㸱nlme㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    collapse

The following object is masked from 㤼㸱package:raster㤼㸲:

    getData

library(rpart)
library(gstat)
library(spatstat)

package 㤼㸱spatstat㤼㸲 was built under R version 4.0.3Registered 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

library(raster)
library(dplyr)

th  <-  as(dirichlet(as.ppp(precip2)), "SpatialPolygons")
crs(th) <- crs(precip2)
crs(boya2) <- 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_defs 

crs(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 

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 SpatialPolygonsDataFrame

th.z     <- over(th, precip2, fn=mean)
th.spdf  <-  SpatialPolygonsDataFrame(th, th.z)
th.clp   <- raster::intersect(boya2,th.spdf)

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)

CRS object has comment, which is lost in output

2. Interpolación ponderada de distancia inversa

Muchos paquetes comparten los mismos nombres de función. Esto puede ser un problema cuando estos paquetes se cargan en una misma sesión de R. Por ejemplo, la función de intersección está disponible en los paquetes base, Spatstat y Ráster, todos los cuales se cargan en esta sesión actual. Para asegurarse de que se selecciona la función adecuada, es una buena idea anteponer el nombre de la función con el nombre del paquete como en raster :: intersect ().

Esta sugerencia se utilizará en el siguiente fragmento de código cuando se llame a la función idw que está disponible tanto en Spatstat como en Gstat.

Tenga en cuenta que la función dirichlet (como la mayoría de las funciones del paquete spatsat) requiere que el objeto de punto esté en formato ppp, de ahí la sintaxis as.ppp (precip2) en línea.

La salida IDW es un ráster. Esto requiere que primero creemos una cuadrícula ráster vacía, luego interpolemos los valores de precipitación en cada celda de la cuadrícula sin muestrear. Se utilizará un valor de potencia IDW de 2 (idp = 2,0).

grd              <- as.data.frame(spsample(precip2, "regular", n=100000))

names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object
proj4string(grd) <- proj4string(precip2)
P.idw <- gstat::idw(rainfall ~ 1, precip2, newdata=grd, idp=2.0)
r       <- raster(P.idw)

r

class      : RasterLayer 
dimensions : 297, 336, 99792  (nrow, ncol, ncell)
resolution : 893.4599, 893.4599  (x, y)
extent     : 935181.8, 1235384, 1006720, 1272077  (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     : 12.85712, 95.76158  (min, max)

boya2

class       : SpatialPolygonsDataFrame 
features    : 123 
extent      : 934932.9, 1235253, 1006592, 1272398  (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  :    ALMEIDA,  15001 
max values  : ZETAQUIRÁ,  15897 

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

r.m 

class      : RasterLayer 
dimensions : 297, 336, 99792  (nrow, ncol, ncell)
resolution : 893.4599, 893.4599  (x, y)
extent     : 935181.8, 1235384, 1006720, 1272077  (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     : 12.85712, 95.76158  (min, max)

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)

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 = "IDW precipitación interpolada en Boyacá 2020 [mm]")

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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_1984 in CRS definition

2.1 Ajuste de 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.

precip2

class       : SpatialPointsDataFrame 
features    : 755 
extent      : 934932.9, 1235253, 1006592, 1272398  (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  :    ALMEIDA,  15001, 12.1801061630249,     12.2 
max values  : ZETAQUIRÁ,  15897, 95.8611221313477,     95.9 

P <- precip2

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

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

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:

sqrt( sum((IDW.out - P$rainfall)^2) / length(P))

[1] 2.573842

2.2 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).

img <- gstat::idw(rainfall~1, P, newdata=grd, idp=2.0)
n   <- length(P)
Zi  <- matrix(nrow = length(img$var1.pred), ncol = n)

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))

  Zi[,i] <- n * img$var1.pred - (n-1) * Z1$var1.pred
}

Zj <- as.matrix(apply(Zi, 1, sum, na.rm=T) / n )

c1 <- apply(Zi,2,'-',Zj)            
c1 <- apply(c1^2, 1, sum, na.rm=T ) 

CI <- sqrt( 1/(n*(n-1)) * c1)

img.sig   <- img
img.sig$v <- CI /img$var1.pred 

r <- raster(img.sig, layer="v")
r.m <- raster::mask(r, boya2)

tm_shape(r.m) + tm_raster(n=7,title="IDW 95% \nintervalo de confianza \n(en mm)") +
  tm_shape(P) + tm_dots(size=0.2) +
  
  tm_legend(legend.outside=TRUE)

2.3 Ajuste polinomial de primer orden

Para ajustar un modelo polinomial de primer orden de la forma precip = intersección + aX + bY

f.1 <- as.formula(rainfall ~ X + Y) 
 

P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]


lm.1 <- lm( f.1, data=P)


dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 

r   <- raster(dat.1st)
r.m <- raster::mask(r, boya2)


tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE, 
            title="polinomio de primer orden \nPrecipitación prevista (en mm)") +
  tm_shape(P) + 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.

2.4 Ajuste polinomial de segundo orden

Para ajustar un modelo polinomial de segundo orden de la forma precip = intersección + aX + bY + dX2 + eY2 + fXY

f.2 <- as.formula(rainfall ~ X + Y + I(X*X)+I(Y*Y) + I(X*Y))


P$X <- coordinates(P)[,1]
P$Y <- coordinates(P)[,2]
 

lm.2 <- lm( f.2, data=P)


dat.2nd <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.2, newdata=grd))) 

r   <- raster(dat.2nd)
r.m <- raster::mask(r, boya2)


tm_shape(r.m) + 
  tm_raster(n=10, palette="RdBu", auto.palette.mapping=FALSE,midpoint = NA,
            title="polinómico de segundo orden \nPrecipitación prevista \n(en mm)") +
  tm_shape(P) + 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.

3. 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 a la función variograma.

f.1 <- as.formula(rainfall ~ X + Y) 


var.smpl <- variogram(f.1, P, cloud = FALSE, cutoff=100000, width=8990)


dat.fit  <- fit.variogram(var.smpl, fit.ranges = TRUE, fit.sills = TRUE,
                          vgm(psill=3, model="Mat", range=150000, nugget=0.0))

dat.fit

plot(var.smpl, dat.fit, xlim=c(0,130000))

3.1 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 residuales 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

f.1 <- as.formula(rainfall ~ X + Y) 

dat.krg <- krige( f.1, P, grd, dat.fit)

r <- raster(dat.krg)
r.m <- raster::mask(r, boya2)

r.m

class      : RasterLayer 
dimensions : 297, 336, 99792  (nrow, ncol, ncell)
resolution : 893.4599, 893.4599  (x, y)
extent     : 935181.8, 1235384, 1006720, 1272077  (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     : 11.93526, 95.82892  (min, max)

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)

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 Boyacá \ en 2020 [mm]")

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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_1984 in CRS definitionSome values were outside the color scale and will be treated as NASome values were outside the color scale and will be treated as NA

3.3 Genere 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, boya2)

tm_shape(r.m) + 
  tm_raster(n=7, palette ="Reds",
            title="Interpolación 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

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

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)

PRESENTACIÓN DE RESULTADOS

1. polígonos de Thiessen

El metodo de interpolacion de Thiessen permite trazar piligonos con una cantidad de puntos los cuales son ponderados posteriormente, en las partes oscuras es donde mas se presenta precipitación en mm, y en las partes claras es donde es menos común

donde se notan mas precipitaciones, de 40 a 100 mm en el departamento, en el Norte, entre el municipio de San Mateo a Cubara hay precipitacion de 40 a 100 mm, en el Sur del departamento entre la frontera con casanare, entre el municipio de Pajarito y Santa Maria, al occidente del departamento entre el municipio de chiquinquira hasta puerto boyacá.

2. Interpolación ponderada de distancia inversa

Con el metodo de interpolacion IDW, tenemos una serie de puntos, se quiere saber cual es la precipitacion en un punto del mapa, el metodo toma los puntos y traza una distancia con la cual mediante una operacion da unos datos representativos, lo cual nos indicaria el valor predicho.

lo que se busca es que la linea roja se encuentre en la linea de refrencia gris para lograr una mejor predicción, se hizo una validacion cruzada, y el margen de error fue bastante pequeño, ya que la linea roja esta bastante cerca de la linea de referencia.Se usó una potencia de 6.0, debido a que es un grado aceptable en el que disminiuye la interacción espacial en función inversa a la distancia.

3. Kriging

La interpolación Kriging es un procedimiento geo estadístico que genera una superficie estimada a partir de un conjunto de puntos dispersos con los valores a analizar.

el semivariograma permite cuantificar la auto correlación ya que grafica la varianza de todos de los datos según la distancia, es decir mide las relaciones estadísticas entre los puntos medidos, permitiendo que relacionar los puntos de muestreo cercano, entre mas cerca estén más relacionadas estaran y tendran una pequeña semivarianza, en cambio si los puntos de muestreo están lejanos estos estaran menos relacionados y por ende tienen una alta semivarianza.

Se grafica un semi-variograma exponencial que permite analizar la información, y se encuentra que los problemas de convergencia de la información no son tan amplios ya que los puntos estan bastante cerca al eje trazado.

este metodo no mostro un area con una superficie uniforme, se logran ver 9 areas ditintas con valores de precipitacion variables, posilemente debido a la faltaa de información

el mapa de Indice de confianza corrobora lo dicho, mostrando valores en casi todo el territorio del departamento de Boyacá de 8 a 12 en el indice de confianza

CONCLUSIONES

La interpolación espacial es el proceso de utilizar puntos con valores conocidos para estimar valores desconocidos en otros puntos. para realizar un mapa de precipitación, se sabe que para el departamento de Boyacá no se encontrarán suficientes estaciones meteorológicas distribuidas uniformemente para cubrir toda la región. La interpolación espacial puede estimar las temperaturas en lugares que no tienen ese dato utilizando lecturas de temperatura conocida en estaciones meteorológicas cercanas

Los metodos mas conliables de intepolación son IDW Y Kriging, la interpolación por polígonos de Thiessen no ofrece datos muy confiables y tiene una alta variacion, por lo que no tiene mucha precisión en sus predicciones.

Referente a los resutados, en geenral se observo que en el departamento de Boyacá, las precipitacion son mas altas en la parte Norte, Occidente y sur, en la parte media del departamento las precipitaciones no menos frecuentes la causa de esta variacion tal vez sea devido a la La Cordillera Oriental la cual atraviesa el departamento de sur a norte, haciendo que el departamento presente terreno quebrados, páramos, valles y altiplanos.

REFERENCIAS

1.Documentación de Qgis2.14[Análisis Espacial (Interpolación)] https://docs.qgis.org/2.14/es/docs/gentle_gis_introduction/spatial_analysis_interpolation.html#

2.Gonçalves João. (2018).Advanced Techniques With Raster Data - Part 3: Exercises.: R-exercises. https://www.r-exercises.com/2018/03/31/advanced-techniques-with-raster-data-part-3-regression-kriging/

3.DEPARTAMENTO DE BOYACÁ, Información general con ciudades y municipios del departamento http://www.colombiamania.com/departamentos/boyaca.html

4.Departamento de Boyacá https://encolombia.com/educacion-cultura/geografia/colombiana/departamentos/boyaca/