Introducción

El paquete sf de R fue diseñado para ofrecer soporte al estándar simple features (sf) —o «objetos simples»—, utilizado para la representación de datos geográficos vectoriales dentro del ecosistema R. El paquete sf se lanzó por primera vez en 2016 y, en pocos años, superó con creces al paquete sp (lanzado en 2005) como la opción más utilizada en R para representar datos geográficos vectoriales. Además de brindar soporte al estándar simple features, sf incluye una amplia gama de funciones para trabajar con datos espaciales vectoriales en R. Si usted está acostumbrado a trabajar con datos espaciales mediante software SIG de escritorio tradicional, como QGIS o ArcGIS, encontrará en el paquete sf —potenciado por los paquetes del tidyverse— equivalentes a las herramientas SIG habituales para el análisis, la consulta, la selección, el geoprocesamiento y la visualización.

Desde la página web de Simple Features para R

«Simple features» (sf), o acceso a «simple features», hace referencia a un estándar formal (ISO 19125-1:2004) que describe cómo pueden representarse en sistemas informáticos los objetos del mundo real, poniendo énfasis en la geometría espacial de dichos objetos. Asimismo, describe cómo es posible almacenar estos objetos en bases de datos y recuperarlos de ellas, además de definir qué operaciones geométricas deben aplicarse a los mismos. Este estándar cuenta con una amplia implementación en bases de datos espaciales (como PostGIS) y sistemas de información geográfica (SIG) comerciales (por ejemplo, ESRI ArcGIS), y constituye la base de datos vectoriales para bibliotecas como GDAL.

Objetivo del cuaderno

Este cuaderno te presentará el paquete sf y el trabajo con datos vectoriales geográficos en el estándar de datos de objetos simples (sf). Este cuaderno solo pretende ser una introducción rápida al paquete sf, ya que muchas de sus funciones no se abordarán aquí.

✨ Requisitos previos ✨

Tener una copia del archivo shapefile MGN2023_MPIO_POLITICO disponible en DANE (y utilizado desde el inicio del curso de Geomática Básica).

Recursos adicionales en la web

https://www.jessesadler.com/post/simple-feature-objects/

https://learning.nceas.ucsb.edu/2023-04-coreR/session_15.html

Descripción general del cuaderno (notebook)

1. Setup

Este notebook debe ubicarse en tu directorio GB2, cerca del directorio donde guardaste tus datos.

Este notebook requiere los siguientes paquetes y funciones de R.

Este cuaderno utiliza las siguientes funciones de sf:

  • st_as_sf· convertir un objeto externo en un objeto sf

  • st_bbox · devolver el cuadro delimitador (bounding box) de una objetos simple o de un conjunto de objetos simples

  • st_geometry · obtener, establecer, reemplazar o renombrar la geometría de un objeto sf

  • st_read · leer entidades simples o capas desde un archivo o una base de datos

  • st_write · escribir un objeto de entidades simples en un archivo o una base de datos

  • https://dplyr.tidyverse.org/ - dplyr es una gramática para la manipulación de datos que ofrece un conjunto coherente de verbos, los cuales ayudan a resolver los desafíos más habituales en la manipulación de datos:

  • mutate() añade nuevas variables que son funciones de variables ya existentes.

  • select() selecciona variables basándose en sus nombres.

  • filter() selecciona casos basándose en sus valores.

  • summarise() reduce múltiples valores a un único resumen.

  • arrange() modifica el orden de las filas.

  • https://ggplot2.tidyverse.org/ ggplot2 es un sistema para crear gráficos de manera declarativa, basado en «La gramática de los gráficos» (The Grammar of Graphics). Tú proporcionas los datos, indicas a ggplot2 cómo asignar las variables a los elementos estéticos y qué primitivas gráficas utilizar; el sistema se encarga del resto de los detalles.

  • https://paleolimbot.github.io/ggspatial/ - Los datos espaciales, sumados a la potencia del marco de trabajo ggplot2, facilitan la creación de mapas.

Puedes instalar las librerías necesarias desde la consola utilizando install.packages(“nombre_de_la_librería”). Si la instalación se realiza correctamente, deberás cargar sf y dplyr para poder utilizar sus funciones:

library(sf)
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

Enlazado con GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() es TRUE

library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Cargando paquete: ‘dplyr’

Los siguientes objetos están enmascarados del paquete ‘stats’:

filter, lag

Los siguientes objetos están enmascarados desde ‘package:base’:

intersect, setdiff, setequal, union

library(ggplot2)
library(ggspatial)

2. Lectura y escritura de datos espaciales

2a. Lectura de datos espaciales con st_read

La función st_read, del paquete sf, lee datos vectoriales espaciales a partir de una variedad de formatos compatibles con GDAL (Geospatial Data Abstraction Library), incluidos ESRI Shapefile (.shp), GeoJSON (.geojson), GeoPackage (.gpkg), KML (.kml) y KMZ (.kmz), archivos CSV con WKT (well-known text) (.csv), entre otros.

st_read simplifica el proceso de lectura de datos espaciales al reconocer automáticamente el formato del archivo y devolver los datos como un objeto de tipo «simple features» (sf).

★ Consejo profesional: https://en.wikipedia.org/wiki/Spatial_reference_system (CRS) está especificado en el archivo, se leerá automáticamente.

#Get the filenames close to your notebook
list.files("D:/MI PC/GB02/Proyecto7/Mun_Cesar.gpkg") 
## character(0)

Ahora puedes especificar la ruta relativa para leer el shapefile:

# read the *Municipios* shapefile as an sf object
Cesar <- st_read("D:/MI PC/GB02/Proyecto7/Mun_Cesar.gpkg")
## Multiple layers are present in data source D:\MI PC\GB02\Proyecto7\Mun_Cesar.gpkg, reading layer `Municipios'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.
## Reading layer `Municipios' from data source 
##   `D:\MI PC\GB02\Proyecto7\Mun_Cesar.gpkg' using driver `GPKG'
## Simple feature collection with 24 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.15633 ymin: 7.683902 xmax: -72.88663 ymax: 10.8524
## Geodetic CRS:  WGS 84
  • Los datos se cargaron mediante el controlador de archivos shape de ESRI.

  • Los datos incluyen 24 objetos (el número de municipios en Cesar) y cuentan con 11 campos de atributos.

Cesar
## Simple feature collection with 24 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.15633 ymin: 7.683902 xmax: -72.88663 ymax: 10.8524
## Geodetic CRS:  WGS 84
## First 10 features:
##    ID_0 ISO   NAME_0 ID_1 NAME_1 ID_2          NAME_2    TYPE_2    ENGTYPE_2
## 1    53 COL Colombia   12  Cesar  445       Aguachica Municipio Municipality
## 2    53 COL Colombia   12  Cesar  446 Agustín Codazzi Municipio Municipality
## 3    53 COL Colombia   12  Cesar  447          Astrea Municipio Municipality
## 4    53 COL Colombia   12  Cesar  448        Becerril Municipio Municipality
## 5    53 COL Colombia   12  Cesar  449        Bosconia Municipio Municipality
## 6    53 COL Colombia   12  Cesar  450     Chimichagua Municipio Municipality
## 7    53 COL Colombia   12  Cesar  451      Chiriguaná Municipio Municipality
## 8    53 COL Colombia   12  Cesar  452        Curumaní Municipio Municipality
## 9    53 COL Colombia   12  Cesar  453        El Copey Municipio Municipality
## 10   53 COL Colombia   12  Cesar  454         El Paso Municipio Municipality
##    NL_NAME_2 VARNAME_2                           geom
## 1       <NA>      <NA> MULTIPOLYGON (((-73.4802 8....
## 2       <NA>      <NA> MULTIPOLYGON (((-72.95625 1...
## 3       <NA>      <NA> MULTIPOLYGON (((-73.8072 9....
## 4       <NA>      <NA> MULTIPOLYGON (((-72.95625 1...
## 5       <NA>      <NA> MULTIPOLYGON (((-73.7618 9....
## 6       <NA>      <NA> MULTIPOLYGON (((-73.6973 9....
## 7       <NA>      <NA> MULTIPOLYGON (((-73.20842 9...
## 8       <NA>      <NA> MULTIPOLYGON (((-73.38644 9...
## 9       <NA>      <NA> MULTIPOLYGON (((-74.0645 9....
## 10      <NA>      <NA> MULTIPOLYGON (((-73.8228 9....
  • El tipo de geometría es multipolígono. ¿Qué es un multipolígono?

Un MULTIPOLYGON (multipolígono) en el contexto de sistemas de información geográfica (SIG) y programación en R (por ejemplo, con el paquete sf) es una geometría vectorial que agrupa varios polígonos independientes dentro de un único objeto o registro.

  • El cuadro delimitador es:

st_bbox · devuelve el límite de un objeto simple o un conjunto de objetos simples.

  • El sistema de referencia de coordenadas (SRC) es MAGNA-SIRGAS. Consulta tu cuaderno físico para conocer los parámetros de este SRC.

Parámetros del Elipsoide (GRS 1980)

MAGNA-SIRGAS utiliza el elipsoide GRS80 (Geodetic Reference System 1980), cuyos valores oficiales son:

• Semieje mayor (a): 6 378 137 m (Radio ecuatorial) • Semieje menor (b): 6 356 752,31414 m (Radio polar) • Aplanamiento (1/f): 298,257222101 • Unidad angular: Grados sexagesimales / Radianes • Meridiano de referencia: Greenwich (0°)

Parámetros del SRC Unificado: Origen Nacional (CTM12)

Mediante resoluciones oficiales del IGAC (como la Resolución 471 de 2020 y 370 de 2021), Colombia adoptó un Origen Único Nacional para la cartografía y el catastro multipropósito (identificado con el código internacional EPSG: 9377):

• Nombre del sistema: MAGNA-SIRGAS 2018 / Origen-Nacional (CTM12) • Tipo de proyección: Transversa de Mercator secante • Elipsoide: GRS 1980 • Falso Este (Coordenada X de origen): 5 000 000 m • Falso Norte (Coordenada Y de origen): 2 000 000 m • Meridiano central (λ₀): \(-73^\circ\,00^\prime\,00^{''} \,(73^\circ\text{ W})\) • Latitud de origen (φ₀): \(4^\circ\,00^\prime\,00^{''} \,(4^\circ\text{ N})\) • Factor de escala (k₀): 0,9992

head(Cesar)
## Simple feature collection with 6 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.15633 ymin: 7.9975 xmax: -72.89888 ymax: 10.25736
## Geodetic CRS:  WGS 84
##   ID_0 ISO   NAME_0 ID_1 NAME_1 ID_2          NAME_2    TYPE_2    ENGTYPE_2
## 1   53 COL Colombia   12  Cesar  445       Aguachica Municipio Municipality
## 2   53 COL Colombia   12  Cesar  446 Agustín Codazzi Municipio Municipality
## 3   53 COL Colombia   12  Cesar  447          Astrea Municipio Municipality
## 4   53 COL Colombia   12  Cesar  448        Becerril Municipio Municipality
## 5   53 COL Colombia   12  Cesar  449        Bosconia Municipio Municipality
## 6   53 COL Colombia   12  Cesar  450     Chimichagua Municipio Municipality
##   NL_NAME_2 VARNAME_2                           geom
## 1      <NA>      <NA> MULTIPOLYGON (((-73.4802 8....
## 2      <NA>      <NA> MULTIPOLYGON (((-72.95625 1...
## 3      <NA>      <NA> MULTIPOLYGON (((-73.8072 9....
## 4      <NA>      <NA> MULTIPOLYGON (((-72.95625 1...
## 5      <NA>      <NA> MULTIPOLYGON (((-73.7618 9....
## 6      <NA>      <NA> MULTIPOLYGON (((-73.6973 9....

2b. Escritura de datos espaciales con st_write

Para exportar datos espaciales, la función st_write() permite guardar un objeto sf en diversos formatos.

Puede especificar el formato de salida añadiendo la extensión de archivo correspondiente o indicando el controlador GDAL adecuado. Se pueden especificar opciones adicionales para controlar cómo se escriben los datos.

★ Consejo: El parámetro append = F permite sobrescribir un archivo existente con el mismo nombre en el directorio de trabajo.

st_write(Cesar, "municipios.gpkg", driver = "GPKG", append = F)
## Deleting layer `municipios' using driver `GPKG'
## Writing layer `municipios' to data source `municipios.gpkg' using driver `GPKG'
## Writing 24 features with 11 fields and geometry type Multi Polygon.
# Check the new geopackage file was created
list.files(pattern="gpkg")
## [1] "cesar_munic.gpkg" "Mun_Cesar.gpkg"   "municipios.gpkg"  "stder_munic.gpkg"
# read the *Municipios* geopackage as an sf object
colombia2 <- st_read("./municipios.gpkg")
## Reading layer `municipios' from data source 
##   `D:\MI PC\GB02\Proyecto7\municipios.gpkg' using driver `GPKG'
## Simple feature collection with 24 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.15633 ymin: 7.683902 xmax: -72.88663 ymax: 10.8524
## Geodetic CRS:  WGS 84

2c. Tipos de datos espaciales compatibles con las funciones st_read y st_write

3. Selección de características

3.1 Selección por atributo

La herramienta “Seleccionar por atributo” en el software SIG permite introducir una consulta SQL para seleccionar las entidades que cumplen los criterios de selección.

Por otro lado, sf y dplyr ofrecen opciones bastante sencillas para realizar operaciones similares. Analicemos una de ellas.

Comprueba el número de entidades seleccionadas. ¿Coincide con la información previa sobre tu departamento?

A continuación, vamos a crear un gráfico muy básico de los municipios de Cesar:

plot(st_geometry(Cesar), col = sf.colors(12, categorical = TRUE), border = 'grey', axes = TRUE, las = 1)
plot(st_geometry(st_centroid(Cesar)), pch = 3, col = 'red', add = TRUE)
## Warning: st_centroid assumes attributes are constant over geometries

Nótese que hemos trazado los límites de los municipios, así como sus centroides.

Ahora, podemos guardar las entidades seleccionadas en un nuevo archivo geopackage para su uso posterior:

st_write(Cesar, "cesar_munic.gpkg", driver = "GPKG", append = F)
## Deleting layer `cesar_munic' using driver `GPKG'
## Writing layer `cesar_munic' to data source `cesar_munic.gpkg' using driver `GPKG'
## Writing 24 features with 11 fields and geometry type Multi Polygon.

5.2 Seleccionar por ubicación

Una operación muy común al trabajar con geodatos es la selección de entidades en función de su relación espacial con otro objeto (capa). En el software SIG, este proceso se suele realizar mediante la herramienta Seleccionar por ubicación. Esta herramienta permite seleccionar entidades según las siguientes relaciones:

Vamos a ilustrar la funcionalidad de selección por ubicación.

Primero, leamos un archivo CSV con ciudades del mundo. Sí, es el mismo archivo que usamos hace varias semanas en este curso. ¿Recuerdas el sistema de coordenadas de referencia (CRS) de este archivo?

cities = read.csv('D:/MI PC/GB02/Proyecto7/Cities.csv') %>% 
         st_as_sf(coords=c("lng","lat"), crs=4326) # remember x=longitude and y=latitude
cities
## Simple feature collection with 105 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -81.7006 ymin: -4.215 xmax: -67.4837 ymax: 12.5847
## Geodetic CRS:  WGS 84
## First 10 features:
##            city  country iso2         admin_name capital population
## 1        Bogotá Colombia   CO             Bogotá primary    7968095
## 2          Cali Colombia   CO    Valle del Cauca   admin    2916790
## 3      Medellín Colombia   CO          Antioquia   admin    2490164
## 4  Barranquilla Colombia   CO          Atlántico   admin    1326588
## 5     Cartagena Colombia   CO            Bolívar   admin     914552
## 6   Bucaramanga Colombia   CO          Santander   admin     893040
## 7       Palermo Colombia   CO              Huila   minor     800000
## 8        Cúcuta Colombia   CO Norte de Santander   admin     748948
## 9       Soledad Colombia   CO          Atlántico   minor     676014
## 10       Soacha Colombia   CO       Cundinamarca   minor     655025
##    population_proper                 geometry
## 1            8034649  POINT (-74.0817 4.6097)
## 2            2471474    POINT (-76.5197 3.44)
## 3            2529403  POINT (-75.5676 6.2502)
## 4            1326588 POINT (-74.8019 10.9833)
## 5             914552 POINT (-75.5144 10.3997)
## 6             614860  POINT (-73.1161 7.1186)
## 7             800000  POINT (-75.4339 2.8883)
## 8             806378  POINT (-72.5047 7.9075)
## 9             676014 POINT (-74.7858 10.9097)
## 10            782647  POINT (-74.2144 4.5781)

Es posible que hayamos notado que el sistema de referencia de coordenadas (CRS) de los dos conjuntos de datos (es decir, Santander y ciudades) no coincide. Vamos a demostrarlo:

st_crs(cities)$epsg
## [1] 4326
st_crs(Cesar)$epsg
## [1] 4326

Antes de utilizar cualquier selección por ubicación, necesitamos transformar el sistema de coordenadas de referencia (SCR) de las ciudades para que coincida con el SCR de Santander. Podemos hacerlo utilizando la función st_transform:

# take note of this function sintax
ncities <- st_transform(cities, crs= st_crs(Cesar))

Ahora, estamos listos para partir:

stder_cities <- ncities[Cesar, , op = st_within]

¿Obtuvimos lo que esperábamos?

Comprobémoslo de nuevo con un gráfico rudimentario:

plot(st_geometry(Cesar), col = sf.colors(12, categorical = TRUE), border = 'grey', axes = TRUE, las = 1)
plot(st_geometry(stder_cities), pch = 20, col = 'red', add = TRUE)

Nótese la variedad de símbolos que se pueden utilizar para representar geometrías de puntos:

4. Trazar con elegancia

4.1 Representación gráfica con ggplot

Usaremos ggplot() y llamaremos a geom_sf(), que le indica a R que estamos trabajando con datos espaciales. Al crear un mapa, también es fundamental considerar el orden de las capas. Si colocamos los datos de las ciudades (stder_cities) antes del polígono de municipios de Cesar (Cesar), las ciudades quedarán debajo del polígono de municipios y no podremos verlas.

ggplot() +
  #Add  municipalities
  geom_sf(data = Cesar) +
  #Add cities layer
  geom_sf(data = stder_cities, aes(color = city, label = city), size = 3) +
  #Add titles
  labs(x = "Longitud", y = "Latitud", title = "Ciudades de Cesar") +
  #Add theme
  theme_bw()

# 4.2 Representación gráfica mediante ggplot y ggspatial

Ahora, agregaremos algunos elementos personalizados a nuestro mapa, como una flecha que indique el norte y una barra de escala, utilizando el paquete ggspatial de R.

ggplot() +
  #Crop Virginia boundary to spatial extent of cities and add Virginia layer
  geom_sf(data = Cesar) +
  #Add cities layer
  geom_sf(data = stder_cities, aes(color = city, label = city), size = 3) +
  #Add scale bar to bottom left from ggspatial
  annotation_scale(location = "tr", height = unit(.25, "cm"), 
                   width = unit(1, "cm"), pad_x = unit(0.3, "in"), 
                   pad_y = unit(0.5, "in")) +
  #Add north arrow to bottom left from ggspatial
  annotation_north_arrow(height = unit(1, "cm"), width = unit(1, "cm"),
                         which_north = "true", location = "tr", 
                         pad_x = unit(0.5, "in"), pad_y = unit(0.05, "in")) +
  #Add titles
  labs(x = "Longitud", y = "Latitud", title = "Ciudades del Cesar") +
  #Add theme
  theme_bw()

Si obtuviste resultados similares a los mostrados anteriormente, ¡felicidades!

5. Entorno informático

sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22000)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggspatial_1.1.11 ggplot2_4.0.3    dplyr_1.2.1      sf_1.1-3        
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       jsonlite_2.0.0     compiler_4.6.1     tidyselect_1.2.1  
##  [5] Rcpp_1.1.2         jquerylib_0.1.4    scales_1.4.0       yaml_2.3.12       
##  [9] fastmap_1.2.0      R6_2.6.1           generics_0.1.4     classInt_0.4-11   
## [13] s2_1.1.13          knitr_1.52         tibble_3.3.1       units_1.0-1       
## [17] DBI_1.3.0          bslib_0.12.0       pillar_1.11.1      RColorBrewer_1.1-3
## [21] rlang_1.3.0        cachem_1.1.0       xfun_0.62          sass_0.4.10       
## [25] S7_0.2.2           otel_0.2.0         cli_3.6.6          withr_3.0.3       
## [29] magrittr_2.0.5     wk_0.9.5           class_7.3-23       digest_0.6.39     
## [33] grid_4.6.1         rstudioapi_0.19.0  lifecycle_1.0.5    vctrs_0.7.3       
## [37] KernSmooth_2.23-26 proxy_0.4-29       evaluate_1.0.5     glue_1.8.1        
## [41] farver_2.1.2       e1071_1.7-17       rmarkdown_2.32     tools_4.6.1       
## [45] pkgconfig_2.0.3    htmltools_0.5.9