10 Introducción al manejo de datos geoespaciales con R
10.1 Trabajo previo
10.1.1 Lecturas
Lovelace, R., Nowosad, J., & Münchow, J. (2019). Geocomputation with R (capítulos 1, 2 y 9). CRC Press. https://geocompr.robinlovelace.net/
10.2 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:
El paquete sf. Ofrece un conjunto de funciones para el manejo de datos vectoriales, de acuerdo con el estándar Simple Features.
El paquete terra. Implementa un conjunto de funciones para el manejo de datos raster. Es una reimplementación del paquete raster.
El paquete tmap. Se utiliza para programar mapas estáticos e interactivos, especialmente mapas temáticos como mapas de coropletas y mapas de burbujas.
Algunos paquetes de graficación estadística, como ggplot2 y plotly, también cuentan con capacidades para visualización de datos geoespaciales.
En CRAN Task View: Analysis of Spatial Data, puede encontrarse un resumen detallado de los paquetes de R con funciones geoespaciales.
10.3 Datos vectoriales
10.3.1 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. postes eléctricos, cabinas telefónicas) o pueden también agruparse para formar geometrías más complejas como líneas (ej. ríos, caminos) o polígonos (ej. fincas, países, provincias). 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.
10.3.2 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, Fiona, GDAL, PostgreSQL/PostGIS, SQLite/SpatiaLite, Oracle Spatial y Microsoft SQL Server, entre muchas otras.
La especificación define 18 tipos de geometrías, de las cuales siete son las más comúnmente utilizadas. Estas últimas se muestran en la Figura 10.1.
10.3.3 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:
- Geospatial Data Abstraction Library (GDAL): es una biblioteca para leer y escribir datos geoespaciales en varios formatos raster y vectoriales. Implementa un único modelo abstracto de datos raster y un único modelo abstracto de datos vectoriales, lo que permite programar aplicaciones geoespaciales sin tener que ocuparse de las particularidades de cada formato que se utilice (GeoTIFF, NetCDF, ESRI Shapefile, GeoJSON, etc.). A pesar de que GDAL está programada en C/C++, cuenta con una interfaz de programación de aplicaciones (API) para varios lenguajes de programación, incluyendo C, C++, Python y Java. Además, ofrece un conjunto de utilitarios de línea de comandos cuyas distribuciones binarias están disponibles para varios sistemas operativos, incluyendo Windows, macOS y Linux.
- Geometry Engine, Open Source (GEOS): es una implementación en C++ de la biblioteca JTS Topology Suite (desarrollada en Java) y que implementa un conjunto de operaciones y predicados geoespaciales (ej. unión, intersección, distancia, área).
- PROJ: es una biblioteca que transforma coordenadas entre diferentes CRS, incluyendo tanto proyecciones cartográficas como transformaciones geodésicas.
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 objetos de una clase también llamada sf, los cuales son data frames que contiene una columna especial para las geometrías. Esta columna se denomina generalmente geom o geometry (aunque pueden tener cualquier otro nombre). El manejo de datos geoespaciales como data frames permite manipularlos con las funciones ya desarrolladas para este tipo de datos y con la misma forma de referenciar las filas (observaciones) y las columnas (variables).
10.3.3.1 Instalación y carga
# Instalación de sf
install.packages("sf")# Carga de sf
library(sf)10.3.3.2 Métodos
La lista de métodos (i.e. funciones) de la clase sf puede obtenerse a través de la función methods():
# Métodos de la clase sf
methods(class = "sf") [1] [ [[<-
[3] [<- $<-
[5] aggregate as.data.frame
[7] cbind coerce
[9] dbDataType dbWriteTable
[11] duplicated identify
[13] initialize merge
[15] plot points
[17] print rbind
[19] show slotsFromS3
[21] st_agr st_agr<-
[23] st_area st_as_s2
[25] st_as_sf st_as_sfc
[27] st_bbox st_boundary
[29] st_break_antimeridian st_buffer
[31] st_cast st_centroid
[33] st_collection_extract st_concave_hull
[35] st_convex_hull st_coordinates
[37] st_crop st_crs
[39] st_crs<- st_difference
[41] st_drop_geometry st_exterior_ring
[43] st_filter st_geometry
[45] st_geometry<- st_inscribed_circle
[47] st_interpolate_aw st_intersection
[49] st_intersects st_is_full
[51] st_is_valid st_is
[53] st_join st_line_merge
[55] st_m_range st_make_valid
[57] st_minimum_rotated_rectangle st_nearest_points
[59] st_node st_normalize
[61] st_point_on_surface st_polygonize
[63] st_precision st_reverse
[65] st_sample st_segmentize
[67] st_set_precision st_shift_longitude
[69] st_simplify st_snap
[71] st_sym_difference st_transform
[73] st_triangulate_constrained st_triangulate
[75] st_union st_voronoi
[77] st_wrap_dateline st_write
[79] st_z_range st_zm
[81] text transform
see '?methods' for accessing help and source codeSeguidamente, se describen y ejemplifican algunos de los métodos básicos de la clase sf.
10.3.3.2.1 st_read() - lectura de datos
El método st_read() lee datos vectoriales de una fuente en formato geoespacial (ej. shapefiles, archivos GeoJSON, bases de datos geoespaciales) y los recupera en un objeto sf.
En el siguiente bloque de código en R, se utiliza el método st_read() para leer un archivo GPKG 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 de una capa vectorial (GPKG) de provincias de Costa Rica
provincias <-
st_read(
"https://github.com/pf0953-programacionr/2024-ii/raw/refs/heads/main/datos/ign/provincias.gpkg",
quiet = TRUE # para evitar el despliegue de mensajes
)st_read() también puede crear objetos sf a partir de archivos de texto. Esta variante se utiliza principalmente cuando el archivo contiene coordenadas correspondientes a geometrías de puntos.
En el siguiente bloque de código, se utiliza st_read() para leer un archivo CSV con registros de presencia de félidos (familia Felidae) de Costa Rica, el cual contiene dos columnas llamadas decimalLongitude y decimalLatitude correspondientes a la longitud decimal y latitud decimal en las que fue observado cada felino. Este archivo proviene de una consulta al portal de datos de la Infraestructura Mundial de Información en Biodiversidad (GBIF).
# Lectura de un archivo CSV con registros de presencia de félidos en Costa Rica
felidos <-
st_read(
"https://raw.githubusercontent.com/pf0953-programacionr/2024-ii/refs/heads/main/datos/gbif/felidos.csv",
options = c(
"X_POSSIBLE_NAMES=decimalLongitude", # columna de longitud decimal
"Y_POSSIBLE_NAMES=decimalLatitude" # columna de latitud decimal
),
quiet = TRUE
)Tanto provincias como felidos son objetos de la clase sf (y además de data.frame).
# Clase del objeto provincias
class(provincias)[1] "sf" "data.frame"# Clase del objeto felidos
class(felidos)[1] "sf" "data.frame"Al escribirse el nombre de un objeto sf en la consola de R, se despliega información general sobre este.
# Información general sobre el objeto provincias
provinciasSimple feature collection with 7 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -87.10193 ymin: 5.499137 xmax: -82.55285 ymax: 11.21964
Geodetic CRS: WGS 84
gml_id gml_id2 OBJECTID CÓDIGO CÓDIGO_PROVINCIA
1 limiteprovincial_5k.1 LIMITE_PROVINCIAL.2715 2715 160103 2
2 limiteprovincial_5k.2 LIMITE_PROVINCIAL.2716 2716 160103 3
3 limiteprovincial_5k.3 LIMITE_PROVINCIAL.2717 2717 160103 5
4 limiteprovincial_5k.4 LIMITE_PROVINCIAL.2718 2718 160103 4
5 limiteprovincial_5k.5 LIMITE_PROVINCIAL.2719 2719 160103 7
6 limiteprovincial_5k.6 LIMITE_PROVINCIAL.2720 2720 160103 6
7 limiteprovincial_5k.7 LIMITE_PROVINCIAL.2721 2721 160103 1
PROVINCIA
1 Alajuela
2 Cartago
3 Guanacaste
4 Heredia
5 Limón
6 Puntarenas
7 San José
ORIGEN_DEL_TOPÓNIMO
1 Se remonta al paraje llamado La Lajuela que por primera vez se cita en los Protocolos de Cartago de 1657
2 Don Juan Vázques de Coronado escogió el sitio en el valle del Guarco para trasladar a la ciudad de Garcimuños, en 1563, bautizando al nuevo asentamiento con el nombre de Cartago
3 En alegoria a un frondoso árbol de Guanacaste ubicado en la intersección de los caminos que se dirigían a Nicoya, Bagaces y Rivas, en lo que hoy día es el parque de Liberia. Esta referencia data del siglo XVIII
4 En correspondiencia al Presidente de la Real Audiencia de Guatemala, Capitán General don Alonso Fernández de Heredia, de la Inmaculada Concepción de Cubujuquí a Villa de Heredia
5 El origen del nombre de la provincia se remonta a 1852, cuando por primera vez se cita en un documento oficial el puerto de Limón
6 En documento de 1720, se menciona la llegada del pirata Chipperton a la zona, en el cual aparece la descripcíon referente a una embarcación pequeña en la Punta de Arena, adoptando con el tiempo el nombre de Puntarenas
7 Se remonta a la creación de la ermita dedicada al Patriarca San José en 1737
VERSIÓN GLOBALID
1 20240703001 {322A624B-14A4-44DB-BA11-F37A656BF296}
2 20240703001 {46E33550-F4C9-436F-ADD7-CE70A6C46EB1}
3 20240703001 {49CD624C-A818-4AF3-9370-DC219973EC03}
4 20240703001 {1CA4FA45-F32C-4088-9031-C48CFC460B34}
5 20240703001 {A062C16C-B115-439D-BB17-B261F64AF776}
6 20240703001 {8D914665-65A7-44B6-9D63-95865A946757}
7 20240703001 {B4447D8B-B85A-4D5E-A1D6-0EA14F13E084}
SHAPE
1 MULTIPOLYGON (((-84.66639 1...
2 MULTIPOLYGON (((-84.03403 9...
3 MULTIPOLYGON (((-85.3575 11...
4 MULTIPOLYGON (((-84.02428 1...
5 MULTIPOLYGON (((-83.7025 10...
6 MULTIPOLYGON (((-84.82622 1...
7 MULTIPOLYGON (((-83.9281 10...10.3.3.2.2 st_crs() y st_transform() - manejo de sistemas de coordenadas
El método st_crs() retorna el CRS de un objeto sf.
# Despliegue del CRS del objeto provincias
st_crs(provincias)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]# Despliegue CRS del objeto felidos
st_crs(felidos)Coordinate Reference System: NAst_crs() también puede asignar un CRS a un objeto sf que no lo tiene.
# Asignación de un CRS al objeto felidos
st_crs(felidos) <- 4326El método st_transform() transforma un objeto sf a un nuevo CRS.
# Transformación del CRS del objeto provincias a WGS84 (EPSG = 4326)
provincias <-
provincias |>
st_transform(4326)10.3.3.2.3 plot() - mapeo
El método plot() grafica objetos sf en un mapa.
# Mapeo de las geometrías del objeto provincias
plot(provincias$SHAPE)
# Mapeo con argumentos adicionales de plot()
plot(
provincias$SHAPE,
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
)
Los argumentos reset y add de plot() permiten generar un mapa con varias capas.
# Primera capa del mapa
plot(
provincias$SHAPE,
extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
main = "Registros de presencia de félidos en Costa Rica",
axes = TRUE,
graticule = TRUE,
reset = FALSE
)
# Segunda capa
plot(felidos$geometry,
add = TRUE,
pch = 16,
col = "orange")
Para conocer los valores del argumento pch, puede consultar R plot pch symbols.
10.3.3.2.4 st_write() - escritura de datos
El método st_write() guarda en el disco un objeto sf en los diferentes formatos vectoriales de GDAL.
# Especificación del directorio de trabajo (debe utilizarse una ruta existente)
setwd("/home/mfvargas")
# Escritura del objeto provincias en formato GeoJSON
provincias |>
st_write("provincias.geojson")
# Escritura del objeto felidos en formato KML
felidos |>
st_write("felidos.kml")10.3.3.2.5 Otros
# Cantidad de filas de un objeto sf
nrow(provincias)[1] 7# Cantidad de columnas de un objeto sf
ncol(provincias)[1] 10# Resumen de la columna de geometría
summary(provincias$SHAPE) MULTIPOLYGON epsg:4326 +proj=long...
7 0 0 10.3.4 Mapeo de objetos sf::sf con otros paquetes
10.3.4.1 tmap
El paquete tmap genera mapas estáticos e interactivos con una sintaxis similar a la que utiliza el paquete ggplot2 para generar gráficos.
10.3.4.1.1 Instalación y carga
De acuerdo con la recomendación que se brinda en el sitio web de tmap, se instala la versión disponible ahí (versión 4) y no la disponible en CRAN (versión 3).
# Instalación de tmap (desde su repositorio en GitHub)
install.packages("remotes")
install_github("r-tmap/tmap")# Carga de tmap
library(tmap)10.3.4.1.2 Ejemplos de uso
El paquete tmap incluye varios conjuntos de datos de ejemplo, como el conjunto de datos World, que contiene información geoespacial y socioeconómica de países del mundo. El siguiente bloque de código muestra un mapa básico con los contornos de los países.
# Carga de datos de ejemplo
data("World")
# Modo estático
tmap_mode("plot")
# Mapa del mundo
tm_shape(World) +
tm_polygons()
La función tm_shape() especifica el objeto espacial que se va a visualizar, mientras que tm_polygons() dibuja los polígonos.
El siguiente bloque de código crea un mapa de coropletas que colorea el mapa de acuerdo con el valor de una variable.
# Modo interactivo
tmap_mode("view")
# Mapa de coropletas por densidad de población
tm_shape(World, name = "Países") +
tm_polygons(
col = "pop_est_dens", # columna para las coropletas
style = "quantile", # método de clasificación
palette = "YlOrRd", # paleta de colores
id = "name", # campo que se muestra al pasar el ratón
popup.vars = c("Área" = "area", # campos de la ventana de pop-up
"Población" = "pop_est",
"Densidad de población" = "pop_est_dens"),
title = "Densidad de población")A continuación, se muestran las capas de provincias de Costa Rica (polígonos) y de registros de presencia de félidos de Costa Rica (puntos).
# Crear mapa de provincias y félidos
mapa_provincias_felidos <-
tm_view(set.view = c(lon = -84.2, lat = 9.6, zoom = 7)) + # centro y zoom inicial
tm_shape(provincias, name = "Provincias") + # capa de provincias
tm_borders() +
tm_shape(felidos, name = "Félidos") + # capa de félidos
tm_dots(
col = "species", # color de los puntos
palette = "Set1", # paleta de colores
title = "Especie", # título de la leyenda
size = 0.05, # tamaño de los puntos
id = "species",
popup.vars = c("Localidad" = "locality",
"Fecha" = "eventDate",
"Fuente" = "institutionCode")
) +
tm_scale_bar(position = c("left", "bottom")) # escala
# Cambiar a modo interactivo
tmap_mode("view")
# Desplegar el mapa
mapa_provincias_felidos