El ecosistema espacial de R

Preparativos

Carga de paquetes

library(sf)
library(rgdal)
library(raster)
library(ggplot2)
library(leaflet)
library(dplyr)

Conjuntos de datos utilizados

Introducción

La comunidad de programadores de R ha desarrollado un conjunto de paquetes para el manejo de datos geoespaciales, tanto en formatos vectoriales como raster. Algunos de los principales de estos paquetes son:

Algunos paquetes de graficación, como ggplot2, también cuentan con algunas capacidades para procesamiento de datos geoespaciales.

Datos vectoriales

El modelo vectorial

El modelo vectorial de datos está basado en puntos localizados en un sistema de referencia de coordenadas (CRS). Los puntos individuales pueden representar objetos independientes (ej. la localización de un poste eléctrico o de una cabina telefónica) o pueden también agruparse para formar geometrías más complejas como líneas o polígonos. Por lo general, los puntos tienen solo dos dimensiones (x, y), a las que se les puede agregar una tercera dimensión z, usualmente correspondiente a la altitud sobre el nivel del mar.

El estándar Simple Features

Simple Features (o Simple Feature Access) es un estándar abierto de la Organización Internacional de Estandarización (ISO) y del Open Geospatial Consortium (OGC) que especifica un modelo común de almacenamiento y acceso para geometrías de dos dimensiones (líneas, polígonos, multilíneas, multipolígonos, etc.). El estándar es implementado por muchas bibliotecas y bases de datos geoespaciales como sf, GDAL, PostgreSQL/PostGIS, SQLite/SpatiaLite, Oracle Spatial y Microsoft SQL Server, entre muchas otras.

La especificación define 17 tipos de geometrías, de las cuales siete son las más comúnmente utilizadas. Estas últimas se muestran en la figura 1.

Figura 1. Tipos de geometrías de Simple Features más usadas. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#vector-data)

El paquete sf

El paquete sf (de Simple Features) de R implementa los modelos de datos de las geometrías de tipo vectorial: puntos, líneas, polígonos, sus versiones múltiples y las colecciones de geometrías. Está basado en bibliotecas de sofware ampliamente utilizadas en aplicaciones geoespaciales:

sf provee acceso, desde un mismo paquete de R, a la funcionalidad de estas tres bibliotecas, proporcionando así una interfaz unificada para leer y escribir datos geoespaciales mediante GDAL, realizar operaciones con geometrías mediante GEOS y efectuar transformaciones entre sistemas de coordenadas mediante PROJ.

En sf, los conjuntos de datos geoespaciales se almacenan en un data frame que contiene una columna especial para las geometrías. Esta columna se denomina generalmente geom o geometry. El manejo de datos geoespaciales como data frames, permite manipularlos con las funciones ya desarrolladas para data frames y con la misma forma de referenciar las filas (observaciones) y las columnas (variables).

Ejemplos de uso de funciones del paquete sf

La función st_read() permite leer una fuenta de datos vectoriales (ej. archivo, base de datos) y recuperarlos en un objeto sf. Este tipo de objetos extiende los data frames con una columna de geometrías.

# Lectura de una capa vectorial mediante st_read()
provincias <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/ign/delimitacion-territorial-administrativa/cr_provincias_simp_wgs84.geojson",
    quiet = TRUE
  )

La función plot() muestra un objeto sf en un mapa.

# Mapeo de un objeto sf mediante plot()
plot(provincias$geometry)

# Mapeo de un objeto sf con argumentos adicionales de plot()
plot(
  provincias$geometry,
  extent = extent(-86,-82.3, 8, 11.3),
  main = "Provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

st_read() puede leer datos de varios formatos.

# Lectura de una archivo CSV con columnas de coordenadas mediante st_read()
manigordos <-
  st_read(
    "https://raw.githubusercontent.com/gf0604-procesamientodatosgeograficos/2021i-datos/main/gbif/leopardus_pardalis-cr-registros.csv",
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude",
      "Y_POSSIBLE_NAMES=decimalLatitude"
    ),
    quiet = TRUE
  )

# Asignación del sistema de coordenadas
st_crs(manigordos) = 4326

# Mapeo
plot(
  manigordos$geometry,
  extent = extent(-86,-82.3, 8, 11.3),
  main = "Registros de Leopardus pardalis (manigordo) en Costa Rica",
  col = "orange",
  axes = TRUE,
  graticule = TRUE
)

Los argumentos reset y add de la función plot() permiten construir un mapa con varias capas.

# Primera capa del mapa
plot(
  provincias$geometry,
  extent = extent(-86,-82.3, 8, 11.3),
  main = "Distribución de Leopardus pardalis (manigordo) en Costa Rica",
  reset = FALSE,
  axes = TRUE,
  graticule = TRUE
)

# Segunda capa
plot(manigordos,
     add = TRUE,
     col = "orange")

La función st_write() guarda en el disco un objeto sf. Pueden utilizarse los formatos vectoriales de GDAL.

# Especificación del directorio de trabajo (debe utilizarse una ruta existente)
setwd("C:/Users/mfvargas/Downloads")

# Se guarda la capa de provincias
provincias %>%
  st_write("provincias.shp")

# Se guarda la capa de registros de manigordos
manigordos %>%
  st_write("manigordos.kml")

Mapeo de objetos sf con otros paquetes

ggplot2

ggplot(data = provincias) +
  geom_sf() +
  geom_point(
    data = manigordos,
    aes(x = decimalLongitude, y = decimalLatitude),
    size = 2,
    col = "orange",
    fill = "orange"
  ) +
  coord_sf(xlim = c(-86, -82.3),
           ylim = c(8, 11.3),
           expand = FALSE) +
  ggtitle("Distribución de Leopardus pardalis (manigordo) en Costa Rica") +
  xlab("Longitud") +
  ylab("Latitud")

Leaflet

leaflet() %>%
  addTiles() %>%
  addPolygons(
    data = provincias,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%
  addCircleMarkers(
    data = manigordos,
    stroke = F,
    radius = 4,
    fillColor = 'orange',
    fillOpacity = 1
  )

Datos raster

El modelo raster

El modelo de datos raster usualmente consiste de un encabezado y de una matriz con celdas (también llamadas pixeles) de un mismo tamaño. El encabezado define el sistema de referencia de coordenadas (CRS), la extensión y el punto de origen de una capa raster. Por lo general, el origen se ubica en la esquina inferior izquierda o en la esquina superior izquierda de la matriz. La extensión se define mediante el número de filas, el número de columnas y el tamaño (resolución) de la celda.

Cada celda tiene una identificación (ID) y almacena un único valor, el cual puede ser numérico o categórico, como se muestra en la figura 2.

Figura 2. El modelo raster: (A) ID de las celdas, (B) valores de las celdas, (C) mapa raster de colores. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#raster-data)

A diferencia del modelo vectorial, el modelo raster no necesita almacenar todas las coordenadas de cada geometría (i.e. las esquinas de las celdas), debido a que la ubicación de cada celda puede calcularse a partir de la información contenida en el encabezado. Esta simplicidad, en conjunto con el álgebra de mapas, permiten que el procesamiento de datos raster sea mucho más eficiente que el procesamiento de datos vectoriales. Por otra parte, el modelo vectorial es mucho más flexible en cuanto a las posibilidades de representación de geometrías y almacenamiento de valores, por medio de múltiples elementos de datos.

Los mapas raster generalmente almacenan fenómenos continuos como elevación, precipitación, temperatura, densidad de población y datos espectrales. También es posible representar mediante raster datos discretos, tales como tipos de suelo o clases de cobertura de la tierra, como se muestra en la figura 3.

Figura 3. Ejemplos de mapas raster continuos y categóricos. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#raster-data)

El paquete raster

El paquete raster proporciona funciones para la lectura, escritura, manipulación, análisis y modelado de datos raster. Por su parte, el paquete rgdal provee enlaces a las bibliotecas GDAL y PROJ.

El paquete raster define tres tipos de datos principales:

Ejemplos de uso de funciones del paquete raster

Para ejemplificar el uso del paquete raster, se accederá a la base de datos climáticos WorldClim, mediante la función getData().

# Especificación del directorio de trabajo (debe utilizarse una ruta existente)
setwd("c:/users/mfvargas/")

# Consulta del directorio de trabajo
getwd()

[1] "c:/users/mfvargas"
# Obtención de la capa de altitud
alt <- getData(
  "worldclim",
  var = "alt",
  res = .5,
  lon = -84,
  lat = 10
)

# Capa de altitud recortada para los límites aproximados de Costa Rica
altitud <- crop(alt, extent(-86, -82.3, 8, 11.3))

# Mapeo
plot(altitud, main = "Altitud recortada para los límites aproximados de Costa Rica")

# Capa de altitud recortada para los límites exactos de Costa Rica
altitud <-
  alt %>%
  crop(provincias) %>%
  mask(provincias)

# Mapeo
plot(altitud, main = "Altitud recortada para los límites exactos de Costa Rica")

# Resumen de información básica de la capa raster
altitud

class      : RasterLayer 
dimensions : 686, 546, 374556  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -87.1, -82.55, 5.5, 11.21667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : alt_23 
values     : -10, 3732  (min, max)
# Clase de la capa raster
class(altitud)

[1] "RasterLayer"
attr(,"package")
[1] "raster"
# Cantidad de capas
nlayers(altitud)

[1] 1
# Obtención de la capa de precipitación
prec <-
  getData(
    "worldclim",
    var = "prec",
    res = .5,
    lon = -84,
    lat = 10
  )

# Capa de precipitación recortada para los límites aproximados de Costa Rica
precipitacion <- crop(prec, extent(-86, -82.3, 8, 11.3))

# Mapeo
plot(precipitacion, main = "Precipitación recortada para los límites aproximados de Costa Rica")

# Capa de precipitación recortada para los límites exactos de Costa Rica
precipitacion <-
  prec %>%
  crop(provincias) %>%
  mask(provincias)

# Mapeo
plot(precipitacion, main = "Precipitación recortada para los límites exactos de Costa Rica")

# Resumen de información básica de la capa raster
precipitacion

class      : RasterBrick 
dimensions : 686, 546, 374556, 12  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -87.1, -82.55, 5.5, 11.21667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : prec1_23, prec2_23, prec3_23, prec4_23, prec5_23, prec6_23, prec7_23, prec8_23, prec9_23, prec10_23, prec11_23, prec12_23 
min values :        0,        0,        0,        9,      116,      179,      124,      156,       90,       140,        79,         7 
max values :      500,      288,      283,      372,      618,      659,      807,      673,      790,       829,       801,       752 
# Clase de la capa raster
class(precipitacion)

[1] "RasterBrick"
attr(,"package")
[1] "raster"
# Cantidad de capas
nlayers(precipitacion)

[1] 12
# Lista e información sobre las capas
unlist(precipitacion)

class      : RasterBrick 
dimensions : 686, 546, 374556, 12  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -87.1, -82.55, 5.5, 11.21667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : prec1_23, prec2_23, prec3_23, prec4_23, prec5_23, prec6_23, prec7_23, prec8_23, prec9_23, prec10_23, prec11_23, prec12_23 
min values :        0,        0,        0,        9,      116,      179,      124,      156,       90,       140,        79,         7 
max values :      500,      288,      283,      372,      618,      659,      807,      673,      790,       829,       801,       752 
# Mapeo de altitud y registros de presencia de manigordos
plot(
  altitud,
  extent = extent(-86,-82.3, 8, 11.3),
  axes = TRUE,
  graticule = TRUE,
  reset = FALSE
)
plot(manigordos, col = "black", add = TRUE)

Mapeo de objetos raster con otros paquetes

Leaflet

# Paleta de colores
pal <- colorNumeric(
  c("#006400", "#FFFF00", "#0000FF"), 
  values(altitud), 
  na.color = "transparent"
)

leaflet() %>%
  addTiles() %>%
  addRasterImage(
    altitud, 
    colors = pal, 
    opacity = 0.8
  ) %>%  
  addPolygons(
    data = provincias,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%
  addCircleMarkers(
    data = manigordos,
    stroke = F,
    radius = 4,
    fillColor = 'orange',
    fillOpacity = 1
  )

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/gf0604-procesamientodatosgeograficos/2021i-leccion-15/, 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 ...".