15 Operaciones con datos espaciales

15.1 Resumen

Las operaciones espaciales para datos vectoriales incluyen creación de subconjuntos espaciales, unión de datos espaciales, agregación de datos espaciales y relaciones de distancia, entre otras. Por su parte, las operaciones espaciales para datos raster incluyen creación de subconjuntos espaciales y álgebra de mapas, entre otras.

15.2 Trabajo previo

15.2.1 Lecturas

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

15.3 Preparativos

15.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

Datos de ejemplo de Lovelace et al.:

# Instalación de paquete de datos de ejemplo de Lovelace et al.
install.packages("spData")
install.packages("spDataLarge", repos = "https://nowosad.r-universe.dev")

# Carga de paquete de datos de ejemplo de Lovelace et al.
library(spData)
library(spDataLarge)

También se utiliza el paquete raster y con sus métodos:

raster::raster(): para convertir objetos SpatRaster a RasterLayer y así poder desplegarlos con leaflet::addRasterImage(), que solo acepta objetos RasterLayer.
raster::aggregate(): para reducir el tamaño de objetos RasterLayer y así poder desplegarlos con leaflet::addRasterImage(), que no despliega objetos muy grandes.

raster no se carga en este capítulo con la función library(). Sus métodos se llaman con la notación paquete::metodo(). raster requiere la instalación del paquete rgdal.

15.3.2 Conjuntos de datos para ejemplos

15.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). Se transforma a WGS84 para combinarlo más fácilmente con otros conjuntos de datos.

# Lectura y visualización de datos geoespaciales de provincias

# Lectura
provincias <-
  st_read(
    dsn = "datos/ign/delimitacion-territorial-administrativa/provincias.geojson",
    quiet = TRUE
  ) %>%
  st_transform(4326) # transformación a WGS84

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

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

15.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). Se transforma a WGS84 para combinarlo más fácilmente con otros conjuntos de datos.

# Lectura y visualización de datos geoespaciales de cantones

# Lectura
cantones <-
  st_read(
    dsn = "datos/ign/delimitacion-territorial-administrativa/cantones.geojson",
    quiet = TRUE
  ) %>%
  st_transform(4326) # transformación a WGS84

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

15.3.2.3 Aeródromos de Costa Rica

Es un archivo GeoJSON con las geometrías de puntos de los aeródromos de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Instituto Geográfico Nacional (IGN). Se transforma a WGS84 para combinarlo más fácilmente con otros conjuntos de datos.

# Lectura y visualización de datos geoespaciales de aeródromos

# Lectura
aerodromos <-
  st_read(
    dsn = "datos/ign/infraestructura/aerodromos.geojson",
    quiet = TRUE
  ) %>%
  st_transform(4326) # transformación a WGS84

# Visualización en un mapa
plot(
  aerodromos$geometry,
  pch = 16,
  main = "Aeródromos de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

15.3.2.4 Cabeceras de cantones y provincias de Costa Rica

Se origina en un filtro a un archivo GeoJSON con las geometrías de puntos de los poblados de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Instituto Geográfico Nacional (IGN). Se transforma a WGS84 para combinarlo más fácilmente con otros conjuntos de datos.

# Lectura y visualización de datos geoespaciales de cabeceras de cantones y provincias

# Lectura
cabeceras_cantones_provincias <-
  st_read(dsn = "datos/ign/nombres-geograficos/poblados.geojson",
          quiet = TRUE) %>%
  filter(
    tipo == "Cabecera de cantón y distrito" |
      tipo == "Cabecera de provincia - cantón y distrito" |
      tipo == "Capital y cabecera de provincia"
  ) %>%
  st_transform(4326) # transformación a WGS84

# Visualización en un mapa
plot(
  cabeceras_cantones_provincias$geometry,
  pch = 16,
  main = "Cabeceras de cantones y provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

15.3.2.5 Vipéridos de Costa Rica

Es un archivo CSV con registros de presencia de la familia Viperidae (serpientes venenosas) 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 aeródromos

# Lectura
viperidos <-
  st_read(
    "datos/gbif/viperidos.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(viperidos) <- 4326

# Visualización en un mapa
plot(
  viperidos$geometry,
  pch = 16,
  main = "Vipéridos de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

15.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"
)

15.4 Introducción

Esta lección brinda una visión general de las operaciones espaciales para datos vectoriales implementadas en el paquete sf, y para datos raster implementadas en el paquete terra.

15.5 Datos vectoriales

Las operaciones espaciales para datos vectoriales incluyen:

Creación de subconjuntos espaciales (spatial subsetting).
Unión de datos espaciales (spatial joining).
Agregación de datos espaciales (spatial aggregation).
Relaciones de distancia.

Seguidamente, se explicará como maneja estas operaciones el paquete sf.

15.6 Manejo de datos espaciales con el paquete sf

15.6.1 Creación de subconjuntos espaciales

Es el proceso de selección de objetos espaciales con base en su relación con otros objetos espaciales. Estas relaciones se expresan como predicados espaciales, los cuales están implementados como métodos de sf.

La creación de subconjuntos espaciales es análoga a la creación de subconjuntos por datos de atributos. Puede realizarse a través de los operadores [ y $ del paquete base de R o por medio de la función filter() de dplyr.

En los dos ejemplos siguientes, se utiliza el método st_within() para filtrar los puntos contenidos en un polígono.

Primero, se utilizan los operadores del paquete base.

# Selección de la provincia de Limón (por atributos)
limon <- provincias[provincias$provincia == "Limón",]

# Selección de los aeródromos ubicados en Limón (espacial)
aerodromos_limon <- aerodromos[limon, , op = st_within]

Mapa leaflet.

# Mapa leaflet
leaflet() %>%
  addTiles() %>% # capa base de OSM
  addPolygons( # capa de provincias (polígonos)
    data = limon,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%  
  addCircleMarkers( # capa de registros de presencia (puntos)
    data = aerodromos_limon,
    stroke = F,
    radius = 4,
    fillColor = 'brown',
    fillOpacity = 1
  )

El mismo resultado se obtiene con las funciones y operadores de Tidyverse.

# Selección de la provincia de Limón (por atributos)
limon <-
  provincias %>%
  filter(provincia == "Limón")

# Selección de los aeródromos ubicados en Limón (espacial)
aerodromos_limon <-
  aerodromos %>%
  filter(st_within(x = ., y = limon, sparse = FALSE))

En el anterior llamado a filter(), la expresión x = . es equivalente a x = aerodromos. Para una explicación sobre el argumento sparse, por favor lea la sección 4.2.2. del libro “Geocomputation with R” de R. Lovelace et. al..

Mapa leaflet.

# Mapa leaflet
leaflet() %>%
  addTiles() %>% # capa base de OSM
  addPolygons( # capa de provincias (polígonos)
    data = limon,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%  
  addCircleMarkers( # capa de registros de presencia (puntos)
    data = aerodromos_limon,
    stroke = F,
    radius = 4,
    fillColor = 'black',
    fillOpacity = 1
  )

Además de st_within(), sf implementa predicados espaciales como, entre otros, st_contains(), st_intersects() y st_disjoint(), entre otros.

15.6.2 Unión de datos espaciales

La unión “no espacial” de dos conjuntos de datos se basa en uno o varios campos (llamados llaves o keys) que están presentes en ambos conjuntos. Las uniones espaciales se basan en un principio similar pero, en lugar de campos, la relación entre los conjuntos se realiza a través de una operación topológica, a veces llamada spatial overlay. Al igual que con los datos de atributos, la unión espacial, ejecutada con el método st_join(), agrega una o varias columnas al conjunto de datos destino (i.e. el argumento x de la función), provenientes del objeto fuente (i.e. el argumento y). La operación topológica que ejecuta por defecto st_join() es st_intersects().

15.6.2.1 Ejemplos

15.6.2.1.1 Vipéridos de Costa Rica

15.6.2.1.1.1 En provincias

En el siguiente ejemplo, se unen los registros de presencia de vipéridos (geometrías de puntos) con la capa de provincias (geometrías de polígonos), para agregar las columnas de código y nombre de provincia al conjunto de registros de presencia.

# Unión de provincias y vipéridos a través st_join()
viperidos <- 
  viperidos %>%
  st_join(provincias[c("cod_provin", "provincia")])

# Despliegue de los datos unidos
viperidos %>%
  st_drop_geometry() %>%
  dplyr::select(species, stateProvince, provincia, locality) %>%
  datatable(options = list(
    pageLength = 10,
    language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
  ))

15.6.2.1.1.2 Registros de presencia

Se cuenta la cantidad de registros de presencia por provincia y se une al data frame provincias.

# Conteo de registros de presencia por código de provincia
registros_viperidos_x_provincia <-
  viperidos %>%
  st_drop_geometry() %>%
  count(cod_provin, name = "registros_viperidos")

# Unión de cantidad de registros de presencia por provincia a provincias
provincias_viperidos <-
  provincias %>%
  left_join(
    registros_viperidos_x_provincia,
    by = "cod_provin",
    copy = FALSE,
    keep = FALSE
  )

Paleta de colores para los mapas de coropletas.

# Paleta de colores para los mapas
colores_provincias_registros_viperidos <-
  colorNumeric(palette = "Reds",
               domain = provincias_viperidos$registros_viperidos,
               na.color = "transparent")

Mapa de coropletas generado con plot().

# Visualización en un mapa generado con plot()
plot(
  provincias_viperidos["registros_viperidos"],
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  col = colores_provincias_registros_viperidos(provincias_viperidos$registros_viperidos),
  main = "Cantidad de registros de presencia de vipéridos en provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

15.6.2.1.1.3 Especies

Se cuenta la cantidad de especies por provincia y se une al data frame provincias.

# Conteo de especies por código de provincia
especies_viperidos_x_provincia <-
  viperidos %>%
  st_drop_geometry() %>%
  filter(taxonRank == "SPECIES" | taxonRank == "SUBSPECIES") %>% # para excluir identificaciones a género o superiores
  group_by(cod_provin) %>%
  summarise(especies_viperidos = n_distinct(species, na.rm = TRUE))

# Unión de cantidad de registros de presencia por provincia a provincias
provincias_viperidos <-
  provincias_viperidos %>%
  left_join(
    especies_viperidos_x_provincia,
    by = "cod_provin",
    copy = FALSE,
    keep = FALSE
  )

Paleta de colores para los mapas de coropletas.

# Paleta de colores para los mapas
colores_provincias_especies_viperidos <-
  colorNumeric(palette = "Blues",
               domain = provincias_viperidos$especies_viperidos,
               na.color = "transparent")

Mapa de coropletas generado con plot().

# Visualización en un mapa generado con plot()
plot(
  provincias_viperidos["especies_viperidos"],
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  col = colores_provincias_especies_viperidos(provincias_viperidos$especies_viperidos),
  main = "Cantidad de especies de vipéridos en provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

15.6.2.1.1.4 Tabla y mapa interactivos de cantidades de especies y registros de presencia

Tabla DT.

# Visualización en formato de tabla
provincias_viperidos %>%
  st_drop_geometry() %>%
  dplyr::select(provincia, especies_viperidos, registros_viperidos) %>%
  arrange(desc(especies_viperidos)) %>%
  datatable(options = list(
    pageLength = 10,
    language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
  ))

Mapa leaflet.

# Mapa leaflet de vipéridos en provincias
leaflet() %>%
  setView(# centro y nivel inicial de acercamiento
    lng = -84.19452,
    lat = 9.572735,
    zoom = 7) %>%
  addTiles(group = "OpenStreetMap") %>% # capa base
  addPolygons(
    data = provincias_viperidos,
    fillColor = ~ colores_provincias_registros_viperidos(provincias_viperidos$registros_viperidos),
    fillOpacity = 0.8,
    color = "black",
    stroke = TRUE,
    weight = 1.0,
    popup = paste(
      paste(
        "<strong>Provincia:</strong>",
        provincias_viperidos$provincia
      ),
      paste(
        "<strong>Registros:</strong>",
        provincias_viperidos$registros_viperidos
      ),
      sep = '<br/>'
    ),
    group = "Cantidad de registros de vipéridos"
  ) %>%
  addPolygons(
    data = provincias_viperidos,
    fillColor = ~ colores_provincias_especies_viperidos(provincias_viperidos$especies_viperidos),
    fillOpacity = 0.8,
    color = "black",
    stroke = TRUE,
    weight = 1.0,
    popup = paste(
      paste(
        "<strong>Provincia:</strong>",
        provincias_viperidos$provincia
      ),
      paste(
        "<strong>Registros:</strong>",
        provincias_viperidos$especies_viperidos
      ),
      sep = '<br/>'
    ),
    group = "Cantidad de especies de vipéridos"
  ) %>%
  addHeatmap(
    data = viperidos,
    lng = ~decimalLongitude,
    lat = ~decimalLatitude,
    radius = 10,
    blur = 20,
    group = "Mapa de calor"
  ) %>%  
  addCircleMarkers(
    data = viperidos,
    stroke = F,
    radius = 3,
    fillColor = 'black',
    fillOpacity = 1,
    popup = paste(
      viperidos$species,
      viperidos$provincia,
      viperidos$eventDate,
      paste0("<a href='", viperidos$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    clusterOptions = markerClusterOptions(),
    group = "Registros de presencia de vipéridos"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c(
      "Cantidad de registros de vipéridos",
      "Cantidad de especies de vipéridos",
      "Mapa de calor",
      "Registros de presencia de vipéridos"
    )
  ) %>%
  addLegend(
    position = "bottomleft",
    pal = colores_provincias_registros_viperidos,
    values = provincias_viperidos$registros_viperidos,
    group = "Cantidad de registros de vipéridos",
    title = "Cantidad de registros"
  ) %>%
  addLegend(
    position = "bottomleft",
    pal = colores_provincias_especies_viperidos,
    values = provincias_viperidos$especies_viperidos,
    group = "Cantidad de especies de vipéridos",
    title = "Cantidad de especies"
  )  %>%
  addResetMapButton() %>%
  addSearchOSM() %>%
  addMouseCoordinates() %>%
  addScaleBar(position = "bottomright", options = scaleBarOptions(imperial = FALSE)) %>%
  addMiniMap(position = "bottomright") %>%
  hideGroup("Cantidad de especies de vipéridos") %>%
  hideGroup("Mapa de calor")

Algunas consultas:

# Conteo de especies de vipéridos en Puntarenas
viperidos %>%
  st_drop_geometry() %>%
  filter(provincia == "Puntarenas") %>%
  filter(taxonRank == "SPECIES" | taxonRank == "SUBSPECIES") %>%
  count(species, name = "registros") %>%
  rename(especie = species) %>%
  arrange(desc(registros))
#>                      especie registros
#> 1             Bothrops asper       238
#> 2      Bothriechis lateralis       182
#> 3     Bothriechis schlegelii        98
#> 4             Crotalus simus        25
#> 5         Porthidium porrasi        13
#> 6     Lachesis melanocephala         6
#> 7  Metlapilcoatlus mexicanus         6
#> 8   Metlapilcoatlus nummifer         6
#> 9   Bothriechis nigroviridis         4
#> 10 Bothriechis supraciliaris         4
#> 11        Porthidium nasutum         3
#> 12    Porthidium ophryomegas         2
#> 13            Bothrops atrox         1
#> 14      Bothrops lanceolatus         1
#> 15             Lachesis muta         1

# Conteo de especies de vipéridos en Limón
viperidos %>%
  st_drop_geometry() %>%
  filter(provincia == "Limón") %>%
  filter(taxonRank == "SPECIES" | taxonRank == "SUBSPECIES") %>%
  count(species, name = "registros") %>%
  rename(especie = species) %>%
  arrange(desc(registros))
#>                    especie registros
#> 1   Bothriechis schlegelii       366
#> 2           Bothrops asper        90
#> 3       Porthidium nasutum        47
#> 4      Lachesis stenophrys        17
#> 5            Lachesis muta         6
#> 6 Bothriechis nigroviridis         4
#> 7    Bothriechis lateralis         3
#> 8 Metlapilcoatlus nummifer         2
#> 9   Porthidium ophryomegas         1

15.6.3 Agregación de datos espaciales

De manera similar al caso de la agregación de atributos, la agregación espacial es una forma de “condensar” o “resumir” datos. Los datos agregados muestran estadísticas de una variable (ej. promedio, suma) en relación con una variable de agrupación. Esto puede lograrse con el método agregate() de sf o con los métodos group_by() y summarise() de dplyr, utilizados en combinación con st_join().

En el siguiente bloque de código, se utiliza summarise() para mostrar el promedio de altitud de los puntos más altos de Nueva Zelanda (NZ) en cada región del país. Los datos de los puntos más altos están en el objeto nz_height y los de las regiones de Nueva Zelanda en el objeto nz.

# Promedio de altitud de puntos más altos para cada región de NZ
nz_regiones_altitud_promedio_puntos_altos <-
  nz_height %>%
  aggregate(by = nz, FUN = mean)

Los resultados se muestran en un mapa leaflet.

# Paleta de colores
colores_nz_regiones_altitud_promedio_puntos_altos <-
  colorNumeric(palette = "Blues",
               domain = nz_regiones_altitud_promedio_puntos_altos$elevation,
               na.color = "transparent")

# Mapa leaflet
leaflet() %>%
  addTiles(group = "OpenStreetMap") %>% # capa base de OSM
  addPolygons(
    # capa de regiones de NZ (polígonos)
    data = st_transform(nz_regiones_altitud_promedio_puntos_altos, 4326),
    color = "black",
    fillColor = ~ colores_nz_regiones_altitud_promedio_puntos_altos(
      nz_regiones_altitud_promedio_puntos_altos$elevation
    ),
    fillOpacity = 0.8,
    stroke = TRUE,
    weight = 1.0,
    popup = paste(
      "<strong>Altitud de la región:</strong>",
      nz_regiones_altitud_promedio_puntos_altos$elevation,
      "m"
    ),
    group = "Regiones de NZ"
  ) %>%
  addCircleMarkers(
    # capa de puntos altos de NZ
    data = st_transform(nz_height, 4326),
    stroke = F,
    radius = 4,
    fillColor = 'brown',
    fillOpacity = 1,
    popup = paste(
      "<strong>Altitud del punto:</strong>",
      nz_height$elevation,
      "m"
    ),    
    group = "Puntos altos"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c(
      "Regiones de NZ",
      "Puntos altos"
    )
  ) %>%
  addLegend(
    position = "bottomleft",
    pal = colores_nz_regiones_altitud_promedio_puntos_altos,
    values = nz_regiones_altitud_promedio_puntos_altos$elevation,
    group = "Regiones de NZ",
    title = "Altitud promedio"
  )

15.6.4 Relaciones de distancia

La distancia entre dos objetos sf se calcula con el método st_distance(). En el siguiente ejemplo, se calcula la distancia entre los puntos correspondientes a varias cabeceras de provincias y cantones de Costa Rica.

# Cálculo de distancias entre cabeceras de cantones y provincias

cat("Distancia entre San José y La Cruz:",
    round(st_distance(
      filter(cabeceras_cantones_provincias, nombre == "San José"),
      filter(cabeceras_cantones_provincias, nombre == "La Cruz")
    ) / 1000, 2),
    "km",
    "\n")
#> Distancia entre San José y La Cruz: 212.12 km

cat("Distancia entre San José y Neily:",
    round(st_distance(
      filter(cabeceras_cantones_provincias, nombre == "San José"),
      filter(cabeceras_cantones_provincias, nombre == "Neily")
    ) / 1000, 2),
    "km",
    "\n")
#> Distancia entre San José y Neily: 189.5 km

cat("Distancia entre San José y Bribrí:",
    round(st_distance(
      filter(cabeceras_cantones_provincias, nombre == "San José"),
      filter(cabeceras_cantones_provincias, nombre == "Bribrí")
    ) / 1000, 2),
    "km",
    "\n")
#> Distancia entre San José y Bribrí: 140.72 km

cat("Distancia entre Heredia y Alajuela:",
    round(st_distance(
      filter(cabeceras_cantones_provincias, nombre == "Heredia"),
      filter(cabeceras_cantones_provincias, nombre == "Alajuela")
    ) / 1000, 2),
    "km",
    "\n")
#> Distancia entre Heredia y Alajuela: 10.81 km

cat("Distancia entre Alajuela y Cartago:",
    round(st_distance(
      filter(cabeceras_cantones_provincias, nombre == "Alajuela"),
      filter(cabeceras_cantones_provincias, nombre == "Cartago")
    ) / 1000, 2),
    "km",
    "\n")
#> Distancia entre Alajuela y Cartago: 36.23 km

Mapa leaflet de cabeceras de cantones y provincias.

# Mapa de cabeceras de cantones y provincias
leaflet() %>%
  addTiles(group = "OpenStreetMap") %>% # capa base de OSM
  addCircleMarkers(
    # capa de puntos altos de NZ
    data = cabeceras_cantones_provincias,
    stroke = F,
    radius = 4,
    fillColor = 'brown',
    fillOpacity = 1,
    popup = cabeceras_cantones_provincias$nombre,
    group = "Cabeceras de cantones y provincias"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c("Cabeceras de cantones y provincias")
  ) %>%
  addScaleBar(position = "bottomright", options = scaleBarOptions(imperial = FALSE))

15.7 Datos raster

Las operaciones espaciales para datos raster incluyen:

Creación de subconjuntos espaciales (spatial subsetting).
Álgebra de mapas (map algebra).

15.8 Manejo de datos espaciales con el paquete terra

En los siguientes ejemplos, se utilizarán los conjuntos de datos elevacion y grano que se crearon en la lección sobre datos de atributos, los cuales están también incluidos en el paquete spData.

# Creación de objetos SpatRaster
elevacion <-  rast(system.file("raster/elev.tif", package = "spData"))
grano <- rast(system.file("raster/grain.tif", package = "spData"))

# Mapeo de los conjuntos de datos de ejemplo

# Elevación
plot(elevacion)


# Tipos de granos
plot(grano)

También se utilizará la capa de altitud de Costa Rica.

15.8.1 Creación de subconjuntos espaciales

En la lección sobre operaciones con atributos, se explicó como recuperar subconjuntos de objetos raster, ya sea mediante su ID o su posición en filas y columnas. Los subconjuntos de objetos raster también pueden obtenerse mediante la especificación de coordenadas o de otros objetos espaciales.

# La función cellFromXY() retorna el ID de la celda correspondiente a una coordenada
id <- cellFromXY(elevacion, xy = matrix(c(0.1, 0.1), ncol = 2))
elevacion[id]
#>   elev
#> 1   16

# El mismo resultado puede obtenerse con terra::extract()
terra::extract(elevacion, data.frame(x = 0.1, y = 0.1))
#>   ID elev
#> 1  1   16

El siguiente bloque de código retorna el valor en la capa de altitud de Costa Rica correspondiente a varios pares de coordenadas (x, y).

# Altitud del punto en (-84, 10)
altitud[cellFromXY(altitud, xy = matrix(c(-84, 10), ncol = 2))]
#>   altitud
#> 1    1451

# Altitud del Cerro Chirripó
altitud[cellFromXY(altitud, xy = matrix(c(-83.488667, 9.484083), ncol = 2))]
#>   altitud
#> 1    3664

# Altitud de la Catedral Metropolitana de San José
altitud[cellFromXY(altitud, xy = matrix(c(-84.078758, 9.932684), ncol = 2))]
#>   altitud
#> 1    1156

También es posible consultar las celdas contenidas en la extensión (i.e. los límites) otro raster.

# Objeto raster para hacer un recorte ("clip") de otro raster
clip <- rast(
  xmin = 0.9,
  xmax = 1.8,
  ymin = -0.45,
  ymax = 0.45,
  resolution = 0.3,
  vals = rep(1, 9)
)

# Celdas de elev contenidas en la extensión de clip
elevacion[clip]
#>      elev
#> [1,]   18
#> [2,]   24

Los métodos explicados anteriormente solo retornan ID y valores de celdas. Con el operador [ y el argumento drop = FALSE pueden retornarse objetos raster.

# Objeto raster creado a partir de un rango de ID de celdas
r1 <- elevacion[1:2, drop = FALSE]
plot(r1)


# Objeto raster creado a partir de posiciones de filas y columnas
r2 <- elevacion[1:3, 1:3, drop = FALSE]
plot(r2)

El método crop() recorta un objeto SpatRaster de acuerdo con el contorno de otro objeto espacial, raster o vectorial. El siguiente bloque de código recorta la capa de altitud de Costa Rica con base en un objeto raster ubicado en el centro del país.

# Creación de un raster ubicado en el centro del país, alrededor de (-84, 10)
clip_centro_cr <-
  rast(
    xmin = -84.10,
    xmax = -83.90,
    ymin = 9.90,
    ymax = 10.10,
    res = 0.10
  )

# Recorte de la capa de altitud con base en el raster del centro del país
altitud_centro_cr <- crop(altitud, clip_centro_cr)
plot(altitud_centro_cr)

Seguidamente, la misma capa de altitud se recorta siguiendo el contorno de una capa vectorial de la clase SpatVector de terra, correspondiente a la provincia de Heredia. El polígono de Heredia debe convertirse a SpatVector mediante el método vect().

# Polígono de la provincia de Heredia
heredia <-
  provincias %>%
  filter(provincia == "Heredia")

# Recorte de la capa de altitud
altitud_heredia <-
  altitud %>%
  crop(vect(heredia)) %>%
  mask(vect(heredia))

plot(altitud_heredia)

El siguiente mapa leaflet muestra la capa de altitud de Costa Rica y los recortes efectuados en esta.

# 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 de Costa Rica"
  ) %>%
  addRasterImage(
    raster::raster(altitud_heredia),
    colors = colores_altitud,
    opacity = 0.8,
    group = "Altitud de Heredia"
  ) %>%
  addRasterImage(
    raster::raster(altitud_centro_cr),
    colors = colores_altitud,
    opacity = 0.8,
    group = "Altitud del centro de Costa Rica"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap"),
    overlayGroups = c(
      "Altitud de Costa Rica",
      "Altitud de Heredia",
      "Altitud del centro de Costa Rica"
    )
  ) %>%
  addLegend(
    title = "Altitud de Costa Rica",
    values = values(altitud),
    pal = colores_altitud,
    position = "bottomleft",
    group = "Altitud de Costa Rica"
  ) %>%
  addLegend(
    title = "Altitud de Heredia",
    values = values(altitud_heredia),
    pal = colores_altitud,
    position = "bottomleft",
    group = "Altitud de Heredia"
  ) %>%
  addLegend(
    title = "Altitud del centro de Costa Rica",
    values = values(altitud_centro_cr),
    pal = colores_altitud,
    position = "bottomright",
    group = "Altitud del centro de Costa Rica"
  )  %>%
  hideGroup("Altitud de Heredia") %>%
  hideGroup("Altitud del centro de Costa Rica")

También es posible obtener subconjuntos raster mediante la aplicación en un objeto raster de una “máscara” (mask) con la misma extensión y resolución, y que contenga valores lógicos o NA. Esta operación puede realizarse con el método mask().

# Creación de una "máscara"
mascara <- elevacion 
values(mascara) <- sample(c(NA, TRUE), 36, replace = TRUE)
plot(mascara)


# Creación de subconjuntos espaciales mediante la máscara
elevacion_mascara <- mask(elevacion, mascara)                   
plot(elevacion_mascara)

15.8.2 Álgebra de mapas

El álgebra de mapas divide las operaciones raster en cuatro clases (Lovelace et al.):

Operaciones locales o de “celda por celda”.
Operaciones focales o de “vecindario” (neighborhood). Generalmente, el valor de salida de cada celda proviene de un bloque de entrada de 3 x 3 celdas.
Operaciones zonales. Son similares a las focales, pero el bloque de entrada puede tener tamaños y formas irregulares.
Operaciones globales. Los valores de salida de cada celda provienen de uno o varios objetos raster completos.

Esta clasificación se basa en la cantidad o forma de las celdas utilizadas por cada pixel durante el procesamiento. Otras clasificaciones pueden estar basadas en el área de aplicación del análisis (ej. terreno, hidrología, teledetección).

15.8.2.1 Operaciones locales

Son operaciones realizadas celda por celda en una o varias capas raster. Por ejemplo, la reclasificación de una capa mediante el método classify().

# Reclasificación de una capa raster
rcl <-  matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
recl <- classify(elevacion, rcl = rcl)
plot(recl)

También operaciones aritméticas, similares a las que se realizan con matrices algebraicas.

# Álgebra raster con operadores aritméticos y de comparación

# Operación aritmética
elevacion_doble <- elevacion + elevacion
plot(elevacion_doble)


# Operación aritmética
elevacion_cuadrado <- elevacion * elevacion
plot(elevacion_cuadrado)


# Expresión lógica de comparación
elevacion_mayor_30 <- elev > 30
plot(elevacion_mayor_30)

El siguiente bloque de código calcula el Índice de vegetación de diferencia normalizada (NDVI, en inglés, Normalized difference vegetation index) en una imagen Landsat. El NDVI se utiliza para estimar la cantidad, calidad y desarrollo de la vegetación con base a la medición de la intensidad de la radiación de ciertas bandas del espectro electromagnético que la vegetación emite o refleja. El valor del NDVI varía entre -1.0 y +1.0.

# Cálculo del NDVI

# Imagen Landsat de 4 bandas (4 = red, 3 = infraroja cercana)
archivo_imagen_landsat = system.file("raster/landsat.tif", package = "spDataLarge")
imagen_landsat = rast(archivo_imagen_landsat)

plot(imagen_landsat)


# Función para calcular el NDVI (red = banda roja, nir = banda infraroja cercana)
ndvi <- function(nir, red) {
  (nir - red) / (nir + red)
}

# Se aplica ndvi() sobre todas las celdas de la capa raster
raster_ndvi = lapp(imagen_landsat[[c(4, 3)]], fun = ndvi)

plot(raster_ndvi)

El resultado se aprecia mejor en un mapa leaflet. Nótese el uso del método raster::aggregate() para reducir el tamaño del objeto raster_ndvi, el cual puede ser muy grande para desplegarse con leaflet::addRasterImage().

# Paleta de colores de NDVI
colores_ndvi <-
  colorNumeric(rev(grDevices::terrain.colors(25)),
               values(raster_ndvi),
               na.color = "transparent")

# Mapa leaflet
leaflet() %>%
  addTiles(group = "OpenStreetMap") %>%
  addRasterImage(
    raster::aggregate(raster::raster(raster_ndvi), fact = 2, fun = mean),
    colors = colores_ndvi,
    opacity = 0.8,
    group = "NDVI"
  ) %>%
  addLayersControl(baseGroups = c("OpenStreetMap"),
                   overlayGroups = c("NDVI")) %>%
  addLegend(
    title = "NDVI",
    values = values(raster_ndvi),
    pal = colores_ndvi,
    position = "bottomleft",
    group = "NDVI"
  )

15.8.2.2 Operaciones focales

En este tipo de operaciones, el valor de salida de cada pixel procesado depende de un bloque compuesto por una celda central y sus vecinas. Este “vecindario” (también llamado kernel, filtro o “ventana móvil”) es típicamente de 3 x 3 celdas, pero puede tomar otras formas y tamaños. Una operación focal aplica una función de agregación (ej. promedio, mínimo, máximo) a todas las celdas del vecindario, utiliza la salida como nuevo valor de la celda central correspondiente y se mueve a la celda siguiente.

# Uso de focal() para obtener el valor mínimo en un vecindario de 3 x 3
minimo_ventana <- focal(elevacion, w = matrix(1, nrow = 3, ncol = 3), fun = min)
plot(minimo_ventana)

Las operaciones focales tienen aplicaciones en áreas como procesamiento de imágenes (ej. remoción de valores extremos) o análisis de terreno (ej. cálculo de pendiente o dirección de flujo).

15.8.2.3 Operaciones zonales

De manera similar a las focales, las operaciones zonales aplicación de agregación a múltiples celdas raster. Sin embargo, en el caso de las zonales, generalmente se usa un raster categórico que define las zonas, a diferencia de la ventana predefinida que se emplea en las focales. Por lo tanto, las celdas que definen el filtro zonal no deben ser necesariamente vecinas.

# Uso de zonal() para encontrar la elevación promedio de cada tipo de grano
zonal(elevacion, grano, fun = "mean") %>%
  as.data.frame()
#>   grain     elev
#> 1  clay 14.80000
#> 2  silt 21.15385
#> 3  sand 18.69231

15.8.2.4 Operaciones globales

Las operaciones globales pueden considerarse un caso particular de las operaciones zonales, en las cuales un raster completo corresponde a una zona. Las operaciones globales más comunes incluyen estadísticas descriptivas para todo un conjunto raster.

# Resumen
terra::summary(altitud)
#> Warning: [summary] used a sample
#>     altitud      
#>  Min.   :   1.0  
#>  1st Qu.:  63.0  
#>  Median : 254.0  
#>  Mean   : 561.3  
#>  3rd Qu.: 852.0  
#>  Max.   :3628.0  
#>  NA's   :83907

# Histograma
terra::hist(altitud)

Otras operaciones globales incluyen cálculos de distancia y rasters de “peso” (weight rasters) (ej. distancia de cada celda a una celda objetivo, costo en combustible de cada celda a una celda objetivo).

14 Operaciones con datos de atributos

16 Operaciones con geometrías