16 Operaciones con geometrías

16.1 Resumen

Las operaciones con geometrías para datos vectoriales incluyen simplificación, creación de centroides, creación de áreas de amortiguamiento (buffers), recortes (clipping) y uniones de geometrías, entre otras. Por su parte, las operaciones con geometrías para datos raster incluyen intersecciones geométricas, agregación y desagregación, entre otras.

16.2 Trabajo previo

16.2.1 Lecturas

Lovelace, R., Nowosad, J., & Münchow, J. (2019). Geocomputation with R (capítulo 5). CRC Press. https://geocompr.robinlovelace.net/

16.3 Preparativos

16.3.1 Carga de paquetes

# Carga de paquetes

library(dplyr) # transformación de datos
library(sf) # manejo de datos vectoriales
library(terra) # manejo de datos raster
library(DT) # tablas interactivas
library(leaflet) # mapas interactivos
library(leaflet.extras) # funciones adicionales de leaflet
library(leafem) # funciones adicionales de leaflet
library(spData) # datos de ejemplo
library(spDataLarge) # datos de ejemplo

Paquete rmapshaper, para la edición y simplificación de geometrías:

# Instalación de rmapshaper
install.packages("rmapshaper")

# Carga de rmapshaper
library(rmapshaper)

16.3.2 Conjuntos de datos para ejemplos

16.3.2.1 Provincias de Costa Rica

Es un archivo GeoJSON con los polígonos de las provincias de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Instituto Geográfico Nacional (IGN).

# Lectura y visualización de datos geoespaciales de provincias

# Lectura
provincias <-
  st_read(
    dsn = "datos/ign/delimitacion-territorial-administrativa/provincias.geojson",
    quiet = TRUE
  )

# Visualización en un mapa
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),
  main = "Provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

16.3.2.2 Cantones de Costa Rica

Es un archivo GeoJSON con los polígonos de los cantones de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Instituto Geográfico Nacional (IGN).

# Lectura y visualización de datos geoespaciales de cantones

# Lectura
cantones <-
  st_read(
    dsn = "datos/ign/delimitacion-territorial-administrativa/cantones.geojson",
    quiet = TRUE
  )

# Transformación
cantones <-
  cantones %>%
  st_transform(5367) %>%
  st_simplify(dTolerance = 100) %>% # simplificación de geometrías
  st_transform(4326)

# Visualización en un mapa
plot(
  cantones$geometry,
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  main = "Cantones de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

16.3.2.3 Áreas silvestres protegidas de Costa Rica

Es un archivo GeoJSON con los polígonos de las áreas silvestres protegidas (ASP) de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Sistema Nacional de Áreas de Conservación (Sinac).

# Lectura y visualización de datos geoespaciales de ASP

# Lectura
asp <-
  st_read(
    dsn = "datos/sinac/areas-silvestres-protegidas.geojson",
    quiet = TRUE
  )

# Visualización en un mapa
plot(
  asp$geometry,
  main = "Áreas silvestres protegidas de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

16.3.2.4 Red vial de Costa Rica

Es un archivo GeoJSON con las líneas de la red vial de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Instituto Geográfico Nacional (IGN).

# Lectura y visualización de datos geoespaciales de la red vial

# Lectura
red_vial <-
  st_read(
    dsn = "datos/ign/infraestructura/redvial.geojson",
    quiet = TRUE
  )

# Visualización en un mapa
plot(
  red_vial$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),
  main = "Red vial de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

16.3.2.5 Mamíferos de Costa Rica

Es un archivo CSV con registros de presencia de la clase Mammalia (mamíferos) de Costa Rica. Este archivo proviene de una consulta al portal de datos de la Infraestructura Mundial de Información en Biodiversidad (GBIF).

# Lectura y visualización de datos geoespaciales de mamíferos

# Lectura
mamiferos <-
  st_read(
    "datos/gbif/mamiferos.csv",
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude", # columna de longitud decimal
      "Y_POSSIBLE_NAMES=decimalLatitude"   # columna de latitud decimal
    ),
    quiet = TRUE
  )

# Asignación del CRS WGS84
st_crs(mamiferos) <- 4326

# Visualización en un mapa
plot(
  mamiferos$geometry,
  pch = 16,
  main = "Mamíferos de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

16.3.2.6 Altitud de Costa Rica

Es un archivo GeoTIFF correspondiente a la altitud de Costa Rica, en resolución de 30 x 30 segundos. Este archivo proviene de WorldClim, un conjunto de capas climáticas disponibles en varias resoluciones espaciales.

# Lectura y visualización de datos geoespaciales de altitud de Costa Rica

# Lectura
altitud <-
  rast(
    "datos/worldclim/altitud.tif"
  )

# Visualización en un mapa
terra::plot(
  altitud,
  main = "Altitud de Costa Rica"
)

16.4 Introducción

Esta lección brinda una visión general de las operaciones con geometrías en datos vectoriales implementadas en el paquete sf y en datos raster implementadas en el paquete terra. Estas operaciones trabajan con la columna de geometrías (ej. geometry, geom) del paquete sf, para el caso de los datos vectoriales, y con la localización geográfica de los pixeles para el caso de los datos raster. En la sección final, se muestran varias operaciones de interacción entre los modelos raster y vectorial.

16.5 Datos vectoriales

Las operaciones con geometrías en datos vectoriales incluyen:

Simplificación.
Centroides.
Áreas de amortiguamiento (buffers).
Recortes (clipping).
Uniones de geometrías.

16.6 Operaciones con geometrías con el paquete sf

Estas operaciones modifican las geometrías de objetos vectoriales.

16.6.1 Simplificación

La simplificación puede realizarse en geometrías de líneas y polígonos. Reduce la cantidad de memoria, disco y ancho de banda que utilizan las geometrías. Para simplificar geometrías, el paquete sf incluye el método st_simplify(), basado en el algoritmo de Douglas-Peucker, el cual recibe el argumento dTolerance para controlar el nivel de generalización de las unidades del mapa. Este argumento se expresa en las unidades de medida del CRS de la capa, por lo que es conveniente utilizar un CRS con unidades de medida de distancias (ej. metros).

El siguiente bloque de código simplifica la capa de provincias, primero sin preservar su topología y luego preservándola.

# Simplificación de la capa de provincias

# Simplificación sin preservación de topología
provincias_simplificado <-
  provincias %>%
  st_simplify(dTolerance = 5000, preserveTopology = FALSE)

# Mapa de la capa de provincias con simplificación y sin preservación de topología
plot(
  provincias_simplificado$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),  
  main = "Provincias simplificadas sin preservación de topología",
  axes = TRUE,
  graticule = TRUE)


# Simplificación con preservación de topología
provincias_simplificado_topologia <-
  provincias %>%
  st_simplify(dTolerance = 5000, preserveTopology = TRUE)

# Mapa de la capa de provincias con simplificación y con preservación de topología
plot(
  provincias_simplificado_topologia$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),  
  main = "Provincias simplificadas con preservación de topología",
  axes = TRUE,
  graticule = TRUE)


# Tamaño de la capa original
object.size(provincias)
#> 12144488 bytes

# Tamaño de la capa simplificada sin preservación de topología
object.size(provincias_simplificado)
#> 18720 bytes

# Tamaño de la capa simplificada con preservación de topología
object.size(provincias_simplificado_topologia)
#> 70608 bytes

La función ms_simplify() del paquete `rmapshaper`` proporciona un método alternativo para la simplificación de geometrías, el cual preserva la topología.

# Simplificación con ms_simplify()
provincias_simplificado_rmapshaper <-
  provincias %>%
  ms_simplify(keep = 0.1, keep_shapes = TRUE)

# Mapa de la capa de provincias con simplificación mediante ms_simplify()
plot(
  provincias_simplificado_rmapshaper$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),  
  main = "Provincias simplificadas con ms_simplify()",
  axes = TRUE,
  graticule = TRUE)


# Tamaño de la capa simplificada con ms_simplify()
object.size(provincias_simplificado_rmapshaper)
#> 1447832 bytes

16.6.2 Centroides

Un centroide es un punto que identifica el centro de un objeto geográfico. Puede calcularse para geometrías de líneas y de polígonos y se utilizan para brindar una representación simplificada de geometrías más complejas. Existen varios métodos para calcularlos.

El paquete sf incluye la función st_centroid(), la cual calcula el centroide geográfico (comúnmente llamado “el centroide”). Es posible que el centroide geográfico se ubique fuera de la geometría “padre” (ej. en el caso de una geometría con forma de anillo). Para evitar este resultado, la función st_point_on_surface() se asegura de que el centroide esté siempre dentro de la geometría “padre”.

El siguiente bloque de código calcula los centroides para Costa Rica, mediante las dos funciones mencionadas.

# Costa Rica y sus centroides calculados con st_centroid() y st_point_on_surface()

# Mapa de provincias
plot(
  st_union(provincias), # unión de los polígonos de provincias
  main = "Centroides de CR: st_centroid (rojo) y st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

# Mapa del centroide calculado con st_centroid()
plot(st_centroid(st_union(provincias)),
     add = TRUE,
     pch = 16,
     col = "red")

# Mapa del centroide calculado con st_point_on_surface()
plot(
  st_point_on_surface(st_union(provincias)),
  add = TRUE,
  pch = 16,
  col = "green")


# Coordenadas del centroide calculado con st_centroid()
# CRTM05
st_coordinates(st_centroid(st_union(provincias)))
#>          X       Y
#> 1 478674.5 1102735
# WGS84
st_coordinates(st_transform(st_centroid(st_union(provincias)), crs = 4326))
#>           X        Y
#> 1 -84.19451 9.972732

# Coordenadas del centroide calculado con st_point_on_surface()
# CRTM05
st_coordinates(st_point_on_surface(st_union(provincias)))
#>          X       Y
#> 1 539373.5 1065147
# WGS84
st_coordinates(st_transform(st_point_on_surface(st_union(provincias)), crs = 4326))
#>           X        Y
#> 1 -83.64124 9.632735

El siguiente bloque de código calcula los centroides de las provincias de Costa Rica, mediante las dos funciones mencionadas.

# Provincias de Costa Rica y sus centroides calculados con st_centroid() y st_point_on_surface()

# Mapa de provincias
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),  
  main = "Centroides de provincias: st_centroid (rojo) y st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

# Mapa de los centroides calculados con st_centroid()
plot(st_centroid(provincias),
     add = TRUE,
     pch = 16,
     col = "red")

# Mapa de los centroides calculados con st_point_on_surface()
plot(
  st_point_on_surface(provincias),
  add = TRUE,
  pch = 16,
  col = "green")

El siguiente bloque de código calcula los centroides para la ruta 32, mediante las dos funciones mencionadas.

# Ruta 32 y sus centroides calculados con st_centroid() y st_point_on_surface()

# Polígonos de San José, Heredia y Limón
sanjose_heredia_limon <-
  provincias %>%
  filter(provincia == "San José" | provincia == "Heredia" | provincia == "Limón")

# Línea de la ruta 32
ruta_32 <-
  red_vial %>%
  filter(num_ruta == "32")

# Mapa de San José, Heredia y Limón
plot(
  sanjose_heredia_limon$geometry,
  main = "Centroides de la ruta 32: st_centroid (rojo) y st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

# Mapa de la ruta 32
plot(
  ruta_32$geometry,
  add = TRUE,
  lwd = 2,
  col = "blue")

# Mapa del centroide calculado con st_centroid()
plot(
  st_centroid(st_union(ruta_32)),
  add = TRUE,
  pch = 16,
  col = "red")

# Mapa del centroide calculado con st_point_on_surface()
plot(
  st_point_on_surface(st_union(ruta_32)),
  add = TRUE,
  pch = 16,
  col = "green")

16.6.3 Áreas de amortiguamiento (buffers)

Los buffers son polígonos creados alrededor de otra geometría, ya sea otro polígono, una línea o un punto. El paquete sf incluye la función st_buffer() para la generación de buffers.

# Buffer alrededor de la ruta 32

# Buffer que rodea la ruta 32
plot(
  st_buffer(st_union(ruta_32), 5000),
  main = "Buffer que rodea la ruta 32",
  axes = TRUE,
  graticule = TRUE)

# Línea de la ruta 32
plot(
  ruta_32$geometry,
  col = "blue",
  add = TRUE
)

16.6.3.1 Ejemplos

16.6.3.1.1 Especies de mamíferos en riesgo de atropello en las cercanías de la ruta 32

Es común el uso de buffers en análisis de datos, para responder preguntas como, por ejemplo, “¿cuántos puntos hay alrededor de una línea?” o “¿cuáles especies pueden encontrarse en las márgenes de un río?”. En este ejemplo, se utiliza un buffer para identificar las especies de mamíferos en riesgo de ser atropellados en las cercanías de la ruta 32.

# Registros de presencia de mamíferos no voladores ubicados alrededor de la ruta 32

# Registros de presencia de mamíferos no voladores
mamiferos_no_voladores <-
  mamiferos %>%
  filter(taxonRank == "SPECIES" | taxonRank == "SUBSPECIES") %>% # para excluir identificaciones a género o superiores
  filter(order != "Chiroptera") # se excluyen los murciélagos

# Línea de la ruta 32
ruta_32 <-
  red_vial %>%
  filter(num_ruta == "32") %>%
  st_transform(4326)

# Buffer de la ruta 32
buffer_ruta_32 <-
  ruta_32 %>%
  st_buffer(dist = 5000) %>%
  st_transform(4326)

# Registros de presencia dentro del buffer
mamiferos_buffer_ruta_32 <-
  st_join(mamiferos_no_voladores, buffer_ruta_32) %>%
  filter(!is.na(codigo))

# Mapa
plot(
  st_union(buffer_ruta_32),
  main = "Mamíferos terrestres alrededor de la ruta 32",
  axes = TRUE,
  graticule = TRUE
)

plot(ruta_32$geometry,
     col = "blue",
     add = TRUE)

plot(
  mamiferos_buffer_ruta_32,
  pch = 16,
  col = "orange",
  add = TRUE
)

Lista de especies y cantidad de registros de presencia:

# Lista de especies
lista_especies <-
  mamiferos_buffer_ruta_32 %>%
  st_drop_geometry() %>%
  filter(!is.na(species) & species != "") %>%
  group_by(species) %>%
  summarise(registros = n()) %>%
  arrange(desc(registros)) %>%
  rename(especie = species)

# Tabla
lista_especies %>%
  datatable(options = list(
    pageLength = 10,
    language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
  ))

Mapa leaflet:

# Mapa de mamíferos cerca de la ruta 32
leaflet() %>%
  addTiles(group = "OpenStreetMap") %>%
  addPolygons(data = st_union(buffer_ruta_32),
              group = "Buffer") %>%
  addHeatmap(
    data = mamiferos_buffer_ruta_32,
    lng = ~ decimalLongitude,
    lat = ~ decimalLatitude,
    radius = 10,
    blur = 20,
    group = "Mapa de calor"
  ) %>%
  addPolylines(data = ruta_32,
               group = "Ruta 32") %>%
  addCircleMarkers(
    data = mamiferos_buffer_ruta_32,
    radius = 1,
    color = "black",
    popup = paste(
      mamiferos_buffer_ruta_32$species,
      paste0(
        "<a href='",
        mamiferos_buffer_ruta_32$occurrenceID,
        "'>Más información</a>"
      ),
      sep = '<br/>'
    ),
    clusterOptions = markerClusterOptions(),
    group = "Registros de presencia"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c("Buffer", "Mapa de calor", "Ruta 32", "Registros de presencia")
  ) %>%
  addScaleBar(position = "bottomright", options = scaleBarOptions(imperial = FALSE))

16.6.4 Recortes (clipping)

El recorte de una geometría con base en la forma de otra puede realizarse con el método st_intersection(), el cual retorna la intersección entre dos geometrías.

En el siguiente bloque de código en R, se recorta la parte del Parque Nacional Braulio Carrillo ubicada en la provincia de San José.

# Recorte de la sección del PN Braulio Carrillo ubicada en San José

# Polígono del PN Braulio Carrillo
braulio_carrillo <-
  asp %>%
  filter(nombre_asp == "Braulio Carrillo")

# Polígono de la provincia de San José
san_jose <-
  provincias %>%
  filter(provincia == "San José")

# Recorte
braulio_carrillo_san_jose <-
  st_intersection(braulio_carrillo, san_jose)
#> Warning: attribute variables are assumed to be spatially
#> constant throughout all geometries

# Mapa
plot(
  braulio_carrillo$geometry,
  main = "Sección del PN Braulio Carrillo ubicada en San José",
  axes = TRUE,
  graticule = TRUE
)

plot(san_jose$geometry,
     add = TRUE)

plot(
  braulio_carrillo_san_jose$geometry,
  col = "orange",
  add = TRUE
)

Mapa leaflet. Por ser muy “pesadas”, las geometrías se simplifican con st_simplify(). Luego se reproyectan con st_transform() para convertirlas al CRS WGS84.

# Mapa leaflet del recorte de la sección del PN Braulio Carrillo ubicada en San José

leaflet() %>%
  addTiles(group = "OpenStreetMap") %>%
  addPolygons(
    data = st_transform(st_simplify(braulio_carrillo, dTolerance = 100), 4326),
    color = "red",
    weight = 1.0,
    group = "PN Braulio Carrillo"
  ) %>%
  addPolygons(
    data = st_transform(st_simplify(san_jose, dTolerance = 100), 4326),
    color = "blue",
    weight = 1.0,
    group = "San José"
  ) %>%
  addPolygons(
    data = st_transform(st_simplify(braulio_carrillo_san_jose, dTolerance = 100), 4326),
    color = "magenta",
    weight = 1.0,
    group = "Intersección"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c("PN Braulio Carrillo", "San José", "Intersección")
  )

16.6.5 Uniones de geometrías

En lecciones anteriores, se ha mostrado como agregar geometrías mediante los métodos aggregate() y summarise(). Internamente, ambos métodos utilizan la función st_union() para combinar las geometrías y disolver sus límites.

En el siguiente bloque de código, se unen los cantones de Cartago en una geometría.

# Unión de los cantones de Cartago

# Cantones de la provincia de Cartago
cantones_cartago <-
  cantones %>%
  filter(provincia == "Cartago")

plot(
  cantones_cartago$geometry,
  main = "Cantones de Cartago",
  axes = TRUE,
  graticule = TRUE)


# Cantones de la provincia de Cartago unificados
cantones_cartago_unificados <- 
  st_union(cantones_cartago)

plot(
  cantones_cartago_unificados,
  main = "Cantones de Cartago unificados",
  axes = TRUE,
  graticule = TRUE)

Mapa leaflet. Las geometrías se reproyectan con st_transform() para convertirlas al CRS WGS84.

# Mapa leaflet de la unión de los cantones de Cartago

leaflet() %>%
  addTiles(group = "OpenStreetMap") %>%
  addPolygons(
    data = st_transform(cantones_cartago_unificados, 4326),
    color = "blue",
    fillColor = "blue",
    opacity = 1.0,
    weight = 1.0,
    group = "Unión de los cantones de Cartago"
  ) %>%
  addPolygons(
    data = st_transform(cantones_cartago, 4326),
    color = "blue",
    fillColor = "trasparent",
    weight = 1.0,
    group = "Cantones de Cartago"
  ) %>%  
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c("Unión de los cantones de Cartago", "Cantones de Cartago")
  )

16.7 Datos raster

Las operaciones con geometrías en datos raster incluyen:

Intersecciones geométricas.
Agregación y desagregación.

16.8 Operaciones con geometrías con el paquete terra

16.8.1 Intersecciones geométricas

Los objetos SpatRaster pueden intersecarse con las funciones intersect() y crop().

# Objeto a recortar
elevacion <- rast(system.file("raster/elev.tif", package = "spData"))

# Objeto según el cual se hace el recorte
clip <- rast(
  xmin = 0.9,
  xmax = 1.8,
  ymin = -0.45,
  ymax = 0.45,
  resolution = 0.3,
  vals = rep(1, 9)
)

# Objeto recortado
elevacion_recortado <- crop(elevacion, clip)

# Metadatos del objeto recortado
elevacion_recortado
#> class       : SpatRaster 
#> dimensions  : 2, 1, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : 1, 1.5, -0.5, 0.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : elev 
#> min value   :   18 
#> max value   :   24

# Mapeo
plot(elevacion)

plot(elevacion_recortado)

16.8.2 Agregación y desagregación

La resolución de un raster puede disminuirse con la función aggregate(), para agregación, o aumentarse con la función disagg(), para desagregación.

En el siguiente bloque de código, se utiliza la función aggregate() para disminuir la resolución del raster de altitud de Costa Rica por un factor de 8, como se especifica con el argumento factor = 8. El argumento fun = mean indica que en el raster agrupado, cada celda será el promedio de las cuatro celdas correspondientes en el raster original.

# Agregación de la capa de altitud de Costa Rica

# Capa de altitud
plot(
  altitud,
  main = "Capa de altitud de Costa Rica",
  axes = TRUE)


# Agregación
altitud_agregada <- 
  altitud %>%
  aggregate(fact = 8, fun = mean)

# Capa agregada
plot(
  altitud_agregada,
  main = "Capa agregada de altitud de Costa Rica",
  axes = TRUE)


# Metadatos de la capa original de altitud
altitud
#> class       : SpatRaster 
#> dimensions  : 687, 546, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : -87.10189, -82.55189, 5.494651, 11.21965  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : altitud.tif 
#> name        : altitud

# Metadatos de la capa agrupada de altitud
altitud_agregada
#> class       : SpatRaster 
#> dimensions  : 86, 69, 1  (nrow, ncol, nlyr)
#> resolution  : 0.06666667, 0.06666667  (x, y)
#> extent      : -87.10189, -82.50189, 5.486317, 11.21965  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        :  altitud 
#> min value   : 6.296875 
#> max value   : 3139.156

La capa de altitud y su capa agregada se muestran en un mapa leaflet.

# Mapa leaflet de agregación de la capa de altitud

# Paleta de colores de altitud de Costa Rica
colores_altitud <-
  colorNumeric(terrain.colors(25),
               values(altitud),
               na.color = "transparent")

# Mapa leaflet
leaflet() %>%
  setView(# centro y nivel inicial de acercamiento
    lng = -84.19452,
    lat = 9.572735,
    zoom = 7) %>%
  addTiles(group = "OpenStreetMap") %>%
  addRasterImage(
    raster::raster(altitud),
    colors = colores_altitud,
    opacity = 0.8,
    group = "Altitud"
  ) %>%
  addRasterImage(
    raster::raster(altitud_agregada),
    colors = colores_altitud,
    opacity = 0.4,
    group = "Altitud agregada"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c(
      "Altitud",
      "Altitud agregada"
    )
  ) %>%
  addLegend(
    title = "Altitud",
    values = values(altitud),
    pal = colores_altitud,
    position = "bottomleft",
    group = "Altitud"
  ) %>%
  addLegend(
    title = "Altitud agregada",
    values = values(altitud_agregada),
    pal = colores_altitud,
    position = "bottomleft",
    group = "Altitud agregada"
  )

La función disagg() genera varias celdas por cada celda del raster original. Debe especificarse un método con el argumento method. Su valor puede ser "", para copiar los valores de la celda de entrada o bilinear, un método de interpolación.

En el siguiente bloque de código, se desagrega el raster de elevación por un factor de dos.

# Desagregación de la capa raster de elevación

# Desagregación de la capa de elevación
elevacion_desagregada <- 
  elevacion %>%
  disagg(fact = 2, method = "bilinear")

plot(elevacion_desagregada,
  main = "Capa desagregada de elevación",
  axes = TRUE
)


# Metadatos de la capa original de elevación
elevacion
#> class       : SpatRaster 
#> dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
#> resolution  : 0.5, 0.5  (x, y)
#> extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : elev.tif 
#> name        : elev 
#> min value   :    1 
#> max value   :   36

# Metadatos de la capa desagregada de elevación
elevacion_desagregada
#> class       : SpatRaster 
#> dimensions  : 12, 12, 1  (nrow, ncol, nlyr)
#> resolution  : 0.25, 0.25  (x, y)
#> extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : elev 
#> min value   :    1 
#> max value   :   36

15 Operaciones con datos espaciales

17 Desarrollo de tableros de control