owned this note
owned this note
Published
Linked with GitHub
# 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
