Descripción general

En esta lección se cubren algunas de las capacidades del lenguaje de programación R para manejo de datos vectoriales. Se presentan los principales estándares, paquetes y algunos ejemplos de mapas elaborados con este tipo de datos.

Objetivos

  1. Conocer el estándar Simple Features de datos vectoriales.
  2. Aprender a instalar y a utilizar el paquete sf.

Trabajo previo

Se recomienda leer el capítulo 2 del libro Geocomputation with R, de Robin Lovelace et al., en el cual está basado este documento.

El modelo de datos vectorial

El modelo de datos vectorial 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)

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.

sf cuenta con varias ventajas con respecto a otros paquetes geoespaciales de R (ej. sp):

  • Rápida lectura y escritura de datos.
  • Rendimiento de graficación mejorado.
  • Manejo de objetos como data frames.
  • Las funciones pueden combinarse con el operador %>%, lo que las combina bien con los paquetes de tidyverse.
  • Los nombres de las funciones son consistentes e intuitivos (comienzan con st_).

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 (ej. summary(), View(), str()) y con la misma forma de referencias las filas (observaciones) y las columnas (variables).

Instalación y carga de sf

Instalación:

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

Carga:

# Carga de sf
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

Documentación:

 # Páginas disponibles
vignette(package = "sf")

# Página introductoria
vignette("sf1")

Instalación de conjuntos de datos utilizados en los ejemplos

Para ejemplos, se utilizarán los conjuntos de datos spData y spDataLarge.

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

# Instalación de spDataLarge
# Si es necesario, instale antes el paquete remotes con install.packages("remotes")
remotes::install_github("Nowosad/spDataLarge")

Carga:

# Carga de spData y spDataLarge
library(spData)
library(spDataLarge)

Ejemplos de uso de sf

En los siguientes ejemplos, se utilizará el conjunto de datos world de spData.

Exploración inicial de datos

Con la función names(), pueden desplegarse los nombres de las columnas de un data frame:

# Despliegue de los nombres de columnas
names(world)
##  [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"     
##  [7] "area_km2"  "pop"       "lifeExp"   "gdpPercap" "geom"
# Estructura del data frame
str(world)
## Classes 'sf', 'tbl_df', 'tbl' and 'data.frame':  177 obs. of  11 variables:
##  $ iso_a2   : chr  "FJ" "TZ" "EH" "CA" ...
##  $ name_long: chr  "Fiji" "Tanzania" "Western Sahara" "Canada" ...
##  $ continent: chr  "Oceania" "Africa" "Africa" "North America" ...
##  $ region_un: chr  "Oceania" "Africa" "Africa" "Americas" ...
##  $ subregion: chr  "Melanesia" "Eastern Africa" "Northern Africa" "Northern America" ...
##  $ type     : chr  "Sovereign country" "Sovereign country" "Indeterminate" "Sovereign country" ...
##  $ area_km2 : num  19290 932746 96271 10036043 9510744 ...
##  $ pop      : num  8.86e+05 5.22e+07 NA 3.55e+07 3.19e+08 ...
##  $ lifeExp  : num  70 64.2 NA 82 78.8 ...
##  $ gdpPercap: num  8222 2402 NA 43079 51922 ...
##  $ geom     :sfc_MULTIPOLYGON of length 177; first list element: List of 3
##   ..$ :List of 1
##   .. ..$ : num [1:8, 1:2] 180 180 179 179 179 ...
##   ..$ :List of 1
##   .. ..$ : num [1:9, 1:2] 178 178 179 179 178 ...
##   ..$ :List of 1
##   .. ..$ : num [1:5, 1:2] -180 -180 -180 -180 -180 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geom"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA
##   ..- attr(*, "names")= chr  "iso_a2" "name_long" "continent" "region_un" ...

Nótese que la columna geom contiene las geometrías correspondientes a los multiplígonos de los países.

sf implementa la función llamada plot() para visualización de datos geográficos:

# Despliegue de datos geoespaciales
plot(world)
## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all

En el ejemplo anterior, se desplegó un mapa por cada variable (i.e. cada columna). plot() permite también desplegar el mapa correspondiente a una variable particular:

# Despliegue de datos geoespaciales de una columna
plot(world["pop"])

Filtrado de datos

Puede separarse una “porción” de un data frame geoespacial de la misma forma en la que se hace con cualquier otro data frame:

# Subconjunto de un data frame
world_muestra <- world[1:2, 1:8]

# Despliegue tabular
world_muestra
## Simple feature collection with 2 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -180 ymin: -18.28799 xmax: 180 ymax: -0.95
## geographic CRS: WGS 84
##   iso_a2 name_long continent region_un      subregion              type
## 1     FJ      Fiji   Oceania   Oceania      Melanesia Sovereign country
## 2     TZ  Tanzania    Africa    Africa Eastern Africa Sovereign country
##    area_km2      pop                           geom
## 1  19289.97   885806 MULTIPOLYGON (((180 -16.067...
## 2 932745.79 52234869 MULTIPOLYGON (((33.90371 -0...
# Despliegue del mapa de población
plot(world_muestra["pop"])

También puede realizarse con filtro con base en el valor de los datos de una columna:

# Países del continente africano
africa_paises = world[world$continent == "Africa", ]

# Mapa de producto interno bruto
plot(africa_paises["gdpPercap"], main="Producto interno bruto per capita")

Confección de mapas

Es posible agregar capas a un mapa mediante el uso de los parámetros reset=FALSE y add=TRUE de la función plot().

# Filtro para crear data frame con países de África
africa_paises <- world[world$continent == "Africa", ]

# La función st_union() realiza una unión de los polígonos
africa_continente <- st_union(africa_paises)

# Primera capa del mapa
plot(world["pop"], reset=FALSE)
# Segunda capa
plot(africa_continente, add=TRUE, col = "red")

Lectura y escritura de datos

Las operaciones de lectura tranforman datos contenidos en diferentes fuentes al formato requerido por el paquete sf. La lectura se lleva a cabo con la función st_read(). A continuación, se presentan ejemplos para tres tipos diferentes de fuentes de datos:

  • Archivos geoespaciales (ej. Shapefile, Geopackage, GeoJSON).
  • Servicios WFS.
  • Archivos CSV.

Lectura de un archivo geoespacial:

# COVID-19 en cantones de Costa Rica
cr_cantones_covid19 <- st_read("https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/covid19/casos/cr/cr-covid19-cantones.geojson")
## Reading layer `cr-covid19-cantones' from data source `https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/covid19/casos/cr/cr-covid19-cantones.geojson' using driver `GeoJSON'
## Simple feature collection with 82 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -87.10189 ymin: 5.498564 xmax: -82.55285 ymax: 11.21964
## geographic CRS: WGS 84
# Mapa de coropletas de casos confirmados de COVID-19
plot(cr_cantones_covid19["confirmados"], main="Casos confirmados de COVID-19 en los cantones de Costa Rica")

Lectura de un Web Feature Service (WFS). Se utiliza una capa publicada por el Instituto Geográfico Nacional (IGN) en la Infraestructura Nacional de Datos Espaciales de Costa Rica (SNIT):

# Dirección base del servicio WFS
url_base <- "http://geos.snitcr.go.cr/be/IGN_5/wfs?"

# Parámetros de la solicitud WFS de la capa:
#   typeName: nombre de la capa
#   outputFormat: formato de salida
solicitud_wfs <- "request=GetFeature&service=WFS&version=2.0.0&typeName=IGN_5:limiteprovincial_5k&outputFormat=application/json"

# URL completo del servicio WFS
cr_provincias_wfs <- paste0(url_base, solicitud_wfs)

# Recuperación de los datos en un data frame
cr_provincias <- st_read(cr_provincias_wfs)
## Reading layer `OGRGeoJSON' from data source `http://geos.snitcr.go.cr/be/IGN_5/wfs?request=GetFeature&service=WFS&version=2.0.0&typeName=IGN_5:limiteprovincial_5k&outputFormat=application/json' using driver `GeoJSON'
## Simple feature collection with 7 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 156152 ymin: 608833.8 xmax: 658879.5 ymax: 1241118
## projected CRS:  CR05 / CRTM05
# Mapeo del data frame
plot(cr_provincias$geometry)

Lectura de un archivo de valores separados por comas (CSV). Se utilizan datos de presencia de especies agrupados por la Infraestructura Mundial de Información en Biodiversidad (GBIF):

# Lectura del archivo CSV
cr_ara_ambiguus <- st_read(
  "https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/biodiversidad/registros-presencia/cr/cr-ara-ambiguus.csv", 
  options=c("X_POSSIBLE_NAMES=decimalLongitude","Y_POSSIBLE_NAMES=decimalLatitude")
)
## options:        X_POSSIBLE_NAMES=decimalLongitude Y_POSSIBLE_NAMES=decimalLatitude 
## Reading layer `cr-ara-ambiguus' from data source `https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/biodiversidad/registros-presencia/cr/cr-ara-ambiguus.csv' using driver `CSV'
## Simple feature collection with 87 features and 50 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -84.75266 ymin: 9.606838 xmax: -82.60557 ymax: 10.95063
## CRS:            NA
# Mapeo del data frame
plot(cr_ara_ambiguus$geometry, col="red")

Mapa de provincias y registros de presencia:

# Se cambia la proyección a WGS84 (EPSG 4326)
cr_provincias_wgs84 = st_transform(cr_provincias, 4326)

# Primera capa del mapa
plot(cr_provincias_wgs84$geometry, 
     main="Distribución de Ara ambiguus en Costa Rica",
     reset=FALSE, 
     axes = TRUE,
     graticule=TRUE
)

# Segunda capa
plot(cr_ara_ambiguus$geometry, add=TRUE, col = "green")