Generation of ENVIREM variables from a set of monthly climatic rasters
Pascal O. Title
2019-11-21
Back to main project site.
We will strive to make available certain sets of precalculated versions of the ENVIREM rasters, but the number of possible climate projections and input datasets is such that users may need to generate these rasters on their own.
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.
- Input rasters no longer need to be renamed - rather, the naming scheme is provided to the R functions.
- The 19 bioclimatic variables are no longer required as inputs.
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 raster package
- the rgeos and maps packages for some polygon work and mapping
- the viridis package for some nice color ramps
- the rgdal package
- GDAL
- monthly rasters for precipitation, max and min temperature, and mean temperature if available
GDAL can be installed in a number of ways. If you have the rgdal R package, then gdal should be installed somewhere on your system. Otherwise, you can install GDAL from source via github, or via homebrew for Mac OS X. For our purposes, we want to make sure that two of GDAL’s programs, gdalinfo and gdal_translate are accessible from within R. To check that these functions will work, we can do the following:
## gdalinfo
## "/usr/local/bin/gdalinfo"
Sys.which('gdal_translate')
## gdal_translate
## "/usr/local/bin/gdal_translate"
If the functions were found, a path will be returned. If they were not found, an empty path will be returned. If these functions are not found in the default search paths, but you know specifically where they are, you can always provide the full path later on.
GDAL is only necessary for processing rasters as tiles. If you are having difficulties locating your GDAL programs, then you can still run the functions in this R package, but with no splitting into tiles.
Acquiring monthly rasters
In this tutorial, we will generate the ENVIREM rasters for a future climate scenario. We downloaded monthly min temperature, max temperature and precipitation variables for the CESM1-BGC model (Lindsay et al. 2014) under the RCP45 emission scenario from CHELSA, for the time period of 2041-2060, at a 30 arcsecond resolution (this particular scenario was randomly selected).
Note that if monthly mean temperature is available, it is preferable to acquire that as well. If not available, the envirem R functions will simply take the mean of min and max temperature, which we do here for convenience.
These rasters, in GeoTiff format, were placed in a directory that we (arbitrarily) named futureClim/global:
list.files("./futureClim/global")
## [1] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_1_2041-2060.tif"
## [2] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_10_2041-2060.tif"
## [3] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_11_2041-2060.tif"
## [4] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_12_2041-2060.tif"
## [5] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_2_2041-2060.tif"
## [6] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_3_2041-2060.tif"
## [7] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_4_2041-2060.tif"
## [8] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_5_2041-2060.tif"
## [9] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_6_2041-2060.tif"
## [10] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_7_2041-2060.tif"
## [11] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_8_2041-2060.tif"
## [12] "CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_9_2041-2060.tif"
## [13] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_1_2041-2060_V1.2.tif"
## [14] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_10_2041-2060_V1.2.tif"
## [15] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_11_2041-2060_V1.2.tif"
## [16] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_12_2041-2060_V1.2.tif"
## [17] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_2_2041-2060_V1.2.tif"
## [18] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_3_2041-2060_V1.2.tif"
## [19] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_4_2041-2060_V1.2.tif"
## [20] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_5_2041-2060_V1.2.tif"
## [21] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_6_2041-2060_V1.2.tif"
## [22] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_7_2041-2060_V1.2.tif"
## [23] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_8_2041-2060_V1.2.tif"
## [24] "CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_9_2041-2060_V1.2.tif"
## [25] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_1_2041-2060_V1.2.tif"
## [26] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_10_2041-2060_V1.2.tif"
## [27] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_11_2041-2060_V1.2.tif"
## [28] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_12_2041-2060_V1.2.tif"
## [29] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_2_2041-2060_V1.2.tif"
## [30] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_3_2041-2060_V1.2.tif"
## [31] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_4_2041-2060_V1.2.tif"
## [32] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_5_2041-2060_V1.2.tif"
## [33] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_6_2041-2060_V1.2.tif"
## [34] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_7_2041-2060_V1.2.tif"
## [35] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_8_2041-2060_V1.2.tif"
## [36] "CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_9_2041-2060_V1.2.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.
## Loading required package: sp
## rgeos version: 0.5-2, (SVN revision 621)
## GEOS runtime version: 3.7.2-CAPI-1.11.2
## Linking to sp version: 1.3-1
## Polygon checking: TRUE
library(maps)
library(viridis)
## Loading required package: viridisLite
poly <- readWKT("POLYGON((
-130 50,
-100 50,
-100 30,
-130 30,
-130 50
))", p4s = CRS("+proj=longlat +datum=WGS84"))
map(database='state', lwd = 0.5)
plot(poly, border = 'blue', xpd = NA, add = TRUE)
Future projection CHELSA climate grids are not masked to terrestrial environments, but one way to do this ourselves is to use a precipitation raster, and convert values < 0 to NA (as precip cannot be below zero mm).
# crop and mask input rasters
inputDir <- "./futureClim/global"
outputDir <- "./futureClim/NorthAmerica"
files <- list.files(inputDir, pattern = '.tif$', full.names = TRUE)
# read in a precipitation raster and crop/mask and keep as a mask layer
## precipitation rasters have the term "_pr_" in the file name
## crop/mask using polygon
grep(files, pattern = '_pr_', value=TRUE)[1]
## [1] "./futureClim/global/CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_1_2041-2060.tif"
precipRaster <- raster(grep(files, pattern = '_pr_', value=TRUE)[1])
precipRaster <- crop(precipRaster, poly)
precipRaster <- mask(precipRaster, poly)
# create a terrestrial/marine mask by setting values < 0 to NA, and values >= 0 to 1
precipRaster[precipRaster < 0] <- NA
precipRaster[!is.na(precipRaster)] <- 1
plot(precipRaster, col = 'blue', legend = FALSE)
title(main = 'terrestrial mask')
Now we have a terrestrial mask that we can apply to the other rasters.
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 future climate layers from CHELSA, we can see that precipitation is in mm, and temperature is in degrees C * 10. Although there is capacity in the envirem R package to tell the functions what the units are, I find it more straightforward to define the units as preferred ahead of time. Therefore, here, we will convert the temperature rasters to degree C. Precipitation rasters can remain in mm.
Therefore, we can do the following to each global raster:
- crop/mask using the polygon defined above
- use the precip mask to set marine cells to
NA
- convert temperature rasters to degrees C
- coarsen the rasters 2x to speed up processing for this demonstration
Process all rasters
# supply some additional flags to writeRaster() to employ compression to achieve smaller file sizes
tifOptions <- c("COMPRESS=DEFLATE", "PREDICTOR=2", "ZLEVEL=6")
for (i in 1:length(files)) {
cat(i, ' ')
r <- raster(files[i])
r <- crop(r, poly)
# apply terrestrial mask
r <- mask(r, precipRaster)
# if temperature, divide values by 10
## recognize temperature files as those containing the terms tasmin or tasmax
if (grepl('tasmin|tasmax', files[i])) {
r <- r / 10
}
r <- aggregate(r, fact = 2)
outfile <- paste0(outputDir, '/', names(r), '.tif')
writeRaster(r, filename = outfile, format = 'GTiff', options = tifOptions, overwrite = TRUE)
}
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
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:
- precipitation: CHELSA_pr_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_##_2041-2060.tif
- min temperature: CHELSA_tasmin_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_##_2041-2060_V1.2.tif
- max temperature: CHELSA_tasmax_mon_CESM1-BGC_rcp45_r1i1p1_g025.nc_##_2041-2060_V1.2.tif
… 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 varnames() returns the naming scheme that has been set. As we have not yet defined one, it lists the defaults:
## Loading required package: palinsol
## Loading required package: gsl
## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo
##
## Please see the vignette at https://envirem.github.io/ENVIREM_tutorial.html for a detailed walk-through of this package.
##
## tmin: tmin_ ##
## tmax: tmax_ ##
## tmean: tmean_ ##
## precip: precip_ ##
## solrad: et_solrad_ ##
##
##
## To change these values, see ?assignNames.
We can now assign our naming scheme
assignNames(tmax = "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_##_2041.2060_V1.2",
tmin = "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_##_2041.2060_V1.2",
precip = "CHELSA_pr_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_##_2041.2060",
solrad = "solar_##_westernUSA"
)
varnames()
##
## tmin: CHELSA_tasmin_mo ## _2041.2060_V1.2
## n_CESM1.BGC_rcp4
## 5_r1i1p1_g025.nc
## _
## tmax: CHELSA_tasmax_mo ## _2041.2060_V1.2
## n_CESM1.BGC_rcp4
## 5_r1i1p1_g025.nc
## _
## tmean: tmean_ ##
## precip: CHELSA_pr_mon_CE ## _2041.2060
## SM1.BGC_rcp45_r1
## i1p1_g025.nc_
## solrad: solar_ ## _westernUSA
##
##
## To change these values, see ?assignNames.
Note that we specified a naming scheme for solar radiation rasters, which will be next in this tutorial.
We can verify that our naming scheme is properly recognizing the appropriate files. This function will print messages if the defined naming scheme does not match exactly with the files that would be expected.
verifyFileStructure('./futureClim/NorthAmerica', returnFileNames = FALSE)
## tmean files are not properly named or missing. Tmean will therefore be calculated.
## solrad files are not properly named.
Here, a message was returned for tmean, which makes sense as we did not include those files. A message was returned for solar radiation, as we defined the naming, but have not yet produced those files. Most important here is that no message was returned for tmin, tmax and precip.
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 to match the raster grid,
- 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 2041-2060, we can define the year as 2050 - 1950 = 100.
# read in a climatic raster for use as a template
rasterTemplate <- raster('./futureClim/NorthAmerica/CHELSA_pr_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_1_2041.2060.tif')
# calculate monthly solar radiation, defined for the year 2050, output to the current directory
ETsolradRasters(rasterTemplate = rasterTemplate, year = 100, outputDir = './futureClim/NorthAmerica', overwrite = TRUE)
## 1 LA04 data loaded
## 2 3 4 5 6 7 8 9 10 11 12
And we can verify the file structure again.
verifyFileStructure('./futureClim/NorthAmerica', returnFileNames = FALSE)
## tmean files are not properly named or missing. Tmean will therefore be calculated.
We are now ready to generate new variables.
How is the envirem R package structured?
The envirem R package may at first be a bit confusing without a roadmap.
Briefly, 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 rasterLayer objects in R.
There are also utility functions, such as assignNames(), verifyFileStructure(), verifyRasterNames(), dataTypeCheck(), etc. Some of these are necessary, and some are convenience functions.
The two main functions in this package are layerCreation() and generateRasters(). Both will produce the envirem variables, but they are intended for different uses.
layerCreation() takes rasterStacks of input rasters and outputs a rasterStack of envirem variables. generateRasters() takes as input a path to a directory of raster files, and outputs raster files, written to disk. If requested, this function will also split the input rasters into tiles to reduce the computational load. Essentially, generateRasters() does some housekeeping, then calls layerCreation() internally, and then does some more housekeeping.
Therefore, generateRasters() is primarily intended for situations where it would be too computationally intensive to read in all rasters, and calculate the new variables.
Additionally, both main functions perform some necessary extra steps, such as calculating the bioclimatic variables that may be required, and calculating mean temperature if it is missing.
With the size of dataset used here for demonstration, we can do both for comparison.
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 layerCreation() or generateRasters().
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.
# Here, we will take the mean of max and min temperature,
# although if mean temp was available, that would be better.
files <- list.files('./futureClim/NorthAmerica', pattern = 'tasmin', full.names = TRUE)
mintempstack <- stack(files)
files <- list.files('./futureClim/NorthAmerica', pattern = 'tasmax', full.names = TRUE)
maxtempstack <- stack(files)
# We need to be sure the months are in chronological order.
names(mintempstack)
## [1] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_1_2041.2060_V1.2"
## [2] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_10_2041.2060_V1.2"
## [3] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_11_2041.2060_V1.2"
## [4] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_12_2041.2060_V1.2"
## [5] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_2_2041.2060_V1.2"
## [6] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_3_2041.2060_V1.2"
## [7] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_4_2041.2060_V1.2"
## [8] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_5_2041.2060_V1.2"
## [9] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_6_2041.2060_V1.2"
## [10] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_7_2041.2060_V1.2"
## [11] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_8_2041.2060_V1.2"
## [12] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_9_2041.2060_V1.2"
## [1] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_1_2041.2060_V1.2"
## [2] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_10_2041.2060_V1.2"
## [3] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_11_2041.2060_V1.2"
## [4] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_12_2041.2060_V1.2"
## [5] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_2_2041.2060_V1.2"
## [6] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_3_2041.2060_V1.2"
## [7] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_4_2041.2060_V1.2"
## [8] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_5_2041.2060_V1.2"
## [9] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_6_2041.2060_V1.2"
## [10] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_7_2041.2060_V1.2"
## [11] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_8_2041.2060_V1.2"
## [12] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_9_2041.2060_V1.2"
# They are not. We will fix this by adding a zero to single digit months, and sorting.
names(mintempstack) <- gsub('_(\\d)_', '_0\\1_', names(mintempstack))
names(maxtempstack) <- gsub('_(\\d)_', '_0\\1_', names(maxtempstack))
mintempstack <- mintempstack[[order(names(mintempstack))]]
maxtempstack <- maxtempstack[[order(names(maxtempstack))]]
names(mintempstack)
## [1] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_01_2041.2060_V1.2"
## [2] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_02_2041.2060_V1.2"
## [3] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_03_2041.2060_V1.2"
## [4] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_04_2041.2060_V1.2"
## [5] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_05_2041.2060_V1.2"
## [6] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_06_2041.2060_V1.2"
## [7] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_07_2041.2060_V1.2"
## [8] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_08_2041.2060_V1.2"
## [9] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_09_2041.2060_V1.2"
## [10] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_10_2041.2060_V1.2"
## [11] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_11_2041.2060_V1.2"
## [12] "CHELSA_tasmin_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_12_2041.2060_V1.2"
## [1] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_01_2041.2060_V1.2"
## [2] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_02_2041.2060_V1.2"
## [3] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_03_2041.2060_V1.2"
## [4] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_04_2041.2060_V1.2"
## [5] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_05_2041.2060_V1.2"
## [6] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_06_2041.2060_V1.2"
## [7] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_07_2041.2060_V1.2"
## [8] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_08_2041.2060_V1.2"
## [9] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_09_2041.2060_V1.2"
## [10] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_10_2041.2060_V1.2"
## [11] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_11_2041.2060_V1.2"
## [12] "CHELSA_tasmax_mon_CESM1.BGC_rcp45_r1i1p1_g025.nc_12_2041.2060_V1.2"
# fixed!
# calculate a proxy for mean temperature
meantempstack <- (mintempstack + maxtempstack) / 2
gdd <- growingDegDays(meantempstack, baseTemp = 5)
plot(gdd, col = inferno(100))
For that to work, we had to do some extra steps to prepare the input data. We could alternatively have the layerCreation() function do this all for us: It will create missing variables like mean temp or bioclimatic variables, and ensure month-order is correct.
chelsaFiles <- list.files('./futureClim/NorthAmerica', pattern = 'CHELSA', full.names = TRUE)
chelsaStack <- stack(chelsaFiles)
solarFiles <- list.files('./futureClim/NorthAmerica', pattern = 'solar', full.names = TRUE)
solarStack <- stack(solarFiles)
# We can verify that variable names will be correctly recognized.
verifyRasterNames(chelsaStack, solradstack = solarStack)
## Names appear to be correct!
gdd <- layerCreation(chelsaStack, solarStack, var = 'growingDegDays5')
## ...splitting rasterstack...
## ...calculating mean temp...
## ...calculating necessary bioclim variables...
## ...growing degree days...
plot(gdd, col = inferno(100))
We could also generate GDD by just specifying the file paths with generateRasters(). This is not the most convenient in this particular situation, but helpful for demonstration purposes:
generateRasters(var = 'growingDegDays5',
prefix = 'test_',
maindir = './futureClim/NorthAmerica',
outputDir = './futureClim/NorthAmerica'
)
## tmean files are not properly named or missing. Tmean will therefore be calculated.
## ...splitting rasterstack...
## ...calculating mean temp...
## ...calculating necessary bioclim variables...
## ...growing degree days...
If this was a high resolution global dataset, we could request that the function split the input variables into tiles. Here, we can request that the function split inputs into 4 tiles, calculate variables, and then merge the tiles together.
generateRasters(var = 'growingDegDays5',
prefix = 'tileTest_',
maindir = './futureClim/NorthAmerica',
outputDir = './futureClim/NorthAmerica',
nTiles = 4
)
## tmean files are not properly named or missing. Tmean will therefore be calculated.
## Splitting rasters into tiles...
## _tile1
## ...splitting rasterstack...
## ...calculating mean temp...
## ...calculating necessary bioclim variables...
## ...growing degree days...
## _tile2
## ...splitting rasterstack...
## ...calculating mean temp...
## ...calculating necessary bioclim variables...
## ...growing degree days...
## _tile3
## ...splitting rasterstack...
## ...calculating mean temp...
## ...calculating necessary bioclim variables...
## ...growing degree days...
## _tile4
## ...splitting rasterstack...
## ...calculating mean temp...
## ...calculating necessary bioclim variables...
## ...growing degree days...
##
## Putting tiles back together...
## Tiles being combined for growingDegDays5.tif...
# read in the resulting raster
result <- raster('./futureClim/NorthAmerica/tileTest_growingDegDays5.tif')
plot(result, col = inferno(100))
Generating multiple variables
You can use layerCreation() or generateRasters() 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 generateRasters(), and we will again split up the work into tiles.
# create output directory
outdir <- './futureClim/NorthAmerica/envirem'
dir.create(outdir)
generateRasters(var = 'all',
maindir = './futureClim/NorthAmerica',
outputDir = outdir,
nTiles = 4
)
## tmean files are not properly named or missing. Tmean will therefore be calculated.
## Splitting rasters into tiles...
## _tile1
## ...splitting rasterstack...
## ...calculating mean temp...
## ...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...
## _tile2
## ...splitting rasterstack...
## ...calculating mean temp...
## ...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...
## _tile3
## ...splitting rasterstack...
## ...calculating mean temp...
## ...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...
## _tile4
## ...splitting rasterstack...
## ...calculating mean temp...
## ...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...
##
## Putting tiles back together...
## Tiles being combined for annualPET.tif...
## Tiles being combined for aridityIndexThornthwaite.tif...
## Tiles being combined for climaticMoistureIndex.tif...
## Tiles being combined for continentality.tif...
## Tiles being combined for embergerQ.tif...
## Tiles being combined for growingDegDays0.tif...
## Tiles being combined for growingDegDays5.tif...
## Tiles being combined for maxTempColdest.tif...
## Tiles being combined for minTempWarmest.tif...
## Tiles being combined for monthCountByTemp10.tif...
## Tiles being combined for PETColdestQuarter.tif...
## Tiles being combined for PETDriestQuarter.tif...
## Tiles being combined for PETseasonality.tif...
## Tiles being combined for PETWarmestQuarter.tif...
## Tiles being combined for PETWettestQuarter.tif...
## Tiles being combined for thermicityIndex.tif...
enviremRasters <- list.files(outdir, pattern = '\\.tif$', full.names = TRUE)
enviremRasters <- stack(enviremRasters)
par(mfrow = c(4, 4), mar = c(0.5, 0.5, 2, 0))
for (i in 1:nlayers(enviremRasters)) {
plot(enviremRasters[[i]], col = inferno(100), box = FALSE, axes = FALSE)
title(main = names(enviremRasters)[i])
}
