# INBO CODING CLUB

26 January 2026

Welcome!

## Share your code snippet

If you want to share your code snippet, copy paste your snippet within a section of three backticks (```):

As an **example**:

library(tidyverse)

(*you can copy paste this example and add your code further down*)

## Yellow sticky notes

No yellow sticky notes online. Put your name + " | " and add a "*" each time you solve a challenge (see below).

## Participants

Name | Challenges
--- | ---
Damiano Oldoni | ***
Nele Mullens |
Pieter Huybrechts |
Hanne Van Gompel |
Anne-Lie Van Praet |
Jolien Buyse |
Suzanna Lettens |**
Lawrence Whatley |
Adriaan | **

## Challenge 1

### Damiano's solution (example)

Copy paste this section to show your solutions.

# dummy code
print("This is how to insert code.")

### Pieter

### Raïsa

#Challenge 1.1
occs <- st_as_sf(occs_df,
                 coords = c("decimalLongitude", "decimalLatitude"),
                 crs = 4326)

#Challenge 1.2
sf::st_layers("data/20260126/20260126_protected_areas_be.gpkg")
#1st layer is multipolygon with known crs, other two have unknown geometry_type and crs.

#Challenge 1.3
pa <- st_read("data/20260126/20260126_protected_areas_be.gpkg", layer = "NaturaSite_polygon")

#Challenge 1.4
pa_bioregion <- st_read("data/20260126/20260126_protected_areas_be.gpkg", layer = "BIOREGION")
pa_habitats <- st_read("data/20260126/20260126_protected_areas_be.gpkg", layer = "HABITATS")
#they are read as simple data frames, not spatial data frames

#Challenge 1.5
municipalities <- st_read(municipalities_folder)

#Challenge 1.6
#unzipping the data put the files in another subfolder.
# I manually moved them out of the subfolder for this to work.
# Alternatively, we could have updated the aae_grid_folder: 
# eea_grid <- st_read(paste0(eea_grid_folder, "/eea_v_3035_10_km_eea-ref-grid-be_p_2013_v02_r00"))
eea_grid <- st_read(eea_grid_folder)

#Challenge 1.7
st_crs(eea_grid) #ETRS89 / ETRS-LAEA 
st_crs(pa) #ETRS89-extended / LAEA Europe
st_crs(occs) #WGS 84
st_crs(municipalities) #WGS 84

#Challenge 1.8
pa_a <- pa %>%
  dplyr::filter(SITETYPE == "A")

### Lawrence

#1
occs <- st_as_sf(occs_df, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
occs
mapview(occs)

#2
protected <- "data/20260126/20260126_protected_areas_be.gpkg"
st_layers(protected)

#3
pa <- st_read(protected, layer = "NaturaSite_polygon")
pa
mapview(pa)

#4
pa_bioregion <- st_read(protected, layer = "BIOREGION")
pa_habitats <- st_read(protected, layer = "HABITATS")

#5
municipalities <- st_read(municipalities_folder)
municipalities
mapview(municipalities)

#6
eea_grid <- st_read("data/20260126/20260126_eea_grid_be/eea_v_3035_10_km_eea-ref-grid-be_p_2013_v02_r00/be_10km.shp")
eea_grid
mapview(eea_grid)

#7
st_crs(occs) #WGS84
st_crs(pa) #ETRS89
st_crs(municipalities) #WGS84
st_crs(eea_grid) #ETRS89

#8
pa_a <- subset(pa, SITETYPE == "A")
pa_a
mapview(pa_a)

### Jolien

occs <- st_as_sf(occs_df, coords = c("decimalLongitude","decimalLatitude"), crs = 4326)

st_layers("data/20260126/20260126_protected_areas_be.gpkg.curltmp", )
pa <- st_read("data/20260126/20260126_protected_areas_be.gpkg.curltmp", "NaturaSite_polygon")
pa_bioregion <- st_read("data/20260126/20260126_protected_areas_be.gpkg.curltmp", "BIOREGION")
pa_habitats <- st_read("data/20260126/20260126_protected_areas_be.gpkg.curltmp", "HABITATS")

municipalities <- read_sf("data/20260126/20260126_municipalities_belgium/georef-belgium-municipality-millesime.shp")
eea_grid <- read_sf("data/20260126/20260126_eea_grid_be/eea_v_3035_10_km_eea-ref-grid-be_p_2013_v02_r00/be_10km.shp")

st_crs(eea_grid) # ETRS89 / ETRS-LAEA 
st_crs(municipalities) # WGS 84
st_crs(pa) # ETRS89-extended / LAEA Europe 
st_crs(occs) # WGS 84

pa_a <- pa %>%
  filter(SITETYPE == "A")
  

### Nele

#1.1
#use decimalLongitude and decimalLatitude to specify coordinates
#start from occs_df
occs <- st_as_sf(occs_df, 
                 coords = c("decimalLongitude", "decimalLatitude"), 
                 crs = 4326)

#1.2
prot_ar <- st_layers("data/20260126/20260126_protected_areas_be.gpkg")
prot_ar$name
prot_ar$geomtype


#1.3
pa <- st_read(dsn = "data/20260126/20260126_protected_areas_be.gpkg", layer = "NaturaSite_polygon")

#1.4
pa_bioregion <- st_read(dsn = "data/20260126/20260126_protected_areas_be.gpkg", layer = "BIOREGION")
pa_habitats <- st_read(dsn = "data/20260126/20260126_protected_areas_be.gpkg", layer = "HABITATS")

#1.5
municipalities <- st_read(municipalities_folder)

#1.6
eea_grid <- st_read(eea_grid_folder)

#1.7
st_crs(occs) # epsg 4326
st_crs(pa) # epsg 3035
st_crs(municipalities) # epsg 4326
st_crs(eea_grid) # epsg 3035

#1.8

pa_a <- pa %>%
  dplyr::filter(SITETYPE == "A")

## Challenge 2

### Raïsa

#Challenge 2.1
occs_3035 <- st_transform(occs, crs = 3035)

#Challenge 2.2
st_write(pa,
         layer = "NaturaSite",
         "data/20260126/pa_occs_3035.gpkg")
st_write(occs_3035,
         layer = "racoon_occs",
         "data/20260126/pa_occs_3035.gpkg")
#check the layers
st_layers("data/20260126/pa_occs_3035.gpkg")

#Challenge 2.3
municipalities_3035 <- st_transform(municipalities,
                                    crs = 3035)

#Challenge 2.4
municipalities_3035_centroids <- st_centroid(municipalities_3035)
#quick visual check
mapview(municipalities_3035) #polygons
mapview(municipalities_3035_centroids) #points

### Nele

#2.1

occs_3035 <- st_transform(occs, 3035)
st_crs(occs_3035)

#2.2
# layers: NaturaSite: protected aread from pa, and raccoon_occs: occs_3035

st_write(pa, layer = "NaturaSite", "data/20260126/pa_occs_3035.gpkg")
st_write(occs_3035, layer = "raccoon_occs", "data/20260126/pa_occs_3035.gpkg")

st_layers("data/20260126/pa_occs_3035.gpkg")

#2.3

municipalities_3035 <- st_transform(municipalities, 3035)
st_crs(municipalities)
st_crs(municipalities_3035)

#2.4

municipalities_3035_centroids <- st_centroid(municipalities_3035)

### Adriaan
occs_3035 <- st_transform(occs, 3035)

st_write(pa, dsn = 'pa_occs_3035.gpkg', layer = 'NaturaSite')
st_write(occs_3035, dsn = 'pa_occs_3035.gpkg', layer = 'raccoon_occs')

municipalities_3035 <- st_transform(municipalities, 3035)

municipalities_3035_centroids <- st_centroid(municipalities_3035)

mapview(municipalities_3035_centroids)

### Jolien

occs_3035 <- st_transform(occs, 3035)
st_crs(occs_3035)

st_write(pa, "pa_occs_3035.gpkg", layer = "NaturaSite")
st_write(occs_3035, "pa_occs_3035.gpkg", layer = "raccoon_occs", append = TRUE)

municipalities_3035 <- st_transform(municipalities, 3035)
st_crs(municipalities_3035)

municipalities_3035_centroids <-  st_centroid(municipalities_3035)

### Lawrence

#1
occs_3035 <- st_transform(occs, crs = 3035)
occs_3035

#2
st_write(pa, layer = "NaturaSite",  "data/20260126/pa_occs_3035.gpkg")
st_write(occs_3035, layer = "raccoon_occs",  "data/20260126/pa_occs_3035.gpkg")
st_layers("data/20260126/pa_occs_3035.gpkg")

#3
municipalities_3035 <- st_transform(municipalities, crs = 3035)
municipalities_3035

#4
municipalities_3035_centroids <- st_centroid(municipalities_3035)
mapview(municipalities_3035_centroids)

## Challenge 3

#Challenge 3.1
cube_sf <- cube %>%
  left_join(eea_grid %>% dplyr::select(CELLCODE)) %>%
  st_as_sf()


#Challenge 3.2
# A list with which pa intersect with each of the municipalities
intersections <- st_intersects(municipalities_3035, pa)
#add number of intersecting pa to municipalities
municipalities_3035 <- municipalities_3035 %>%
  mutate(n_pa = sapply(intersections,
                       FUN = function(x) {length(x)}))

#Challenge 3.3
no_pa_municipalities_3035 <- municipalities_3035 %>%
  dplyr::filter(n_pa == 0)
nrow(no_pa_municipalities_3035) #there are 114 municipalities without pa

#Challenge 3.4
no_pa_municipalities_3035 %>%
  mutate(area = st_area(no_pa_municipalities_3035)) %>%
  arrange(-area) %>%
  head(10)

#Challenge 3.5
intersections <- st_intersects(municipalities_3035, occs_3035)
#I'm unsure whether we need the number of occurences or the number of individuals...
#I'm going with the number of individuals since it would otherwise be similar to 3.2 :-)
#Note that I'm ignoring occurences where "individualCount" is not filled in (it's NA in 129 rows)

municipalities_3035 <- municipalities_3035 %>%
  mutate(n_occs = sapply(intersections,
                         FUN = function(x) {
                           if (length(x) > 0) {
                             sum(occs_3035[x, "individualCount"] %>%
                                   st_drop_geometry(), na.rm = T)}
                           else {
                             0
                             }})
  )

#Challenge 3.6
municipalities_3035 %>%
  arrange(-n_occs) %>%
  head(10)

#Challenge 3.7
#Note that there are overlapping pa which could lead to double counting of the area.
#We use st_cast(st_union()) to simplify the pa first
intersecting_pa <- st_intersection(municipalities_3035, st_cast(st_union(pa),"POLYGON"))
area_pa <- intersecting_pa %>%
  mutate(area = st_area(intersecting_pa)) %>%
  st_drop_geometry() %>%
  group_by(mun_code) %>%
  summarize(protected_area_m2 = sum(area))
municipalities_3035 <- municipalities_3035 %>%
  left_join(area_pa)
mapview(municipalities_3035, zcol = "protected_area_m2")

#Challenge 3.8
municipalities_3035 <- municipalities_3035 %>%
  mutate(area = st_area(municipalities_3035),
         protected_percentage = units::drop_units(protected_area_m2 / area))
mapview(municipalities_3035, zcol = "protected_percentage")
#still, the percentage is sometimes higher than 100% for Hechtel-Eksel:-( 
municipalities_3035 %>%
  dplyr::filter(!is.na(protected_percentage) & protected_percentage > 1) %>%
  arrange(-protected_percentage)