library(tidyverse)
library(sf)
library(terra)
library(mapSpain)
options(scipen = 999)18 Datos espaciotemporales
18.1 Introducción
La información espacio-temporal es clave en muchas disciplinas, como la geografía, la planificación urbana, la climatología o la meteorología, entre otras. Ello hace necesario disponer de un formato de datos que permita una estructura multidimensional. Además es importante que ese formato tenga un alto grado de compatibilidad de intercambio y pueda almacenar un elevado número de datos. Estas características llevaron al desarrollo del estándar abierto NetCDF (Network Common Data Form). El formato NetCDF es un estándar abierto de intercambio de datos científicos multidimensionales que se utiliza con datos de observaciones o modelos, principalmente en disciplinas como la climatología, la meteorología y la oceanografía. La convención NetCDF es gestionada por Unidata. Se trata de un formato espacio-temporal con una cuadrícula regular o irregular. La estructura multidimensional en forma de array permite usar no sólo datos espacio-temporales, sino también multivariables. Las características generales del NetCDF se refieren al uso de un sistema de coordenadas n-dimensional, de múltiples variables y de una rejilla (malla o grid) regular o irregular. Además se incluyen metadatos que describen los contenidos. La extensión del formato NetCDF es “.nc”.
En esta práctica trabajaremos con datos de sequía de España, en formato NetCDF y con una resolución de 1 km. Para medir la sequía, existen diferentes indicadores, siendo el índice SPEI (Standardized Precipitation-Evapotranspiration Index) uno de los más utilizados. El Consejo Superior de Investigaciones Científicas (CSIC) pone a disposición de la comunidad científica datos desde 1960. Nosotros haremos uso de un subconjunto: el año 2017 del SPEI de 12 meses, que puedes descargar de aquí.
18.2 Librerías
El manejo de datos en formato NetCDF es posible a través de varios paquetes de forma directa o indirecta. Destaca el paquete ncdf4(Pierce, 2023) específicamente diseñado, del que hacen uso también otros paquetes, aunque no lo veamos. El manejo con ncdf4 es algo complejo, particularmente por la necesidad de gestionar la memoria RAM cuando tratamos grandes conjuntos de datos o también por la forma de manejar la clase array. Otro paquete muy potente es terra(Hijmans, 2023), que utilizamos cuando trabajamos con datos raster, y que permite usar sus funciones también para el manejo del formato NetCDF.
18.3 Preparando los datos
Una vez que hemos descargado el conjunto de datos (fichero spei12_2017.nc), lo importamos al entorno R mediante la función rast(). En realidad esta función no importa los datos, sino que establece el acceso a tales datos (así evitamos saturar la memoria de nuestro equipo).
spei <- terra::rast("data/spei12.nc")
speiclass : SpatRaster
dimensions : 834, 1115, 3024 (nrow, ncol, nlyr)
resolution : 1100, 1100 (x, y)
extent : -80950, 1145550, 3979450, 4896850 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs
source : spei12.nc
varname : spei
names : spei_1, spei_2, spei_3, spei_4, spei_5, spei_6, ...
unit : standardized values, standardized values, standardized values, standardized values, standardized values, standardized values, ...
time (days) : 1961-01-01 to 2023-12-23 Vemos en los metadatos el número de capas (layers) disponibles. El índice SPEI-12 está calculado semanalmente con 4 semanas por mes.
Si nos fijamos en los metadatos, falta la definición del sistema de coordenadas. Asignamos el código EPSG:25830 (ETRS89/UTM 30N):
# definimos el sistema de coordenadas
terra::crs(spei) <- "EPSG:25830"Seleccionamos las capas de 2017:
fechas <- as.Date(terra::time(spei))
capas <- which(year(fechas) == 2017)
spei <- spei[[capas]]Graficamos las primeras semanas:
plot(spei)
Atención
En caso de que el mapa se mostrara invertido, podemos hacer uso de la función terra::flip(). Por ejemplo, si España aparece “boca abajo”, la sentencia sería terra::flip(spei, direction="vertical").
Existen diferentes funciones para acceder a los metadatos, como las fechas, los nombres de las capas o los nombres de las variables. Recordemos que los archivos NetCDF también pueden contener varias variables…
# fechas
t <- time(spei)
head(t)[1] "2017-01-01" "2017-01-09" "2017-01-16" "2017-01-23" "2017-02-01"
[6] "2017-02-09"# nombres de capas
names(spei) %>% head()[1] "spei_2689" "spei_2690" "spei_2691" "spei_2692" "spei_2693"
[6] "spei_2694"# nombres de variables
varnames(spei)[1] "spei"18.4 Extracción de series temporales
Una posibilidad que permiten los datos en formato NetCDF es la extracción de series temporales, a partir de puntos o de áreas.
A continuación, crearemos las series temporales del SPEI-12 para la ciudad de Valencia y el promedio de toda la Comunitat Valenciana:
# coordenadas de Valencia
val <- sf::st_point(c(-0.431551, 39.4078969)) %>%
sf::st_sfc(crs = 4326) %>%
sf::st_as_sf() %>%
sf::st_transform(25830)
valSimple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 721138.4 ymin: 4365191 xmax: 721138.4 ymax: 4365191
Projected CRS: ETRS89 / UTM zone 30N
x
1 POINT (721138.4 4365191)El paquete terra sólo acepta su propia clase vectorial: SpatVector. Por eso es necesario convertir el punto de clase sf con la función vect(). Para extraer la serie temporal empleamos la función extract(). Los datos extraídos los encontramos en forma de una tabla, cada fila es un elemento de los datos vectoriales y cada columna una capa. En nuestro caso sólo hay una fila, correspondiente a la ciudad de Valencia.
# extraer la serie temporal
spei_val <- terra::extract(spei, terra::vect(val))
# dimensiones
dim(spei_val)[1] 1 49spei_val[,1:5] ID spei_2689 spei_2690 spei_2691 spei_2692
1 1 0.09 0.12 0.7 0.72# creamos un data.frame con fechas y valores
spei_val <- tibble(date = t, val = unlist(spei_val)[-1])
head(spei_val)# A tibble: 6 × 2
date val
<date> <dbl>
1 2017-01-01 0.0900
2 2017-01-09 0.120
3 2017-01-16 0.700
4 2017-01-23 0.720
5 2017-02-01 0.690
6 2017-02-09 0.760 El promedio de la Comunidad Valenciana lo obtenemos usando la geometría de polígono e indicando el tipo de función con la que queremos resumir el área. Como ya hemos visto en otras ocasiones, la función esp_get_ccaa() del paquete mapSpain es muy útil a la hora de importar límites administrativos españoles, en este caso a nivel autonómico. En la extracción es importante que pasemos el argumento na.rm = TRUE de la función mean() para excluir píxeles sin valor.
# límites de la Comunitat Valenciana
cval <- esp_get_ccaa("Comunitat Valenciana") %>%
st_transform(25830)
# extraemos los valores medios del SPEI-12
spei_cval <- terra::extract(spei, terra::vect(cval), fun = "mean", na.rm = TRUE)
# dimensiones
dim(spei_cval)[1] 1 49spei_cval[,1:5] ID spei_2689 spei_2690 spei_2691 spei_2692
1 1 -0.01527766 -0.02279898 0.5954508 0.6369707# añadimos los nuevos valores a nuestro data.frame
spei_val <- mutate(spei_val, cval = unlist(spei_cval)[-1])En el siguiente paso transformamos la tabla a formato largo (long format) con pivot_longer(), fusionando el valor del índice SPEI de Valencia y de la Comunitat Valenciana. Además añadimos una columna con la interpretación del índice y cambiamos las etiquetas.
spei_val <- tidyr::pivot_longer(data = spei_val,
cols = 2:3,
names_to = "reg",
values_to = "spei") %>%
dplyr::mutate(sign = case_when(spei < -0.5 ~ "sequía",
spei > 0.5 ~ "húmedo",
TRUE ~ "normal"),
date = as_date(date),
reg = factor(reg,
c("val", "cval"),
c("Valencia", "Comunitat Valenciana")))Ahora construimos un gráfico en el que comparamos el SPEI-12 de la ciudad de Valencia con el promedio de la Comunitat Valenciana. La función geom_rect() nos ayuda a dibujar diferentes rectángulos de fondo para marcar la sequía, episodio normal o húmedo.
ggplot(spei_val) +
geom_rect(aes(xmin = min(date), xmax = max(date),
ymin = -0.5, ymax = 0.5),
fill = "#41ab5d") +
geom_rect(aes(xmin = min(date), xmax = max(date),
ymin = -1, ymax = -0.5),
fill = "#ffffcc") +
geom_rect(aes(xmin = min(date), xmax = max(date),
ymin = -1.5, ymax = -1),
fill = "#F3641D") +
geom_hline(yintercept = 0, size = 0.75, color = "#296C3B") +
geom_line(aes(date, spei, linetype = reg), size = 1, alpha = .7) +
scale_x_date(date_breaks = "month", date_labels = "%b") +
scale_y_continuous(limits = c(-1.5, 1.5)) +
labs(linetype = "", y = "SPEI-12", x = "2017") +
coord_cartesian(expand = FALSE) +
theme_minimal() +
theme(legend.position = "top",
legend.direction = "horizontal",
legend.text = element_text(size = 12),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank())
18.5 Mapa de sequía de España
Antes de crear un mapa de la severidad de la sequía en 2017, primero debemos hacer algunas modificaciones en el dataset. Con la función subset() obtenemos una o varias capas como subconjunto. En este caso seleccionamos la última para poder ver el estado de sequía de todo el año.
En el siguiente paso reemplazamos todos los valores mayores de -0,5 con NA (se considera sequía cuando el índice SPEI está por debajo de -0,5).
La clase raster no es directamente compatible con ggplot, por eso, lo convertimos en una tabla xyz con longitud, latitud y la variable de interés. Cuando hacemos la misma conversión de varias capas, cada columna representa una capa.
Finalmente, renombramos la columna del índice y añadimos una nueva variable (columna) con los distintos grados de severidad de sequía.
# extraemos capa(s) con su índice
spei_anual <- subset(spei, 48)
# sustituimos valores de no-sequía con NA
spei_anual[spei_anual > -0.5] <- NA
# convertimos nuestro raster en una tabla de xyz
spei_df <- as.data.frame(spei_anual, xy = TRUE)
head(spei_df) x y spei_2736
38096 123100 4858900 -1.48
39195 105500 4857800 -1.59
39197 107700 4857800 -1.40
39211 123100 4857800 -1.47
39212 124200 4857800 -1.50
40310 105500 4856700 -1.63# cambiamos el nombre de la variable
names(spei_df)[3] <- "spei"
# categorizamos el índice y fijamos el orden del factor
spei_df <- mutate(spei_df, spei_cat = case_when(spei > -0.9 ~ "leve",
spei > -1.5 & spei < -0.9 ~ "moderada",
spei > -2 & spei <= -1.5 ~ "severa",
TRUE ~ "extrema") %>%
forcats::fct_relevel(c("leve", "moderada", "severa", "extrema")))Creamos un mapa raster con la geometría geom_tile() indicando longitud, latitud y el color de los píxeles con nuestra variable categorizada.
ccaa <- esp_get_ccaa() %>%
filter(!ine.ccaa.name %in% c("Canarias", "Ceuta", "Melilla")) %>%
st_transform(25830)
ggplot(spei_df) +
geom_tile(aes(x , y, fill = spei_cat)) +
geom_sf(data = ccaa, fill = NA, size = .1, color = "white", alpha = .4) +
scale_fill_manual(values = c("#ffffcc", "#F3641D", "#DE2929", "#8B1A1A"),
na.value = NA) +
guides(fill = guide_legend(keywidth = 2,
keyheight = .3,
label.position = "bottom",
title.position = "top")) +
coord_sf() +
labs(fill = "SEQUÍA 2017") +
theme_void() +
theme(legend.position = "top",
legend.justification = 0.5,
plot.background = element_rect(fill = "black", color = NA),
legend.title = element_text(color = "white", size = 20, hjust = .5),
legend.text = element_text(color = "white"),
plot.margin = margin(t = 10))
18.6 Mapa de sequía de la Comunitat Valenciana
En este último ejemplo seleccionamos la situación de sequía a 12 meses vista, es decir, al principio y al final del año. La función principal es crop(), que recorta a la extensión de un objeto espacial (en nuestro caso la Comunitat Valenciana). Después aplicamos la función mask(), que enmascara todos aquellos píxeles dentro de los límites dejando en NA los demás.
# subconjunto primera y ultima semana 2017
spei_sub <- subset(spei, c(1, 48))
# recortamos y enmascaramos la Comunitat Valenciana
spei_cval <- terra::crop(spei_sub, cval) %>%
terra::mask(vect(cval))
# convertimos los datos a xyz
spei_df_cval <- as.data.frame(spei_cval, xy = TRUE)
# renombramos las dos capas
names(spei_df_cval)[3:4] <- c("Enero", "Diciembre")
# pasamos al formato de tabla larga fusionando ambos meses
spei_df_cval <- pivot_longer(spei_df_cval, 3:4,
names_to = "mes",
values_to = "spei") %>%
mutate(mes = forcats::fct_relevel(mes, c("Enero", "Diciembre")))Estos dos mapas los hacemos de la misma forma que el de España. La diferencia principal es que usamos el índice SPEI directamente como variable continua. Además, para crear dos mapas con una fila añadimos la función facet_grid(). Por último, el índice muestra valores negativos y positivos, por tanto, es necesario una gama divergente de colores. Con el objetivo de centrar el punto medio en 0 debemos reescalar con ayuda de la función rescale() del paquete scales.
ggplot(spei_df_cval) +
geom_tile(aes(x , y, fill = spei)) +
geom_sf(data = cval, fill = NA, size = .1, color = "white", alpha = .4) +
scale_fill_distiller(palette = "RdYlGn", direction = 1,
values = scales::rescale(c(-2.1, 0, 0.9)),
breaks = seq(-2, 1, .5)) +
guides(fill = guide_colorbar(barwidth = 8, barheight = .3, label.position = "bottom")) +
facet_grid(. ~ mes) +
coord_sf() +
labs(fill = "SPEI-12", title = "Comunitat Valenciana") +
theme_void() +
theme(legend.position = "top",
legend.justification = 0.5,
legend.title = element_text(color = "white", vjust = 1.1),
strip.text = element_text(color = "white"),
plot.background = element_rect(fill = "black", color = NA),
plot.title = element_text(color = "white", size = 20, hjust = .5, vjust = 2.5,
margin = margin(b = 10, t = 10)),
legend.text = element_text(color = "white"),
plot.margin = margin(10, 10, 10, 10))
18.7 Más posibilidades
Es posible agrupar las diferentes capas aplicando una función. Usando los meses de cada semana del SPEI-12 podemos calcular el promedio mensual en 2017. Para ello hacemos uso de la función tapp() que a su vez aplica sobre índices otra función. Es importante que el grupo o bien sea un factor o el índice de cada capa. Las funciones tapp() y app() tienen un argumento para procesar en paralelo usando más de un núcleo (core).
# meses como factor
mo <- lubridate::month(t, label = TRUE)
mo [1] ene ene ene ene feb feb feb feb mar mar mar mar abr abr abr abr
[17] may may may may jun jun jun jun jul jul jul jul ago ago ago ago
[33] sep sep sep sep oct oct oct oct nov nov nov nov dic dic dic dic
12 Levels: ene < feb < mar < abr < may < jun < jul < ... < dic# promedio por mes
spei_mo <- terra::tapp(x = spei,
index = mo,
fun = mean)
spei_moclass : SpatRaster
dimensions : 834, 1115, 12 (nrow, ncol, nlyr)
resolution : 1100, 1100 (x, y)
extent : -80950, 1145550, 3979450, 4896850 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830)
source(s) : memory
names : ene, feb, mar, abr, may, jun, ...
min values : -1.2800, -1.4675, -2.2400, -2.6500, -2.5775, -2.4675, ...
max values : 1.3875, 1.9175, 1.7475, 1.8375, 1.7500, 1.7000, ... plot(spei_mo)
La función mean() directamente usada sobre un objeto de clase SpatRaster multidimensional devuelve el promedio por celda. El mismo resultado lo podemos obtener con la función app() que aplica cualquier función.
El número de capas resultante depende de la función, por ejemplo, al aplicar range() el resultado son dos capas, una del valor mínimo y otra del máximo.
Por último, la función global() resume con la función indicada cada capa en forma de una tabla.
# promedio sobre capas
spei_mean <- terra::mean(spei)
spei_meanclass : SpatRaster
dimensions : 834, 1115, 1 (nrow, ncol, nlyr)
resolution : 1100, 1100 (x, y)
extent : -80950, 1145550, 3979450, 4896850 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830)
source(s) : memory
name : mean
min value : -2.127083
max value : 1.568542 plot(spei_mean)
# alternativa
spei_max <- terra::app(x = spei, fun = max)
spei_maxclass : SpatRaster
dimensions : 834, 1115, 1 (nrow, ncol, nlyr)
resolution : 1100, 1100 (x, y)
extent : -80950, 1145550, 3979450, 4896850 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830)
source(s) : memory
name : max
min value : -1.06
max value : 2.02 spei_min <- terra::app(x = spei, fun = min)
spei_minclass : SpatRaster
dimensions : 834, 1115, 1 (nrow, ncol, nlyr)
resolution : 1100, 1100 (x, y)
extent : -80950, 1145550, 3979450, 4896850 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830)
source(s) : memory
name : min
min value : -3.33
max value : 0.29 spei_range <- terra::app(x = spei, fun = range)
names(spei_range) <- c("min", "max")
spei_rangeclass : SpatRaster
dimensions : 834, 1115, 2 (nrow, ncol, nlyr)
resolution : 1100, 1100 (x, y)
extent : -80950, 1145550, 3979450, 4896850 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830)
source(s) : memory
names : min, max
min values : -3.33, -1.06
max values : 0.29, 2.02 plot(spei_range)
# resumen estadístico por capa
terra::global(x = spei,
fun = "mean",
na.rm = TRUE) %>%
head() mean
spei_2689 -0.03389126
spei_2690 -0.17395742
spei_2691 -0.13228593
spei_2692 -0.07536089
spei_2693 0.06718260
spei_2694 -0.03461822