# INBO CODING CLUB

25 October 2022

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 |
Anja Leyman |**
Lynn Pallemaerts | **
Emma Cartuyvels|
Thomas Scheppers |
Hans Van Calster |
Marijke Thoonen |
Teun Everts |
Adriaan Seynaeve |
Amber Mertens |**
Raïsa Carmen |
Zoé Hermans
An Leyssen |
Wouter Depaepe |
Leen Govaere
Natalie Beenaerts
Wim Kuypers
Heleen Deroo
Helee

## Challenge 1
Raisa:

# Create a geospatial data.frame called spatial_ludwigia_df starting from ludwigia_df. Note that GBIF data are stored using WGS 84. Hint: find first which numeric code is associated with WGS84 coordinate reference system and use the cheatsheet.
spatial_ludwigia_df <- st_as_sf(ludwigia_df,
                                coords = c("decimalLongitude", "decimalLatitude"),
                                crs = "EPSG:4326")


# 2. How many layers does the  geospatial file `20221025_protected_areas.gpkg`
# contain?
lyrs <- sf::st_layers("./data/20221025/20221025_protected_areas.gpkg")
length(lyrs$name)


# 3. Import the layer `ps_hbtrl`: call it `prot_areas`
prot_areas <- st_read("./data/20221025/20221025_protected_areas.gpkg",
        layer = lyrs$name[1])


# 4. What is the CRS declared by user? Does it coincide with the real Geographic
# Coordinate Reference System (GEOCRS)?
st_crs(prot_areas)


# 5. Do `prot_areas` and `spatial_ludwigia_df` have the same CRS?
st_crs(prot_areas) == st_crs(spatial_ludwigia_df)


# 6. Read the Belgian provinces rds file as `be_provinces_sp` (the code is
# given!). What is the class of this variable? From which package? How to
# transform it to a sf object?

be_provinces_sp <- read_rds(
  file = "./data/20221025/20221025_be_provinces_sp.rds"
)
summary(be_provinces_sp )
typeof(be_provinces_sp)
be_provinces_sf <- st_as_sf(be_provinces_sp)
st_crs(prot_areas) == st_crs(be_provinces_sf)

# 7. Extract the Flemish provinces.
provinces_flanders <- be_provinces_sf %>% 
  filter(TX_RGN_DESCR_NL == "Vlaams Gewest")

Hans

library(tidyverse)
library(sf)

ludwigia_df <- read_tsv("./data/20221025/20221025_ludwigia_grandiflora.txt")

spatial_ludwigia_df <- ludwigia_df %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)


gpkg_file <- file.path("data/20221025/20221025_protected_areas.gpkg")



# 1. Transform both `prot_areas` and `spatial_ludwigia_df` to [European
# Terrestrial Reference System 1989](https://epsg.io/3035) (EPSG: 3035), the
# coordinate reference system used at EU level

View(spatial_ludwigia_df)
View(prot_areas)


# 2. How many layers does the  geospatial file `20221025_protected_areas.gpkg`
# contain?

st_layers(gpkg_file)

length(st_layers(gpkg_file)$name)

# 3. Import the layer `ps_hbtrl`: call it `prot_areas`

prot_areas <- read_sf(gpkg_file, layer = "ps_hbtrl")
plot(st_geometry(prot_areas))

# 4. What is the CRS declared by user? Does it coincide with the real Geographic
# Coordinate Reference System (GEOCRS)?

prot_areas

# 5. Do `prot_areas` and `spatial_ludwigia_df` have the same CRS?

waldo::compare(
  st_crs(prot_areas),
  st_crs(spatial_ludwigia_df))

st_crs(prot_areas) == st_crs(spatial_ludwigia_df)

# 6. Read the Belgian provinces rds file as `be_provinces_sp` (the code is
# given!). What is the class of this variable? From which package? How to
# transform it to a sf object?

be_provinces_sp <- read_rds(
  file = "./data/20221025/20221025_be_provinces_sp.rds"
)

class(be_provinces_sp)

be_provinces_sf <- st_as_sf(be_provinces_sp)


# 7. Extract the Flemish provinces.

be_prov_fla <- be_provinces_sf %>%
  filter(TX_RGN_DESCR_NL == "Vlaams Gewest")

plot(st_geometry(be_prov_fla))

Lynn:

# 1. Create a geospatial data.frame called spatial_ludwigia_df starting from
# ludwigia_df. Note that GBIF data are stored using WGS 84. Hint: find first
# which numeric code is associated with WGS84 coordinate reference system and
# use the cheatsheet.

ludwigia_df <- read_tsv("./data/20221025/20221025_ludwigia_grandiflora.txt")
spatial_ludwigia_df <- st_as_sf(ludwigia_df,
                                coords = c("decimalLongitude", "decimalLatitude"),
                                crs = 4326)

# 2. How many layers does the geospatial file 20221025_protected_areas.gpkg
# contain?

st_layers("./data/20221025/20221025_protected_areas.gpkg") # maar dit toont maar de eerste layer

# 3. Import the layer ps_hbtrl: call it prot_areas

prot_areas <- st_read("./data/20221025/20221025_protected_areas.gpkg",
                      layer = "ps_hbtrl")

# 4. What is the CRS declared by user? Does it coincide with the real Geographic
# Coordinate Reference System (GEOCRS)?

prot_areas
st_crs(prot_areas)


# 5. Do prot_areas and spatial_ludwigia_df have the same CRS?

st_crs(spatial_ludwigia_df) == st_crs(prot_areas)

# 6. Read the Belgian provinces rds file as `be_provinces_sp` (the code is
# given!). What is the class of this variable? From which package? How to
# transform it to a sf object?

be_provinces_sp <- read_rds("./data/20221025/20221025_be_provinces_sp.rds")
class(be_provinces_sp)
st_crs(be_provinces_sp)
be_provinces_sp <- st_as_sf(be_provinces_sp,
                            crs = 6326)

# 7. Extract the Flemish provinces.

fl_provinces_sp <- be_provinces_sp %>%
  filter(CD_RGN_REFNIS == 2000 )

## Challenge 2

Lynn:

# 1. Transform both `prot_areas` and `spatial_ludwigia_df` to [European
# Terrestrial Reference System 1989](https://epsg.io/3035), the coordinate
# reference system used at EU level

ludwigia_3035 <- st_transform(spatial_ludwigia_df,
                              crs = 3035)
prot_areas_3035 <- st_transform(prot_areas,
                                crs = 3035)

# 2. Write the transformed data as a geopackage file called
# `prot_areas_and_ludwigia_3035.gpkg` with two layers: the first called
# `prot_areas`, containing the protected areas, the second layer,
# `ludwigia_obs`, containing the observations of water primrose

st_write(prot_areas_3035,
         layer = "prot_areas", 
         "./data/20221025/prot_areas_and_ludwigia_3035.gpkg",
         append = F)
st_write(ludwigia_3035,
         layer = "ludwigia_obs",
         "./data/20221025/prot_areas_and_ludwigia_3035.gpkg",
         append = T)

# 3. Due to spatial uncertainty (gridded data, GPS uncertainty, etc.) GBIF
# observations should not be idealized as points in space, but as circles.
# Create such circles using the values stored in column
# `coordinateUncertaintyInMeters` of `spatial_ludwigia_df_3035`

st_write(prot_areas_3035,
         layer = "prot_areas", 
         "./data/20221025/prot_areas_and_ludwigia_3035.gpkg",
         append = F)
st_write(ludwigia_3035,
         layer = "ludwigia_obs",
         "./data/20221025/prot_areas_and_ludwigia_3035.gpkg",
         append = T)

Hans:

# 1. Transform both `prot_areas` and `spatial_ludwigia_df` to [European
# Terrestrial Reference System 1989](https://epsg.io/3035), the coordinate
# reference system used at EU level

prot_areas <- st_transform(prot_areas, crs = 3035)
spatial_ludwigia_df <- st_transform(spatial_ludwigia_df, crs = 3035)


# 2. Write the transformed data as a geopackage file called
# `prot_areas_and_ludwigia_3035.gpkg` with two layers: the first called
# `prot_areas`, containing the protected areas, the second layer,
# `ludwigia_obs`, containing the observations of water primrose



st_write(prot_areas, "prot_areas_and_ludwigia_3035.gpkg", layer = "prot_areas")
st_write(spatial_ludwigia_df, "prot_areas_and_ludwigia_3035.gpkg",
         layer = "ludwigia_obs")

st_layers("prot_areas_and_ludwigia_3035.gpkg")

# 3. Due to spatial uncertainty (gridded data, GPS uncertainty, etc.) GBIF
# observations should not be idealized as points in space, but as circles.
# Create such circles using the values store in column
# `coordinateUncertaintyInMeters` of `spatial_ludwigia_df_3035`


spatial_ludwigia_df %>%
  st_buffer(dist = .$coordinateUncertaintyInMeters)

## Challenge 3

Hans:

# Using data in CRS 3035:
# 1. Find which observations, as points, are _contained_ in each protected area?


st_contains(x = prot_areas, y = spatial_ludwigia_df)

# You can use `st_intersects()` as well: same result. However, st_contains gives
# more the idea of what you are doing




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


st_intersects(prot_areas, spatial_ludwigia_buffer)

# 3. So, how many observations in the protected area called `"Demervallei"`?

st_join(spatial_ludwigia_buffer, prot_areas, left = FALSE) %>%
  filter(NAAM == "Demervallei") %>%
  count(NAAM)

# 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)

# Using data in CRS 3035:
# 1. Find which observations, as points, are _contained_ in each protected area?
prot_area_ob <- st_contains(x = prot_areas_etrs,
                            y = spatial_ludwigia_df_etrs)



# You can use `st_intersects()` as well: same result. However, st_contains gives
# more the idea of what you are doing
prot_area_ob2 <- st_intersects(y = prot_areas_etrs,
              x = spatial_ludwigia_df_etrs) #the other way around works >< contains
prot_area_ob3 <- st_intersects(x = prot_areas_etrs,
              y = spatial_ludwigia_df_etrs) #the same as before




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

prot_area_ob4 <- st_intersects(x = prot_areas_etrs,
                               y = circles)


# 3. So, how many observations in the protected area called `"Demervallei"`?
demervallei <- which(prot_areas$NAAM == "Demervallei")
length(prot_area_ob4[[demervallei]])

# 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(st_make_valid(prot_areas))

## Bonus challenge