BFN_data.R

# Only needed to be run one time to get bfn data
# can be adjusted to load data for other layers (nationalparks etc.)
library(xml2)
library(terra)
library(sf)
library(tidyverse)
library(gdalUtilities)

url <- "https://geodienste.bfn.de/ogc/wfs/schutzgebiet?REQUEST=GetCapabilities&SERVICE=WFS&VERSION=2.0.0"
wms_capabilities <- read_xml(url)

layers <- st_layers(url)

data_nasc <- read_sf(url, layer= "bfn_sch_Schutzgebiet:Naturschutzgebiete")

data_nasc <- data_nasc %>% mutate(value=1)

ggplot(data=data_nasc) +
  #geom_sf(data=germany_geo) +
  geom_sf(fill="darkgreen", col="darkgreen") 

#ggsave("Figures/bfn_map.png", width = 5, height = 4)


# Rasterize data  #

# Issue with conversion use function to convert data first before vectorizing
ensure_multipolygons <- function(X) {
  tmp1 <- tempfile(fileext = ".gpkg")
  tmp2 <- tempfile(fileext = ".gpkg")
  st_write(X, tmp1)
  ogr2ogr(tmp1, tmp2, f = "GPKG", nlt = "MULTIPOLYGON")
  Y <- st_read(tmp2)
  st_sf(st_drop_geometry(X), geom = st_geometry(Y))
}

data_nasc_con <- ensure_multipolygons(data_nasc)

data_nasc_con <- st_transform(data_nasc_con, crs="epsg:3035")

vec_nasc <- vect(data_nasc_con)

empty_raster <- rast(vec_nasc, resolution(1,1))
res(empty_raster) <- c(100, 100)

raster_nasc <- rasterize(vec_nasc, empty_raster, "value")

plot(raster_nasc, col="darkgreen")

saveRDS(raster_nasc, "data/gis/protected_bfn_raster.rds")

# Check if conversion and storage worked #
raster_bfn <- readRDS("data/gis/protected_bfn_raster.rds") # works