2.2 Introducción al paquete sf

El modelo de geometrías de características simples (o rasgos simples) es un estándar (ISO 19125) desarrollado por el Open Geospatial Consortium (OGC) para formas geográficas vectoriales, que ha sido adoptado por gran cantidad de software geográfico (entre otros por GeoJSON, ArcGIS, QGIS, PostGIS, MySQL Spatial Extensions, Microsoft SQL Server…). Como ya se comentó, este tipo de datos espaciales está implementado en R en el paquete sf.

Los objetos principales, del tipo sf, son extensiones de data.frame (o tibble) y como mínimo contienen una columna denominada simple feature geometry list column que contiene la geometría de cada observación (se trata de una columna tipo list). Cada fila, incluyendo la geometría y otras posibles variables (denominados atributos de la geometría), se considera una característica simple (SF).

library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
nc <- nc[c(5, 9:15)]
nc
## Simple feature collection with 100 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 10 features:
##           NAME BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
## 1         Ashe  1091     1      10  1364     0      19
## 2    Alleghany   487     0      10   542     3      12
## 3        Surry  3188     5     208  3616     6     260
## 4    Currituck   508     1     123   830     2     145
## 5  Northampton  1421     9    1066  1606     3    1197
## 6     Hertford  1452     7     954  1838     5    1237
## 7       Camden   286     0     115   350     2     139
## 8        Gates   420     0     254   594     2     371
## 9       Warren   968     4     748  1190     2     844
## 10      Stokes  1612     1     160  2038     5     176
##                          geometry
## 1  MULTIPOLYGON (((-81.47276 3...
## 2  MULTIPOLYGON (((-81.23989 3...
## 3  MULTIPOLYGON (((-80.45634 3...
## 4  MULTIPOLYGON (((-76.00897 3...
## 5  MULTIPOLYGON (((-77.21767 3...
## 6  MULTIPOLYGON (((-76.74506 3...
## 7  MULTIPOLYGON (((-76.00897 3...
## 8  MULTIPOLYGON (((-76.56251 3...
## 9  MULTIPOLYGON (((-78.30876 3...
## 10 MULTIPOLYGON (((-80.02567 3...
str(nc)
## Classes 'sf' and 'data.frame':   100 obs. of  8 variables:
##  $ NAME    : chr  "Ashe" "Alleghany" "Surry" "Currituck" ...
##  $ BIR74   : num  1091 487 3188 508 1421 ...
##  $ SID74   : num  1 0 5 1 9 7 0 0 4 1 ...
##  $ NWBIR74 : num  10 10 208 123 1066 ...
##  $ BIR79   : num  1364 542 3616 830 1606 ...
##  $ SID79   : num  0 3 6 2 3 5 2 2 2 5 ...
##  $ NWBIR79 : num  19 12 260 145 1197 ...
##  $ geometry:sfc_MULTIPOLYGON of length 100; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:7] "NAME" "BIR74" "SID74" "NWBIR74" ...

El nombre de la columna de geometrías está almacenado en el atributo "sf_column" del objeto y se puede acceder a ella mediante la función st_geometry() (además de poder emplear los procedimientos habituales para acceder a los componentes de un data.frame). Esta columna es un objeto de tipo sfc (simple feature geometry list column), descritos más adelante.

# geom_name <- attr(nc, "sf_column")
# nc[, geom_name]; nc[[geom_name]]
# nc$geometry
st_geometry(nc)
## Geometry set for 100 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 5 geometries:
## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...
## MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...
## MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...
## MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...
## MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...

En este paquete, todos los métodos y funciones que operan sobre datos espaciales comienzan por st_ (spatial type; siguiendo la implementación de PostGIS):

methods(class="sf")
##  [1] $<-                   [                     [[<-                 
##  [4] aggregate             as.data.frame         cbind                
##  [7] coerce                dbDataType            dbWriteTable         
## [10] filter                identify              initialize           
## [13] merge                 plot                  print                
## [16] rbind                 show                  slotsFromS3          
## [19] st_agr                st_agr<-              st_area              
## [22] st_as_s2              st_as_sf              st_bbox              
## [25] st_boundary           st_buffer             st_cast              
## [28] st_centroid           st_collection_extract st_convex_hull       
## [31] st_coordinates        st_crop               st_crs               
## [34] st_crs<-              st_difference         st_filter            
## [37] st_geometry           st_geometry<-         st_inscribed_circle  
## [40] st_interpolate_aw     st_intersection       st_intersects        
## [43] st_is                 st_is_valid           st_join              
## [46] st_line_merge         st_m_range            st_make_valid        
## [49] st_nearest_points     st_node               st_normalize         
## [52] st_point_on_surface   st_polygonize         st_precision         
## [55] st_reverse            st_sample             st_segmentize        
## [58] st_set_precision      st_shift_longitude    st_simplify          
## [61] st_snap               st_sym_difference     st_transform         
## [64] st_triangulate        st_union              st_voronoi           
## [67] st_wrap_dateline      st_write              st_z_range           
## [70] st_zm                 transform            
## see '?methods' for accessing help and source code

Los objetos geométricos básicos son del tipo sfg (simple feature geometry) que contienen la geometría de una única característica definida a partir de puntos en dos (XY), tres (XYZ, XYM) o cuatro dimensiones (XYZM). Admite los 17 tipos de geometrías simples del estándar, pero de forma completa los 7 tipos básicos:

Tipo Description Creación
POINT, MULTIPOINT Punto o conjunto de puntos st_point(), st_multipoint()
LINESTRING, MULTILINESTRING Línea o conjunto de líneas st_linestring(), st_multilinestring()
POLYGON, MULTIPOLYGON Polígono6 o conjunto de polígonos st_polygon(), st_multipolygon()
GEOMETRYCOLLECTION Conjunto de geometrías de los tipos anteriores st_geometrycollection()

Las geometrías se imprimen empleando la representación well-known text (WKT) del estándar (se exportan empleando la representación well-known binary, WKB).

nc$geometry[[1]]
## MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))

Los objetos básicos sfg (normalmente del mismo tipo) se pueden combinar en un objeto sfc (simple feature geometry list column) mediante la función st_sfg(). Estos objetos pueden incorporar un sistema de referencia de coordenadas (por defecto NA_crs_), descritos en la Sección 2.2.1. Posteriormente se puede crear un objeto sf mediante la función st_sf().

p1 <- st_point(c(-8.395835, 43.37087))
p2 <- st_point(c(-7.555851, 43.01208))
p3 <- st_point(c(-7.864641, 42.34001))
p4 <- st_point(c(-8.648053, 42.43362))
sfc <- st_sfc(list(p1, p2, p3, p4))
cprov <- st_sf(names = c('Coruña (A)', 'Lugo', 'Ourense', 'Pontevedra'),
    geom = sfc)
cprov
## Simple feature collection with 4 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -8.648053 ymin: 42.34001 xmax: -7.555851 ymax: 43.37087
## CRS:           NA
##        names                       geom
## 1 Coruña (A) POINT (-8.395835 43.37087)
## 2       Lugo POINT (-7.555851 43.01208)
## 3    Ourense POINT (-7.864641 42.34001)
## 4 Pontevedra POINT (-8.648053 42.43362)

Esta forma de proceder puede resultar de interés cuando se construyen geometrías tipo líneas o polígonos, pero en el caso de datos puntuales (las observaciones habituales en geoestadística), resulta mucho más cómodo emplear un data.frame que incluya las coordenadas en columnas y convertirlo a un objeto sf mediante la función st_as_sf().


Ejercicio 2.1 (Creación de una columna de geometrías) Crear una geometría (un objeto sfc) formada por: dos puntos en las posiciones (1,5) y (5,5), una línea entre los puntos (1,1) y (5,1), y un polígono, con vértices {(0,0), (6,0), (6,6), (0,6), (0,0)} y con un agujero con vértices {(2,2), (2,4), (4,4), (4,2), (2,2)} (NOTA: consultar la ayuda ?st, puede resultar cómodo emplear matrix(... , ncol = 2, byrow = TRUE)).


Como ejemplo consideraremos el conjunto de datos meuse del paquete sp que contiene concentraciones de metales pesados, junto con otras variables del terreno, en una zona de inundación del río Meuse (cerca de Stein, Holanda)7 (ver Figura 2.1).

data(meuse, package="sp")
str(meuse)
## 'data.frame':    155 obs. of  14 variables:
##  $ x      : num  181072 181025 181165 181298 181307 ...
##  $ y      : num  333611 333558 333537 333484 333330 ...
##  $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##  $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
##  $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
##  $ zinc   : num  1022 1141 640 257 269 ...
##  $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
##  $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
##  $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##  $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##  $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##  $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##  $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...
# ?meuse 
# Sistema de coordenadas Rijksdriehoek (RDH) (Netherlands topographical) 
# https://epsg.io/28992 # EPSG:28992
meuse_sf <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992, agr = "constant")

# Rio Meuse 
data(meuse.riv, package="sp")
str(meuse.riv)
##  num [1:176, 1:2] 182004 182137 182252 182315 182332 ...
meuse_riv <- st_sfc(st_polygon(list(meuse.riv)), crs = 28992)

# Rejilla
data(meuse.grid, package="sp")
str(meuse.grid)
## 'data.frame':    3103 obs. of  7 variables:
##  $ x     : num  181180 181140 181180 181220 181100 ...
##  $ y     : num  333740 333700 333700 333700 333660 ...
##  $ part.a: num  1 1 1 1 1 1 1 1 1 1 ...
##  $ part.b: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ dist  : num  0 0 0.0122 0.0435 0 ...
##  $ soil  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
meuse_grid <- st_as_sf(meuse.grid, coords = c("x", "y"), 
                       crs = 28992, agr = "constant")
# Almacenar
# save(meuse_sf, meuse_riv, meuse_grid, file = "datos/st_meuse.RData")

# Representar
plot(meuse_sf["zinc"], pch = 16, cex = 1.5, main = "",
     breaks = "quantile", key.pos = 4, reset = FALSE)
plot(meuse_riv, col = "lightblue", add = TRUE)
plot(st_geometry(meuse_grid), pch = 3, cex = 0.2, col = "lightgray", add = TRUE)
Concentración de zinc (ppm) en el entorno del río Meuse (datos `sp::meuse`).

Figura 2.1: Concentración de zinc (ppm) en el entorno del río Meuse (datos sp::meuse).


Ejercicio 2.2 (Creación y representación de datos espaciales) Cargar los datos del acuífero Wolfcamp (aquifer.RData), generar el correspondiente objeto sf y representarlo mostrando los ejes.


2.2.1 Sistemas de referencia de coordenadas

El sistema de referencia de coordenadas (CRS) especifica la correspondencia entre valores de las coordenadas y puntos concretos en la superficie de la Tierra (o del espacio), y resulta fundamental cuando se combinan datos espaciales. En general se consideran dos tipos de CRS:

  • Geodésico: las coordenadas en tres dimensiones (latitud, longitud y altura) se basan en un elipsoide de referencia (global o local) que sirve como aproximación del globo terrestre (se tiene en cuenta que no es una esfera perfecta e incluso que puede haber variaciones locales). Este elipsoide, junto con información adicional sobre como interpretar las coordenadas (incluyendo el orden y el origen), define el denominado datum. Normalmente se asume que las coordenadas son en la superficie terrestre y solo se consideran:

    • latitud: ángulo entre el plano ecuatorial y la línea que une la posición con el centro de la Tierra. Varía desde -90 (polo sur) hasta 90 (polo norte). Un grado equivale aproximadamente a 110.6 km. Los paralelos son las líneas en la superficie terrestre correspondientes a la misma latitud (siendo el 0 el ecuador).
    • longitud: ángulo (paralelo al plano ecuatorial) entre un meridiano de referencia (arco máximo que une los polos pasando por una determinado punto, normalmente el observatorio de Greenwich) y la línea que une la posición con el centro de la Tierra. Varía desde -180 (oeste) hasta 180 (este). Un grado en el ecuador equivale a aproximadamente a 111.3 km. Los meridianos son las líneas en la superficie terrestre correspondientes a la misma longitud (siendo el 0 el meridiano de Greenwich y -180 o 180 el correspondiente antimeridiano).

    La rejilla correspondiente a un conjunto de paralelos y meridianos se denomina gratícula (ver st_graticule()).

    Coordenadas geográficas en la superficie terrestre (Fuente Wikimedia Commons).

    Figura 2.2: Coordenadas geográficas en la superficie terrestre (Fuente Wikimedia Commons).

    Uno de los CRS más empleados es el WGS84 (World Geodetic System 1984) en el que se basa el Sistema de Posicionamiento Global (GPS).

  • Proyectado (cartesiano): sistema (local) en dos dimensiones que facilita algún tipo de cálculo (normalmente distancias o áreas). Por ejemplo el UTM (Universal Transverse Mercator), que emplea coordenadas en metros respecto a una cuadrícula de referencia (se divide la tierra en 60 husos de longitud, numerados, y 20 bandas de latitud, etiquetadas con letras; por ejemplo Galicia se encuentra en la cuadricula 29T). Se define relacionando estas coordenadas cartesianas con coordenadas geodésicas con un determinado datum.

En sf se emplea la librería PRØJ para definir el CRS y convertir coordenadas en distintos sistemas8. Para obtener o establecer el CRS se puede emplear la función st_crs(). Se puede especificar mediante una cadena de texto que admita GDAL (por ejemplo "WGS84", que se corresponde con el World Geodetic System 1984), que típicamente es de la forma ESTÁNDAR:CÓDIGO (también puede ser una cadena de texto PROJ.4). El estándar más empleado es el EPSG (European Petroleum Survey Group), y es que que da por hecho el paquete sf cuando se especifica el CRS mediante un número. También admite el estándar OGC WKT (Open Geospatial Consortium well-known text) que es el que emplea internamente, pero resulta complicado manejar en la práctica.

st_crs("WGS84")
## Coordinate Reference System:
##   User input: WGS84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     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]],
##     ID["EPSG",4326]]
all.equal(st_crs(4326), st_crs("EPSG:4326"), st_crs("WGS84"))
## [1] TRUE
st_crs(nc)
## Coordinate Reference System:
##   User input: NAD27 
##   wkt:
## GEOGCRS["NAD27",
##     DATUM["North American Datum 1927",
##         ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4267]]

En spatialreference.org se puede obtener información detallada sobre una gran cantidad de proyecciones (y permite realizar búsquedas). También puede ser de utilidad epsg.io o este listado con detalles de los parámetros.

El CRS ideal dependerá del tipo de problema y de la zona cubierta por los datos (ver e.g Lovelace et al, 2021, Sección 6.3, para más información). En general en estadística espacial nos interesará trabajar con coordenadas proyectadas, de forma que tenga sentido emplear la distancia euclídea (algo que puede ser poco o nada razonable si se trabaja con coordenadas geodésicas en una zona muy amplia del globo o cerca de los polos). En el caso de coordenadas sin proyectar (latitud/longitud) puede ser preferible trabajar con distancias ortodrómicas (longitud del arco del círculo máximo que une los puntos, great circle distances)9. Es importante destacar que cambiar el CRS no reproyecta los datos, hay que emplear st_transform() para hacerlo, como se describe en la Sección 2.4.

Finalmente hay que insistir también en que el campo de aplicación de la estadística espacial no se restringe al análisis de datos geográficos (por ejemplo nos puede interesar analizar el desgaste en la pared de un crisol empleado en fundición) y en estos casos los CRS geográficos carecen de sentido. De todos modos habrá que emplear un sistema de coordenadas que permita calcular algún tipo de salto o distancia entre puntos (aunque siempre se pueden considerar coordenadas espaciales tres dimensiones con la distancia euclídea).

2.2.2 Integración con el ecosistema tidyverse

El paquete sf es compatible con tidyverse y proporciona métodos para interactuar con los paquetes dplyr, tidyr y ggplot2.

Algunos de los métodos de interés para manipular datos espaciales con el paquete dplyr son:

  • filter(), select(), mutate(), summarise(..., do_union = TRUE, is_coverage = FALSE), group_by(), ungroup(), etc.
  • inner_join(), left_join(), right_join(), full_join(), semi_join(), anti_join(), st_join().
  • st_drop_geometry(), st_set_crs().

Para detalles ver la referencia.

En el caso del paquete ggplot2 se puede consultar la referencia y el tutorial Drawing beautiful maps programmatically with R, sf and ggplot2:

Por ejemplo, se puede generar un gráfico similar al de la Figura 1.2 (porcentaje de incremento de las defunciones en el año 2020 respecto al 2019 en las CCAA españolas; datos provisionales INE), con el siguiente código:

library(dplyr)
library(mapSpain) 
mortalidad <- read.csv2("datos/mortalidad.csv")
CCAA_sf <- esp_get_ccaa() %>% left_join(mortalidad) %>%
  mutate(incremento = 100*(mort.2020 - mort.2019)/mort.2019)

library(ggplot2)
ggplot(CCAA_sf) +
  geom_sf(aes(fill = incremento), color = "grey70") +
  scale_fill_gradientn(colors = hcl.colors(10, "Blues", rev = TRUE)) +
  geom_sf_label(aes(label = paste0(round(incremento, 1), "%")), alpha = 0.5) +
  geom_sf(data = esp_get_can_box(), color = "grey70") + 
  theme_void()
Ejemplo de gráfico generado empleando los paquetes `dplyr` y `ggplot2`.

Figura 2.3: Ejemplo de gráfico generado empleando los paquetes dplyr y ggplot2.

Sin embargo, en este libro se supone que no se está familiarizado con estas herramientas y se evitará su uso (aunque pueden resultar más cómodas después de su aprendizaje). Para una introducción a dplyr, ver por ejemplo la viñeta Introduction to dplyr, el Capítulo 5 del libro R for Data Science o el Capítulo 4 de los apuntes Prácticas de Tecnologías de Gestión y Manipulación de Datos.

No obstante, en ciertas ocasiones emplearemos el operador pipe %>% (tubería, redirección) por comodidad. Este operador permite canalizar la salida de una función a la entrada de otra. Por ejemplo segundo(primero(datos)) se traduce en datos %>% primero %>% segundo (facilitando la lectura de expresiones de izquierda a derecha).


  1. Secuencia de puntos que forma un anillo cerrado, que no se interseca; el primero anillo definen el anillo exterior, anillos posteriores definen agujeros. Según la norma, los puntos del anillo exterior deben especificarse en sentido contrario a las agujas del reloj y los de los agujeros en sentido de las agujas del reloj.↩︎

  2. Empleado en la viñeta del paquete gstat con el paquete sp.↩︎

  3. El paquete sf admite las últimas versiones PROJ 5 y 6, incluyendo el formato WKT-2 de 2019, mientras que el paquete sp está diseñado para cadenas de texto PROJ.4 que se recomiendan abandonar (las últimas versiones permiten añadir una cadena WKT2 como comment).↩︎

  4. Algo que ya hace de forma automática el paquete gstat.↩︎