Density map of wolf in the Iberian Peninsula

1 Introduction

Today we will see how to create a map of the density of wolves’ observations in the Iberian Peninsula during the last decade (2014-2023).

We will use data from the Global Biodiversity Information Facility (GBIF) through the {geodata} package, and we will create a map of the Kernel Density of the observations using {ggplot2}.

2 Loading packages

We will use the following packages:

## Load pacman
library(pacman)

## Load rest of the packages
p_load(
  geodata, ggnewscale, giscoR, rayshader, sf, terra, tidyterra, tidyverse
)

{geodata} (Hijmans et al. 2024): we will use it to download data from the GBIF database, and the Digital Elevation Model (DEM).
{ggnewscale} (Campitelli 2024): to create a new scale in {ggplot2} of the same aesthetic.
{giscoR} (Hernangómez 2024): to download the boundaries of the study area.
{rayshader} (Morgan-Wall 2024): to create a 3D map.
{sf} (Pebesma and Bivand 2023): to manipulate vectorial data.
{terra} (Hijmans 2024): to manipulate raster data.
{tidyterra} (Hernangómez 2023): to use a nice colour palette for the DEM.
{tidyverse} (Wickham et al. 2019): to manipulate data in general.

3 Load the data

In this exercise we need to download three sources of data:

Study area
Wolf data
Digital Elevation Model

3.1 Study area

The study area for this exercise is the Iberian Peninsula, which includes the main land of Spain and Portugal without islands. To extract this area, we need first to get the boundaries of Spain and Portugal, eliminate the islands, and combine both countries. So, one step at a time. We start downloading the country of Portugal, and casting the MULTIPOLYGON to POLYGON:

## Get Portugal
portugal_sf <- gisco_get_countries(
  resolution = "01",
  country    = "Portugal"
) |> 
  st_cast("POLYGON")

## Visualize
plot(st_geometry(portugal_sf))

To filter the main land, we can calculate the area of each polygon, and slice the polygon with the maximum value of area:

## Get main land
portugal_sf <- portugal_sf |> 
  mutate(
    area = st_area(portugal_sf)
  ) |> 
  slice_max(area)

## Visualize
plot(st_geometry(portugal_sf))

Now, we can do the same for Spain:

## Get Spain
spain_sf <- gisco_get_countries(
  resolution = "01",
  country    = "Spain"
) |> 
  st_cast("POLYGON")

## Get main land
spain_sf <- spain_sf |> 
  mutate(
    area = st_area(spain_sf)
  ) |> 
  slice_max(area)

## Visualize
plot(st_geometry(spain_sf))

Up to this point, we have two different {sf} objects with the main land of Portugal Figure 2 and Spain Figure 3. The final step, is to combine them to form only one polygon:

## Make union
iberia_sf <- st_union(
  spain_sf,
  portugal_sf
) |> 
  st_transform(25830)

## Visualize
plot(st_geometry(iberia_sf))

Note that we also transformed the Coordinates Reference System to a projected system (EPSG 25830).

3.2 Wolf data

In this exercise, we will study the distribution of the gray wolf (Canis lupus; Figure 5) in the Iberian Peninsula:

As a proxy of the distribution, we will use data from GBIF for the period 2014-2023. This database contains information about observations of any species in the world. We will use the sp_occurrence() function from the {geodata} package, and we will add two arguments to specify the countries (ES for Spain, PT for Portugal) and the period (2014-2023).

Tip

To specify a period, separate two years by a comma.

## Download wolf data
## - Years: 2014-2023
wolf_lst <- map(
  .x = c("ES", "PT"),
  .f = \(x) sp_occurrence(
    genus   = "Canis",
    species = "lupus",
    args    = c(
      paste0("country=", x),
      "year=2014,2023"
    )
  )
)

Note that we use the purrr::map() function to iterate over two elements (ES and PO) which are used in the country argument.

3.3 Elevation

The elevation data can be downloaded using the {geodata} package. We use the elevation_global() function to download the DEM for all the world at 0.5 minutes of a degree spatial resolution. Immediately, we use the terra::crop() function to crop and mask only the Iberian Peninsula (note that the data must be transformed before cropping to the same CRS), and then we transform the CRS to the same projected CRS as the Iberian Peninsula.

## Download global DEM, crop it to Iberia, and project it
dem_sr <- elevation_global(res = .5, path = tempdir()) |> 
  crop(
    st_transform(iberia_sf, 4326),
    mask = TRUE
  ) |> 
  project("EPSG:25830")
## Print DEM
print(dem_sr)

class       : SpatRaster 
dimensions  : 1142, 1492, 1  (nrow, ncol, nlyr)
resolution  : 774.7848, 774.7848  (x, y)
extent      : -86157.96, 1069821, 3984505, 4869310  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 30N (EPSG:25830) 
source(s)   : memory
name        : wc2.1_30s_elev 
min value   :      -4.587787 
max value   :    3250.524902

4 Prepare the data

We did already most of the data preparation steps previously, but there are a couple of things that I wanted to leave for this section.

4.1 Wolf data

We are using a projected CRS, but the coordinates data that we have in wolf_lst are in the WGS 84. Therefore, we need to transform this list to a sf object, transform the coordinates, and convert it back to a tibble with the coordinates in two columns. We could use the geom_sf(), but it might give us some issues with rayshader so we will use geom_point().

Another thing to note is that we have data from the entire countries, but we want only from the Iberian Peninsula. In this sense, we can use st_intersection() to filter only the observations in the Iberian Peninsula.

## Bind data, convert to sf, transform and select points within Iberia
wolf_sf <- list_rbind(wolf_lst) |> 
  st_as_sf(
    coords = c("lon", "lat"),
    crs    = 4326
  ) |> 
  st_transform(25830) |> 
  st_intersection(iberia_sf) |> 
  select(country) 

## Visualize
plot(st_geometry(wolf_sf))
plot(st_geometry(iberia_sf), add = TRUE)

Figure 6: Gray wolf observations in the Iberian Peninsula

## Get data as tibble with coordinates
wolf_tbl <- wolf_sf |> 
  mutate(
    x = st_coordinates(wolf_sf)[, 1],
    y = st_coordinates(wolf_sf)[, 2]
  ) |> 
  as_tibble()

## Data head
head(wolf_tbl)

# A tibble: 6 × 4
  country           geometry       x        y
  <chr>          <POINT [m]>   <dbl>    <dbl>
1 Spain   (602078.7 4335512) 602079. 4335512.
2 Spain   (519809.4 4060480) 519809. 4060480.
3 Spain   (193188.3 4115153) 193188. 4115153.
4 Spain   (428693.5 4560619) 428693. 4560619.
5 Spain   (279953.7 4815493) 279954. 4815493.
6 Spain     (227613 4157558) 227613. 4157558.

4.2 Elevation

We can use the {tidyterra} package to add the DEM to a {ggplot2} object. However, this would give us some problems when translating to the {rayshader} package. Because of this, we will use the geom_raster() function instead, and for that, we need a tibble with the coordinates and the values:

## Convert to tibble 
dem_tbl <- dem_sr |>   
  as_tibble(xy = TRUE) |>   
  na.omit() |>    
  rename(elevation = wc2.1_30s_elev)
## Print head 
head(dem_tbl)

# A tibble: 6 × 3
        x        y elevation
    <dbl>    <dbl>     <dbl>
1 122647. 4858850.      12  
2 107926. 4858075.      51  
3 123421. 4858075.      91.7
4 124196. 4858075.     111  
5 105601. 4857300.     173. 
6 106376. 4857300.     166.

In the previous code, we eliminated the NA values (pixels in the bounding box outside the boundaries of the Iberian Peninsula), and renamed the elevation value to a more intuitive name.

5 Visualization

5.1 2D visualization

We will start creating a 2D visualization using {ggplot2}. You can see a better explanation of the code in the video, but one thing to note is that I am using the new_scale_fill() function because we have two different geometries (geom_raster and stat_density_2d) using the fill aesthetic. Therefore, if we want to modify them differently, we need to treat them in different scale_fill_* functions. The {ggnewscale} package help us to do this very easily, which is basically creating the first layer and its scale, then adding the new_scale_* function, and then adding the next geometry and its scale.

## Create ggplot2 map
map <- wolf_tbl |> 
  ggplot() +
  geom_raster(
    data = dem_tbl,
    aes(x, y, fill = elevation)
  ) +
  scale_fill_gradientn(
    colours = whitebox.colors(20),
    guide   = guide_colourbar(
      title          = "Elevation (m)",
      title.hjust    = .5,
      position       = "inside",
      barwidth       = unit(4, "mm"),
      barheight      = unit(4, "cm")
    )
  ) +
  new_scale_fill() +
  stat_density_2d(
    aes(x = x, y = y, fill = after_stat(level)), 
    geom = "polygon", 
    alpha = .1,
    bins = 25,
    show.legend = FALSE
  ) +
  geom_point(
    aes(x, y),
    alpha = .3
  ) +
  scale_fill_gradientn(
    colours = c("darkblue", "blue", "green", "yellow", "red"),
    # colours = pal
  ) +
  labs(
    title   = "Where are wolves seen the most in the Iberian Peninsula?",
    caption = "Author: Adrián Cidre | Data source: GBIF (2013-2024)"
  ) +
  theme_void(
    base_size = 10
  ) +
  theme(
    text = element_text(color = "snow"),
    plot.title = element_text(
      family = "Merriweather",
      face   = "bold",
      size   = 12,
      hjust  = .5
    ),
    plot.caption = element_text(
      hjust = .5
    ),
    legend.position.inside = c(.85, .3)
  )
## Visualize
map

Figure 7: Density map of the distribution of the gray wolf

5.2 3D visualization

The {rayshader} package has a convenient function to transform a {ggplot2} object to a 3D graph. To maintain the dimensions of the Iberian Peninsula, we create two variables with the number of rows and columns of the DEM.

## Dimensions of the DEM
h <- nrow(dem_sr) 
w <- ncol(dem_sr) 
## Print
print(h)

[1] 1142

print(w)

[1] 1492

Now, we are ready to transform our map to 3D:

plot_gg(
  ggobj  = map,
  width  = w / 200,
  height = h / 200,
  scale  = 50,
  solid  = FALSE,
  zoom   = .8,
  phi    = 85,
  theta  = 0,
  shadow = FALSE
)

What this function does is to create a map with the specified with and height in inches, and the scale factor to create the 3D effect. It will open an rgl window that looks like Figure 8.

The final step is to render and export the map in high quality:

render_highquality(
  filename          = "07_wolf_density/wolf_density.png",
  width             = w,
  height            = h,
  preview           = TRUE,
  light             = TRUE,
  environment_light = "data/env_lights/misty_farm_road_4k.hdr",
  intensity_env     = 1,
  interactive       = FALSE
)

During the rendering you will see Figure 9. At the top of the figure there’s a red loading bar that shows you how the rendering process is going.

Figure 9: Rendering window of render_highquality

Caution

Sometimes this rendering function will give weird rendering results. This might be happening because we are using more than 1 fill aesthetic, and it might produce some errors.

It is not recommended to create 3D maps with more than 1 fill aesthetic, although I can also work as we can see in this example.

After some minutes, we can see the final result:

References

Campitelli, Elio. 2024. “Ggnewscale: Multiple Fill and Colour Scales in ’Ggplot2’.” https://CRAN.R-project.org/package=ggnewscale.

Hernangómez, Diego. 2023. “Using the Tidyverse with Terra Objects: The Tidyterra Package” 8: 5751. https://doi.org/10.21105/joss.05751.

———. 2024. “giscoR: Download Map Data from GISCO API - Eurostat.” https://doi.org/10.5281/zenodo.4317946.

Hijmans, Robert J. 2024. “Terra: Spatial Data Analysis.” https://CRAN.R-project.org/package=terra.

Hijmans, Robert J., Márcia Barbosa, Aniruddha Ghosh, and Alex Mandel. 2024. “Geodata: Download Geographic Data.” https://CRAN.R-project.org/package=geodata.

Morgan-Wall, Tyler. 2024. “Rayshader: Create Maps and Visualize Data in 2D and 3D.” https://CRAN.R-project.org/package=rayshader.

Pebesma, Edzer, and Roger Bivand. 2023. “Spatial Data Science: With Applications in r.” https://doi.org/10.1201/9780429459016.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse” 4: 1686. https://doi.org/10.21105/joss.01686.