Pinus sylvestris Distribution in Europe: Climate Scenarios Mapped with R

Create a map of the distribution of tree species under climate change scenarios
forestry
ggplot2
R
spatial
Author

Adrián Cidre

Published

June 17, 2024

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.

Caution

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:

remotes::install_github("Cidree/forestdata")

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
psylvestris_2005_sr <- fd_forest_eutrees4f(
  species  = "Pinus sylvestris",
  period   = 2005,
  type     = "bin",
  distrib  = "pot"
)

## Load Pinus sylvestris data for 2095 (RCP 4.5)
psylvestris_2095_rcp45_sr <- fd_forest_eutrees4f(
  species  = "Pinus sylvestris",
  period   = 2095,
  scenario = "rcp45"
)

## Load Pinus sylvestris data for 2095 (RCP 8.5)
psylvestris_2095_rcp85_sr <- fd_forest_eutrees4f(
  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.

Tip

The list of possible species is available through forestdata:

forestdata::eutrees4f_species
 [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:

reclassify_rcp <- function(raster, classes) {
  
  ## Reclassify 1 for 2
  raster_class <- ifel(
    raster == 1, 2, raster
  )
  ## Sum to current (2005) distribution
  raster_class <- as.factor(raster_class + psylvestris_2005_sr)
  ## 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
ps_classes <- c("Absent", "Will disappear", "Will appear",  "Present and stable")
ps_colors  <- c("#BE92A2", "#96031A", "#6DA34D", "#CEEDDB")
## Apply function
psylvestris_rcp45_sr <- reclassify_rcp(psylvestris_2095_rcp45_sr, ps_classes)
psylvestris_rcp85_sr <- reclassify_rcp(psylvestris_2095_rcp85_sr, ps_classes)

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.

Note

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:

map_rcp <- function(data, title = "(a) RCP 4.5") {
  
  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:

rcp45_gg <- map_rcp(psylvestris_rcp45_sr)
rcp85_gg <- map_rcp(psylvestris_rcp85_sr, title = "(b) RCP 8.5")

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
absent_txt   <- str_glue("<b style = 'color: {ps_colors[1]};'>almost absent</b>")
decrease_txt <- str_glue("<b style = 'color: {ps_colors[2]};'>decrease</b>")
increase_txt <- str_glue("<b style = 'color: {ps_colors[3]};'>shift</b>")
present_txt  <- str_glue("<b style = 'color: {ps_colors[4]};'>present</b>")

## Plot title
title_txt <- str_glue(
  "*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.

References

Mauri, Achille, Marco Girardello, Giovanni Strona, et al. 2022. “EU-Trees4F, a Dataset on the Future Distribution of European Tree Species.” Scientific Data 9 (1): 37. https://doi.org/10.1038/s41597-022-01128-5.