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.
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
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]]
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)
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()
The Belgian regions are available by package
BelgiumMaps.StatBel
:
data(BE_ADMIN_REGION)
BE_ADMIN_REGION <- st_as_sf(BE_ADMIN_REGION)
In this section we add three columns to know whether a protected area belongs to one or more Belgian regions.
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)
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)
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)
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)
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()
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))