Spanish Burned Area during 2006-2015 by Municipality

Creating a 3D visualization of the burned area in Spain by municipality during the period 2006-2015 using {rayshader}.
R
forestry
wildfires
ggplot2
rayshader
Author

Adrián Cidre

Published

December 31, 2023

1 Introduction

In this post, we will explore the distribution of the burned area in Spain by municipality creating a 3D visualization using the awesome package rayshader 😎.

Honestly, I am not fan of 3D plots in general, because they are usually horrible and deceiving when trying to communicate something. However, with rayshader we can create spikes maps which can somehow look great and communicating the message properly.

Important

If you can, avoid always to create a 3D plot (3D pie charts (even 2D, please), 3D bar charts…).

The data was obtained from the MITECO (Ministerio para la Transición Ecológica y el Reto Demográfico), and it can be found in this url.

I will use the data for the decade 2006-2015, which is the latest available data for this purpose.

2 Loading packages

As usually, we will start loading all the necessary packages by using the pacman package.

library(pacman)

p_load(
  ## Data manipulation
  tidyverse, readxl,
  
  ## Spatial data manipulation
  sf, 
  
  ## Download data
  mapSpain,
  
  ## Visualization
  rayshader
)

We will use the tidyverse for data wrangling and visualization, the readxl package for reading our data in Excel format, the sf package for vectorial data manipulation, the package mapSpain for downloading some convenient spatial data, and finally rayshader for creating 3D visualizations.

3 Load the data

The data is the same I used in the previous post, but this time we will work with the municipalities instead of autonomous communities, and with burned area instead number of wildfires. Nevertheless, the first step will be to load the data, and fix the table names so we can all understand what’s going on with this dataset:

## Load data
wildfires_tbl <- read_excel("wildfires_2006_2015.xlsx", sheet = 3) 

## Set names
wildfires_tbl <- wildfires_tbl %>% 
  select(
    mun = NOMBRE, code = CODIGOINE, woodland_area = ARBOLADO,
    treeless_area = NOARBOLADO, total_area = TOTAL
  )

glimpse(wildfires_tbl)
Rows: 4,376
Columns: 5
$ mun           <chr> "A Arnoia", "A Baña", "A Bola", "A Cañiza", "A Capela", …
$ code          <dbl> 32003, 15007, 32014, 36009, 15018, 15030, 36017, 27018, …
$ woodland_area <dbl> 30, 615, 16, 978, 394, 18, 262, 600, 8, 803, 30, 364, 60…
$ treeless_area <dbl> 8, 82, 61, 1260, 165, 104, 791, 641, 17, 3395, 6, 1558, …
$ total_area    <dbl> 38, 697, 77, 2238, 559, 122, 1053, 1241, 25, 4198, 36, 1…

I selected a total of 5 columns: the first two columns corresponds to the name and ID of the municipalities, and the remaining columns refer to the burned area (woodland, treeless, and total burned area).

4 Prepare data

The data preparation is not very demanding to be honest. We have to download the map of Spain by municipalities, and for doing so, we have the convenient function esp_get_munic() from the mapSpain package:

## Get data
esp_mun_sf <- esp_get_munic() %>% 
  filter(!(ine.ccaa.name %in% c("Melilla", "Ceuta"))) %>% 
  mutate(code = as.numeric(LAU_CODE))

## Explore data
head(esp_mun_sf)
Simple feature collection with 6 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -3.14019 ymin: 36.73845 xmax: -2.05701 ymax: 37.54576
Geodetic CRS:  ETRS89
    codauto ine.ccaa.name cpro ine.prov.name cmun      name LAU_CODE
382      01     Andalucía   04       Almería  001      Abla    04001
379      01     Andalucía   04       Almería  002  Abrucena    04002
374      01     Andalucía   04       Almería  003      Adra    04003
375      01     Andalucía   04       Almería  004 Albánchez    04004
358      01     Andalucía   04       Almería  005 Alboloduy    04005
373      01     Andalucía   04       Almería  006     Albox    04006
                          geometry code
382 POLYGON ((-2.77744 37.23836... 4001
379 POLYGON ((-2.88984 37.09213... 4002
374 POLYGON ((-2.93161 36.75079... 4003
375 POLYGON ((-2.13138 37.29959... 4004
358 POLYGON ((-2.70077 37.09674... 4005
373 POLYGON ((-2.15335 37.54576... 4006

Here, we got the Spanish map as a simple feature object from the sf package, we eliminated Ceuta and Melilla (no information about wildfires), and converted the LAU_CODE column from character to numeric.

The last step before the visualization is to join the wildfires table with the municipalities, so we have the information of the wildfires stored in a spatial object. This is pretty straightforward by using the left_join function:

esp_mun_wildfires_sf <- esp_mun_sf %>% 
  left_join(
    wildfires_tbl,
    by = join_by(code == code)
  )

This will generate many NA values because in the wildfires_tbl we have a total of 4376, while in esp_mun_sf we have a total of 8129. This is because not every municipality of Spain had suffered wildfires. These NA values can be then changed by the value of 0:

esp_mun_wildfires_sf <- esp_mun_wildfires_sf %>% 
  mutate(
    total_area = if_else(is.na(total_area), 0, total_area)
  )

5 2D visualization

In this post we will learn how to turn a basic ggplot2 plot into a basic 3D visualization using rayshader. Therefore, the first step is to create a ggplot2 object:

## Create ggplot2 object
esp_mun_wildfires_gg <- esp_mun_wildfires_sf %>% 
  ggplot() +
  # Geometries
  geom_sf(aes(fill = total_area), color = NA) +
  # Scales
  scale_fill_distiller(palette = "Reds", direction = 1) +
  # Labels
  labs(title    = "Burned Area in Spain during 2006-2015",
       caption  = "@2023 Adrián Cidre (https://adrian-cidre.com)\nData: MITECO",
       fill     = "Area (ha)") +
  # Theme
  theme_void(base_family = "Rockwell") +
  theme(
    plot.title = element_text(
      size   = 14,
      face   = "bold", 
      hjust  = .95, 
      colour = "indianred4",
      margin = margin(b = 23, t = 5)
    ),
    plot.caption      = element_text(
      size   = 8, 
      colour = "indianred4",
      hjust  = .5
    ),
    legend.title      = element_text(size = 8, colour = "indianred4"),
    legend.text       = element_text(size = 8, colour = "indianred4"),
    legend.position   = "bottom",
    legend.key.width  = unit(1, "cm"),
    legend.key.height = unit(0.3, "cm"),
    plot.background   = element_rect(fill = "#B0C4DE", colour = NA)
  )

## Display plot
esp_mun_wildfires_gg
Figure 1: Distribution of total burned area in Spain by Municipality during the period 2006-2015

6 3D Visualization

Once we have this visualization, we can use the rayshader::plot_gg function. This function has many arguments to control the appearance of our map. In this tutorial I will show you the basic ones:

  • Multicore = TRUE: it will use multiple cores to compute the shadow matrix. This is computationally demanding, so this can improve the performance.

  • Width and height: they control the size of the ggplot2 object.

  • Scale: the vertical scaling. The higher this number, the taller the spikes.

  • Shadow intensity: a number between 0 and 1, representing the intensity of the shadows. When set to 0 the shadows are completely black, while 1 means no shadows.

  • Sunangle: the angle from where the sun would generate the shadows. The default 315, means that the sun is shining from the North West.

  • Offset edges: when TRUE, adds some space between polygons.

  • Window size: the size of the RGL device (preview window) when executing the function.

  • Zoom: amount of zoom to the plot.

  • Phi (\(\phi\)): azimuth angle (Figure 2).

  • Theta (\(\theta\)): rotation around z-axis (Figure 2).

Figure 2: Phi and theta angles

With all these parameters set, we can execute the function plot_gg() which will open a new window to visualize the rendered plot:

plot_gg(
  esp_mun_wildfires_gg,
  multicore        = TRUE,
  width            = 5,
  height           = 5,
  scale            = 150,
  shadow_intensity = .6,
  sunangle         = 315,
  offset_edges     = FALSE,
  windowsize       = c(600, 600),
  zoom             = .5, 
  phi              = 35, 
  theta            = -30
)

After we set the parameters that makes us happy, we can export the last rendered object by using the rayshader::render_snapshot() function:

render_snapshot(filename = "wildfires_ha_muni.png", clear = FALSE)

The result can be seen in Figure 3. As we saw in my previous post, the area that is most affected by the wildfires is by far the northwestern area of Spain, both in number of fires and burned hectares.

Figure 3: Burned area in Spain during 2006-2016 by municipality. The tallest spikes represent higher burned area.