Computando tendencias en una rejilla climática con Tidyverse

Jul 9, 2020 6 min de lectura R, Tidyverse

Lectura de un NetCDF al estilo Tidyverse!

Actualmente es posible leer un fichero NetCDF al estilo Tidyverse mediante el package creado por Elio Campitelli llamado metR. Además esta librería contiene muchisimas funciones útiles para la representación y el análisis de datos meteorológicos y climáticos.

En esta ocasión trabajaremos con una rejilla diaria de 5x5 km de temperatura máxima de las Islas Baleares (España) (Serrano-Notivoli 2019). Para la lectura usaremos un par de funciones del package metR: GlanceNetCDF, para echarle un vistazo al nombre de la variable, dimensiones, etc, y ReadNetCDF, para la lectura del fichero NetCDF. Lógicamente haremos uso de la librería tidyverse, la cual nos permitirá organizar los datos y realizar algun cálculo sencillo (recuerda que tidyverse es una colección de packages!), así como las librerias lubridate, para el manejo de fechas, y la librería trend para calcular la significación estadística de la tendencia de nuestros datos. Además usaremos la librería pals para escoger la paleta de colores cuando representemos los resultados gráficamente.

library(tidyverse)
library(lubridate)
library(metR)
library(trend)
library(pals) # paletas de colores

Como queremos trabjar los datos del NetCDF en formato Tidyverse, vamos a especificar que el output de la función ReadNetCDF sea un data.frame:

# Descargamos el fichero de Temperatura máxima (Tx) de Baleares 

# https://digital.csic.es/handle/10261/177655

# Extraemos información del NetCDF (nombre variable, dimensiones, etc)
file = "tmax_bal.nc"

variable <- GlanceNetCDF(file = file)
variable
## ----- Variables ----- 
## tx:
##     Maximum temperature in Celsius degrees
##     Dimensions: lon by lat by Time
## 
## 
## ----- Dimensions ----- 
##   lon: 63 values from 0.970530815018679 to 4.55413081501868 degrees_east
##   lat: 48 values from 38.3196292450381 to 40.4252292450381 degrees_north
##   Time: 16071 values from 1971-01-02 to 2015-01-01

# Leemos el NetCDF
tx <- ReadNetCDF(file = file,
                 vars = c(var = names(variable$vars[1])),# extracción del nombre de la variable.
                 out = "data.frame") # objeto de salida
tx # voilà.
##                 Time      lat       lon var
##        1: 1971-01-02 40.42523 0.9705308  NA
##        2: 1971-01-02 40.42523 1.0283308  NA
##        3: 1971-01-02 40.42523 1.0861308  NA
##        4: 1971-01-02 40.42523 1.1439308  NA
##        5: 1971-01-02 40.42523 1.2017308  NA
##       ---                                  
## 48598700: 2015-01-01 38.31963 4.3229308  NA
## 48598701: 2015-01-01 38.31963 4.3807308  NA
## 48598702: 2015-01-01 38.31963 4.4385308  NA
## 48598703: 2015-01-01 38.31963 4.4963308  NA
## 48598704: 2015-01-01 38.31963 4.5541308  NA

Gestión y manejo de los datos

A continuación el objetivo es agregar los datos a escala anual, eso es, pasar del dato diario al promedio anual. Para ello tenemos que jugar con varias funciones de la librería dplyr, cargada en el tidyverse. Primero, sin embargo, hay que filtrar los valores NA:

tx <- tx %>% filter(!is.na(var)) %>% as.tibble()
## Warning: `as.tibble()` is deprecated as of tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
tx
## # A tibble: 4,114,176 x 4
##    Time                  lat   lon   var
##    <dttm>              <dbl> <dbl> <dbl>
##  1 1971-01-02 00:00:00  40.1  3.92  7.64
##  2 1971-01-02 00:00:00  40.1  3.98  7.69
##  3 1971-01-02 00:00:00  40.1  4.03  7.44
##  4 1971-01-02 00:00:00  40.1  4.09  8.29
##  5 1971-01-02 00:00:00  40.1  4.15  8.15
##  6 1971-01-02 00:00:00  40.1  4.21  8.35
##  7 1971-01-02 00:00:00  40.0  3.80  8.10
##  8 1971-01-02 00:00:00  40.0  3.86  8.04
##  9 1971-01-02 00:00:00  40.0  3.92  7.63
## 10 1971-01-02 00:00:00  40.0  3.98  7.70
## # ... with 4,114,166 more rows

Una vez filtrados los NA, procedemos a realizar el agregado anual en cada punto de grid:

tx_annual <- tx %>%
  mutate(Time = as_date(Time)) %>% #convertimos las fechas en formato lubridate
  mutate(year = year(Time)) %>% # Creamos la columna de años para realizar agregado
  group_by(year, lon, lat) %>% # agrupamos por año, longitud y latitud
  summarise(annual_tx = mean(var)) %>% # Ahora solo falta realizar el promedio
  ungroup() # una vez terminamos, desagrupamos
## `summarise()` regrouping output by 'year', 'lon' (override with `.groups` argument)

tx_annual
## # A tibble: 11,520 x 4
##     year   lon   lat annual_tx
##    <dbl> <dbl> <dbl>     <dbl>
##  1  1971  1.20  38.9      20.8
##  2  1971  1.20  38.9      20.9
##  3  1971  1.26  38.9      20.2
##  4  1971  1.26  38.9      20.3
##  5  1971  1.26  38.9      20.8
##  6  1971  1.32  38.9      20.7
##  7  1971  1.32  38.9      20.6
##  8  1971  1.32  38.9      20.7
##  9  1971  1.32  39.0      21.0
## 10  1971  1.32  39.0      20.3
## # ... with 11,510 more rows

Como podemos observar, se muestra la temperatura máxima promedio para cada año y punto de grid.

Cálculo de la tendencia y de su significación estadística

De este modo, ya podemos calcular la tendencia anual y su significación estadística con la librería trend. Para ello emplearemos la función sens.slope para el cálculo de la pendiente de la tendencia (qué magnitud de temperatura aumenta/disminuye en un año); y la función mk.test (Test de Mann Kendall), para verificar si la tendencia es estadísticamente significativa o si realmente es fruto de la variabilidad de la propia serie termométrica. Además, al tratarse de un test no paramétrico, nuestros datos no tendrán por que responder a una distribución gausiana.

tx_trend <- tx_annual %>%
  group_by(lon,lat) %>% # Agrupamos por lon y lat para realizar el cálculo en cada celda
  summarise(slope = sens.slope(annual_tx)$estimates *10,
            sign = mk.test(annual_tx)$p.value)
## `summarise()` regrouping output by 'lon' (override with `.groups` argument)

tx_trend
## # A tibble: 256 x 4
## # Groups:   lon [38]
##      lon   lat slope   sign
##    <dbl> <dbl> <dbl>  <dbl>
##  1  1.20  38.9 0.122 0.156 
##  2  1.20  38.9 0.122 0.156 
##  3  1.26  38.9 0.110 0.115 
##  4  1.26  38.9 0.112 0.120 
##  5  1.26  38.9 0.147 0.0834
##  6  1.32  38.9 0.136 0.0983
##  7  1.32  38.9 0.133 0.111 
##  8  1.32  38.9 0.146 0.0869
##  9  1.32  39.0 0.150 0.0645
## 10  1.32  39.0 0.133 0.0983
## # ... with 246 more rows

Ya tenemos el valor del pendiente Sen, que nos indica la tasa de cambio cada 10 años (hemos multiplicado x 10), así como la significación de la tendencia indicada a través del p-value. Ahora simplemente realizaremos la visualización de los resultados indicando aquellas celdas que muestran una significación estadística al aumento o disminución de la temperatura máxima anual.

ggplot()+
  geom_tile(data = tx_trend, aes(x = lon, y = lat, fill = slope))+
  scale_fill_gradientn(colors = rev(pals::linearlhot(100)), name = "ºC/10y", limits = c(0.1,0.45)) +
  geom_point(data = filter(tx_trend, sign < 0.01),aes(x = lon, y = lat, color = "Tendencia sign. \n p-value <0.01"),
             size = 0.4, show.legend = T) +
  scale_color_manual(values = c("black"), name = "")+
  coord_fixed(1.3)+
  xlab("Longitud") + ylab("Latitud")+
  labs(title = "Tendencia decadal de la temperatura máxima en las Islas Baleares", 
       subtitle = "(1971-2014)",
       caption = "Base de datos: STEAD (Serrano-Notivoli,2019)")+
  theme_bw() +
  guides(fill = guide_colourbar(barwidth = 9, barheight = 0.5, title.position="right"))+
  theme(legend.position = "bottom")