En este tutorial trabajaremos con datos climáticos descargados de la pagina web de CRU (Climatic Researc Unit), se realizará la descarga, y tratamiento de los datos espaciales con el fin de realizar gráficas de tendencia a lo largo de 35 años (1980-2015); todo lo anterior lo realizaremos en el lenguaje de programación R. Las variables disponibles para descarga son las siguientes:

– cld: cobertura de nubes (%)
– dtr: rango temperatura diurno (grados centigrados)
– fdf: frost day frequency (días)
– pet: evapotranspiración potencial (mm/day)
– pre: precipitation (mm por mes)
– rhm: humedad relativa (porcentaje)
– ssh: sunshine duración (hours)
– tmp: temperatura promedio (grados centigrados)
– tmn: promedio de temperatura minima (grados centigrados)
– tmx: promedio de temperatura máxima (grados centigrados)
– presión de vapor (hectopascal hPa)
– wet: frecuencia de días húmedos (días)
– wnd: velocidad del viento (metros por segundo)

En este ejercicio usaremos las variables de precipitación, temperatura máxima, promedio y mínima. Para la realización de este ejercicio utilizaremos información que descargaremos durante el mismo tutorial. Se recomienda tener un nivel básico de programación para poder entender este tutorial, sin embargo, se intentará explicar las lineas a utilizar.

A continuación, listaremos los pasos a realizar explicando cada uno de ellos.

1. Cómo primer paso esta la instalación de los paquetes a utilizar en el sitio web.

# Load libraries
if(!require(raster)){install.packages('raster'); library(raster)} else {library(raster)}
if(!require(raster)){install.packages('rgdal'); library(raster)} else {library(rgdal)}
if(!require(tidyverse)){install.packages('tidyverse'); library(tidyverse)} else {library(tidyverse)}
if(!require(utils)){install.packages('utils'); library(utils)} else {library(utils)}

if(!require(gtools)){install.packages('gtools'); library(gtools)} else {library(gtools)}
if(!require(velox)){install.packages('velox'); library(velox)} else {library(velox)}
if(!require(stringr)){install.packages('stringr'); library(stringr)} else {library(stringr)}

2. Luego procedemos a configurar las variables y los años a descargar, para ello creamos los objetos vars, yr.start y yr.end; además establecemos la ruta base.

path_base <- 'https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.01/cruts.1709081022.v4.01/'
vars <- c('pre', 'tmn', 'tmx', 'tmp')
yr.start <-  c(1981, 1991, 2001, 2011)
yr.end <- c(1990, 2000, 2010, 2016)

3. Ahora procedemos a crear dos ciclos, haciendo uso del comando lapply, para descargar tanto las 4 variables como los distintos periodos de tiempo, de manera automática.

lapply(1:length(vars), function(v){
  print(vars[v])
  lapply(1:length(yr.start), function(x){
  print(yr.start[x])
  download.file(url = paste0(path_base, vars[v], '/cru_ts4.01.', yr.start[x], '.', yr.end[x], '.', vars[v], '.dat.nc.gz'),
                destfile = paste0('_raster/_gz/', vars[v], '_', yr.start[x], '_', yr.end[x], '.nc.gz'),
                mode = 'wb')
  })
})

4. Como siguiente paso procedemos a descomprimir los archivo .gz obtenidos, para ello hacemos uso de las siguientes lineas. En el parametro destname establecemos la ruta donde queremos guardar nuestro conjunto de datos.

Todos los archivos resultantes los guardaremos en la carpeta nc (parámetro destname de la línea en la que usamos el comando gunzip).

gzfiles <- list.files('../_download/_raster/_gz', full.names = T, pattern = '.gz$')
nms <- basename(gzfiles) %>% gsub('.gz', '', .)
lapply(1:length(gzfiles), function(x){
  print(gzfiles[x])
  R.utils::gunzip(filename = gzfiles[x], destname = paste0('../_download/_raster/_nc/', nms[x]))
})

5. Ahora procedemos a listar los objetos que hemos descomprimido y luego las decadas de dichos archivos, estos elementos serán los parámetros a utilizar en la función que hará la extracción por máscara para Colombia de los archiovs, y luego la escritura de los mismos en una carpeta de nuestro espacio de trabajo.

ncfiles <- list.files('../_download/_raster/_nc', full.names = T, pattern = '.nc$')
decs <- basename(ncfiles) %>% str_sub(., start = 5, end = 13) %>% mixedsort() %>% unique()
col <- getData('GADM', country = 'COL', level = 0)
extNC <- function(vr, dec){
  ncfile <- grep(vr, ncfiles, value = T) %>% grep(dec, ., value = T)
  stk <- brick(ncfile)
  stk.col <- raster::crop(stk, col) %>% raster::mask(., col)
  nms <- names(stk.col) %>% gsub('\\.', '_', .) %>% str_sub(string = ., start = 2, end = 8)
  lyrs <- unstack(stk.col)
  Map('writeRaster', x = lyrs, filename = paste0('../_download/_raster/_tif/', vr, '_', nms, '.tif'))
  print('Done!')
} 

Aquí la aplicación de la función anteriormente creada. Este ciclo intera primero en la variable y luego en la decena, así primero extraera todos los períodos de tiempo y luego pasará a la siguiente variable, hasta completar las cuatro variables descargadas.

lapply(1:length(vars), function(v){
  print(vars[v])
  lapply(1:length(decs), function(d){
    print(decs[d])
    extNC(vr = vars[v], dec = decs[d])
  })
})

6. En la función anterior logramos obtener, entonces, los archivos de precipitación, temperatura máxima, media y mínima para Colombia. Ahora plotearemos un archivo de ellos, y después convertiremos los archivos raster en tablas para de esta manera hacer gráficas de tendencia para las variables mencionadas anteriormente.

Lo primero que haremos en este paso será listar los archivos resultantes del pasos anterior, y luego a partir de un for crearemos 4 stacks, uno para cada variable, tal como se muestra a continuación:

fls <- list.files('../_download/_raster/_tif', full.names = T, pattern = '.tif$') %>% mixedsort()
for(i in 1:length(vars)){
  eval(parse(text = paste0(vars[i], '<- grep(vars[', i, '], fls, value = T) %>% stack()')))
}

plot(pre[[1]])

7. Lo siguiente a realizar corresponde a la conversión de los archivos raster a tabla, para poder así crear las gráficas de tendencias. Aquí crearemos la función que realiza dicha conversión.

Lo primero que se realiza en la función es la conversión del stack a un objeto velox, esto pues es mucho más rapido las funciones de dicha libreria, luego extraemos las coordenadas de los raster, y luego con SpatialPoints convertimos las coordeanas en un shapefile, el siguiente paso es la extracción de los valores para todos y cada uno de los raster, luego y quizá estas son las lineas más complejas de entender, es la conversión de la matrix resultante en un data frame, aquí se realiza la asignación de un ID para cada píxel, así como también la creación de columnas que representan el año y el mes.

rstToTbl <- function(rst){
  print(rst)
  #rst <- pre
  vlx <- velox(rst)
  crd <- vlx$getCoordinates()
  spn <- SpatialPoints(coords = crd)
  vls <- vlx$extract_points(sp = spn)
  rsl <- vls %>% 
    tbl_df() %>%
    .[complete.cases(.),] %>%
    mutate(ID = 1:nrow(.)) %>% 
    setNames(c(names(rst), 'ID')) %>%
    gather(var, value, -ID) %>%
    mutate(year = str_sub(var, 5, 8),
           mnth = str_sub(var, 10, 11))
  print('Done')
  return(rsl)
}

Ahora aplicamos la función a las variables, para ello primero creamos una lista a partir de los objetos: pre, tmx, tmp, y tmn. Y seguido de ello con el ciclco del lapply aplicamos la función a la lista, el resultado será el objeto fnl (de final) que contendrá los cuatro tablas con la información climática de Colombia.

stks <- list(pre, tmx, tmp, tmn)
fnl  <- lapply(1:length(stks), function(k) rstToTbl(rst = stks[[k]]))

8. Ahora haremos uso de funciones de la libreria ggplot para crear una gráfica para un píxel en particular. Sí Ud. lo desea puede crear un ciclo (lapply o for, por ejemplo) para crear una gráfica para cada píxel.

En este paso lo primero que haremos será un resumen anual de la precipitación, es decir, haremos la suma de la precipitación para los 12 meses de cada año, esto para un píxel en particular (ID = 300); luego crearemos el gráfico y estimaremos una gráfica de tendencia, vale aclarar que existen muchas maneras de realizarlo, este es solamente uno de ellas.

x <- fnl[[1]] %>%
  filter(ID == 300)
x <- x %>% 
  group_by(year) %>%
  summarize(sum = sum(value)) %>%
  mutate(year = as.numeric(year))
gg <- ggplot(data = x, aes(x = year, y = sum)) +
  geom_line() +
  geom_smooth(method = 'lm') + 
  xlab('year') +
  ylab('Prec (mm)') +
  theme_bw()
ggsave(plot = gg, filename = 'prec_line.png', units = 'cm', width = 15, height = 13, dpi = 150)  

El resultado debería ser algo como esto:

Sí has llegado hasta gracias por leer el tutorial, y por favor, no olvides compartirlo para que así este sea difundido de la mejor manera.