In EvolEcolGroup/pastclim: Manipulate Time Series of Climate Reconstructions

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

THIS IS WORK IN PROGRESS, DO NOT USE FOR ANY ANALYSIS

Note that we will use the draft functions for downscaling, which can be accessed by using the prefix pastclim:::. The ::: accesses internal functions. You will also get a message reminding you that these functions are not final and thus should not be used for real analysis.

Delta downscaling a dataset in pastclim

The terra package can downscale reconstructions using bilinear interpolation. Other R packages also offer various approaches to downscale rasters.

For palaeoclimate reconstructions, the delta method has been shown to be very effective (Beyer et al, REF), but sea level changes require some care on how to apply such a method. pastclim includes functions to use the delta method for downscaling. We will focus on South East Asia, as this area was greatly affected by sea level changes over the last glaciation. Specifically, we will look at the area around Borneo (for real applications, we would recommend using a bigger extent in areas of large changes in land extent, as interpolating over a small extent can lead to greater artefacts; for this example, we keep the extent small to reduce computational time).

An example for one variable

When downscaling, it is generally a better idea to downscale the monthly variables, and then recompute the BIOCLIM variables. We will use the monthly variables from the Beyer2020 dataset, and focus on the time steps that are also available in the Example dataset (i.e. a subset of the full Beyer2020 dataset, to reduce computational time):

library(terra)
library(pastclim)
sea_ext<- terra::ext(110, 120, -5, 5)
tavg_vars <- c(paste0("temperature_0",1:9),paste0("temperature_",10:12))
time_steps <- get_time_steps(dataset = "Example")

The series for this region can be generated with the following command (note that you will need the appropriate variables for the Beyer2020 dataset downloaded on your computer, use download_dataset if needed):

tavg_series <- region_series(bio_variables =tavg_vars,
                             time_bp =  time_steps,
                             dataset = "Beyer2020",
                              ext = sea_ext)

tavg_series <- terra::sds(system.file("extdata/delta/tavg_series.nc",
                                       package="pastclim"))

Downscaling is performed one variable at a time. We will start with temperature in January. So, we first need to extract the SpatRaster of model low resolution data from the SpatRasterDataset:

tavg_model_lres_rast <- tavg_series$temperature_01
tavg_model_lres_rast

And we can now plot it:

plot(tavg_model_lres_rast, main = time_bp(tavg_model_lres_rast))

We can see how that the reconstructions are rather coarse (the Beyer2020 dataset uses 0.5x0.5 degree cells). We now need a set of high resolutions observations for the variable of interest that we will use to generate the delta raster used to downscale reconstructions. We will use data from WorldClim2 at 10 minute resolution (but other datasets such as CHELSA would be equally suitable). For WorldClim, there are some convenient functions to download and then load variables:

pastclim:::download_worldclim("tavg", 10)

Once the variable is downloaded, we can load it at any time with:

tavg_obs_hres_all<- pastclim:::load_worldclim("tavg", 10)

tavg_obs_hres_all<- terra::rast(system.file("extdata/delta/tavg_obs_hres_all.nc",
                                            package="pastclim"))

We want to crop these reconstructions to the extent of interest

tavg_obs_hres_all <- terra::crop(tavg_obs_hres_all, sea_ext)
# extract the January raster
tavg_obs_hres_rast <- tavg_obs_hres_all[[1]]
plot(tavg_obs_hres_rast)

We need to make sure that the extent of the modern observations is the same as the extent of the model reconstructions:

ext(tavg_obs_hres_rast)==ext(tavg_model_lres_rast)

If that was not the case, we would use terra::crop to match the extents.

We also need a high resolution global relief map (i.e. integrating both topographic and bathymetric values) to reconstruct past coastlines following sea level change. The relief raster will need to have the same extent and resolution as the high resolution observations. We can download one (based on ETOPO2022) with:

relief_rast <- pastclim:::download_relief(tavg_obs_hres_rast)

relief_rast <- terra::rast(system.file("extdata/delta/relief.nc",
                                       package="pastclim"))

We can quickly confirm that the resulting relief raster has the same extent as the model reconstructions, as well as the same resolution of the high resolution climate observations:

ext(relief_rast) == ext(tavg_model_lres_rast)

ext(relief_rast) == ext(tavg_obs_hres_rast)
ncol(relief_rast) == ncol(tavg_obs_hres_rast)
nrow(relief_rast) == nrow(tavg_obs_hres_rast)

We can now generate a high resolution land mask for the periods of interest. By default, we use the sea level reconstructions from Spratt et al 2016, but a different reference can be used by setting sea levels for each time step (see the man page for make_land_mask for details):

high_res_mask <- pastclim:::make_land_mask(relief_rast = relief_rast, 
                                time_bp = time_bp(tavg_model_lres_rast))
plot(high_res_mask, main=time_bp(high_res_mask))

We can now compute a delta raster and use it to downscale the model reconstructions:

delta_rast<-pastclim:::delta_compute(x=tavg_model_lres_rast, ref_time = 0, 
                                     obs = tavg_obs_hres_rast)
model_downscaled <- pastclim:::delta_downscale (x = tavg_model_lres_rast, 
                                                delta_rast = delta_rast,
                                                x_landmask_high = high_res_mask)
model_downscaled

Let's inspect the resulting data:

plot(model_downscaled, main = time_bp(model_downscaled))

And, as a reminder, the original reconstructions (note that the colour scales are not the same! terra chooses a scale for each time step based on the time specific range):

plot(tavg_model_lres_rast, main = time_bp(tavg_model_lres_rast))

Computing the bioclim variables

To compute the bioclim variables, we need to repeat the procedure above for temperature and precipitation for all months. Let us start with temperature. We loop over each month, create a SpatRaster of downscaled temperature, add it to a list, and finally convert the list into a SpatRasterDataset

tavg_downscaled_list<-list()
for (i in 1:12){
  delta_rast<-pastclim:::delta_compute(x=tavg_series[[i]], ref_time = 0, 
                                       obs = tavg_obs_hres_all[[i]])
  tavg_downscaled_list[[i]] <- pastclim:::delta_downscale (x = tavg_series[[i]], 
                                                  delta_rast = delta_rast,
                                                  x_landmask_high = high_res_mask)
}
tavg_downscaled <- terra::sds(tavg_downscaled_list)

Quickly inspect the resulting dataset:

tavg_downscaled

As expected, we have 12 months (subdatasets), each with 5 time steps.

We now need to create a series for precipitation:

prec_vars <- c(paste0("precipitation_0",1:9),paste0("precipitation_",10:12))
prec_series <- region_series(bio_variables = prec_vars,
                             time_bp =  time_steps,
                             dataset = "Beyer2020",
                             ext = sea_ext)

prec_vars <- c(paste0("precipitation_0",1:9),paste0("precipitation_",10:12))
prec_series <- terra::sds(system.file("extdata/delta/prec_series.nc",
                           package="pastclim"))

Get some high resolution observations:

pastclim:::download_worldclim("prec",10)
prec_obs_hres_all <- pastclim:::load_worldclim("prec",10)
prec_obs_hres_all <- terra::crop(prec_obs_hres_all, sea_ext)

prec_obs_hres_all <-terra::rast(system.file("extdata/delta/prec_obs_hres_all.nc",
                                            package = "pastclim"))

And finally downscale precipitation:

prec_downscaled_list<-list()
for (i in 1:12){
  delta_rast<-pastclim:::delta_compute(x=prec_series[[i]], ref_time = 0, 
                                       obs = prec_obs_hres_all[[i]])
  prec_downscaled_list[[i]] <- pastclim:::delta_downscale (x = prec_series[[i]], 
                                                  delta_rast = delta_rast,
                                                  x_landmask_high = high_res_mask)
}
prec_downscaled <- terra::sds(prec_downscaled_list)

We are now ready to compute the bioclim variables:

bioclim_downscaled<-bioclim_vars(tavg = tavg_downscaled, prec = prec_downscaled)

Let's inspect the object:

bioclim_downscaled

And plot the first variable (bio01):

plot(bioclim_downscaled[[1]], main = time_bp(bioclim_downscaled[[1]]))

We can now save the downscaled sds to a netcdf file:

terra::writeCDF(bioclim_downscaled,paste0(tempdir(),"/EA_bioclim_downscaled.nc"), overwrite=TRUE)

And then use it as a custom dataset for any function in pastclim. Let's extract a region series for three variables:

custom_data <- region_series(bio_variables =c("bio01","bio04","bio19"),
                             dataset = "custom",
                             path_to_nc = paste0(tempdir(),"/EA_bioclim_downscaled.nc"))

We can quickly inspect the resulting sds object:

custom_data

And plot it:

plot(custom_data$bio01,main=time_bp(custom_data$bio01))

EvolEcolGroup/pastclim documentation built on April 5, 2025, 9:22 a.m.