Vamos a utilizar R como una aplicación GIS (Sistemas de información geográfica); es decir, lo vamos a utilizar para almacenar, manipular, analizar y presentar datos geográficos o geoespaciales. GIS es un área enorme, por ejemplo, aquí tienes un listado de herramientas y paquetes para trabajar con datos espaciales en R. Hay muchísimas!!!, es un área de conocimiento en sí misma, así que sólo voy a intentar hacer una pequeña introducción a GIS con R.
Durante este tutorial usaremos distintos data.frames con geometrías de:
Venga, vamos a hacer mapas con ellos!!! Sí, enseguida!! A todos nos gustan los mapas y hacer mapas con R es easy … bueno, no tan fácil, pero vamos a ello!!
Más que mapas lo que vamos a hacer es visualizar datos geoespaciales. Unos datos no se pueden considerar geoespaciales si no vienen acompañados de un CRS (Coordinate Reference System) de forma que las aplicaciones GIS puedan manipular y mostrar correctamente esos datos. Los CRS utilizan un modelo matemático para conectar los datos, las coordenadas de un objeto, con su posición física en la superficie terrestre. Un sistema de referencia de coordenadas geográficas permite que todos los puntos de la tierra puedan ser descritos como un conjuntos de coordenadas (como latitud, longitud y altitud). Diversos sistemas se utilizan para representar la tierra tridimensional en un plano o mapa bidimensional.
Los datos geoespaciales representan lugares geográficos de la Tierra y para poder representar la posición de un objeto en la Tierra necesitamos dos cosas:
las coordenadas del objeto
un sistema de referencia (CRS) para conectar esos datos con la superficie terrestre.
Un CRS consiste en modelo geométrico elipsoide de la superficie terrestre y un datum que identifica el origen y orientación de los ejes y también las unidades de medida.
Como sabéis la Tierra no es esférica y esto complica la utilización de un único CRS, de forma que se han creado unos cuantos para tratar de representar la superficie terrestre de forma adecuada.
Además, hay dos tipos de CRS: CRS geográficos y proyectados. En los sistema geográficos, el marco de referencia es el globo terráqueo y cuantos grados norte o este, está el objeto; mientras que en los sistemas proyectados, el marco de referencia es una representación plana de la Tierra (o de parte de ella). Si te interesa el tema, aquí y aquí tienes dos posts divulgativos.
Como la Tierra es esférica, para hacer un mapa bidimensional (o plano) de ella, se tiene que usar algún tipo de proyección, por lo que las áreas se mostrarán más o menos distorsionadas. Por ejemplo:
p1 <- ggplot(CCAA, aes(fill = ine_ccaa) ) +
geom_sf() +
scale_fill_viridis_d(guide = "none") +
labs(title = "Coord. geográficas")
p2 <- ggplot(CCAA, aes(fill = ine_ccaa) ) +
geom_sf() + coord_sf(crs = 2154, datum = sf::st_crs(2154)) +
scale_fill_viridis_d(guide = "none") +
labs(title = "Lambert 93")
p1 + p2
Para poder representar adecuadamente datos espaciales, estos deben contener la información de su CRS y algunas veces hay que convertir entre CRS’s; de forma que la estructura tradicional de datos en R (los data.frames) no es adecuada para almacenar datos espaciales1. Piensa además que los datos espaciales más simples (los puntos) serían fáciles de almacenar en dos columnas de un data.frame, pero piensa en objetos espaciales más complejos como lineas, polígonos, polígonos con huecos etc… hace que para almacenar datos espaciales en R se utilicen estructuras de datos más complejas que un data.frame.
Además, para complicar un poco más las cosas, hay dos tipos fundamentales de datos geográficos: raster y vectoriales. En el curso solo vamos a ver (un poquito) los datos espaciales vectoriales.
Los datos ráster se almacenan como una cuadrícula de valores que se representan en un mapa como píxeles. Cada valor de píxel representa un área en la superficie de la Tierra. Los datos ráster son datos pixelados (o cuadriculados) en los que cada píxel está asociado a una ubicación geográfica específica. El valor de un píxel puede ser continuo (por ejemplo, elevación) o categórico (por ejemplo, uso del suelo). Un ráster geoespacial solo es diferente de una foto digital en que está acompañado de información espacial que conecta los datos a una ubicación en particular. Esto incluye la extensión y el tamaño de celda del ráster, el número de filas y columnas y su sistema de referencia de coordenadas (o CRS)
Las estructuras de datos vectoriales representan posiciones específicas en la superficie de la Tierra.
sp versus sfEn R hay 2 formas alternativas de almacenar (y de manipular) datos
geográficos vectoriales: el paquete sp y el paquete
sf. Las principales diferencias entre los paquetes sp y sf
se derivan de la manera en que las clases sp y sf almacenan información
sobre CRS, distinguen entre puntos, líneas y polígonos, y cómo conectan
estos datos espaciales con datos no espaciales.
Cuando se usan estos paquetes, tanto sp como
sf, proporcionan enlaces a potentes bibliotecas SIG de
código abierto, como GDAL (Geospatial Data
Abstraction Library: es una biblioteca de software para la lectura y
escritura de formatos de datos geoespaciales, publicada bajo la MIT
License por la fundación geoespacial de código abierto (Open Source
Geospatial Foundation), GEOS
(Geometry Engine Open Source: para operaciones de topología en
geometrías) y PROJ (PROJ es una
biblioteca que proporciona métodos para transformar entre sistemas de
referencia de coordenadas diferentes y operaciones de proyección). Por
eso cuando cargamos el paquete sf nos aparece el siguiente
mensaje en la consola:
Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1
Muchos paquetes R para GIS utilizan la metodología/arquitectura
proporcionada por sp, pero por varias razones la
metodología implementada por el paquete el paquete sf ya es
prácticamente el estándar.
Las principales razones son 4:
El paquete sf implementa el estándar abierto “simple
feautures” para la representación de datos geográficos vectoriales en R.
Simple features es un estándar, un conjunto común de reglas, para
representar las geometrías espaciales de objetos del mundo real en un
ordenador, en las que se define desde cómo se almacenan los datos en el
ordenador hasta qué operaciones geométricas se pueden aplicar.
sf utiliza para contener/almacenar los datos
espaciales una estructura muy similar a un dataframe: los spatial
dataframes o sf2 , que son muy parecidos a un data frame
salvo que incluyen una columna adicional, generalmente llamada
geometry que almacena la información geográfica, los datos
vectoriales espaciales, ya sean estos puntos, lineas, polígonos
etc…
Muchas de las funciones del paquete sf comienzan con
st_; por ejemplo st_nearest_points(), lo que
facilita su uso.
las funciones de sf se integran en el
tidyverse workflow. De hecho la estructura de datos que
utiliza el paquete sf para almacenar datos geográficos son
los objetos de tipo `sf una extensión de los data.frames,
de forma que pueden ser manipulados con dplyr y graficados
con ggplot2. Los objetos sf pueden almacenar
datos normales (cuantitativos, cualitativos, texto) pero almacenan datos
espaciales en una columna especial llamada habitualmente
geometry con clase especial (sfc), y allí se
almacena la información geográfica: el CRS, coordenadas y el tipo de
objeto geométrico. La clase sfc tiene siete subclases para denotar los
distintos tipos de objetos geométricos que se derivan del estándar
sf.
filter():#- elegimos 4 CC.AA
zz_prov_chosen <- c("Andalucía", "Galicia", "Aragón", "Canarias", "Illes Balears")
zz_prov_1 <- Provincias %>% filter(ine_ccaa.n %in% zz_prov_chosen)
zz_prov_2 <- Provincias %>% filter(!(ine_ccaa.n %in% zz_prov_chosen))
p1 <- ggplot(zz_prov_1, aes(fill = ine_prov) ) +
geom_sf() +
scale_fill_viridis_d(guide = "none") +
labs(title = "Coord. geográficas")
p2 <- ggplot(zz_prov_2, aes(fill = ine_prov) ) +
geom_sf() +
scale_fill_viridis_d(guide = "none") +
labs(title = "Coord. geográficas")
p1 + p2
summarize()#- creamos las geometrías de las CC.AA por agrupación de las geometrías provinciales
CCAA_2 <- Provincias %>% group_by(ine_ccaa, ine_ccaa.n) %>% summarize()
p1 <- ggplot(Provincias, aes(fill = ine_prov) ) +
geom_sf() +
scale_fill_viridis_d(guide = "none") +
labs(title = "Coord. geográficas")
p2 <- ggplot(CCAA_2, aes(fill = ine_ccaa) ) +
geom_sf() +
scale_fill_discrete(guide = "none") +
labs(title = "Coord. geográficas")
p1 + p2
Bueno, pues este es toda la teoría del tutorial. Vamos a hacer gráficos!!
Primero vamos cargar un conjunto de datos espaciales del paquete rnaturalearth.
Bueno, en realidad ya lo hemos cargado antes.
#- antes usaba las worlg geometies del paquete rnaturalearth, pero ...
# world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
#- pero ahora prefiero utilizar las geometrias del paquete tmap,
#- las guarde en noviembre de 2021. Se cargaban así:
library(tmap)
data(World) #- hacemos accesibles 1 df
world <- World ; rm(World)
Veamos que hay en world
class(world) #- es un sf, pero tb es un data.frame
#> [1] "sf" "data.frame"
names(world) #- fijate que la ultima columna se llama geometry y almacena los datos espaciales
#> [1] "iso_a3" "name" "sovereignt" "continent" "area"
#> [6] "pop_est" "pop_est_dens" "economy" "income_grp" "gdp_cap_est"
#> [11] "life_exp" "well_being" "footprint" "inequality" "HPI"
#> [16] "geometry"
Los gráficos los vamos a hacer, de momento, con ggplot2.
Es tan fácil como usar geom_sf(). Las siguientes 3
expresiones producen el mismo gráfico:
ggplot(data = world) + geom_sf()
ggplot(data = world, aes(geometry = geometry)) + geom_sf()
ggplot() + geom_sf(data = world, aes(geometry = geometry))
Como el mapa que hemos hecho, en realidad es un gráfico ggplot, pues ya sabemos cambiar/tunear algunas de las características del gráfico. Por ejemplo:
ggplot(data = world) + geom_sf(color = "black", fill = "lightgreen", lwd = 0.2)
Podemos poner títulos y leyendas:
p <- ggplot(data = world) + geom_sf() +
labs(title = "Gráfico 1: Mapa del mundo",
caption = "Datos provenientes de tmap")
p
Podemos hacer mapas de COROPLETAS: las distintas regiones se colorean de forma que muestren los valores de una variable estadística. Por ejemplo, una coropleta con la variable población(pob_est)
p + geom_sf(aes(fill = pop_est))
p + geom_sf(aes(fill = pop_est)) + scale_fill_distiller()
p + geom_sf(aes(fill = pop_est)) + scale_fill_viridis_c(option = "plasma", trans = "sqrt")
Para cambiar la proyección de un objeto “sf” a otro CRS usamos la
función st_transform(), o alternativamente podemos
especificar directamente el CRS en el gráfico con
coord_sf():
#- world2 <- st_transform(world, "+proj=laea +y_0=0 +lon_0=155 +lat_0=-90 +ellps=WGS84 +no_defs")
p + coord_sf(crs = "+proj=laea +y_0=0 +lon_0=155 +lat_0=-90 +ellps=WGS84 +no_defs") +
labs(subtitle = "Lambert Azimuthal Equal Area projection")
Vamos a jugar un poco con el CRS con la función
coord_sf(). Hay que saber que ggplot usará el CRS de la
primera capa que contenga uno. Si no hubiese ningún CRS se usa
WGS84 (latitude/longitude, the reference system used in
GPS).
Se puede cambiar el CRS de varias maneras:
usando un
**string PROJ** válido. Por ejemplo este string“+proj=laea
+lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m
+no_defs”` es el European-centric ETRS89 Lambert Azimuthal Equal-Area
projection)
usando un código SRID (Spatial Reference System Identifier)
Usando un código EPGS (European Petroleum Survey Group)3
Los siguientes tres gráficos hacen lo mismo. Representan los datos
del spatial.data.frame world usando the ETRS89 Lambert
Azimuthal Equal-Area projection, cuyo código EPSG es 3035. Aquí
tienes algunos códigos de proyecciones:
p + coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ") #- usamos PROJ string
p + coord_sf(crs = st_crs(3035)) #- usamos SRID identifier
p + coord_sf(crs = "+init=epsg:3035") #- usamos EPSG code (mejor no usarlo xq GDAL ha deprecated esta sintaxis)
p + coord_sf(crs = st_crs(3035)) + #- usamos SRID identifier
labs(subtitle = "(European-centric ETRS89 \nLambert Azimuthal Equal-Area projection)")
Volvemos a la proyección nativa del conjunto de datos
world y hacemos un zoom en nuestro mapa con
coord_sf(). El zoom lo hacemos pare centrarnos en España y
alrededores. Para hacer el zoom en realidad ajustamos los límites de los
ejes del gráfico
p_esp <- p + coord_sf(xlim = c(-15.00, 15.00), ylim = c(21, 47.44), expand = TRUE)
p_esp
Con el paquete ggspatial podemos poner elementos de
mapas como una Scale bar con annotation_scale() y
una North arrow con annotation_north_arrow():
p_esp +
ggspatial::annotation_scale(location = "bl", width_hint = 0.5) +
ggspatial::annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering)
El mapa se ha hecho con ggplot2, así que podemos añadir,
por ejemplo, texto con geom_text(); pero antes de ello,
vamos a crear/calcular centroides de cada país para ver donde situamos
el texto en el mapa.
world_points <- st_centroid(world, of_largest_polygon = TRUE) #- sustituye la geometry por el centroide
De esta forma, con st_centroid() podemos calcular los
centroides de cada geometría, de cada país, pero Los datos del centroide
(longitud, latitud) se almacenan en la columna geometry; es
decir sustituyen o eliminan las geometrías anteriores que eran los
límites o fronteras de los países.
Si queremos que las geometrías originales y los centroides estén en un mismo data.frame, hemos de hacer:
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry, of_largest_polygon = TRUE))) #- ahora el centroide pasa a estar en 2 columnas, llamadas X e Y, donde figuran la longitud y latitud del centroide
Lo que hacemos es fusionar los los dos data.frame’s, y el centroide pasará a estar en dos nuevas columnas llamadas X e Y, sin sustituir las geometrías de los países.
Para ver esas dos nuevas columnas (X e Y), vamos a hacer un “pequeño truco”, vamos a quitar la columna de geometrías. Esto facilita el visionado de los datos, pero lógicamente cuando vayamos a hacer el gráfico, realmente no podremos eliminar la geometry porque si no, los datos dejarían de ser datos espaciales.
#- truco para ver mejor los datos espaciales: quitar la geometría
names(world_points)
#> [1] "iso_a3" "name" "sovereignt" "continent" "area"
#> [6] "pop_est" "pop_est_dens" "economy" "income_grp" "gdp_cap_est"
#> [11] "life_exp" "well_being" "footprint" "inequality" "HPI"
#> [16] "X" "Y" "geometry"
zz <- world_points %>% select(name, iso_a3, X, Y) %>%
st_set_geometry(NULL) #- por si quiero ver mejor los datos del df
Hacemos el mapa con los nombres de los países con
geom_text() y el nombre del Mar Mediterráneo con
annotate()
p_esp +
geom_text(data = world_points,
aes(x = X, y = Y, label = name),
color = "darkblue", fontface = "bold", check_overlap = TRUE, size = 2.5) +
annotate(geom = "text", x = 5.5, y = 38,
label = "Mar Mediterráneo", fontface = "italic", color = "grey22", size = 4)
Vamos a tunear el mapa. Por ejemplo, vamos a poner el mar (en
realidad al background del theme) del color
aliceblue, y a los continentes les pondremos un color
terroso con fill = "antiquewhite". Arreglamos un poco las
grid-lines dentro del theme, y le ponemos la estrellita. No
nos va a salir bonito-bonito porque los nombres de los países se van a
superponer. Podríamos quitar nombres, pero otro día:
p <- p +
geom_sf(fill = "antiquewhite") +
geom_text(data = world_points,
aes(x = X, y = Y, label = name),
color = "darkblue", fontface = "bold",
check_overlap = TRUE, size = 2) +
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"),
pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
theme(panel.grid.major = element_line(color = gray(.5), linetype = "dashed", size = 0.5),
panel.background = element_rect(fill = "aliceblue"))
p