2.4 Operaciones con datos espaciales

A continuación se describe una selección de las herramientas disponibles para datos espaciales. Para un listado completo de las funciones implementadas en el paquete sf se puede consultar la referencia (o la “chuleta”, aunque puede contener algunos errores).

Puede que algunas herramientas (o recursos) admitan únicamente objetos Spatial* del paquete sp, aunque siempre se pueden emplear las funciones para convertir tipos de objetos:

  • st_as_sf(x, ...): convierte x a un objeto sf (por ejemplo objetos Spatial*).
  • as(x, "Spatial"): convierte x a un objeto Spatial*.

2.4.1 Importación y exportación de datos espaciales

El paquete sf permite importar y exportar una gran cantidad de formatos de datos espaciales, almacenados en ficheros o en bases de datos, mediante las funciones st_read() y st_write(). Como se mostró al principio de la Sección 2.2, estas funciones deducen el formato automáticamente a partir de la extensión del archivo (por ejemplo .shp para ESRI Shapefile) o a partir del prefijo (por ejemplo PG: para PostGIS/PostgreSQL):

dir <- system.file("shape", package="sf")
list.files(dir, pattern="^[nc]")
## [1] "nc.dbf" "nc.prj" "nc.shp" "nc.shx"
# ESRI Shapefile, consta de por lo menos de 3 ficheros, el principal .shp
file <- paste0(dir, "/nc.shp")
file
## [1] "C:/Program Files/R/R-4.1.1/library/sf/shape/nc.shp"
nc_sf <- st_read(file)
## Reading layer `nc' from data source 
##   `C:\Program Files\R\R-4.1.1\library\sf\shape\nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27

Se admiten los formatos de datos vectoriales soportados por GDAL (que emplea internamente), se puede obtener un listado con la función st_drivers():

drivers <- st_drivers()
str(drivers)
## 'data.frame':    89 obs. of  7 variables:
##  $ name     : chr  "ESRIC" "FITS" "PCIDSK" "netCDF" ...
##  $ long_name: chr  "Esri Compact Cache" "Flexible Image Transport System" "PCIDSK Database File" "Network Common Data Format" ...
##  $ write    : logi  FALSE TRUE TRUE TRUE TRUE TRUE ...
##  $ copy     : logi  FALSE FALSE FALSE TRUE TRUE TRUE ...
##  $ is_raster: logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ is_vector: logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ vsi      : logi  TRUE FALSE TRUE FALSE TRUE TRUE ...

Figura 2.6: Listado de drivers en la instalación (local) de GDAL`.

Además, se han desarrollado una gran cantidad de paquetes de R que permiten acceder directamente desde R a datos espaciales. Muchos incluyen conjuntos de datos espaciales y otros implementan interfaces a bases de datos espaciales o geoportales disponibles en Internet. Algunos de ellos son los siguientes:

library(osmdata) 
# Cuidado: descarga mucha información
# https://nominatim.openstreetmap.org/ui/search.html
# https://wiki.openstreetmap.org/wiki/Map_features
osm_coru <- opq('A Coruña') %>%
    add_osm_feature(key = 'highway') %>%
    osmdata_sf() 
plot(st_geometry(osm_coru$osm_lines), main = "", 
     xlim = c(-8.45, -8.38), ylim = c(43.32, 43.39))
Representación de las carreteras, calles y caminos en A Coruña (generado con el paquete `osmdata`).

Figura 2.7: Representación de las carreteras, calles y caminos en A Coruña (generado con el paquete osmdata).

También están disponibles una gran cantidad de páginas web y geoportales desde donde es posible descargar datos espaciales (algo que se puede hacer directamente desde R). Algunas de ellas son:

Muchos de los archivos de datos están en formato NetCDF (Network Common Data Form) y se pueden importar a R con el paquete ncdf4.

2.4.2 Operaciones con geometrías

Operaciones unarias (operan sobre un único conjunto de geometrías simples, el primer argumento) con resultado geométrico:

  • st_geometry(): devuelve (o establece) la columna sfc de un objeto sf.
  • st_transform(x, crs, ...): transforma o convierte las coordenadas de x a un nuevo sistema de referencia.
  • st_cast(x, to, ...): cambia la geometría x a otro tipo de geometría.
  • st_centroid(): devuelve los centroides de las geometrías.
  • st_buffer(): crea un buffer en torno a la geometría o a cada geometría.
  • st_boundary(): devuelve la frontera de la geometría.
  • st_convex_hull(): crea el envoltorio convexo de un conjunto de puntos.
  • st_voronoi(): crea una teselación de Voronoi.
  • st_make_grid(x, cellsize, offset, n, what = c("polygons", "corners", "centers")): genera una rejilla rectangular (o exagonal) de geometrías (what) que cubre los límites de x.

Como ya se comentó en la Sección 2.2.1, nos puede interesar transformar las coordenadas a un nuevo sistema de referencia (algo necesario para poder combinar conjuntos de datos espaciales con distintos CRS). Por ejemplo podemos utilizar la proyección de Mollweide para representar datos globales (en este caso estimaciones de la población de países; Figura 2.8 derecha).

library(rnaturalearth) 
par_old <- par(mfrow = c(1, 2), mar = c(bottom = 0, left = 0, top = 0, right = 0))
# NOTA: plot.sf() con escala no es compatible con mfrow 
world_pop <- ne_countries(returnclass = "sf")["pop_est"]
plot(world_pop, logz = TRUE, main = "", key.pos = NULL, reset = FALSE)
grat <- st_graticule(crs=st_crs("WGS84"), lon = seq(-180, 180, by = 20), lat = seq(-90, 90, by = 10))
plot(grat[1], col = 'darkgray', add = TRUE)
# https://spatialreference.org/ref/esri/54009/
world_pop2 <- st_transform(world_pop, "ESRI:54009") 
plot(world_pop2, logz = TRUE, main = "", key.pos = NULL, reset = FALSE)
grat <- st_graticule(world_pop2, lon = seq(-180, 180, by = 20), lat = seq(-90, 90, by = 10))
plot(grat[1], col = 'darkgray', add = TRUE)
Mapa de la población estimada por paises (en escala logarítmica), datos sin proyectar (izquierda) y con proyección de Mollweide (derecha).

Figura 2.8: Mapa de la población estimada por paises (en escala logarítmica), datos sin proyectar (izquierda) y con proyección de Mollweide (derecha). 

par(par_old)

Operaciones binarias (operan sobre dos conjuntos de geometrías simples) con resultado geométrico:

  • st_union(x, y, ..., by_feature): une varias geometrías.
  • st_intersection(x, y, ...): intersección de pares de geometrías.
  • st_crop(x, y, ..., xmin, ymin, xmax, ymax): intersección con rectángulo delimitador o especificado.
  • st_difference(x, y, ...): diferencia de pares de geometrías.
  • st_sym_difference(x, y, ...): diferencia simétrica (xor) de pares de geometrías.
  • st_nearest_points(x, y, ...): obtiene los puntos más cercanos entre pares de geometrías.

Operaciones unarias con resultado numérico o lógico:

  • st_coordinates(x): devuelve una matriz con las coordenadas.
  • st_bbox(obj, ...): devuelve los límites del conjunto de geometrías.
  • st_area(x, ...): devuelve el área de polígonos.
  • st_length(x, ...): devuelve la longitud de líneas.
  • st_is(x, type): verifica si la geometría es de un determinado tipo o conjunto de clases.

Operaciones binarias con resultado numérico o lógico:

  • st_distance(x, y, ..., by_element, which): devuelve la matriz de distancias mínimas entre geometrías.
  • st_nearest_feature(x, y): devuelve el índice de la geometría de y más cercana a cada geometría de x.
  • st_intersects(x, y, ...): determina si las geometrías se solapan o tocan.
  • st_disjoint(x, y, ...): determina si las geometrías no se solapan o tocan.
  • st_touches(x, y, ...): determina si las geometrías se tocan.
  • st_overlaps(x, y, ...): determina si las geometrías se solapan, pero no están completamente contenidas la una en la otra.
  • st_crosses(x, y, ...): determina si las geometrías se cruzan, pero no se tocan.
  • st_within(x, y, ...): determina si x está en y.
  • st_contains(x, y, ...): determina si y está en x.
  • st_covers(x, y, ...): determina si todos los puntos de y están dentro de x.
  • st_covered_by(x, y, ...): determina si todos los puntos de x están dentro de y.
  • st_equals(x, y, ...): determina si x es geométricamente igual a y.
  • st_equals_exact(x, y, par, ...): determina si x es igual a y con cierta tolerancia.

El resultado de las operaciones lógicas es una matriz dispersa (de clase sgbp, sparse geometry binary predicate), que se puede convertir a una matriz densa con as.matrix().

Ejemplo 2.1 (Creación de una rejilla de predicción) Continuando con los datos del Ejercicio 2.2, para crear un objeto con las posiciones de predicción, podríamos generar un buffer (st_buffer()) de radio 40 en torno a las posiciones de observación y a partir de él crear una rejilla vectorial (st_make_grid(..., what = "centers")) de dimensiones 50 por 50 e intersecarla con el buffer. 
load("datos/aquifer.RData")
aquifer$head <- aquifer$head/100 # en cientos de pies
aquifer_sf <- st_as_sf(aquifer, coords = c("lon", "lat"), agr = "constant")
buffer <- aquifer_sf %>% st_geometry() %>%  st_buffer(40)
grid <- buffer %>% st_make_grid(n = c(50, 50), what = "centers") %>% 
  st_intersection(buffer)
plot(buffer)
plot(grid, pch = 3, cex = 0.5, add = TRUE)
Rejilla en torno a las posiciones de los datos de `aquifer`.

Figura 2.9: Rejilla en torno a las posiciones de los datos de aquifer.

Sin embargo, en lugar de emplear una rejilla sf, puede resultar preferible (por ejemplo para la representación gráfica) emplear una rejilla stars

library(stars)
grid <- buffer %>%  st_as_stars(nx = 50, ny = 50) %>% st_crop(buffer)
idw <- gstat::idw(formula = head ~ 1, locations = aquifer_sf, newdata = grid)
## [inverse distance weighted interpolation]
plot(idw["var1.pred"], col = sf.colors(64), main = "")
Interpolación por IDW (Inverse Distance Weighting) de los datos del acuífero Wolfcamp.

Figura 2.10: Interpolación por IDW (Inverse Distance Weighting) de los datos del acuífero Wolfcamp. 

# Error gstat::idw, cambia las coordenadas del objeto stars
# summary(st_coordinates(grid))
# summary(st_coordinates(idw))
# Posible solución: añadir el resultado a `grid` y emplearlo en lugar de `idw`
# grid$var1.pred <- idw$var1.pred
# plot(grid["var1.pred"], col = sf.colors(64), axes = TRUE, main = "")