# INBO CODING CLUB

27 June 2023

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 |
Pieter Huybrechts | ***
Anja Leyman | **
Leen Govaere |
Nele Mullens | **
Jonas Bertels | **
Pieter Verschelde | **
Adriaan Seynaeve|

## Challenge 1

### Pieter H's solution

# 1
lepidoptera_df <-
  read_tsv("./data/20230627/20230627_lepidoptera_iNaturalist_2023.csv")

## st_as_sf() will try to convert an object to an sf, sometimes you need to
## provide some extra information
lepidoptera_sf <-
  st_as_sf(lepidoptera_df,
    coords = c("decimalLongitude", "decimalLatitude"), crs = 4326
  )
# 2
st_layers("data/20230627/20230627_protected_areas.gpkg")
# 3
prot_areas <-
  st_read("data/20230627/20230627_protected_areas.gpkg", layer = "ps_hbtrl")
## geopackages are really handy, they can contain a lot of extra information in
## just one file!

# 4
st_crs(prot_areas)
# 5
## not exactly the same
waldo::compare(
  st_crs(prot_areas),
  st_crs(lepidoptera_sf)
)
## but close enough really, only the user provided string is different
st_crs(prot_areas) == st_crs(lepidoptera_sf)
# 6
provinces_be <- readRDS(file = "./data/20230627/20230627_provinces_be.rds")

class(provinces_be)
st_as_sf(provinces_be)
st_as_sf(provinces_be) %>% mapview()

# 7
provinces_fl <-
  st_as_sf(provinces_be) %>%
  dplyr::filter(TX_RGN_DESCR_NL == "Vlaams Gewest")
mapview(provinces_fl)

## Challenge 2

### Pieter H's Solution

# 1

## Transform both prot_areas and lepidoptera to the European Terrestrial
## Reference System 1989 (EPSG: 3035)

lepidoptera_3035 <- st_transform(lepidoptera_sf, 3035)
## let's check if it worked
st_crs(lepidoptera_3035)
mapview(lepidoptera_3035)

prot_areas_3035 <- st_transform(prot_areas, 3035)
## let's check if it worked
st_crs(prot_areas_3035)
mapview(prot_areas_3035)
# 2

## Write the transformed data as a geopackage file called
## prot_areas_and_lepidoptera_3035.gpkg with two layers: the first called
## prot_areas, containing the protected areas, the second layer,
## lepidoptera_obs, containing the observations of lepidoptera

st_write(
  prot_areas_3035,
  dsn = file.path("data", "20230627", "prot_areas_and_lepidoptera_3035.gpkg"),
  layer = "prot_areas"
)

st_write(
  lepidoptera_3035,
  dsn = file.path("data", "20230627", "prot_areas_and_lepidoptera_3035.gpkg"),
  layer = "lepidoptera_obs"
)

# 3

## Due to spatial uncertainty (gridded data, GPS uncertainty, etc.) observations
## should not be idealized as points in space, but as circles. Create such
## circles using the values stored in column coordinateUncertaintyInMeters of
## lepidoptera_3035. Call this sf object lepidoptera_3035_circles.


## You could probably get away with not specifying the units
lepidoptera_3035_circles <-
  st_buffer(
    lepidoptera_3035,
    dist = units::set_units(lepidoptera_3035$coordinateUncertaintyInMeters, m)
  )

### Nele's Solution

# 2.1
lepidoptera_3035 <- st_transform(lepidoptera_sf, epsg = 3035)
prot_areas_3035 <- st_transform(prot_areas, epsg = 3035)

# 2.2

st_write(prot_areas_3035, 
         dsn = "data/20230627/prot_areas_and_lepidoptera_3035.gpkg",
         layer = "prot_areas")

st_write(lepidoptera_3035, 
         dsn = "data/20230627/prot_areas_and_lepidoptera_3035.gpkg",
         layer = "lepidoptera_obs")

# 2.3
lepidoptera_3035_circles <-
  st_buffer(
    lepidoptera_3035,
    dist = lepidoptera_3035$coordinateUncertaintyInMeters, endCapStyle="ROUND"
  )

## Challenge 3

### Pieter H's Solution

# 1

## Which observations in lepidoptera_3035 (points), are contained in which
## protected area?
st_contains(prot_areas_3035, lepidoptera_3035)

lepidoptera_3035 %>%
  mutate(in_prot_area =
           purrr::map_dbl(
             st_within(lepidoptera_3035,prot_areas_3035),
             ~ifelse(is.null(.x),NA,.x)
             )
         ) %>%
  filter(!is.na(in_prot_area))
# 2

## But we should maybe better check which observations, as circles (!),
## intersect each protected area. How to do it?

st_intersects(lepidoptera_3035_circles, prot_areas_3035)

# as a filtered tibble
lepidoptera_3035[
  !purrr::map_lgl(
    st_intersects(lepidoptera_3035_circles, prot_areas_3035),
    rlang::is_empty),]
## check results
mapview(
  list(intersects =
         lepidoptera_3035[!purrr::map_lgl(
           st_intersects(lepidoptera_3035_circles, prot_areas_3035),
           rlang::is_empty),
           ],
       prot_areas = prot_areas_3035),
  col.regions = c("red", "green")
  )

# 3

## So, how many observations as circles intersect the Sonian Forest
## ("Zoniënwoud")?

## Method 1: spatial filtering
sonian_forest <- filter(prot_areas_3035, NAAM == "Zoniënwoud")


st_filter(lepidoptera_3035_circles,
          sonian_forest,
          .predicate = st_intersects) # you could also set st_contains for example


## Method 2: by using dplyr filtering, way more complicated

lepidoptera_in_prot_area_circl <-
  lepidoptera_3035 %>%
  mutate(in_prot_area_circl =
           purrr::map_dbl(
             st_intersects(lepidoptera_3035_circles, prot_areas_3035),
             ~ifelse(is.null(.x),NA,.x)
           )
  ) %>%
  filter(!is.na(in_prot_area_circl))


## what is the number of Zoniënwoud?
filter(mutate(prot_areas_3035, id = row_number()), NAAM == "Zoniënwoud") %>%
  pull(id)

## now we can filter on this id, it's probably better to use dplyr joins, and to
## filter on the name directly

filter(lepidoptera_in_prot_area_circl,
       in_prot_area_circl == 29)

## Let's have a look:


mapview(list(lepidoptera = filter(lepidoptera_in_prot_area_circl,
                                  in_prot_area_circl == 29),
             sonian_forest = filter(mutate(prot_areas_3035,id = row_number()),
                                    NAAM == "Zoniënwoud")),
        col.regions = c("springgreen", "slateblue"))

## add the uncertaintly, you could filter this one down as well..
mapview(list(lepidoptera = filter(lepidoptera_in_prot_area_circl,
                                  in_prot_area_circl == 29),
             sonian_forest = filter(mutate(prot_areas_3035,id = row_number()),
                                    NAAM == "Zoniënwoud"),
             uncert = lepidoptera_3035_circles),
        col.regions = c("cyan2", "wheat4", "red"))

# 4

## Sometimes it's interesting to calculate the centroid of a polygon, e.g. for
## visualizations. Easy by using sf function st_centroids(). However, you get an
## error while calculating the centroids of prot_areas. What does it mean? How
## to solve the problem?

st_centroid(prot_areas)
## because they aren't all valid geometries!
st_is_valid(prot_areas)

## but reprojecting has seemed to solve it
mapview(list(centroid = st_centroid(prot_areas_3035),
             area = prot_areas_3035),
        zcol = "NAAM",
        legend = FALSE
        )

## or you could use a function to magically make geometries valid
st_make_valid(prot_areas) %>%
  st_is_valid()

# 5
provinces_be_sf <- st_as_sf(provinces_be)
## How to get the portion of protected areas located in province Antwerp?
st_intersection(
  st_make_valid(prot_areas),
  filter(provinces_be_sf, TX_PROV_DESCR_NL == "Provincie Antwerpen")
) %>%
  mapview(zcol = "NAAM",
          legend = FALSE)
# 6

## How to get the union of all the protected areas? How such union is expressed
## in sf?

## st_combine returns a single, combined geometry, with no resolved boundaries; returned geometries may well be invalid.
## If y is missing, st_union(x) returns a single geometry with resolved boundaries, else the geometries for all unioned pairs of x[i] and y[j].
st_union(prot_areas_3035) %>% mapview()

## but is invalid
st_combine(prot_areas_3035) %>% mapview()
st_combine(prot_areas_3035) %>% st_is_valid()`

## Bonus challenge