Operaciones con geometrías

Este documento está estructurado con base en el capítulo 5 de Lovelace, R., Nowosad, J. & Muenchow, J. (2020). Geocomputation with R.

Trabajo previo

Lea el capítulo 5 de Lovelace, R., Nowosad, J. & Muenchow, J. (2020). Geocomputation with R.

Preparativos

Carga de paquetes

# Paquete para manipulación de datos
library(dplyr)

# Paquete para manejo de datos vectoriales
library(sf)

# Paquete para manejo de datos raster
library(terra)

# Paquete para simplificación y edición de geometrías
library(rmapshaper)

# Paquetes con datos geoespaciales para ejemplos
library(spData)
library(spDataLarge)

Conjuntos de datos

Provincias de Costa Rica (Instituto Geográfico Nacional)
# Carga de datos
provincias <-
  st_read(
    "https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/main/datos/ign/delimitacion-territorial-administrativa/provincias-simplificadas_100m.geojson",
    quiet = TRUE
  )

# Mapeo
plot(
  provincias$geometry,
  main = "Provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

Cantones de Costa Rica (Instituto Geográfico Nacional)
# Carga de datos
cantones <-
  st_read(
    "https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/main/datos/ign/delimitacion-territorial-administrativa/cantones-simplificadas_100m.geojson",
    quiet = TRUE
  )

# Mapeo
plot(
  cantones$geometry,
  main = "Cantones de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

Áreas silvestres protegidas (Sistema Nacional de Áreas de Conservación)
# Carga de datos
asp <-
  st_read(
    "https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/main/datos/sinac/areas-silvestres-protegidas-simplificadas_100m.geojson",
    quiet = TRUE
  )

# Mapeo
plot(
  asp$geometry,
  main = "Áreas silvestres protegidas (ASP) de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

Red vial (Instituto Geográfico Nacional)
# Carga de datos
red_vial <-
  st_read(
    "https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/main/datos/ign/infraestructura/redvial-simplificadas_500m.geojson",
    quiet = TRUE
  )

# Mapeo
plot(
  provincias$geometry,
  main = "Red vial de Costa Rica",
  axes = TRUE,
  graticule = TRUE,
  reset = FALSE
)
plot(
  add = TRUE,
  red_vial$geometry,
  col = "blue"
)

Altitud (WorldClim) 
# Carga de datos
altitud <-
  rast(
    "/vsicurl/https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/master/datos/worldclim/altitud.tif"
  )

# Mapeo
plot(altitud)

Mamíferos (datos agrupados por la Infraestructura Mundial de Información en Biodiversidad)
# Carga de datos
mamiferos <-
  st_read(
    "/vsicurl/https://raw.githubusercontent.com/tpb728O-programaciongeoespacialr/2021ii/main/datos/gbif/mammalia-cr-registros.csv",
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude",
      "Y_POSSIBLE_NAMES=decimalLatitude"
    ),    
    quiet = TRUE
  )

# Asignación de CRS
st_crs(mamiferos) = 4326

# Mapeo
plot(
  provincias$geometry,
  main = "Mamíferos de Costa Rica",
  axes = TRUE,
  graticule = TRUE,
  reset = FALSE
)
plot(
  add = TRUE,
  mamiferos$geometry,
  pch = 16,
  col = "brown"
)

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.

Datos vectoriales

Las operaciones con geometrías en datos vectoriales incluyen:

Operaciones con geometrías con el paquete sf

Estas operaciones modifican las geometrías de objetos vectoriales (sf).

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

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

# Simplificación sin preservación de topología
provincias_simp <-
  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_simp$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_simp_topo <-
  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_simp_topo$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)

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

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

68280 bytes

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

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

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

# Tamaño de la capa simplificada con rmapshaper::ms_simplify()
object.size(provincias_simp)

36328 bytes

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 una con forma de rosca). 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”.

# Costa Rica y sus centroides calculados con st_centroid() y st_point_on_surface()
plot(
  st_union(provincias),
  main = "Centroides de CR: st_centroid (rojo) y st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

plot(st_centroid(st_union(provincias)),
     add = TRUE,
     pch = 16,
     col = "red")

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 478661.3 1102731
# WGS84
st_coordinates(st_transform(st_centroid(st_union(provincias)), crs = 4326))

          X        Y
1 -84.19463 9.972706
# Coordenadas del centroide calculado con st_point_on_surface()
# CRTM05
st_coordinates(st_point_on_surface(st_union(provincias)))

         X       Y
1 539364.2 1065091
# WGS84
st_coordinates(st_transform(st_point_on_surface(st_union(provincias)), crs = 4326))

          X        Y
1 -83.64133 9.632234
# Provincias de Costa Rica y sus centroides calculados con st_centroid() y st_point_on_surface()
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)

plot(st_centroid(provincias),
     add = TRUE,
     pch = 16,
     col = "red")

plot(
  st_point_on_surface(provincias),
  add = TRUE,
  pch = 16,
  col = "green")

# Ruta 32 y sus centroides calculados con st_centroid() y st_point_on_surface()
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = 280000, xmax = 660000, ymin = 880000, ymax= 1250000)),  
  main = "Centroides de la ruta 32: st_centroid (rojo) y st_point_on_surface (verde)",
  axes = TRUE,
  graticule = TRUE)

ruta_32 <-
  red_vial %>%
  filter(num_ruta == "32")

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

plot(
  st_centroid(st_union(ruta_32)),
  add = TRUE,
  pch = 16,
  col = "red")

plot(
  st_point_on_surface(st_union(ruta_32)),
  add = TRUE,
  pch = 16,
  col = "green")

Á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 que rodea la ruta 32
plot(
  st_buffer(st_union(ruta_32), 5000),
  main = "Buffer que rodea la ruta 32",
  axes = TRUE,
  graticule = TRUE)

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

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?”.

# Registros de presencia de mamíferos no voladores ubicados alrededor de la ruta 32
mamiferos <-
  mamiferos %>%
  filter(order != "Chiroptera") %>% # se excluye el orden de los murciélagos
  st_transform(crs = 5367)

buffer_ruta_32 <-
  ruta_32 %>%
  st_buffer(dist = 5000)

mamiferos_buffer_ruta_32 <-
  st_join(mamiferos, buffer_ruta_32) %>%
  filter(!is.na(codigo))

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
)

# 10 especies con mayor cantidad de registros
mamiferos_buffer_ruta_32 %>% 
  st_drop_geometry() %>%
  filter(!is.na(species) & species != "") %>%
  group_by(species) %>% 
  summarise(registros = n()) %>%
  arrange(desc(registros)) %>%
  slice(1:10)

# A tibble: 10 × 2
   species                  registros
   <chr>                        <int>
 1 Sciurus variegatoides           83
 2 Bradypus variegatus             81
 3 Heteromys desmarestianus        48
 4 Choloepus hoffmanni             46
 5 Peromyscus mexicanus            38
 6 Melanomys caliginosus           30
 7 Alouatta palliata               20
 8 Orthogeomys cherriei            20
 9 Didelphis marsupialis           16
10 Cebus capucinus                 15

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.

# Recorte de la sección del Parque Internacional La Amistad (PILA) ubicada en Puntarenas
pila <-
  asp %>%
  filter(nombre_asp == "Internacional La Amistad")

puntarenas_limon <-
  provincias %>%
  filter(provincia == "Puntarenas" | provincia == "Limón")

# Mapa de Puntarenas, Limón y el PILA
plot(
  puntarenas_limon$geometry,
  main = "Puntarenas, Limón y el PILA",
  extent = st_bbox(c(xmin = 350000, xmax = 660000, ymin = 880000, ymax= 1250000)),  
  axes = TRUE,
  graticule = TRUE)

plot(
  pila$geometry,
  border = "green",
  add = TRUE)

puntarenas <-
  provincias %>%
  filter(provincia == "Puntarenas")

# Recorte de la sección del PILA ubicada en Puntarenas
pila_puntarenas <- st_intersection(pila, puntarenas)

# Mapa de la sección recortada
plot(
  pila_puntarenas$geometry,
  main = "Sección del PILA ubicada en Puntarenas",
  col = "red",
  axes = TRUE,
  graticule = TRUE)

plot(
  puntarenas_limon$geometry,
  add = TRUE
)

Uniones de geometrías

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

# Cantones de la provincia de San José
cantones_sanjose <-
  cantones %>%
  filter(provincia == "San José")

plot(
  cantones_sanjose$geometry,
  main = "Cantones de San José",
  axes = TRUE,
  graticule = TRUE)  

# Cantones de la provincia de San José unificados
cantones_sanjose_unificados <- 
  st_union(cantones_sanjose)

plot(
  cantones_sanjose_unificados,
  main = "Cantones de San José unificados",
  axes = TRUE,
  graticule = TRUE)

Datos raster

Las operaciones con geometrías en datos raster incluyen:

Operaciones con geometrías con el paquete terra

Intersecciones geométricas

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

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

clip <- rast(
  xmin = 0.9,
  xmax = 1.8,
  ymin = -0.45,
  ymax = 0.45,
  resolution = 0.3,
  vals = rep(1, 9)
)

# Objeto recortado
elev_clipped <- crop(elev, clip)

# Metadatos del objeto recortado
elev_clipped

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. : +proj=longlat +datum=WGS84 +no_defs 
source      : memory 
name        : elev 
min value   :   18 
max value   :   24 
# Mapeo
plot(elev)

plot(elev_clipped)

plot(elev)
plot(elev_clipped, add = TRUE)

Agregación y desagregación

La resolución de un raster puede disminuirse con la función aggregate() o aumentarse con la función disaggregate().

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.

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

# Agrupación de la capa de altitud de Costa Rica
altitud_agregada <- 
  altitud %>%
  aggregate(fact = 8, fun = mean)

# Capa de altitud agrupada de Costa Rica
plot(
  altitud_agregada,
  main = "Capa agrupada 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. : +proj=longlat +datum=WGS84 +no_defs 
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. : +proj=longlat +datum=WGS84 +no_defs 
source      : memory 
name        :  altitud 
min value   : 6.296875 
max value   : 3139.156 
# Tamaño de la capa original de altitud
object.size(altitud)

1304 bytes
# Tamaño de la capa agrupada de altitud
object.size(altitud_agregada)

1304 bytes

La función disaggregate() 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, que es 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 de elevación
elev_desagregada <- 
  elev %>%
  disaggregate(fact = 2, method = "bilinear")

plot(elev_desagregada)

# Metadatos de la capa original de elevación
elev

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. : +proj=longlat +datum=WGS84 +no_defs 
source      : elev.tif 
name        : elev 
min value   :    1 
max value   :   36 
# Metadatos de la capa desagrupada de elevación
elev_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. : +proj=longlat +datum=WGS84 +no_defs 
source      : memory 
name        : elev 
min value   :    1 
max value   :   36 
# Tamaño de la capa original de elevación
object.size(elev)

1304 bytes
# Metadatos de la capa desagrupada de elevación
object.size(elev_desagregada)

1304 bytes

Interacciones raster-vector

Recorte (cropping) de datos raster

Los métodos crop() y mask() pueden utilizarse para recortar (crop) un objeto raster con base en el contorno de un objeto vectorial.

# Polígono de la provincia de Limón
limon <-
  provincias %>%
  filter(provincia == "Limón") %>%
  st_transform(4326)

# crop() + mask() a capa de altitud
altitud_limon <-
  altitud %>%
  crop(vect(limon)) %>%
  mask(vect(limon))

plot(
  altitud_limon,
  main = "crop() + mask() a capa de altitud de Costa Rica",
  axes = TRUE)

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. Source code is available at https://github.com/tpb728O-programaciongeoespacialr/2021ii/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".