2 Simple Features

El paquete sf(Pebesma, 2018, 2022), abreviatura de Simple Features, es una herramienta moderna y poderosa en el lenguaje de programación R diseñada específicamente para el manejo y análisis de datos geoespaciales. Este paquete proporciona un enfoque simple y coherente para trabajar con datos espaciales, permitiendo a los usuarios de R manipular, visualizar y analizar datos espaciales de manera eficiente y efectiva.

Desarrollado para reemplazar y ampliar las capacidades de paquetes más antiguos como sp(Pebesma y Bivand, 2023) y rgdal(Bivand et al., 2017), sf ofrece una integración completa con el sistema de objetos de R y trabaja de manera armoniosa con paquetes populares en el ecosistema de R como dplyr(Wickham, François, et al., 2023), tidyr(Wickham, Vaughan, et al., 2023) y ggplot2(Wickham, 2016). Esto facilita la realización de tareas comunes de manipulación de datos, como la selección, filtración y agregación, directamente en conjuntos de datos espaciales.

El paquete sf soporta múltiples tipos de geometrías espaciales, incluyendo puntos, líneas, polígonos y sus variantes combinadas, lo que lo hace adecuado para una amplia gama de aplicaciones geoespaciales. La estructura de datos de sf está alineada con el estándar Simple Features de la Open Geospatial Consortium, lo que garantiza la compatibilidad con una amplia variedad de herramientas y formatos de datos espaciales.

Las capacidades clave de sf incluyen:

Lectura y escritura eficiente: sf puede leer y escribir una variedad de formatos de archivo geoespacial comunes, incluyendo Shapefile, GeoJSON, KML y muchos otros, a través de la integración con la biblioteca GDAL.
Manipulación de datos espaciales: Permite realizar operaciones espaciales como el cálculo de distancias, áreas, intersecciones, uniones y otros tipos de análisis geométrico.
Visualización de datos espaciales: Integración con ggplot2 para la creación de visualizaciones de datos espaciales avanzadas y de alta calidad.

2.1 Instalación del paquete sf

Del mismo modo que hacemos con otros paquetes, para poder utilizarlo hay que instalarlo y cargarlo.

install.packages("sf")

library(sf)

Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE

Advertencia

Cuando cargamos la librería, nos aparece un mensaje informativo. Hablaremos de esto más adelante pero, ¿a qué crees que hace referencia?

Recordemos una forma de comprobar qué versión del paquete tenemos instalada:

packageVersion("sf")

[1] '1.0.19'

2.2 Cómo crear un objeto sf

Los objetos sf se forman como integración de unas piezas o subclases que heredan las características de la norma UNE ISO19125. Básicamente un objeto sf se descompone en dos objetos: uno de tipo sfc (simple feature column) y otro de tipo data.frame. A su vez, un objeto sfc se origina a partir de un objeto sfg (simple feature geometry), pudiéndole añadir un CRS (sistema de coordenadas de referencia). Los objetos sfg pueden ser de tipo POINT, LINESTRING, etc. Veamos más detenidamente cada tipo de objeto, comenzando por los más básicos hasta llegar a sf.

2.2.1 sfg

La clase sfg almacena la geometría con la información gráfica del objeto: las coordenadas, dimensiones, tipo de geometría… Hay 17 tipos posibles, pero los 7 más utilizados son:
- POINT: un punto simple (un vector x,y).
- MULTIPOINT: múltiples puntos (una matriz cada fila un punto).
- LINESTRING: secuencia de 2 o mas puntos conectados por lineas rectas.
- MULTILINESTRING: múltiples lineas (lista de matrices).
- POLYGON: un polígono cerrado (lista de matrices).
- MULTIPOLYGON: múltiples polígonos.
- GEOMETRYCOLLECTION: cualquier combinación de los anteriores.

Para crear cualquiera de estas geometrías básicas la librería sf contiene una serie de funciones creadoras: st_point(), st_multipoint(), st_linestring(), st_polygon(). Veamos algunos ejemplos.

2.2.1.1 Point

sf::st_point(c(5, 2))                 # XY point

POINT (5 2)

sf::st_point(c(5, 2, 3))              # XYZ point

POINT Z (5 2 3)

sf::st_point(c(5, 2, 1), dim = "XYM") # XYM point

POINT M (5 2 1)

sf::st_point(c(5, 2, 3, 1))           # XYZM point

POINT ZM (5 2 3 1)

Los resultados muestran que los tipos de puntos XY (coordenadas 2D), XYZ (coordenadas 3D) y XYZM (3D con una variable adicional) se crean a partir de vectores de longitud 2, 3 y 4, respectivamente. El tipo XYM debe especificarse utilizando el argumento dim (que es la abreviatura de dimensión).

2.2.1.2 Multipoint

Para generar objetos tipo multipoint o linestring, trabajaremos con matrices, en lugar de con vectores. La función rbind() facilita esta tarea.

punto1 <- sf::st_point(c(5,2))
punto2 <- sf::st_point(c(4,3))
punto1

POINT (5 2)

class(punto1)

[1] "XY"    "POINT" "sfg"

plot(punto1, axes = T, xlab = "X", ylab = "Y", pch=19)

matriz_puntos <- rbind(punto1,punto2)
matriz_puntos

       [,1] [,2]
punto1    5    2
punto2    4    3

class(matriz_puntos)

[1] "matrix" "array"

multipunto <- sf::st_multipoint(matriz_puntos)
multipunto

MULTIPOINT ((5 2), (4 3))

class(multipunto)

[1] "XY"         "MULTIPOINT" "sfg"

par(mar = c(2,2,1,2))
plot(multipunto, axes = T, xlab = "X", ylab = "Y", pch=19)

2.2.1.3 Linestring

linea <- sf::st_linestring(matriz_puntos)
linea

LINESTRING (5 2, 4 3)

class(linea)

[1] "XY"         "LINESTRING" "sfg"

par(mar = c(2,2,1,2))
plot(linea, axes = T, xlab = "X", ylab = "Y")

2.2.1.4 Multilinestring

Para generar objetos de tipo multiniestring, polygon, mulipolygon o geometrycollection trabajaremos con listas.

lineas_list <- list(rbind(c(1, 5), c(5, 4), c(3, 2), c(2, 4), c(0, 4)),
                    rbind(c(0, 3), c(1, 1)))
lineas <- st_multilinestring(lineas_list)
class(lineas)

[1] "XY"              "MULTILINESTRING" "sfg"

par(mar = c(2,2,1,2))
plot(lineas, axes = T, xlab = "X", ylab = "Y")

2.2.1.5 Polygon

pol1_list <- list(rbind(c(1, 5), c(5, 4), c(3, 2), c(2, 4), c(1, 5)))
pol1 <- sf::st_polygon(pol1_list)
pol1

POLYGON ((1 5, 5 4, 3 2, 2 4, 1 5))

class(pol1)

[1] "XY"      "POLYGON" "sfg"

par(mar = c(2,2,1,2))
plot(pol1, axes = T, xlab = "X", ylab = "Y")

Pregunta

Este polígono tiene 4 vértices, pero lo hemos construido a partir de 5 puntos… ¿Qué ha pasado? ¿Qué pasa si quitamos el quinto punto?

2.2.1.6 Multipolygon

multipolygon_list <- list(list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
                          list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2))))
multipoly <- st_multipolygon(multipolygon_list)
multipoly

MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5)), ((0 2, 1 2, 1 3, 0 3, 0 2)))

class(multipoly)

[1] "XY"           "MULTIPOLYGON" "sfg"

par(mar = c(2,2,1,2))
plot(multipoly, axes = T, xlab = "X", ylab = "Y")

2.2.1.7 Geometrycollection

geom_collection_list <- list(sf::st_multipoint(matriz_puntos),
                             sf::st_linestring(linea))
geom_collection <- sf::st_geometrycollection(geom_collection_list)
geom_collection

GEOMETRYCOLLECTION (MULTIPOINT ((5 2), (4 3)), LINESTRING (5 2, 4 3))

class(geom_collection)

[1] "XY"                 "GEOMETRYCOLLECTION" "sfg"

par(mar = c(2,2,1,2))
plot(geom_collection, axes = T, xlab = "X", ylab = "Y")

2.2.2 sfc

Como hemos visto en el bloque anterior, un objeto sfg contiene solo una geometría simple (es, por tanto, la unidad más básica de una geometría en el paquete sf). Una sfc (simple feature column) es una lista de objetos sfg (una especie de contenedor que agrupa una o más geometrías sfg), y que además puede contener información sobre el sistema de referencia de coordenadas en uso. Para combinar dos (o más) geometrías simples en un único objeto, podemos usar la función st_sfc().

conjunto <- sf::st_sfc(punto1, punto2, linea, lineas, pol1, multipoly)
conjunto

Geometry set for 6 features 
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 0 ymin: 1 xmax: 5 ymax: 5
CRS:           NA
First 5 geometries:

POINT (5 2)

POINT (4 3)

LINESTRING (5 2, 4 3)

MULTILINESTRING ((1 5, 5 4, 3 2, 2 4, 0 4), (0 ...

POLYGON ((1 5, 5 4, 3 2, 2 4, 1 5))

class(conjunto)

[1] "sfc_GEOMETRY" "sfc"

par(mar = c(2,2,1,2))
plot(conjunto, axes = T, xlab = "X", ylab = "Y")

Ejercicio

Crea la siguiente figura:

2.2.3 sfc + CRS

Pregunta

Observemos el objeto creado… ¿Hay algo que te llame la atención?

Para ver qué CRS tiene asignado un objeto, usaremos la función st_crs:

sf::st_crs(geom_col)

Coordinate Reference System: NA

Como dijimos anteriormente, es posible asignar un sistema de coordenadas de referencia al objeto sfc. De hecho, suele ser lo habitual.

Para asignar un CRS a un objeto que aun no tiene asignado ninguno, también se usa la función st_crs:

sf::st_crs(geom_col) <- 4326

Para modificar el CRS de un objeto que ya tiene uno asignado, se usa la función st_transform():

geom_col2 <- sf::st_transform(geom_col, crs = 4258)
sf::st_crs(geom_col2)

Coordinate Reference System:
  User input: EPSG:4258 
  wkt:
GEOGCRS["ETRS89",
    ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
        MEMBER["European Terrestrial Reference Frame 1989"],
        MEMBER["European Terrestrial Reference Frame 1990"],
        MEMBER["European Terrestrial Reference Frame 1991"],
        MEMBER["European Terrestrial Reference Frame 1992"],
        MEMBER["European Terrestrial Reference Frame 1993"],
        MEMBER["European Terrestrial Reference Frame 1994"],
        MEMBER["European Terrestrial Reference Frame 1996"],
        MEMBER["European Terrestrial Reference Frame 1997"],
        MEMBER["European Terrestrial Reference Frame 2000"],
        MEMBER["European Terrestrial Reference Frame 2005"],
        MEMBER["European Terrestrial Reference Frame 2014"],
        MEMBER["European Terrestrial Reference Frame 2020"],
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[0.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]],
    USAGE[
        SCOPE["Spatial referencing."],
        AREA["Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Moldova; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain; Sweden; Switzerland; United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
        BBOX[32.88,-16.1,84.73,40.18]],
    ID["EPSG",4258]]

Pregunta

¿Qué CRS conoces?

2.2.4 sf

Siguiendo el diagrama inicial, nos queda un último paso para crear un objeto sf. Para ello, lo único que faltaría sería añadir un data.frame al objeto sfc que ya tenemos.

Pregunta

¿Qué información le añadirías?

data <- data.frame(nombre = c("Valencia","Madrid","Barcelona","Ebro","Guadiana","España"),
                   tipo = c(rep("ciudad",3), rep("rio",2), "pais"))
data

objeto_sf <- sf::st_sf(data, geom_col)
class(objeto_sf)

[1] "sf"         "data.frame"

Ahora ya tenemos nuestro objeto sf, el cual podemos representar gráficamente con ggplot2:

library(ggplot2)

ggplot() +
  geom_sf(data = objeto_sf)

Pregunta

¿Qué problemas ves en el gráfico?

Digamos que el objeto sf que hemos creado “mezcla” información (datos) de diversa índole. Lo recomendable es tener un objeto para cada tipo de datos (puntos, líneas, polígonos) y luego usar lo que necesitemos.

library(dplyr)
sf1 <- objeto_sf %>% filter(tipo == "ciudad")
sf2 <- objeto_sf %>% filter(tipo == "rio")
sf3 <- objeto_sf %>% filter(tipo == "pais")

Teniendo la información estructurada correctamente, podemos representarla mejor:

ggplot() +
  geom_sf(data = sf3, fill="white") +
  geom_sf(data = sf2, color="blue") +
  geom_sf(data = sf1, shape=19, color="red")

2.2.5 Descomponer un objeto sf

Ya hemos visto los pasos para crear un objeto sf pero, ¿cómo podemos extraer los datos (data.frame)? ¿y las geometrías?

Para lo primero usaremos la función st_set_geometry(). Para lo segundo, la función st_geometry():

datos <- sf::st_set_geometry(objeto_sf, NULL)
class(datos)

[1] "data.frame"

datos

geom <- sf::st_geometry(objeto_sf)
class(geom)

[1] "sfc_GEOMETRY" "sfc"

2.2.6 Resumen

En este apartado hemos visto las siguientes funciones del paquete sf:

st_point
st_multipoint
st_linestring
st_multilinestring
st_polygon
st_multipolygon
st_geometrycollection
st_sfc
st_crs
st_transform
st_sf
st_set_geometry
st_geometry
st_drop_geometry

Ejercicio

¿Para qué sirve cada una de estas funciones?

2.3 Crecimiento urbano

La Dirección General del Catastro dispone de información espacial de toda la edificación de España, a excepción del País Vasco y Navarra. El conjunto de datos con el que vamos a trabajar forma parte de la implantación de INSPIRE, la Infraestructura de Información Espacial en Europa. Podemos encontrar más información aquí. Podríamos utilizar los enlaces (urls) disponibles en la web del catastro en formato ATOM (un formato de redifusión de tipo RSS). No obstante, para agilizar esta parte, trabajaremos con un fichero rds disponible aquí.

Como siempre, comenzamos cargando las librerías necesarias:

library(dplyr)
library(ggplot2)
library(lubridate)
library(stringr)
library(sf)

Ejercicio

¿Para qué sirve cada una de estas librerías?

Importamos el fichero con los datos (cartografía de edificios y otras variables de interés):

buildings <- readRDS("edificios_valencia.rds")
class(buildings)

[1] "sf"         "data.frame"

dplyr::glimpse(buildings)

Rows: 36,328
Columns: 27
$ gml_id                                  <chr> "ES.SDGC.BU.000100…
$ lowerCorner                             <chr> "725589.5805 43771…
$ upperCorner                             <chr> "725605.5005 43771…
$ beginLifespanVersion                    <chr> "2008-10-20T00:00:…
$ conditionOfConstruction                 <chr> "functional", "fun…
$ beginning                               <chr> "1890-01-01T00:00:…
$ end                                     <chr> "1890-01-01T00:00:…
$ endLifespanVersion                      <chr> NA, NA, NA, NA, NA…
$ informationSystem                       <chr> "https://www1.sede…
$ reference                               <chr> "000100100YJ27F", …
$ localId                                 <chr> "000100100YJ27F", …
$ namespace                               <chr> "ES.SDGC.BU", "ES.…
$ horizontalGeometryEstimatedAccuracy     <dbl> 0.1, 0.1, 0.1, 0.1…
$ horizontalGeometryEstimatedAccuracy_uom <chr> "m", "m", "m", "m"…
$ horizontalGeometryReference             <chr> "footPrint", "foot…
$ referenceGeometry                       <lgl> TRUE, TRUE, TRUE, …
$ currentUse                              <chr> "1_residential", "…
$ numberOfBuildingUnits                   <int> 1, 1, 1, 1, 1, 2, …
$ numberOfDwellings                       <int> 1, 0, 0, 1, 1, 0, …
$ numberOfFloorsAboveGround               <chr> NA, NA, NA, NA, NA…
$ documentLink                            <chr> "http://ovc.catast…
$ format                                  <chr> "jpeg", "jpeg", "j…
$ sourceStatus                            <chr> "NotOfficial", "No…
$ officialAreaReference                   <chr> "grossFloorArea", …
$ value                                   <int> 201, 934, 35, 615,…
$ value_uom                               <chr> "m2", "m2", "m2", …
$ geometry                                <MULTIPOLYGON [m]> MULTI…

2.3.1 Preparar los datos

No nos vamos a detener demasiado en este punto (aunque podríamos hacerlo…). Por el momento, únicamente vamos a convertir la variable beginning (fecha de construcción del edificio) en clase Date. Esta variable contiene algunas fechas en formato --01-01, es decir, no tienen ningún año asignado. Por eso, reemplazamos el primer - por 0000.

buildings <- buildings %>%
  mutate(beginning2 = stringr::str_replace(beginning, "^-", "0000")) %>%
  mutate(beginning3 = lubridate::ymd_hms(beginning2)) %>%
  mutate(beginning4 = lubridate::as_date(beginning3))

sum(is.na(buildings$beginning4))

[1] 6

Pregunta

¿Cómo podríamos mejorar el código anterior?

Descartamos las observaciones con valores ausentes para la variable de interés:

buildings <- buildings %>% filter(!is.na(beginning4))

2.3.2 Gráfico de distribución

Antes de crear el mapa de la edad del edificado, lo que reflejará el crecimiento urbano, haremos un gráfico de distribución de la fecha de construcción de los edificios. De este modo, podremos identificar fácilmente los períodos de expansión urbana. Usaremos el paquete ggplot2(Wickham et al., 2022) con la geometría de geom_density() para este objetivo.

options(scipen = 999) # evita notación científica

buildings %>%
  ggplot(aes(beginning4)) + 
  geom_density(color = "#F06720", fill = "#E6DB61", alpha = 0.5) +
  scale_x_date(breaks = seq(ymd("1800-01-01"),ymd("2020-01-01"), by = "20 years"),
               labels = scales::date_format("%Y")) +
  coord_cartesian(xlim = c(ymd("1800-01-01"),ymd("2020-01-01"))) +
  theme_minimal() +
  labs(y = "", x = "", title = "Evolución del desarrollo urbano")

2.3.3 Discretización

Para poder visualizar bien las diferentes épocas de crecimiento, categorizamos el año de construcción de los edificios en varios grupos. También es posible modificar el número de clases o bien el método aplicado (p.ej. jenks, fisher, etc.). Para más detalles ver la ayuda (?classIntervals).

n <- 11
breaks <- classInt::classIntervals(var = year(buildings$beginning4),
                                   n = n,
                                   style = "quantile")

labels <- names(print(breaks, under = "<", over = ">", cutlabels = FALSE))

style: quantile
     < 1908 1908 - 1922 1922 - 1930 1930 - 1942 1942 - 1955 
       3289        3280        2140        4425        3002 
1955 - 1963 1963 - 1968 1968 - 1974 1974 - 1983 1983 - 1999 
       3579        2952        3608        3356        3288 
     > 1999 
       3403

buildings <- mutate(buildings, 
                    yr_cl = cut(year(beginning4), 
                                breaks$brks, 
                                labels = labels, 
                                include.lowest = TRUE))

2.3.4 Mapa

Antes de representar los años de construcción de los edificios de Valencia en un mapa, vamos a seleccionar únicamente los edificios situados a cierta distancia del centro de la ciudad. Para ello usamos algunas funciones del paquete sf (las veremos con detalle más adelante):

centro_valencia <- sf::st_as_sf(data.frame(lon = -0.37680, lat = 39.46979), 
                                coords = c("lon", "lat"), 
                                crs = 4326)

centro_valencia <- sf::st_transform(centro_valencia, sf::st_crs(buildings))

buffer <- sf::st_buffer(centro_valencia, 3000)

buildings2 <- sf::st_intersection(buildings,buffer)

Con esto, ya estamos en disposición de crear nuestro primer mapa:

ggplot(data = buildings2) +
  geom_sf(aes(fill = yr_cl), colour = "transparent") +
  labs(title = "VALENCIA", fill = "") +
  guides(fill = guide_legend(keywidth = .7, keyheight = 1.6)) +
  scale_fill_manual(values = RColorBrewer::brewer.pal(n,"Spectral")) +
  #RColorBrewer::display.brewer.all()
  theme_void() +
  theme(panel.background = element_rect(fill = "black"),
        plot.background = element_rect(fill = "black"),
        legend.justification = .5,
        legend.text = element_text(colour = "white", size = 8),
        plot.title = element_text(colour = "white", hjust = .5,
                                  size = 50, fontface="bold", margin = margin(t = 25)),
        plot.title.position = "plot",
        plot.margin = margin(r = 40, l = 40))

2.3.5 Agradecimientos

Esta sección se ha basado en el trabajo realizado por Dominic Royé y Roberto Serrano-Notivoli.