library(raster)
## Warning: package 'raster' was built under R version 4.3.1
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.3.1
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
library(RNetCDF)
library(ncdf4)
library(sp)
library(sf)
## Warning: package 'sf' was built under R version 4.3.1
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(R.matlab)
## R.matlab v3.7.0 (2022-08-25 21:52:34 UTC) successfully loaded. See ?R.matlab for help.
## 
## Attaching package: 'R.matlab'
## The following objects are masked from 'package:base':
## 
##     getOption, isOpen
library(easyNCDF)
library(raster)
library(rgdal)
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-6, (SVN revision 1201)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: C:/Users/javel/AppData/Local/R/win-library/4.3/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: C:/Users/javel/AppData/Local/R/win-library/4.3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(plotKML)
## Warning: package 'plotKML' was built under R version 4.3.1
## plotKML version 0.8-3 (2022-06-05)
## URL: https://github.com/Envirometrix/plotKML/
library(lattice)
library(rasterVis)

En este ejercicio veremos cómo se implementa el método delta para generar escenarios futuros de temperatura y precipitación. La metodología fue explicada en las primeras clases por el Mtro. Oscar Calderón Bustamante. También está descrita en Navarro-Racines et al. (2020; https://www.nature.com/articles/s41597-019-0343-8).

# cargar un template de raster
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  # geographical, datum WGS84
r <- raster("bio_01.tif")
projection(r)<-crs.geo

Primero vamos a hacer la extracción de datos históricos y del periodo de referencia y posteriormente el cálculo de los deltas. Se cargan los datos de temperatura mínima, máxima y precipitación mensual de un ensamble de modelos (192 realizaciones) para el SSP585. Los datos los pueden bajar de Cimate Explorer (https://climexp.knmi.nl/start.cgi). La temperatura deben convertirla a grados celsius y los datos de precipitación a unidades de mm/mes.

# temperatura mínima
tasmin.f <- brick("tasmin_mon_ens_ssp585_192_ave.nc")
tasmin.f <- tasmin.f-273.15

# Extraer los datos del periodo de referencia climatológico
tasmin.hist <- tasmin.f[[1441:1812]]
t2 <- split(seq_len(372), rep(1:12, length.out=372))
tasmin.hist.avg <- lapply(1:12, function(x) tasmin.hist[[t2[[x]]]])
tasmin.hist.avg <- lapply(1:12, function(x) mean(tasmin.hist.avg[[x]]))
tasmin.hist.avg <- stack(tasmin.hist.avg)
meses <- seq(as.Date("2000-01-01"), as.Date("2000-12-01"), by = "months")
names(tasmin.hist.avg) <- meses
p1 <- tasmin.f[[1813:3012]]
delta_tasmin_ens_ssp585_192 <- p1 - tasmin.hist.avg
dates <- seq(as.Date("2001-01-01"), as.Date("2100-12-01"), by = "months")
names(delta_tasmin_ens_ssp585_192) <- dates
plot(delta_tasmin_ens_ssp585_192[[1000]])

# temperatura máxima mensual
tasmax.f <- brick("tasmax_mon_ens_ssp585_192_ave.nc")
tasmax.f <- tasmax.f-274.15

tasmax.hist <- tasmax.f[[1441:1812]]
t2 <- split(seq_len(372), rep(1:12, length.out=372))
tasmax.hist.avg <- lapply(1:12, function(x) tasmax.hist[[t2[[x]]]])
tasmax.hist.avg <- lapply(1:12, function(x) mean(tasmax.hist.avg[[x]]))
tasmax.hist.avg <- stack(tasmax.hist.avg)
meses <- seq(as.Date("2000-01-01"), as.Date("2000-12-01"), by = "months")
names(tasmax.hist.avg) <- meses
p1 <- tasmax.f[[1813:3012]]
delta_tasmax_ens_ssp585_192 <- p1 - tasmax.hist.avg
dates <- seq(as.Date("2001-01-01"), as.Date("2100-12-01"), by = "months")
names(delta_tasmax_ens_ssp585_192) <- dates
plot(delta_tasmax_ens_ssp585_192[[1000]])

# precipitación
pr.f <- brick("pr_mon_ens_ssp585_192_ave.nc")
pr.f <- pr.f*(30*24*60*60)
pr.hist <- pr.f[[1441:1812]]
t2 <- split(seq_len(372), rep(1:12, length.out=372))
pr.hist.avg <- lapply(1:12, function(x) pr.hist[[t2[[x]]]])
pr.hist.avg <- lapply(1:12, function(x) mean(pr.hist.avg[[x]]))
pr.hist.avg <- stack(pr.hist.avg)
meses <- seq(as.Date("2000-01-01"), as.Date("2000-12-01"), by = "months")
names(pr.hist.avg) <- meses
p1 <- pr.f[[1813:3012]]
delta_pr_ens_ssp585_192 <- (p1 - pr.hist.avg)/pr.hist.avg
dates <- seq(as.Date("2001-01-01"), as.Date("2100-12-01"), by = "months")
names(delta_pr_ens_ssp585_192) <- dates
plot(delta_pr_ens_ssp585_192[[100]])

Rotar los rasters

delta_tasmax_ens_ssp585_192 <- terra::rotate(delta_tasmax_ens_ssp585_192, left=TRUE)
delta_tasmin_ens_ssp585_192 <- terra::rotate(delta_tasmin_ens_ssp585_192, left=TRUE)
delta_pr_ens_ssp585_192 <- terra::rotate(delta_pr_ens_ssp585_192, left=TRUE)

Cálculo de Delta para temperatura promedio anual

p <- lapply(1:1200, function(x) mean(delta_tasmax_ens_ssp585_192[[x]], delta_tasmin_ens_ssp585_192[[x]]))
tasmean_ssp585 <- stack(p)
dates <- seq(as.Date("2001-01-01"), as.Date("2100-12-01"), by = "months")
names(tasmean_ssp585) <- dates
t1 <- split(seq_len(1200), rep(1:100, each=12))
t2 <- lapply(1:100, function(x) tasmean_ssp585[[t1[[x]]]])
delta_bio1_ssp585 <- lapply(1:100, function(x) mean(t2[[x]]))
delta_bio1_ssp585 <- stack(delta_bio1_ssp585)
delta_bio1_ssp585
## class      : RasterStack 
## dimensions : 144, 192, 27648, 100  (nrow, ncol, ncell, nlayers)
## resolution : 1.875, 1.25  (x, y)
## extent     : -179.0625, 180.9375, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :      layer.1,      layer.2,      layer.3,      layer.4,      layer.5,      layer.6,      layer.7,      layer.8,      layer.9,     layer.10,     layer.11,     layer.12,     layer.13,     layer.14,     layer.15, ... 
## min values :  0.028977258,  0.043724390,  0.027354700,  0.035543544, -0.003240774, -0.044730862,  0.036270548,  0.071346314,  0.038030026,  0.057134025,  0.120097659,  0.148231511,  0.122626428,  0.123610917,  0.157369010, ... 
## max values :     2.262512,     2.544159,     2.898864,     3.334587,     3.480233,     3.877360,     4.174153,     4.499359,     4.450591,     4.582206,     4.808756,     4.922228,     5.119691,     5.327759,     5.151492, ...

Cálculo de Delta para precipitación anual

t1 <- split(seq_len(1200), rep(1:100, each=12))
t2 <- lapply(1:100, function(x) delta_pr_ens_ssp585_192[[t1[[x]]]])
delta_bio12_ssp585 <- lapply(1:100, function(x) mean(t2[[x]]))
delta_bio12_ssp585 <- stack(delta_bio12_ssp585)

Exploramos los deltas para bio1 y bio12

delta_bio1_ssp585
## class      : RasterStack 
## dimensions : 144, 192, 27648, 100  (nrow, ncol, ncell, nlayers)
## resolution : 1.875, 1.25  (x, y)
## extent     : -179.0625, 180.9375, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :      layer.1,      layer.2,      layer.3,      layer.4,      layer.5,      layer.6,      layer.7,      layer.8,      layer.9,     layer.10,     layer.11,     layer.12,     layer.13,     layer.14,     layer.15, ... 
## min values :  0.028977258,  0.043724390,  0.027354700,  0.035543544, -0.003240774, -0.044730862,  0.036270548,  0.071346314,  0.038030026,  0.057134025,  0.120097659,  0.148231511,  0.122626428,  0.123610917,  0.157369010, ... 
## max values :     2.262512,     2.544159,     2.898864,     3.334587,     3.480233,     3.877360,     4.174153,     4.499359,     4.450591,     4.582206,     4.808756,     4.922228,     5.119691,     5.327759,     5.151492, ...
delta_bio12_ssp585
## class      : RasterStack 
## dimensions : 144, 192, 27648, 100  (nrow, ncol, ncell, nlayers)
## resolution : 1.875, 1.25  (x, y)
## extent     : -179.0625, 180.9375, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ... 
## min values : -0.1287373, -0.1750321, -0.1327746, -0.1325774, -0.1300850, -0.1017301, -0.1469527, -0.1140420, -0.1097782, -0.1241216, -0.1367964, -0.1681148, -0.1452880, -0.1302457, -0.1374241, ... 
## max values :  0.4218842,  0.4902218,  0.3677308,  0.4336316,  0.3922701,  0.4383883,  0.5456928,  0.6159500,  0.3678838,  0.4691615,  0.4231700,  0.4288572,  0.3792099,  0.4142094,  0.5729208, ...

Ahora necesitamos ajustar los datos a la resolución espacial de la climatología actual. En este caso tenemos un raster de ~1°x1° de temperatura promedio anual (bio1) generado a partir de datos del WorldClim.

delta_bio1_ssp585.r <- projectRaster(from=delta_bio1_ssp585, to=r, crs=crs(r), method="bilinear")
delta_bio1_ssp585.r
## class      : RasterBrick 
## dimensions : 180, 360, 64800, 100  (nrow, ncol, ncell, nlayers)
## resolution : 1, 1  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : r_tmp_2023-10-30_183147.682356_25412_16904.grd 
## names      :       layer.1,       layer.2,       layer.3,       layer.4,       layer.5,       layer.6,       layer.7,       layer.8,       layer.9,      layer.10,      layer.11,      layer.12,      layer.13,      layer.14,      layer.15, ... 
## min values :  0.0313846005,  0.0461961208,  0.0279455474,  0.0361689325,  0.0001373057, -0.0425215474,  0.0377466023,  0.0756863605,  0.0399598896,  0.0618616620,  0.1247502250,  0.1487914975,  0.1241111506,  0.1251624320,  0.1626065266, ... 
## max values :      2.216503,      2.507488,      2.873116,      3.327629,      3.446888,      3.816464,      4.112271,      4.479281,      4.439846,      4.559672,      4.774508,      4.894906,      5.078867,      5.272557,      5.108799, ...
plot(delta_bio1_ssp585.r[[1]])

delta_bio12_ssp585.r <- projectRaster(from=delta_bio12_ssp585, to=r, crs=crs(r), method="bilinear")
delta_bio12_ssp585.r
## class      : RasterBrick 
## dimensions : 180, 360, 64800, 100  (nrow, ncol, ncell, nlayers)
## resolution : 1, 1  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : r_tmp_2023-10-30_183150.091265_25412_69981.grd 
## names      :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ... 
## min values : -0.1263510, -0.1641282, -0.1273414, -0.1309281, -0.1245210, -0.1000774, -0.1418430, -0.1077873, -0.1058037, -0.1227938, -0.1276665, -0.1671196, -0.1440533, -0.1275869, -0.1353514, ... 
## max values :  0.4016378,  0.4633300,  0.3635910,  0.4238174,  0.3786027,  0.3883792,  0.5009215,  0.5216511,  0.3428648,  0.4615619,  0.3747924,  0.3827242,  0.3782315,  0.3882016,  0.5370588, ...
plot(delta_bio12_ssp585.r[[1]])

Ahora vamos a agregar los deltas a la climatología de WorldClim para crear escenarios futuros de temperatura promedio anual (bio1) y precipitación anual (bio12).

Primero cargamos las dos variables de worldclim

bio1 <- raster("bio_01.tif")
bio12 <- raster("bio_12.tif")

Luego adicionamos los deltas a la climatología.

bio1_ssp585 <- bio1+delta_bio1_ssp585.r
bio12_ssp585 <- (1+delta_bio12_ssp585.r)*bio12

Vamos a graficar para el año 2100

plot(bio1_ssp585[[100]])

plot(bio12_ssp585[[100]])

Tarea: ¿Cómo generarían los escenarios futuros pero para los datos mensuales de temperatura y precipitación?

Pista: Con lo aprendido en las sesiones anteriores fácilmente lo podrían implementar.

Elaborado por Julián A. Velasco Material para el curso “Herramientas para el análisis de impactos de cambio climático en sistemas naturales y humanos”.