Computant tendències en una graella climàtica amb Tidyverse

Lectura d’un NetCDF a l’estil Tidyverse!

Actualment és possible llegir un fitxer NetCDF a l’estil Tidyverse mitjançant el package creat per Elio Campitelli anomenat metR. A més, aquesta llibreria conté moltíssimes funcions útils per a la representació i l’anàlisi de dades meteorològiques i climàtiques. En aquesta ocasió treballarem amb una graella diària de 5x5 km de temperatura màxima de les Illes Balears (Serrano-Notivoli 2019). Per a la seva lectura farem ús d’un parell de funcions del package metR:GlanceNetCDF, ens permetrà fer una ullada al nom de la variable, dimensions, etc, i ReadNetCDF, realitzarà la lectura de l’arxiu NetCDF. Lògicament farem ús de la llibreria tidyverse, la qual ens permetrà organitzar les dades i realitzar algun càlcul senzill (recorda que el tidyverse és una col·lecció de packages!), Així com les llibreries lubridate, per a la gestió de dates, i la llibreria trend per calcular la significació estadística de la tendència. A més farem servir la llibreria pals per escollir la paleta de colors quan representem els resultats gràficament.

library(tidyverse)
library(lubridate)
library(metR)
library(trend)
library(pals) # paletes de colors

Al voler treballar amb les dades del NetCDF en format Tidyverse, especificarem que el output de la funció ReadNetCDF sigui un data.frame:

# Descarreguem el fitxer de Temperatura màxima (Tx) de les Balears 

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

# Extraiem informació de l'NetCDF (nom variable, dimensions, 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

# Llegim el NetCDF
tx <- ReadNetCDF(file = file,
                 vars = c(var = names(variable$vars[1])),# extracció del nom de la variable.
                 out = "data.frame") # objecte de sortida
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ó y manipulació de les dades

A continuació l’objectiu és agregar les dades a escala anual, això és, passar de la dada diària a la mitjana anual. Per a tal objectiu hem de jugar amb diverses funcions de la llibreria dplyr, carregada altidyverse. Primer, però, cal filtrar els valors 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

Un cop filtrats els NA, procedim a calcular les mitjanes anuals a cada punt de grid:

tx_annual <- tx %>%
  mutate(Time = as_date(Time)) %>% # convertim les dates en format lubridate
  mutate(year = year(Time)) %>% # Creem la columna d'anys per realitzar l'agregat
  group_by(year, lon, lat) %>% # agrupem per any, longitud i latitud
  summarise(annual_tx = mean(var)) %>% # Ara només falta fer la mitjana
  ungroup() # desagrupem
## `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

Com podem observar, es mostra la temperatura màxima mitjana per a cada any i punt de grid.

Càlcul de la tendència i de la seva significació estadística

D’aquesta manera, ja podem calcular la tendència anual i la seva significació estadística mitjançant la llibreria trend. Emprarem la funció sens.slope per al càlcul del pendent de la tendència (quina magnitud de temperatura augmenta / disminueix en un any); i la funció mk.test (Test de Mann Kendall), per verificar si la tendència és estadísticament significativa o si realment és fruit de la variabilitat de la pròpia sèrie termomètrica. A més, al tractar-se d’un test no paramètric, les nostres dades no hauran d’ajustar-se a una distribució gaussiana necessàriament.

tx_trend <- tx_annual %>%
  group_by(lon,lat) %>% # Agrupem per lon i lat per a realitzar el càlcul en cada cel·la
  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

Ja tenim el valor del pendent de Sen, que ens indica la taxa de canvi cada 10 anys (hem multiplicat x 10), així com la significació de la tendència indicada a través del p-value. Ara simplement realitzarem la visualització dels resultats indicant aquelles cel·les que mostren una significació estadística a l’augment o disminució 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 = "Tendència 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 = "Tendència decadal de la temperatura màxima a les Illes Balears", 
       subtitle = "(1971-2014)",
       caption = "Base de dades: STEAD (Serrano-Notivoli,2019)")+
  theme_bw() +
  guides(fill = guide_colourbar(barwidth = 9, barheight = 0.5, title.position="right"))+
  theme(legend.position = "bottom")

Per a qualsevol qüestió, no dubtis en deixar un comentari.

Salutacions!