Get cells of Belgian EEA reference grid containing protected areas

1 Goal

Given the EEA reference grid of Belgium at 1km resolution and the Belgian Natura2000 protected areas, this document assesses the cells covering the protected areas.

2 Setup

library(tidyverse)  # To do datascience
library(tidylog) # To provide feedback on dplyr functions
library(here) # To find files
library(sf) # To work with geospatial data
library(BelgiumMaps.StatBel)  # To load Belgian administrative boundaries
library(INBOtheme)  # To apply INBO style to plots
library(leaflet)  # To make interactive maps

3 Read input data

3.1 Import EEA reference grid

We import UTM grid data of Belgium at 1 by 1 km resolution. All grids have CRS (Coordinate Reference System) equal to Belge Lambert 1972:

be_grid <- st_read(here::here("data", "external", "utm1_bel"))

## Reading layer `be_1km' from data source 
##   `C:\Users\soria_delva\Documents\GitHub\indicators\data\external\utm1_bel' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 51726 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 3768000 ymin: 2926000 xmax: 4080000 ymax: 3237000
## Projected CRS: ETRS89 / ETRS-LAEA

be_grid_crs <- st_crs(be_grid)
be_grid_crs

## Coordinate Reference System:
##   User input: ETRS89 / ETRS-LAEA 
##   wkt:
## PROJCRS["ETRS89 / ETRS-LAEA",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["unnamed",
##         METHOD["Lambert Azimuthal Equal Area",
##             ID["EPSG",9820]],
##         PARAMETER["Latitude of natural origin",52,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",10,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",4321000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",3210000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["x",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["y",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     ID["EPSG",3035]]

3.2 Read Belgian protected areas Natura2000

GIS data of protected areas as defined in Natura2000 can be found on related webpage of European Environment Agency. We downloaded the GIS data as shapefile in folder data/external/Natura2000_end2019_Shapefile:

folder_path <- here::here(
  "data",
  "external",
  "Natura2000_end2019_Shapefile")

dir.create(folder_path)


path_natura2000 <- here::here(
  "data",
  "external",
  "Natura2000_end2019_Shapefile",
  "Natura2000_end2021_rev1.gpkg"
)

if (!file.exists(path_natura2000)) {
  file.create(path_natura2000, recursive = TRUE)
  natura2000_file <- file(path_natura2000, "wb")
  response <- 
    httr2::request("https://sdi.eea.europa.eu/datashare/s/fAL7bpzZGHjiFL4/download") %>% 
    httr2::req_url_query(
      path = "/GPKG",
      files = "Natura2000_end2021_rev1.gpkg"
    ) %>% 
    httr2::req_perform()
  # write to disk
  httr2::resp_body_raw(response) %>% 
    writeBin(con = natura2000_file)
  close(natura2000_file)
}

We read the shapefile:

protected_areas <- st_read(path_natura2000)

## Multiple layers are present in data source C:\Users\soria_delva\Documents\GitHub\indicators\data\external\Natura2000_end2019_Shapefile\Natura2000_end2021_rev1.gpkg, reading layer `NaturaSite_polygon'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Reading layer `NaturaSite_polygon' from data source 
##   `C:\Users\soria_delva\Documents\GitHub\indicators\data\external\Natura2000_end2019_Shapefile\Natura2000_end2021_rev1.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 27020 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 746551.8 ymin: 905053.4 xmax: 6506127 ymax: 5298655
## Projected CRS: ETRS89-extended / LAEA Europe

protected_areas_crs <- st_crs(protected_areas)
protected_areas_crs

## Coordinate Reference System:
##   User input: ETRS89-extended / LAEA Europe 
##   wkt:
## PROJCRS["ETRS89-extended / LAEA Europe",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["Europe Equal Area 2001",
##         METHOD["Lambert Azimuthal Equal Area",
##             ID["EPSG",9820]],
##         PARAMETER["Latitude of natural origin",52,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",10,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",4321000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",3210000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (Y)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (X)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Statistical analysis."],
##         AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
##         BBOX[24.6,-35.58,84.73,44.83]],
##     ID["EPSG",3035]]

Convert CRS of protected areas:

protected_areas <- st_transform(protected_areas, crs = be_grid_crs)

3.2.1 Select Belgian protected areas

We are interested in protected areas in Belgium:

protected_areas <-
  protected_areas %>%
  tidylog::filter(MS == "BE")

Summary by site type:

protected_areas %>%
  as_tibble() %>%
  tidylog::group_by(SITETYPE) %>%
  tidylog::count()

3.3 Read Belgian regions

The Belgian regions are available by package BelgiumMaps.StatBel:

data(BE_ADMIN_REGION)
BE_ADMIN_REGION <- st_as_sf(BE_ADMIN_REGION)

3.4 Add regional information to Belgian protected areas

In this section we add three columns to know whether a protected area belongs to one or more Belgian regions.

3.4.1 Flemish protected areas

Add column flanders to indicate whether the protected area is a Flemish protected area (TRUE/FALSE):

protected_areas_flanders <- tibble(
  SITECODE = c(
    "BE2200043",
    "BE2200032",
    "BE2100017",
    "BE2300006",
    "BE2500121",
    "BE2100024",
    "BE2400008",
    "BE2200030",
    "BE2200037",
    "BE2200041",
    "BE2200042",
    "BE2200038",
    "BE2101437",
    "BE2219312",
    "BE2500831",
    "BE2100323",
    "BE2100045",
    "BE2100020",
    "BE2200028",
    "BE2400009",
    "BE2200033",
    "BE2300044",
    "BE2400012",
    "BE2200626",
    "BE2200727",
    "BE2200029",
    "BE2400010",
    "BE2223316",
    "BE2301235",
    "BE2301336",
    "BE2300007",
    "BE2422315",
    "BE2524317",
    "BE2100026",
    "BE2100040",
    "BE2200031",
    "BE2500003",
    "BE2200034",
    "BE2100424",
    "BE2400014",
    "BE2101538",
    "BE2200039",
    "BE2500932",
    "BE2500001",
    "BE2218311",
    "BE2101639",
    "BE2217310",
    "BE2300222",
    "BE2400011",
    "BE2500002",
    "BE2220313",
    "BE2200035",
    "BE2100015",
    "BE2300005",
    "BE2301134",
    "BE2100016",
    "BE2100019",
    "BE2200036",
    "BE2500004",
    "BE2200525",
    "BE2501033",
    "BE2221314"
  )) %>%
  tidylog::mutate(flanders = TRUE,
         wallonia = FALSE,
         brussels = FALSE)

3.4.2 Brussels protected areas

Add column brussels to indicate whether the protected area is located in Brussels region (TRUE/FALSE):

protected_areas_brussels <- tibble(
  SITECODE = c("BE1000001",
               "BE1000002",
               "BE1000003")) %>%
  tidylog::mutate(flanders = FALSE,
         wallonia = FALSE,
         brussels = TRUE)

3.4.3 Protected areas under federal administration

Areas not belonging to any region as they are under federal administration:

protected_areas_federal <- tibble(
  SITECODE = c("BEMNZ0001",
               "BEMNZ0002",
               "BEMNZ0003",
               "BEMNZ0004",
               "BEMNZ0005")) %>%
  tidylog::mutate(flanders = FALSE,
         wallonia = FALSE,
         brussels = FALSE)

3.4.4 Wallonia protected areas

All the other protected areas are in Wallonia (wallonia = TRUE):

protected_areas_wallonia <- tibble(
  SITECODE = as.character(protected_areas$SITECODE[
    !protected_areas$SITECODE %in% c(protected_areas_flanders$SITECODE,
                                     protected_areas_brussels$SITECODE,
                                     protected_areas_federal$SITECODE)
  ])) %>%
  tidylog::mutate(flanders = FALSE,
         wallonia = TRUE,
         brussels = FALSE)

3.4.5 Add regional information

We can now merge the regional information

protected_areas_regional_info <-
  bind_rows(protected_areas_flanders,
            protected_areas_wallonia,
            protected_areas_brussels,
            protected_areas_federal) %>%
  tidylog::mutate(SITECODE = as.factor(SITECODE))

to add it to protected_areas:

protected_areas <- 
  protected_areas %>%
  tidylog::left_join(protected_areas_regional_info,
            by = "SITECODE")

Preview:

protected_areas %>%
  as_tibble() %>%
  tidylog::select(SITECODE, SITENAME, SITETYPE, flanders, wallonia, brussels) %>%
  head()

3.4.6 Distribution of Belgian protected areas

Number of protected areas per type:

protected_areas %>%
  st_drop_geometry() %>%
  tidylog::group_by(SITETYPE) %>%
  tidylog::count()

In Flanders:

protected_areas %>%
  tidylog::filter(flanders == TRUE) %>%
  st_drop_geometry() %>%
  tidylog::group_by(SITETYPE) %>%
  tidylog::count()

# Transform to wgs84
protected_areas_wgs84 <- 
  protected_areas %>%
          st_transform(crs = 4326)
# Flemish protected areas
protected_areas_wgs84_flanders <- 
  protected_areas_wgs84 %>%
  tidylog::filter(flanders == TRUE)

prot_area_fl_palette <- colorFactor(
  topo.colors(nrow(protected_areas_flanders)),
  protected_areas_flanders$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_flanders %>%
                tidylog::filter(SITETYPE == "A"),
              fillColor = ~prot_area_fl_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type A") %>%
  addPolygons(data = protected_areas_wgs84_flanders %>%
                tidylog::filter(SITETYPE == "B"),
              fillColor = ~prot_area_fl_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addLayersControl(
    overlayGroups = c("type A", "type B"),
    options = layersControlOptions(collapsed = FALSE))

In Wallonia:

protected_areas %>%
  tidylog::filter(wallonia == TRUE) %>%
  st_drop_geometry() %>%
  tidylog::group_by(SITETYPE) %>%
  tidylog::count()

protected_areas_wgs84_wallonia <- 
  protected_areas_wgs84 %>%
  tidylog::filter(wallonia == TRUE)

prot_area_wallonia_palette <- colorFactor(
  topo.colors(nrow(protected_areas_wallonia)),
  protected_areas_wallonia$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_wallonia %>%
                tidylog::filter(SITETYPE == "A"),
              fillColor = ~prot_area_wallonia_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type A") %>%
  addPolygons(data = protected_areas_wgs84_wallonia %>%
                tidylog::filter(SITETYPE == "B"),
              fillColor = ~prot_area_wallonia_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addPolygons(data = protected_areas_wgs84_wallonia %>%
                tidylog::filter(SITETYPE == "C"),
              fillColor = ~prot_area_wallonia_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type C") %>%
  addLayersControl(
    overlayGroups = c("type A", "type B", "type C"),
    options = layersControlOptions(collapsed = FALSE))

In Brussels:

protected_areas %>%
  tidylog::filter(brussels == TRUE) %>%
  st_drop_geometry() %>%
  tidylog::group_by(SITETYPE) %>%
  tidylog::count()

protected_areas_wgs84_brussels <- 
  protected_areas_wgs84 %>%
  tidylog::filter(brussels == TRUE)

prot_area_brussels_palette <- colorFactor(
  topo.colors(nrow(protected_areas_brussels)),
  protected_areas_brussels$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_brussels %>%
                tidylog::filter(SITETYPE == "B"),
              fillColor = ~prot_area_brussels_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addLayersControl(
    overlayGroups = c("type B"),
    options = layersControlOptions(collapsed = FALSE))

Areas under federal administration

protected_areas %>%
  tidylog::filter(flanders == FALSE,
         wallonia == FALSE,
         brussels == FALSE) %>%
  st_drop_geometry() %>%
  tidylog::group_by(SITETYPE) %>%
  tidylog::count()

# Protected areas under federal administration
protected_areas_wgs84_federal <- 
  protected_areas_wgs84 %>%
  tidylog::filter(flanders == FALSE,
         wallonia == FALSE,
         brussels == FALSE)

prot_area_palette_federal <- colorFactor(
  topo.colors(nrow(protected_areas_federal)),
  protected_areas_federal$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_federal %>%
                tidylog::filter(SITETYPE == "A"),
              fillColor = ~prot_area_palette_federal(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type A") %>%
  addPolygons(data = protected_areas_wgs84_federal %>%
                tidylog::filter(SITETYPE == "B"),
              fillColor = ~prot_area_palette_federal(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addLayersControl(
    overlayGroups = c("type A", "type B"),
    options = layersControlOptions(collapsed = FALSE))

4 Define Bird and Habitat directive areas

4.1 Bird directive areas

Bird directive areas (SPA areas) are defined as the areas with SITETYPE: A and C.

protected_areas_bird <-
  protected_areas %>%
  tidylog::filter(SITETYPE %in% c("A", "C"))

4.2 Habitat directive areas

Habitat directive areas (SAC areas) are defined as the protected areas with SITETYPE values B and C.

protected_areas_habitat <-
  protected_areas %>%
  tidylog::filter(SITETYPE %in% c("B", "C"))

5 Intersect grid of Belgium with protected areas

5.1 Interset grid with protected areas

Intersect cells of the grid and protected areas:

grid_intersects_areas <- st_intersects(be_grid, protected_areas,
  sparse = FALSE
)
grid_intersects_areas <- as_tibble(grid_intersects_areas,
  .name_repair = "universal"
)
# Assign column names based on protected areas code
prot_areas_sitecode <- as.character(protected_areas$SITECODE)
names(grid_intersects_areas) <- prot_areas_sitecode
# Add 1x1km grid cellcode as new column
grid_intersects_areas <-
  grid_intersects_areas %>%
 tidylog::mutate(CELLCODE = be_grid$CELLCODE) %>%
  tidylog::select(CELLCODE, everything())

5.2 Natura2000 protected areas

Select cells which intersect any of the protected areas:

cells_protected <-
  grid_intersects_areas %>%
  tidylog::filter_at(vars(starts_with("BE")), any_vars(. == TRUE))

5.3 Get grid cells for each protected area

For each protected area we can get the list of intersecting cells:

cells_of_prot_areas <-
  map_dfr(
    prot_areas_sitecode,
    function(pa_cell_code) {
      site_type <-
        protected_areas %>%
        dplyr::filter(SITECODE == pa_cell_code) %>%
        dplyr::pull(SITETYPE)
      cells <-
        cells_protected %>%
        # avoid long logs by using dplyr functions instead of tidylog
        dplyr::select(CELLCODE, !!sym(pa_cell_code)) %>%
        dplyr::filter(!!sym(pa_cell_code) == TRUE) %>%
        dplyr::mutate(SITECODE = pa_cell_code,
                      SITETYPE = site_type) %>%
        dplyr::select(SITECODE, SITETYPE, CELLCODE)
    }
  )

For example, grid cells intersecting protected area BE1000003:

cells_of_prot_areas %>%
  tidylog::filter(SITECODE == "BE1000003")

Distribution of the number of cells covering the protected areas:

cells_of_prot_areas %>%
  tidylog::group_by(SITECODE) %>%
  tidylog::count() %>%
  arrange(desc(n)) %>%
  ggplot() +
  geom_histogram(aes(x = n), bins = 30) +
  scale_x_log10() +
  xlab("Number of cells (km2)") +
  ylab("Number of protected areas")

5.4 Natura2000

We can now define whether a cell of the reference grid intersects a protected area as defined by Natura2000. We add this information in column natura2000:

cellcode_protected_natura2000 <-
  cells_protected %>%
  as_tibble() %>%
  pull(CELLCODE)

be_grid_membership_protected_areas <-
  be_grid %>%
  tidylog::mutate(natura2000 = ifelse(CELLCODE %in% cellcode_protected_natura2000,
    TRUE,
    FALSE
  ))

5.4.1 Bird directive areas

We add also a column, spa, assessing whether the cell intersects a SPA area:

cells_protected_spa <-
  cells_protected %>%
  tidylog::select(
    CELLCODE,
    one_of(as.character(protected_areas_bird$SITECODE))
  ) %>%
  tidylog::filter_at(vars(starts_with("BE")), any_vars(. == TRUE)) %>%
  as_tibble() %>%
  pull(CELLCODE)

be_grid_membership_protected_areas <-
  be_grid_membership_protected_areas %>%
  tidylog::mutate(spa = ifelse(CELLCODE %in% cells_protected_spa,
    TRUE,
    FALSE
  ))

5.4.2 Habitat directive areas

We add also a column, habitat, assessing whether the cell intersects a habitat directive area:

cells_protected_hab <-
  cells_protected %>%
  tidylog::select(
    CELLCODE,
    one_of(as.character(protected_areas_habitat$SITECODE))
  ) %>%
  tidylog::filter_at(vars(starts_with("BE")), any_vars(. == TRUE)) %>%
  as_tibble() %>%
  pull(CELLCODE)

be_grid_membership_protected_areas <-
  be_grid_membership_protected_areas %>%
  tidylog::mutate(habitat = ifelse(CELLCODE %in% cells_protected_hab,
    TRUE,
    FALSE
  ))

5.5 Summary

How many cells contain/intersect protected areas?

be_grid_membership_protected_areas <-
  be_grid_membership_protected_areas %>%
  as_tibble()
be_grid_membership_protected_areas %>%
  tidylog::group_by(natura2000) %>%
  tidylog::count()

How many cells contain/intersect SPA areas?

be_grid_membership_protected_areas %>%
  tidylog::group_by(spa) %>%
  tidylog::count()

How many cells contain/intersect habitat directive areas?

be_grid_membership_protected_areas %>%
  tidylog::group_by(habitat) %>%
  tidylog::count()

5.6 Save data

We save be_grid_membership_protected_areas as tab-separated file (tsv):

write_tsv(be_grid_membership_protected_areas,
  path = here::here(
    "data",
    "interim",
    "intersect_EEA_ref_grid_protected_areas.tsv"
  ),
  na = ""
)

We also save cells_of_prot_areas in the same format:

write_tsv(cells_of_prot_areas,
  path = here::here(
    "data",
    "interim",
    "EEA_ref_grid_cells_covering_protected_areas.tsv"
  ),
  na = ""
)

We save the Belgian protected areas as a geopackage as it will be used elsewhere:

protected_areas %>%
  st_write(
    here::here("data",
               "output",
               "Belgian_Natura2000_protected_areas.gpkg"),
    delete_dsn =TRUE
)

## Deleting source `C:/Users/soria_delva/Documents/GitHub/indicators/data/output/Belgian_Natura2000_protected_areas.gpkg' using driver `GPKG'
## Writing layer `Belgian_Natura2000_protected_areas' to data source 
##   `C:/Users/soria_delva/Documents/GitHub/indicators/data/output/Belgian_Natura2000_protected_areas.gpkg' using driver `GPKG'
## Writing 310 features with 8 fields and geometry type Multi Polygon.

We are also interested in saving basic information about protected areas:

protected_areas_metadata <-
  protected_areas %>%
  st_drop_geometry() %>%
  tidylog::select(SITECODE,
         SITENAME,
         SITETYPE,
         flanders,
         wallonia,
         brussels)
protected_areas_metadata

We save these metadata:

protected_areas_metadata %>%
  write_tsv(
    path = here::here(
      "data",
      "interim",
      "protected_areas_metadata.tsv"
    ),
    na = ""
  )