# INBO CODING CLUB

26 October 2021

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

## Participants

Put your name + " | " and add a "*" each time you solve a challenge (see below).
Example below, Damiano solved 2 challenges

Name | Challenges
--- | ---
Damiano Oldoni | **
Lynn Pallemaerts |**
Loïc van Doorn |**
Joost Vanoverbeke |
Hans Van Calster | ***
Emma Cartuyvels|*
Stien Heremans |*
Marijke Thoonen |*
Amber Mertens |***
Camille Van Eupen |**
Jasmijn Hillaert|*
Adriaan Seynaeve|
Jeroen Speybroeck|**
Geert Spanoghe |
Jolien Goossens |***

## Challenge 1

Hans

impatiens_df <- read_tsv("./data/20211026/20211026_impatiens_glandulifera.txt")


#' Transform impatiens_df to a geospatial data.frame using sf package. Note that
#' GBIF data are stored using WGS 84. Hint: find which numeric code is
#' associated with WGS84 coordinate reference system.

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


#' How many layers does the  geospatial file 20211026_protected_areas1.gpkg contain?

st_layers(dsn = "./data/20211026/20211026_protected_areas.gpkg")


#' Read the layer `ps_hbtrl`: call it `prot_areas`
prot_areas <- read_sf("./data/20211026/20211026_protected_areas.gpkg",
                      layer = "ps_hbtrl")

prot_areas

#' What is the coordinate reference system declared by user? Does it coincide
#' with the real Geographic Coordinate Reference System (GEOCRS)?

st_crs(prot_areas)

#' Check that the CRS of `prot_areas` and `spatial_impatiens_df` are the same

waldo::compare(st_crs(prot_areas), st_crs(impatiens_sf))
# apart from the input name, they are the same

Jolien

#' Transform impatiens_df to a geospatial data.frame using sf package. Note that
#' GBIF data are stored using WGS 84. Hint: find which numeric code is
#' associated with WGS84 coordinate reference system.

impatiens_sp = st_as_sf(impatiens_df, 
                        coords = c("decimalLongitude", "decimalLatitude"),
                        crs = 4326)

#' How many layers does the  geospatial file 20211026_protected_areas1.gpkg contain?
st_layers(dsn = "./data/20211026/20211026_protected_areas.gpkg")

#' Read the layer `ps_hbtrl`: call it `prot_areas`

prot_areas = st_read("./data/20211026/20211026_protected_areas.gpkg", layer = "ps_hbtrl")

#' What is the coordinate reference system declared by user? Does it coincide
#' with the real Geographic Coordinate Reference System (GEOCRS)?

st_crs(prot_areas)

#' Check that the CRS of `prot_areas` and `spatial_impatiens_df` are the same

st_crs(prot_areas) == st_crs(impatiens_sp)

Emma:

#' Transform impatiens_df to a geospatial data.frame using sf package. Note that
#' GBIF data are stored using WGS 84. Hint: find which numeric code is
#' associated with WGS84 coordinate reference system.
spatial_impatiens_df <- st_as_sf(impatiens_df,
                                 crs = "EPSG:4326",
                                 coords = c("decimalLongitude", "decimalLatitude"))


#' How many layers does the  geospatial file 20211026_protected_areas.gpkg contain?
st_layers("./data/20211026/20211026_protected_areas.gpkg")


#' Read the layer `ps_hbtrl`: call it `prot_areas`
prot_areas <- st_read("./data/20211026/20211026_protected_areas.gpkg", "ps_hbtrl")


#' What is the coordinate reference system declared by user? Does it coincide
#' with the real Geographic Coordinate Reference System (GEOCRS)?
st_crs(prot_areas)


#' Check that the CRS of `prot_areas` and `spatial_impatiens_df` are the same
st_crs(prot_areas)
st_crs(spatial_impatiens_df)

Jeroen S

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

library(tidyverse)
library(sf)

impatiens_df <- read_tsv("20211026_impatiens_glandulifera.txt")


#' Transform impatiens_df to a geospatial data.frame using sf package. Note that
#' GBIF data are stored using WGS 84. Hint: find which numeric code is
#' associated with WGS84 coordinate reference system.
imp_sf <- st_as_sf(impatiens_df, coords = c("decimalLongitude", "decimalLatitude"))
st_crs(imp_sf)<- 4326
st_crs(imp_sf)


#' How many layers does the  geospatial file 20211026_protected_areas1.gpkg contain?
st_layers("20211026_protected_areas.gpkg")


#' Read the layer `ps_hbtrl`: call it `prot_areas`
prot_areas <- st_read("20211026_protected_areas.gpkg", layer = "ps_hbtrl")


#' What is the coordinate reference system declared by user? Does it coincide
#' with the real Geographic Coordinate Reference System (GEOCRS)?
   
#' Check that the CRS of `prot_areas` and `spatial_impatiens_df` are the same
ifelse(st_crs(prot_areas) == st_crs(imp_sf),"same","not the same")

Camille

library(tidyverse)
library(sf)

impatiens_df <- read_tsv("./data/20211026/20211026_impatiens_glandulifera.txt")


#' Transform impatiens_df to a geospatial data.frame using sf package. Note that
#' GBIF data are stored using WGS 84. Hint: find which numeric code is
#' associated with WGS84 coordinate reference system.

impatiens_sd <- st_as_sf(impatiens_df, coords = 17:16, crs = 4326)


#' How many layers does the  geospatial file 20211026_protected_areas.gpkg contain?

st_layers('./data/20211026/20211026_protected_areas.gpkg')

#' Read the layer `ps_hbtrl`: call it `prot_areas`

prot_areas <- st_read('./data/20211026/20211026_protected_areas.gpkg', 
                      layer = 'ps_hbtrl')

#' What is the coordinate reference system declared by user? Does it coincide
#' with the real Geographic Coordinate Reference System (GEOCRS)?

st_crs(prot_areas)

#' Check that the CRS of `prot_areas` and `spatial_impatiens_df` are the same

st_crs(prot_areas) == st_crs(impatiens_sd)

Amber

#' Transform impatiens_df to a geospatial data.frame using sf package. Note that
#' GBIF data are stored using WGS 84. Hint: find which numeric code is
#' associated with WGS84 coordinate reference system.

impatiens_sf <- st_as_sf(impatiens_df, coords = c("decimalLongitude","decimalLatitude"), crs = 4326)


#' How many layers does the  geospatial file 20211026_protected_areas.gpkg contain?

st_read("./data/20211026/20211026_protected_areas.gpkg")
st_layers("./data/20211026/20211026_protected_areas.gpkg")
# there are two layers: ps_hbtrl, ps_hbtrl_deel

#' Read the layer `ps_hbtrl`: call it `prot_areas`

prot_areas <- st_read("./data/20211026/20211026_protected_areas.gpkg", "ps_hbtrl")


#' What is the coordinate reference system declared by user? Does it coincide
#' with the real Geographic Coordinate Reference System (GEOCRS)?

# CRS declared is WGS84
st_crs(prot_areas)
# ID = 4326, so yes

#' Check that the CRS of `prot_areas` and `spatial_impatiens_df` are the same

if (st_crs(prot_areas) == st_crs(impatiens_sf)) {
  print("crs is the same")
} else {
  print("crs not the same")
}

## Challenge 2

Hans

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

prot_areas <- prot_areas %>%
  st_transform(crs = 3035)
spatial_impatiens_df_3035 <- impatiens_sf %>%
  st_transform(crs = 3035)

plot(prot_areas %>% st_geometry())

plot(spatial_impatiens_df_3035 %>% st_geometry())

#' Write the transformed data as a geopackage file called
#' `prot_areas_and_impatiens_trs1989.gpkg` with two layers: the first called
#' `prot_areas`, containing the protected areas, the second layer,
#' `impatiens_obs`, containing the observations of Himalayan balsam
prot_areas %>%
  st_write(dsn = "./temp/prot_areas_and_impatiens_trs1989.gpkg",
           layer = "prot_areas")

spatial_impatiens_df_3035 %>%
  st_write(dsn = "./temp/prot_areas_and_impatiens_trs1989.gpkg",
           layer = "impatiens_obs")


#' 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 store in column `coordinateUncertaintyInMeters` for
#' `spatial_impatiens_df_3035`

spatial_impatiens_df_3035 %>%
  st_buffer(dist = .$coordinateUncertaintyInMeters) %>%
  st_geometry() %>%
  plot()

Jolien

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

prot_areas_etrs = st_transform(prot_areas, crs = 3035)
impatiens_etrs = st_transform(impatiens_sp, crs = 3035)


#' Write the transformed data as a geopackage file called
#' `prot_areas_and_impatiens_trs1989.gpkg` with two layers: the first called
#' `prot_areas`, containing the protected areas, the second layer,
#' `impatiens_obs`, containing the observations of Himalayan balsam

st_write(prot_areas_etrs, layer = "prot_areas", 
         dsn =  "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg")
st_write(impatiens_etrs, layer = "impatiens_obs", 
         dsn =  "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg")

st_layers(dsn =  "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg")

#' 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 store in column `coordinateUncertaintyInMeters` for
#' `spatial_impatiens_df_3035`
st_crs(impatiens_etrs)$units # m -> ok
impatiens_etrs_buffer = st_buffer(impatiens_etrs, impatiens_etrs$coordinateUncertaintyInMeters)
impatiens_etrs_buffer %>% st_geometry() %>% plot()

Lynn

prot_areas_3035 <- st_transform(prot_areas, crs=3035)
impatiens_3035 <- st_transform(impatiens_sf, crs=3035)

prot_areas_3035 %>%
  st_write(dsn = "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg",
           layer = "prot_areas")
impatiens_3035 %>%
  st_write(dsn = "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg",
           layer = "impatiens_obs")

impatiens_3035 %>%
  st_buffer(dist = .$coordinateUncertaintyInMeters) %>%
  # st_geometry() %>%
  ggplot() +
  geom_sf()

Amber

## CHALLENGE 2

#' Transform both `prot_areas` and `spatial_impatiens_df` to [European
#' Terrestrial Reference System 1989](https://epsg.io/3035), the coordinate
#' reference system used at EU level
impatiens_sf_89 <- st_transform(impatiens_sf, 3035)
st_crs(impatiens_sf_89)

prot_areas_89 <- st_transform(prot_areas, 3035)
st_crs(prot_areas_89)


#' Write the transformed data as a geopackage file called
#' `prot_areas_and_impatiens_trs1989.gpkg` with two layers: the first called
#' `prot_areas`, containing the protected areas, the second layer,
#' `impatiens_obs`, containing the observations of Himalayan balsam

st_write(prot_areas_89, "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg", "prot_areas")
st_write(impatiens_sf_89, "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg", "impatiens_obs", append = TRUE)

st_layers("./data/20211026/prot_areas_and_impatiens_trs1989.gpkg")

#' 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 store in column `coordinateUncertaintyInMeters` for
#' `spatial_impatiens_df_3035`
impatiens_sf_89_circles <- st_buffer(impatiens_sf_89, dist = impatiens_sf_89$coordinateUncertaintyInMeters)

Camille

## CHALLENGE 2

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

prot_areas_ETR <- prot_areas %>% st_transform(3035)
impatiens_sd_ETR <- impatiens_sd %>% st_transform(3035)

#' Write the transformed data as a geopackage file called
#' `prot_areas_and_impatiens_trs1989.gpkg` with two layers: the first called
#' `prot_areas`, containing the protected areas, the second layer,
#' `impatiens_obs`, containing the observations of Himalayan balsam

st_write(prot_areas_ETR, dsn = "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg", 
         layer = 'prot_areas', delete_dsn = T)
st_write(impatiens_sd_ETR, dsn = "./data/20211026/prot_areas_and_impatiens_trs1989.gpkg", 
         layer = 'impatiens_obs')

st_layers("./data/20211026/prot_areas_and_impatiens_trs1989.gpkg")

## Q: How to delete a layer from an existing geopackage?

#' 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 store in column `coordinateUncertaintyInMeters` for
#' `spatial_impatiens_df_3035`

impatiens_circles <- impatiens_sd_ETR %>% 
  st_buffer(dist = impatiens_sd_ETR$coordinateUncertaintyInMeters)

plot(st_geometry(impatiens_circles))
plot(st_geometry(impatiens_sd_ETR), col = 'red',pch = 16,cex = 0.2, add = T)

## INTERMEZZO

## Challenge 3

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

spatial_impatiens_df_3035 %>%
  st_filter(prot_areas,
          join = st_contains)

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

impatiens_in_protected <- spatial_impatiens_df_3035 %>%
  st_buffer(dist = .$coordinateUncertaintyInMeters) %>%
  st_join(prot_areas,
          left = FALSE,
          join = st_intersects)



#' So, how many observations in the protected area "Bos- en heidegebieden ten
#' oosten van Antwerpen"?

impatiens_in_protected %>%
  filter(NAAM == "Bos- en heidegebieden ten oosten van Antwerpen") %>%
  count()



#' 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 issue?

prot_areas_centroids <- prot_areas %>% st_centroid()

Jolien

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

protarcont = st_contains(prot_areas_etrs, impatiens_etrs)
plot(prot_areas_etrs)
prot_areas_etrs %>% st_geometry() %>% plot()

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

protarsect = st_intersects(prot_areas_etrs, impatiens_etrs)

#' So, how many observations in the protected area "Bos- en heidegebieden ten
#' oosten van Antwerpen"?

prot_areas_etrs %>% 
  filter(NAAM == "Bos- en heidegebieden ten oosten van Antwerpen") %>% 
  st_intersects(impatiens_etrs) 

#' 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 issue?


prot_areas_etrs %>% st_centroid() %>% st_geometry() %>%  plot()
prot_areas_etrs %>% st_geometry() %>%  plot()

Amber

## CHALLENGE 3

#' Using data in CRS 3035:
#' 1. Find which observations, as points, are _contained_ in each protected area?
obs_in_prot_areas <- st_contains(prot_areas_89, impatiens_sf_89)


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

obs_intersect_prot_areas <- st_intersects(prot_areas_89, impatiens_sf_89_circles)


#' So, how many observations in the protected area "Bos- en heidegebieden ten
#' oosten van Antwerpen"?
row_nr <- which(prot_areas_89$NAAM == "Bos- en heidegebieden ten oosten van Antwerpen")
intersection_3 <- st_intersects(prot_areas_89, impatiens_sf_89_circles)[row_nr]
# result is a list of 1 element
nr_obs <- length(intersection_3[[1]])
# 4 observations found


#' 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 issue?

st_centroid(prot_areas_89) # this works
st_centroid(prot_areas) # this does not work because 2 features with invalid geometry
prot_areas_valid <- st_make_valid(prot_areas)
st_centroid(prot_areas_valid)

## Bonus challenge
Amber

#' 2. How to get only the observations, as circles, **totally** contained in
#' protected areas? Hint: check the cheat sheet
obs_completely_within_prot_areas <- st_covered_by(impatiens_sf_89_circles, prot_areas_89)