::install_github("Cidree/forestdata") remotes
1 Introduction
In this post, we will create a map of the change of the potential distribution of Pinus sylvestris under two climate scenarios (RCP 4.5 and RCP 8.5) in Europe.
We will use the EU-Trees4F dataset (Mauri et al. 2022), a dataset with the current (2005) and future (2035, 2065, 2095) potential distribution of 67 tree species in Europe at 10 km spatial resolution. The dataset is accessible through an R package that I am developing, and it’s accessible through GitHub.
The objective of this exercise is to compare the current (2005) and the future (2095) potential distribution of Pinus sylvestris, and determine the areas where it will maintain its distribution or absence, disappear, or expand.
2 Loading packages
First of all, I will announce that I am developing a package called {forestdata}
which gives you access to several sources to easily load forest and land cover data into R.
Note that the package is still in early development and it might be prone to errors.
In order to install the package, you need to have installed the {remotes}
package in your system, and then you can use:
After that, we can load the rest of the packages for this exercise:
library(pacman)
p_load(
## Core
tidyverse,
## Donwload data
forestdata,
## Spatial data manipulation
terra,
## Visualization
ggtext, patchwork, tidyterra )
We will use the {terra}
package to manipulate the data that we download with {forestdata}
, the {ggtext}
package to use markdown and HTML in our plots, the {patchwork}
package for assembling the plots, and {tidyterra}
which makes it easy to plot SpatRasters
in {ggplot2}
.
3 Load the data
The data that we need to create the maps are:
Current Pinus sylvestris potential distribution as for 2005
Future Pinus sylvestris potential distribution as for 2095 under RCP 4.5, a climate change scenario of intermediate emissions
Future Pinus sylvestris potential distribution as for 2095 under RCP 8.5, a climate change scenario o high emissions
To do so, we can use the fd_forest_eutrees4f
function of the {forestdata}
package as follows:
## Load Pinus sylvestris data for 2005
<- fd_forest_eutrees4f(
psylvestris_2005_sr species = "Pinus sylvestris",
period = 2005,
type = "bin",
distrib = "pot"
)
## Load Pinus sylvestris data for 2095 (RCP 4.5)
<- fd_forest_eutrees4f(
psylvestris_2095_rcp45_sr species = "Pinus sylvestris",
period = 2095,
scenario = "rcp45"
)
## Load Pinus sylvestris data for 2095 (RCP 8.5)
<- fd_forest_eutrees4f(
psylvestris_2095_rcp85_sr species = "Pinus sylvestris",
period = 2095,
scenario = "rcp85"
)
The arguments that we used are:
species
: Latin name of the species.period
: period of the data (2005, 2035, 2065, 2095).type
: the are three different types of maps. In this case we choose binary maps (1 = presence; 0 = absence). This is the option by default.distrib
: type of distribution. In this case, we chose the potential distribution which is set by default.scenario
: one of “rcp45” or “rcp85”. When period is 2005, this argument is ignored.
The list of possible species is available through {forestdata}
:
::eutrees4f_species forestdata
[1] "Abies alba" "Acer campestre"
[3] "Acer opalus" "Acer platanoides"
[5] "Acer pseudoplatanus" "Alnus glutinosa"
[7] "Alnus incana" "Arbutus unedo"
[9] "Aria edulis" "Betula pendula"
[11] "Betula pubescens" "Borkhausenia intermedia"
[13] "Carpinus betulus" "Carpinus orientalis"
[15] "Castanea sativa" "Celtis australis"
[17] "Ceratonia siliqua" "Cormus domestica"
[19] "Corylus avellana" "Cupressus sempervirens"
[21] "Fagus sylvatica" "Fraxinus angustifolia"
[23] "Fraxinus excelsior" "Fraxinus ornus"
[25] "Juglans regia" "Juniperus thurifera"
[27] "Larix decidua" "Laurus nobilis"
[29] "Malus sylvestris" "Olea europaea"
[31] "Ostrya carpinifolia" "Picea abies"
[33] "Pinus brutia" "Pinus cembra"
[35] "Pinus halepensis" "Pinus nigra"
[37] "Pinus pinaster" "Pinus pinea"
[39] "Pinus sylvestris" "Pistacia lentiscus"
[41] "Pistacia terebinthus" "Populus alba"
[43] "Populus nigra" "Populus tremula"
[45] "Prunus avium" "Prunus padus"
[47] "Pyrus communis" "Quercus cerris"
[49] "Quercus coccifera" "Quercus faginea"
[51] "Quercus frainetto" "Quercus ilex"
[53] "Quercus petraea" "Quercus pubescens"
[55] "Quercus pyrenaica" "Quercus robur"
[57] "Quercus suber" "Robinia pseudoacacia"
[59] "Salix alba" "Sorbus aucuparia"
[61] "Taxus baccata" "Tilia cordata"
[63] "Tilia platyphyllos" "Torminalis glaberrima"
[65] "Ulmus glabra" "Ulmus laevis"
[67] "Ulmus minor"
We can see how one of these rasters look like:
psylvestris_2095_rcp85_sr
class : SpatRaster
dimensions : 410, 400, 1 (nrow, ncol, nlyr)
resolution : 10000, 10000 (x, y)
extent : 2600000, 6600000, 1400000, 5500000 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035)
source : Pinus_sylvestris_ens-clim_rcp85_fut2095_bin_pot.tif
name : Pinus_sylvestris_ens-clim_rcp85_fut2095_bin_pot
min value : 0
max value : 1
4 Prepare the data
The three rasters that we downloaded have two values:
1 = presence of Pinus sylvestris
0 = absence of Pinus sylvestris
We need to combine the current raster with the other two in a way that we can create the four classes:
Absent: areas where P. sylvestris is absent, and it will still be absent
Will disappear: areas where P. sylvestris is present, but it won’t be in the future models
Will appear: areas where P. sylvestris is absent, but it will be present in the future models
Present and stable: areas where P. sylvestris is present, and it will still be present
To achieve this, I created the next function:
<- function(raster, classes) {
reclassify_rcp
## Reclassify 1 for 2
<- ifel(
raster_class == 1, 2, raster
raster
)## Sum to current (2005) distribution
<- as.factor(raster_class + psylvestris_2005_sr)
raster_class ## Rename levels
levels(raster_class)[[1]][, 2] <- classes
raster_class }
This function will classify the input rasters (RCP 4.5 and RCP 8.5) into 0 for absence and 2 for presence. Then, we sum up the psylvestris_2005_sr
so we have the following values:
0 + 0: absent
0 + 2: will appear
1 + 0: will disappear
1 + 2: present
Finally, we convert it to a categorical raster and use a vector of classes (also a vector of colors) to rename it. Let’s do this:
## Classes and colors for each class
<- c("Absent", "Will disappear", "Will appear", "Present and stable")
ps_classes <- c("#BE92A2", "#96031A", "#6DA34D", "#CEEDDB")
ps_colors ## Apply function
<- reclassify_rcp(psylvestris_2095_rcp45_sr, ps_classes)
psylvestris_rcp45_sr <- reclassify_rcp(psylvestris_2095_rcp85_sr, ps_classes) psylvestris_rcp85_sr
To verify that it worked correctly, let’s create a quick plot:
with(
par(mfrow = c(2, 1)), {
plot(psylvestris_rcp45_sr, main = "RCP 4.5", col = ps_colors)
plot(psylvestris_rcp85_sr, main = "RCP 8.5", col = ps_colors)
} )
The distribution of Pinus sylvestris seems correct to me, so let’s go to the next step.
Take into account that the spatial resolution is 10 km, and we are mapping the potential distribution.
5 Final map
To create the final map, I will create a function, since we need to apply it to two different objects:
<- function(data, title = "(a) RCP 4.5") {
map_rcp
ggplot() +
geom_spatraster(
data = data,
show.legend = FALSE
+
) scale_fill_manual(
values = ps_colors,
labels = ps_classes,
na.value = NA,
na.translate = FALSE
+
) labs(
title = title
+
) theme_void(base_family = "Roboto") +
theme(
plot.title = element_text(
face = "bold", hjust = .5, color = "snow"
)
)
}
Next, we apply the function to both rasters:
<- map_rcp(psylvestris_rcp45_sr)
rcp45_gg <- map_rcp(psylvestris_rcp85_sr, title = "(b) RCP 8.5") rcp85_gg
We will assemble both maps using {patchwork}
, but first I will create an appealing title for it. Instead of using a legend, we can use colors in our title to identify the different classes. There’s a new R package called {marquee}
which makes it easier without using HTML, but I found some limitations for long titles when creating line breaks and modifying the line height, so for the moment I will stick with {ggtext}
:
## Wrappers for ggtext
<- str_glue("<b style = 'color: {ps_colors[1]};'>almost absent</b>")
absent_txt <- str_glue("<b style = 'color: {ps_colors[2]};'>decrease</b>")
decrease_txt <- str_glue("<b style = 'color: {ps_colors[3]};'>shift</b>")
increase_txt <- str_glue("<b style = 'color: {ps_colors[4]};'>present</b>")
present_txt
## Plot title
<- str_glue(
title_txt "*Pinus sylvestris*, {absent_txt} in southern Europe and {present_txt} in Northern and Central Europe, is<br>
projected to {increase_txt} its potential distribution northward and {decrease_txt} in Central Europe under two<br>
climatic scenarios by 2095"
)
There, I create one label or each class using each of the colors defined before, and put the together in the title. Finally, we can assemble every piece of out plot with {patchwork}
:
+
rcp45_gg +
rcp85_gg plot_annotation(
title = title_txt,
caption = "Author: Adrián Cidre | https://adrian-cidre.com | Data source: EU-Trees4F",
theme = theme(
text = element_text(colour = "snow"),
plot.title = element_markdown(
family = "Merriweather",
face = "bold",
lineheight = 1.2,
margin = margin(t = 5, l = 5)
),plot.caption = element_text(
hjust = .5,
family = "Roboto"
),plot.background = element_rect(
fill = "gray10",
colour = "gray10"
)
) )
6 Conclusions
Creating a map of the future distribution of European tree species in R is very easy thanks to the EU-Trees4F dataset, and the {forestdata}
package.
If you found this post interesting and helpful, please share it with others or create your own, share it and tag me in linkedin.