Introduction

This is a proof of concept script intended to demonstrate and document the process of calculating the species sensitivity categories from the input data, including modelling outputs (predicted potential habitat), important point layers (Eurobats points, Golden Eagle nests, etc.), polygon layers (protected areas, EU Natura 2000 sites, etc.).

The data included is for demonstration purposes only and includes preliminary model outputs for six species of bat and a made up point dataset standing in for Eurobats sites.

Setup

Required packages

  • tidyverse 2.0.0
  • terra 1.7-78
  • sf 1.0-16

CRS definition

The coordinate reference system used for the analysis is the Croatian Terrestrial Reference System (HTRS96, EPSG:3765).

# install and load required packages
# install.packages(c("tidyverse", "terra", "sf", "geodata"))

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(terra)
## terra 1.7.78
## 
## Attaching package: 'terra'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
library(geodata)

# set working directory

setwd("C:/Users/mate.zec/Documents/Code/hrvatska-raa/")

# define crs

wgs84 <- crs("EPSG:4326")
htrs96 <- crs("EPSG:3765")

Reference grid definition

The reference grid for this analysis is defined as a 100-meter resolution raster in the HTRS96 CRS with coordinates divisible by 100 and extending 500 meters in each direction beyond the land borders of Croatia.

The Croatian land border was obtained from GADM (level 0) and projected / transformed to HTRS96:

# get country borders from GADM, transform to HTRS96 reference system

cro_border_wgs <- gadm(country = "HRV", level = 0, path = "input/geodata/")
cro_border <- cro_border_wgs %>% 
  project(y = htrs96)

plot(cro_border)

An extent for the raster is defined from the borders of this polygon, rounded to the nearest 100 m and extended by 500 m in every direction:

# get exact extent from border vector, round to outer 100 m coordinates
cro_extent_exact <- ext(cro_border)
cro_extent_rounded <- ext(c(round(cro_extent_exact$xmin, digits = -2) - 500,
                            round(cro_extent_exact$xmax, digits = -2) + 500,
                            round(cro_extent_exact$ymin, digits = -2) - 500,
                            round(cro_extent_exact$ymax, digits = -2) + 500))

# define empty reference grid / raster from rounded extent in HTRS96 
# reference system at a 100 m x 100 m spatial resolution

ref_grid <- rast(extent = cro_extent_rounded, 
                 resolution = c(100, 100), 
                 crs = htrs96)
bkg <- rasterize(cro_border, ref_grid, field = 0)
{
  plot(bkg)
  lines(cro_border)
}

Test: Point data

A randomly sampled dataset is standing in for the Eurobats sites in this notebook, defined using the spatSample() function:

# test: points to buffers, rasterizing buffers

pts <- spatSample(cro_border, size = 50, method = "random")
{ 
  plot(cro_border)
  points(pts)
}

These points are used to define buffers. In this case, two generously spaced buffers of 5 and 15 km were defined, to make them stand out on the map, but an arbitrary number of buffer zones can be defined and used this way.

buffer <- c("inner" = 5000, "outer" = 15000)

{plot(cro_border)
  points(pts, col = "red")
  lines(buffer(pts, buffer[['inner']]), col = "blue")
  lines(buffer(pts, buffer[['outer']]), col = "green")}

When rasterizing, the function to convert from vector attributes (the sensitivity score) can be chosen - in this case, the sensitivity scores for inner and outer buffers are set to 4 and 3, respectively, and the maximum value is taken if features overlap when rasterizing:

sensit <- c("inner" = 4, "outer" = 3)

rast_buf <- max(bkg, 
                rasterize(buffer(pts, buffer[['inner']]), 
                          ref_grid, 
                          field = sensit[['inner']], background = 0, fun = "max"),
                rasterize(buffer(pts, buffer[['outer']]), 
                          ref_grid, 
                          field = sensit[['outer']], background = 0, fun = "max"))
{ 
  plot(rast_buf)
  points(pts)
  lines(cro_border)
}

Test: model outputs

For one species (Hypsugo savii in the example), model output raster (30 m precision using the European Terrestrial Reference Grid) was reprojected to the reference raster. Note that the nearest neighbor shouldn’t be used for reprojection and resampling of continuous values, but is done here to speed up the analysis:

hsav <- rast("input/models/bats/Hypsugo savii/sdm_estimates.tif") %>% 
  project(y = htrs96, method = "near")
## |---------|---------|---------|---------|=========================================                                          
hsav_res <- resample(hsav, ref_grid, method = "near")
## |---------|---------|---------|---------|=========================================                                          
buf <- rasterize(buffer(pts, 15000), 
                 ref_grid, 
                 field = 0.4, background = 0, fun = "max") %>% 
  mask(bkg)

plot(hsav_res$X1 > 0.3 | buf > 0.1)

Test: aggregating outputs for different species

For this test, multiple model outputs were loaded, reprojected, resampled and combined into a single spatRaster object with bands/layers named after the species names. In this case, the continuous probability layer is resampled using the bilinear method, while the probability classes are resampled using the nearest neighbor method. This will be unnecessary in the actual sensitivity calculation because the model outputs will already be on the same HTRS96 100 meter resolution grid:

## [1] "Barbastella barbastellus"
## [1] "Eptesicus serotinus"
## [1] "Hypsugo savii"
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          [1] "Rhinolophus ferrumequinum"
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          [1] "Rhinolophus hipposideros"
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          [1] "Tadarida teniotis"
## |---------|---------|---------|---------|=========================================                                          |---------|---------|---------|---------|=========================================                                          
## |---------|---------|---------|---------|=========================================