Generation of ENVIREM variables from a set of monthly climatic rasters
Pascal Title
2023-08-23
Rasters of ENVIREM variables were made available when this package was originally published in 2017. Those products were based on Worldclim v1.4 monthly grids. We anticipate preparing an updated precomputed dataset in the near future. However, as the number of possible climate projections and input datasets is large, we have written this R package to enable users to generate these rasters on their own, based on any input dataset.
In this vignette we will demonstrate how to generate the ENVIREM variables from a set of monthly temperature and precipitation variables. Some of the steps presented here are purely organizational, and worked well for the authors, but could be handled differently.
Major differences from previous versions of this vignette
The envirem R package has been updated in a few important ways.
- Solar radiation rasters are now generated by functions in the R package (version 1.4).
- Input rasters no longer need to be renamed - rather, the naming scheme is provided to the R functions (version 1.4).
- The 19 bioclimatic variables are no longer required as inputs (version 2.0).
- Removal of dependencies
rasterandrgdaland switch toterrafunctionality (version 3.0). - Restructuring and simplification of the R package (version 3.0).
Goal
The goal in this tutorial is to generate the climatic ENVIREM variables for western North America for a future climate projection.
Set up
In order to accomplish this, we will need:
- the envirem package
- the terra package
- the sf and rnaturalearth packages for some polygon work and mapping
- the viridis package for some nice color ramps
- monthly rasters for precipitation, max and min temperature, and mean temperature if available
Acquiring monthly rasters
In this tutorial, we will generate the ENVIREM rasters for current climate, based on the CHELSA-CLIM dataset v2.1. We downloaded monthly min/max/mean temperature and precipitation variables for the present at a 30 arcsecond resolution.
Note that if monthly mean temperature is not available (as is often the case), the envirem package will calculate it internally as the mean of the min and max temperature.
These rasters, in GeoTiff format, were placed in a directory that we (arbitrarily) named chelsaClim/global:
list.files("./chelsaClim/global")## [1] "CHELSA_pr_01_1981-2010_V.2.1.tif"
## [2] "CHELSA_pr_02_1981-2010_V.2.1.tif"
## [3] "CHELSA_pr_03_1981-2010_V.2.1.tif"
## [4] "CHELSA_pr_04_1981-2010_V.2.1.tif"
## [5] "CHELSA_pr_05_1981-2010_V.2.1.tif"
## [6] "CHELSA_pr_06_1981-2010_V.2.1.tif"
## [7] "CHELSA_pr_07_1981-2010_V.2.1.tif"
## [8] "CHELSA_pr_08_1981-2010_V.2.1.tif"
## [9] "CHELSA_pr_09_1981-2010_V.2.1.tif"
## [10] "CHELSA_pr_10_1981-2010_V.2.1.tif"
## [11] "CHELSA_pr_11_1981-2010_V.2.1.tif"
## [12] "CHELSA_pr_12_1981-2010_V.2.1.tif"
## [13] "CHELSA_tas_01_1981-2010_V.2.1.tif"
## [14] "CHELSA_tas_02_1981-2010_V.2.1.tif"
## [15] "CHELSA_tas_03_1981-2010_V.2.1.tif"
## [16] "CHELSA_tas_04_1981-2010_V.2.1.tif"
## [17] "CHELSA_tas_05_1981-2010_V.2.1.tif"
## [18] "CHELSA_tas_06_1981-2010_V.2.1.tif"
## [19] "CHELSA_tas_07_1981-2010_V.2.1.tif"
## [20] "CHELSA_tas_08_1981-2010_V.2.1.tif"
## [21] "CHELSA_tas_09_1981-2010_V.2.1.tif"
## [22] "CHELSA_tas_10_1981-2010_V.2.1.tif"
## [23] "CHELSA_tas_11_1981-2010_V.2.1.tif"
## [24] "CHELSA_tas_12_1981-2010_V.2.1.tif"
## [25] "CHELSA_tasmax_01_1981-2010_V.2.1.tif"
## [26] "CHELSA_tasmax_02_1981-2010_V.2.1.tif"
## [27] "CHELSA_tasmax_03_1981-2010_V.2.1.tif"
## [28] "CHELSA_tasmax_04_1981-2010_V.2.1.tif"
## [29] "CHELSA_tasmax_05_1981-2010_V.2.1.tif"
## [30] "CHELSA_tasmax_06_1981-2010_V.2.1.tif"
## [31] "CHELSA_tasmax_07_1981-2010_V.2.1.tif"
## [32] "CHELSA_tasmax_08_1981-2010_V.2.1.tif"
## [33] "CHELSA_tasmax_09_1981-2010_V.2.1.tif"
## [34] "CHELSA_tasmax_10_1981-2010_V.2.1.tif"
## [35] "CHELSA_tasmax_11_1981-2010_V.2.1.tif"
## [36] "CHELSA_tasmax_12_1981-2010_V.2.1.tif"
## [37] "CHELSA_tasmin_01_1981-2010_V.2.1.tif"
## [38] "CHELSA_tasmin_02_1981-2010_V.2.1.tif"
## [39] "CHELSA_tasmin_03_1981-2010_V.2.1.tif"
## [40] "CHELSA_tasmin_04_1981-2010_V.2.1.tif"
## [41] "CHELSA_tasmin_05_1981-2010_V.2.1.tif"
## [42] "CHELSA_tasmin_06_1981-2010_V.2.1.tif"
## [43] "CHELSA_tasmin_07_1981-2010_V.2.1.tif"
## [44] "CHELSA_tasmin_08_1981-2010_V.2.1.tif"
## [45] "CHELSA_tasmin_09_1981-2010_V.2.1.tif"
## [46] "CHELSA_tasmin_10_1981-2010_V.2.1.tif"
## [47] "CHELSA_tasmin_11_1981-2010_V.2.1.tif"
## [48] "CHELSA_tasmin_12_1981-2010_V.2.1.tif"Clip rasters to western USA
As we chose the western US in this particular demonstration, we will clip our input rasters to a relevant extent, which is much less computationally intensive than generating the ENVIREM variables at a global extent, and clipping those later.
We will first define an extent with a polygon. Here we have already
identified coordinates that, when linked, create a polygon that
encompasses western North America. We will then read in each input
raster, crop/mask it and write to disk, to a directory called
NorthAmerica. The file names will be the same as the
original files, although this could be set up differently if
desired.
The specifics of how this step is done is entirely up to the user, and not unique to the envirem package.
library(terra)## terra 1.7.39library(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(rnaturalearth)## 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)## Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`library(viridis)## Loading required package: viridisLitepolyWKT <- "POLYGON((
-130 50,
-100 50,
-100 30,
-130 30,
-130 50
))"
# convert wkt string into a sf polygon
# crs = 4326 because these are longitude/latitude coordinates
poly <- st_as_sfc(polyWKT, crs = 4326)
poly## Geometry set for 1 feature
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -130 ymin: 30 xmax: -100 ymax: 50
## Geodetic CRS: WGS 84## POLYGON ((-130 50, -100 50, -100 30, -130 30, -...land <- ne_download(scale = 50, type = 'land', category = 'physical', returnclass = 'sf')
land <- st_geometry(land)
plot(poly, border = 'blue', xpd = NA)
plot(land, border = NA, col = gray(0.9, alpha = 0.5), add = TRUE, xpd = NA)
Climate products are often not masked to terrestrial environments, so we may want to do so ourselves. One way to do this would be with the landmasses polygon we downloaded above. We will also crop our input rasters to the extent of our bounding box.
# crop and mask input rasters
inputDir <- "./chelsaClim/global"
files <- list.files(inputDir, pattern = '.tif$', full.names = TRUE)
# As all rasters that we downloaded share the same reference grid, we can read in and process them all as a stack
clim <- rast(files)
# crop to poly and mask with land
clim <- crop(clim, poly)
clim <- mask(clim, vect(land))
# We will coarsen these rasters for the purposes of this tutorial.
clim <- aggregate(clim, 4)
# check an example
plot(clim[[1]])
Now all monthly rasters are clipped to the desired extent, and masked to terrestrial environments.
Units
When working with climatic datasets for different time periods, it is imperative to have units that are in agreement. For instance, it is not uncommon to have data for the present in degrees C for temperature and in mm for precipitation, but to have temperature in degrees C * 10 for the past or future. In the case of climate layers from CHELSA, the documentation indicates that current climate temperature rasters are in degrees C, but that monthly temperature rasters for the LGM, for instance, are in degrees C / 10. Although there is capacity in the envirem R package to tell the functions what the units are, it is more straightforward to define the units as preferred ahead of time. Easiest might be to just have all monthly temperature rasters in deg C and precipitation rasters in mm, across all datasets you plan to work with. This is the stage at which this should happen.
We will now write our cropped and masked rasters to file. We will create an output folder called westernNA.
outputDir <- "./chelsaClim/westernNA"
if (!dir.exists(outputDir)) {
dir.create(outputDir)
}
# define file names
fn <- paste0(outputDir, '/', names(clim), '.tif')
writeRaster(clim, filename = fn, overwrite = TRUE)Raster naming scheme
The only requirement for the monthly climate raster file names is that they follow a standardized and predictable format: The different variables should have distinct names (precipitation vs min temperature vs max temperature, etc), and the numbers 1:12 should be present to indicate which month each raster corresponds to.
In the files downloaded and listed above, we can see that the naming scheme is as follows:
names(clim)## [1] "CHELSA_pr_01_1981-2010_V.2.1" "CHELSA_pr_02_1981-2010_V.2.1"
## [3] "CHELSA_pr_03_1981-2010_V.2.1" "CHELSA_pr_04_1981-2010_V.2.1"
## [5] "CHELSA_pr_05_1981-2010_V.2.1" "CHELSA_pr_06_1981-2010_V.2.1"
## [7] "CHELSA_pr_07_1981-2010_V.2.1" "CHELSA_pr_08_1981-2010_V.2.1"
## [9] "CHELSA_pr_09_1981-2010_V.2.1" "CHELSA_pr_10_1981-2010_V.2.1"
## [11] "CHELSA_pr_11_1981-2010_V.2.1" "CHELSA_pr_12_1981-2010_V.2.1"
## [13] "CHELSA_tas_01_1981-2010_V.2.1" "CHELSA_tas_02_1981-2010_V.2.1"
## [15] "CHELSA_tas_03_1981-2010_V.2.1" "CHELSA_tas_04_1981-2010_V.2.1"
## [17] "CHELSA_tas_05_1981-2010_V.2.1" "CHELSA_tas_06_1981-2010_V.2.1"
## [19] "CHELSA_tas_07_1981-2010_V.2.1" "CHELSA_tas_08_1981-2010_V.2.1"
## [21] "CHELSA_tas_09_1981-2010_V.2.1" "CHELSA_tas_10_1981-2010_V.2.1"
## [23] "CHELSA_tas_11_1981-2010_V.2.1" "CHELSA_tas_12_1981-2010_V.2.1"
## [25] "CHELSA_tasmax_01_1981-2010_V.2.1" "CHELSA_tasmax_02_1981-2010_V.2.1"
## [27] "CHELSA_tasmax_03_1981-2010_V.2.1" "CHELSA_tasmax_04_1981-2010_V.2.1"
## [29] "CHELSA_tasmax_05_1981-2010_V.2.1" "CHELSA_tasmax_06_1981-2010_V.2.1"
## [31] "CHELSA_tasmax_07_1981-2010_V.2.1" "CHELSA_tasmax_08_1981-2010_V.2.1"
## [33] "CHELSA_tasmax_09_1981-2010_V.2.1" "CHELSA_tasmax_10_1981-2010_V.2.1"
## [35] "CHELSA_tasmax_11_1981-2010_V.2.1" "CHELSA_tasmax_12_1981-2010_V.2.1"
## [37] "CHELSA_tasmin_01_1981-2010_V.2.1" "CHELSA_tasmin_02_1981-2010_V.2.1"
## [39] "CHELSA_tasmin_03_1981-2010_V.2.1" "CHELSA_tasmin_04_1981-2010_V.2.1"
## [41] "CHELSA_tasmin_05_1981-2010_V.2.1" "CHELSA_tasmin_06_1981-2010_V.2.1"
## [43] "CHELSA_tasmin_07_1981-2010_V.2.1" "CHELSA_tasmin_08_1981-2010_V.2.1"
## [45] "CHELSA_tasmin_09_1981-2010_V.2.1" "CHELSA_tasmin_10_1981-2010_V.2.1"
## [47] "CHELSA_tasmin_11_1981-2010_V.2.1" "CHELSA_tasmin_12_1981-2010_V.2.1"- precipitation: CHELSA_pr_##_1981-2010_V.2.1
- min temperature: CHELSA_tasmin_##_1981-2010_V.2.1
- max temperature: CHELSA_tasmax_##_1981-2010_V.2.1
- mean temperature: CHELSA_tas_##_1981-2010_V.2.1
… where ## denotes the month number. The ##
symbol will be used to denote month, therefore it should not be present
elsewhere in the file names.
In the envirem R package, we specify the naming scheme via the
assignNames() function, which stores these names in a
special R environment. This environment is then queried by the various
envirem R functions. This only needs to be done once per R session – it
is not necessary to do this for each call to particular functions.
The function namingScheme() returns the naming scheme
that has been set. As we have not yet defined one, it lists the
defaults:
library(envirem)## Loading required package: palinsol##
## Please see the vignette at https://ptitle.github.io/envirem/ for a detailed walk-through of this package.namingScheme()##
##
## variable prefix index suffix
## --------- ----------- ------- -------
## tmin: tmin_ ##
## tmax: tmax_ ##
## tmean: tmean_ ##
## precip: precip_ ##
## solrad: et_solrad_ ##
## pet: PET_ ####
## To change these values, see ?assignNames.This shows us the naming scheme for temperature and precipitation rasters, but also for solar radiation and potential evapotranspiration. We’ll get to the latter two later.
We can now assign our naming scheme
assignNames(tmax = "CHELSA_tasmax_##_1981-2010_V.2.1",
tmin = "CHELSA_mintas_##_1981-2010_V.2.1",
tmean = "CHELSA_tas_##_1981-2010_V.2.1",
precip = "CHELSA_pr_##_1981-2010_V.2.1",
)
namingScheme()##
##
## variable prefix index suffix
## --------- --------------- ------- -----------------
## tmin: CHELSA_mintas_ ## _1981-2010_V.2.1
## tmax: CHELSA_tasmax_ ## _1981-2010_V.2.1
## tmean: CHELSA_tas_ ## _1981-2010_V.2.1
## precip: CHELSA_pr_ ## _1981-2010_V.2.1
## solrad: et_solrad_ ##
## pet: PET_ ####
## To change these values, see ?assignNames.We can verify that our naming scheme is properly recognizing the appropriate rasters. This function will print messages if the defined naming scheme does not match exactly with the rasters that would be expected.
verifyRasterNames(clim)## The following rasters are missing or are not recognized: CHELSA_mintas_01_1981-2010_V.2.1, CHELSA_mintas_02_1981-2010_V.2.1, CHELSA_mintas_03_1981-2010_V.2.1, CHELSA_mintas_04_1981-2010_V.2.1, CHELSA_mintas_05_1981-2010_V.2.1, CHELSA_mintas_06_1981-2010_V.2.1, CHELSA_mintas_07_1981-2010_V.2.1, CHELSA_mintas_08_1981-2010_V.2.1, CHELSA_mintas_09_1981-2010_V.2.1, CHELSA_mintas_10_1981-2010_V.2.1, CHELSA_mintas_11_1981-2010_V.2.1, CHELSA_mintas_12_1981-2010_V.2.1
## Ensure that you have defined the proper naming scheme. See ?assignNames.Here, the function could not find certain rasters that it deemed missing. In fact, what happened is we incorrectly defined the naming scheme for min temp above. Let’s fix that.
assignNames(tmin = "CHELSA_tasmin_##_1981-2010_V.2.1")
namingScheme()##
##
## variable prefix index suffix
## --------- --------------- ------- -----------------
## tmin: CHELSA_tasmin_ ## _1981-2010_V.2.1
## tmax: CHELSA_tasmax_ ## _1981-2010_V.2.1
## tmean: CHELSA_tas_ ## _1981-2010_V.2.1
## precip: CHELSA_pr_ ## _1981-2010_V.2.1
## solrad: et_solrad_ ##
## pet: PET_ ####
## To change these values, see ?assignNames.verifyRasterNames(clim)## Names appear to be correct!Now we are all set.
Solar radiation rasters
Previously, extraterrestrial solar radiation needed to be acquired elsewhere, but these can now be generated by the envirem R package, thanks to the palinsol R package.
To generate monthly extraterrestrial solar radiation, we need to provide:
- a climatic raster, which will be used as a grid template,
- a year, which will be used to determine the angle of the sun,
- an output directory where the rasters will be written. The names of the output rasters will be defined by the naming scheme defined previously.
In the palinsol R package, year = 0 corresponds to the
year 1950. Therefore, as this is a climate projection for the years
1981-2010, we can define the year as 1995 - 1950 = 45.
# we will use clim[[1]] as our template
# calculate monthly solar radiation, defined for the year 1995, output to the output directory
ETsolradRasters(rasterTemplate = clim[[1]], year = 45, outputDir = "./chelsaClim/westernNA", overwrite = TRUE)## 1 2 3 4 5 6 7 8 9 10 11 12# let's read them in
solrad <- rast(list.files('./chelsaClim/westernNA', pattern = 'solrad', full.names = TRUE))And we can verify the file structure again.
verifyRasterNames(masterstack = clim, solradstack = solrad)## Names appear to be correct!We are now ready to generate new variables.
How is the envirem R package structured?
In the envirem R package, there are a number of functions to produce
individual climatic variables, such as embergerQ(),
monthlyPET(), thermicityIndex(), etc. These
can be used on their own if the necessary input rasters are available as
SpatRaster objects in R.
The function generateEnvirem() will take as input
monthly temperature/precipitation/solar radiation rasters, and will
output ENVIREM variables. This function will also calculate other
variables that may be needed.
Calling functions for individual variables
A number of envirem variables require as input some of the 19
bioclimatic variables. The user can either acquire those bioclimatic
variables elsewhere, calculate them themselves (for example, by using
the biovars() function in the dismo
package), or have the envirem package calculate them via
generateEnvirem().
For example, if we were interested in calculating growing degree days, with a base temperature of 5 deg C, we can directly call the specific function:
# Inputs needed will be monthly mean temperature only.
meanTempVar <- grep('tas_', names(clim), value = TRUE)
meanTempVar## [1] "CHELSA_tas_01_1981-2010_V.2.1" "CHELSA_tas_02_1981-2010_V.2.1"
## [3] "CHELSA_tas_03_1981-2010_V.2.1" "CHELSA_tas_04_1981-2010_V.2.1"
## [5] "CHELSA_tas_05_1981-2010_V.2.1" "CHELSA_tas_06_1981-2010_V.2.1"
## [7] "CHELSA_tas_07_1981-2010_V.2.1" "CHELSA_tas_08_1981-2010_V.2.1"
## [9] "CHELSA_tas_09_1981-2010_V.2.1" "CHELSA_tas_10_1981-2010_V.2.1"
## [11] "CHELSA_tas_11_1981-2010_V.2.1" "CHELSA_tas_12_1981-2010_V.2.1"gdd <- growingDegDays(clim[[meanTempVar]], baseTemp = 5)
plot(gdd, col = inferno(100))
Some other variables, such as the climatic moisture index, require as
inputs bioclim 12 (annual precipitation) and annual PET (total annual
potential evapotranspiration). If the user has these, then the function
can be run directly with those inputs. Otherwise, the function
generateEnvirem() will generate these intermediate products
internally, based on the monthly rasters we are providing.
cmi <- generateEnvirem(masterstack = clim, solradstack = solrad, var = 'climaticMoistureIndex')## ...splitting rasterstack...## ...calculating necessary bioclim variables...## ...annual PET...## ...climatic moisture index...plot(cmi, col = inferno(100))
Generating multiple variables
You can use generateEnvirem() to produce any number of
the envirem variables in one call. If you want to create all available
variables, you can specify var = 'all'.
Here, we will create all variables, using
generateEnvirem(). We will then write them to file.
# create output directory
outdir <- './chelsaClim/westernNA_envirem'
if (!dir.exists(outdir)) {
dir.create(outdir)
}
allGrids <- generateEnvirem(masterstack = clim, solradstack = solrad, var = 'all')## ...splitting rasterstack...## ...calculating necessary bioclim variables...## ...temp extremes...## ...growing degree days...## ...month count by deg...## ...continentality index...## ...thermicity index...## ...emberger's Q...## ...PET extremes...## ...annual PET...## ...PET seasonality...## ...climatic moisture index...## ...Thornthwaite aridity index...par(mfrow = c(5, 4))
for (i in 1:nlyr(allGrids)) {
plot(allGrids[[i]], col = inferno(100), box = FALSE, axes = FALSE)
title(main = names(allGrids)[i], cex = 0.65)
}
filenames <- paste0(outdir, '/', names(allGrids), '.tif')
writeRaster(allGrids, filenames, overwrite = TRUE)