# INBO CODING CLUB

28 May 2024

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 | ***
Raïsa Carmen |***
Rhea Maesele | *
Adriaan Seynaeve |**
Lies Teunen |*
Sanne govaert | ***
Suzanna Lettens |**
Lotte Pohl | ***
Siebe Indestege | ***
Jo-Hannes Nowé | ***
Soria Delva| ***

## Challenge 1

### Damiano's solution (example)

Copy paste this section to show your solutions.

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

### Sanne's solution

### 1.1
lepidoptera <- st_as_sf(
  lepidoptera_df,
  coords = c("decimalLongitude", "decimalLatitude"),
  crs = 4326
  )

### 1.2
st_layers(pa_filepath)

### 1.3
pa <- st_read(pa_filepath, layer = "NaturaSite_polygon")

### 1.4
pa_bioregion <- st_read(pa_filepath, layer = "BIOREGION")
pa_habitats <- st_read(pa_filepath, layer = "HABITATS")
class(pa_bioregion)
class(pa_habitats)

### 1.5
st_crs(pa)

### 1.6
utm_10 <- st_read(utm_10_filepath)

### 1.7
st_crs(pa) == st_crs(lepidoptera)
st_crs(pa) == st_crs(utm_10)

### 1.8
pa_b <- pa %>% 
  filter(SITETYPE == "B")
unique(pa_b$SITETYPE)

### Siebe's Solution

# 1.1
lepidoptera_sf <- st_as_sf(lepidoptera_df, 
                           coords = c("decimalLongitude", "decimalLatitude"), 
                           crs = 4326)

# 1.2
st_layers(pa_filepath) # First has a geometry_type and specified CRS

# 1.3
pa <- st_read(dsn = pa_filepath, layer = "NaturaSite_polygon")

# 1.4
pa_bioregion <- st_read(dsn = pa_filepath, layer = "BIOREGION")
pa_habitats <- st_read(dsn = pa_filepath, layer = "HABITATS")

# 1.5
st_crs(pa)

# 1.6
utm_10 <- st_read(dsn = utm_10_filepath, layer = "20240528_utm10_belgium")

# 1.7
st_crs(pa) == st_crs(lepidoptera_sf)
st_crs(pa) == st_crs(utm_10)

# 1.8
pa_b <- pa %>% filter(SITETYPE == "B")

### Lotte's solution

# 1. geospatial_df 
lepidoptera <- 
  lepidoptera_df %>%
    sf::st_as_sf(coords = c('decimalLongitude', 'decimalLatitude'), crs = 'EPSG:4326')

# 2. layers
pa_layers <- sf::st_layers(dsn = pa_filepath)
pa_layers

# 3. load layer NaturaSite_polygon
pa <- sf::st_read(dsn = pa_filepath, layer = 'NaturaSite_polygon')

# 4. other layers
pa_bioregion <- sf::st_read(dsn = pa_filepath, layer = 'BIOREGION')
pa_habitats <- sf::st_read(dsn = pa_filepath, layer = 'HABITATS')

# 5. CRS of pa
pa_crs <- sf::st_crs(pa)
pa_crs

# 6. import EEA reference grid
utm_10 <- sf::st_read(dsn = utm_10_filepath)

# 7. check CRS

sf::st_crs(pa) == sf::st_crs(lepidoptera)

sf::st_crs(pa) == sf::st_crs(utm_10)

# 8. extract pa_b

pa_b <- 
  pa %>%
    dplyr::filter(SITETYPE == 'B')

### Jo-Hannes' solution

#1.
lepidoptera <- st_as_sf(lepidoptera_df, coords=c("decimalLongitude","decimalLatitude"),crs=4326)

#2.
layers <- st_layers(pa_filepath)
pa_layer1 <- st_read(pa_filepath,layer=layers[[1]][1])
pa_layer2 <- st_read(pa_filepath,layer=layers[[1]][2])
pa_layer2 <- st_read(pa_filepath,layer=layers[[1]][3])

#Only the first layer contains geometry

#3.
pa <- st_read(pa_filepath,layer="NaturaSite_polygon")

#4.
pa_bioregion <- st_read(pa_filepath,layer="BIOREGION")
pa_habitats <- st_read(pa_filepath,layer="HABITATS")
#These are not spatial dataframes

#5. 
crs <- st_crs(pa)

#6.
utm_10 <- st_read(utm_10_filepath)

#7.
st_crs(pa)[1]
st_crs(lepidoptera)[1]
st_crs(utm_10)[1]
#pa and utm_10 have the same crs, but lepidoptera has a different one

#8. 
pa_b <- pa[pa$SITETYPE == "B",]

## Challenge 2
### Soria's solution

#2.1
lepidoptera_3035<-st_transform(lepidoptera, crs=3035)

#2.2
st_write(pa,"./data/20240528/pa_and_lepidoptera_3035.gpkg", "pa")
st_write(lepidoptera_3035, "./data/20240528/pa_and_lepidoptera_3035.gpkg", "lepidoptera_obs", append = TRUE)

#2.3
lepidoptera_3035_circles<-st_buffer(lepidoptera_3035, dist=lepidoptera_3035$coordinateUncertaintyInMeters)
plot(lepidoptera_3035_circles)

### Jo-Hannes' solution

#1.
lepidoptera_3035 <- st_transform(lepidoptera,crs="EPSG:3035")

#2.
st_write(pa, dsn="data/pa_and_lepidoptera_3035.gpkg", layer='pa')
st_write(lepidoptera_3035, dsn="data/pa_and_lepidoptera_3035.gpkg", layer='lepidoptera_obs')

#3. 

lepidoptera_3035_circles <- st_buffer(lepidoptera_3035$geometry,dist = lepidoptera$coordinateUncertaintyInMeters)

### Siebe's solution

### 2.1
lepidoptera_3035 <- st_transform(lepidoptera_sf, crs = 3035)

### 2.2
st_write(pa, "./data/20240528/pa_and_lepidoptera_3035.gpkg", "pa")
st_write(lepidoptera_3035, "./data/20240528/pa_and_lepidoptera_3035.gpkg", "lepidoptera_obs", append = TRUE)
st_layers("./data/20240528/pa_and_lepidoptera_3035.gpkg")

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

### Sanne's solution

### 2.1
lepidoptera_3035 <- st_transform(lepidoptera,ncrs = 3035)

### 2.2
st_write(pa, dsn = "./data/20240528/pa_and_lepidoptera_3035.gpkg", layer = "pa")
st_write(lepidoptera_3035, dsn = "./data/20240528/pa_and_lepidoptera_3035.gpkg", layer = "lepidoptera_obs")

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

### Lotte's solution

# 1. transform crs of df
lepidoptera_3035 <- lepidoptera %>% sf::st_transform(crs = 'EPSG: 3035')

# 2. write pa_and_lepidoptera_3035.gpkg

pa_lep_path <- "./data/20240528/pa_and_lepidoptera_3035.gpkg"

sf::st_write(obj = pa, dsn = pa_lep_path, layer = 'pa')
sf::st_write(obj = lepidoptera_3035, dsn = pa_lep_path, layer = 'lepidoptera_obs', append = T)

# 3. make circles of observations

lepidoptera_3035_circles <- sf::st_buffer(x = lepidoptera_3035, dist = lepidoptera_3035$coordinateUncertaintyInMeters)
lepidoptera_3035_circles %>% mapview()

## Challenge 3

### Sanne's solution

### 3.1
st_contains(pa, lepidoptera_3035)

### 3.2
st_intersects(pa, lepidoptera_3035_circles)

### 3.3
st_centroid(pa)
st_area(pa)

### 3.4
st_covered_by(pa, lepidoptera_3035_circles)

### 3.5
union_pa <- st_union(pa)
class(union_pa)

### 3.6
st_intersects(pa, utm_10)

### Jo-Hannes' solution

#1. 
contain <-st_contains(pa,lepidoptera_3035)
contained_area <-cbind(pa$SITENAME,contain)


#2. 
intersects <- st_intersects(pa,lepidoptera_3035_circles)
intersect_area <- cbind(pa$SITENAME,intersects)

#3.
pa_centroid <- st_centroid(pa)
mapview(pa_centroid)
pa_area <- st_area(pa)

#4. 
covered <- st_covered_by(lepidoptera_3035_circles,pa)
lepidoptera_covered <- cbind(pa$SITENAME,covered)

#5. 
pa_union <- st_union(pa)
mapview(pa_union)
pa_union


#6. 
intersect_utm_pa <- st_intersects(pa,utm_10)
intersection_utm_pa <- cbind(pa$SITENAME,covered)

### Lotte's solution

# 1. which obs in lep_3035 are CONTAINED in which protected area?
lep_3035_contains_pa <- sf::st_contains(pa, lepidoptera_3035)

## question: how to transform it into a df or tibble nicely, to be able to access the rownum of the points that are within each pa?

# 2. intersection
lep_3035_intersection_pa <- sf::st_intersects(pa, lepidoptera_3035_circles)

# 3. centroids and areas
pa_centroids <- sf::st_centroid(pa)

pa_area <- sf::st_area(pa) %>% units::set_units(km^2)
pa_area

# 4. totally contained circles
lep_3035_covered_by_pa <- sf::st_covered_by(pa, lepidoptera_3035_circles)
## question: the dfs 'lep_3035_covered_by_pa' and 'lep_3035_contains_pa' differ a lot. I would have expected that the df 'lep_3035_covered_by_pa' is 
## basically the same as 'lep_3035_contains_pa' with a few less entries since not all big circles are fully within a pa.

# 5. union
pa_union <- sf::st_union(pa)
class(pa_union)

## weird, viz works in leaflet but not in mapview (even not when crs is transformed before plotting)
pa_union %>% sf::st_transform(crs = 4326) %>% mapview()

leaflet() %>%
  addTiles %>%
  addPolygons(data = pa_union %>% sf::st_transform(crs = 4326))

# 6. utm_10 grid cells intersect pa

utm_10_intersects_pa <- sf::st_intersects(pa, utm_10)

### Siebe's solution

### 3.1
contain <- st_contains(pa,lepidoptera_3035)
contained_pa <- cbind(pa$SITENAME,cont)

### 3.2
intersect <- st_intersects(pa,lepidoptera_3035_circles)
intersect_pa <- cbind(pa$SITENAME,intersect)

### 3.3
st_centroid(pa)
st_area(pa)

### 3.4
contain_prop <- st_contains_properly(pa,lepidoptera_3035_circles)
contained_prop_pa <- cbind(pa$SITENAME,contain_prop)

### 3.5
union_pa <- st_union(pa)
class(union_pa)

### 3.6
intersect_utm <- st_intersects(pa,utm_10)
intersect_utm_pa <- cbind(pa$SITENAME,intersect_utm)