owned this note
owned this note
Published
Linked with GitHub
# INBO CODING CLUB
24 April 2025
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 | ***
Nele Mullens | ***
Jorre |
Sanne Govaert | *
Pieter Huybrechts | ***
Robin Daelemans | ***
Falk Mielke | .*
Amber Mertens | **
Lawrence Whatley | *
## Challenge 1
### Sanne's solution
```r
heracleum_df <- readr::read_tsv(
"data/20250424/20250424_heracleum_BE.tsv"
)
# 1. import shapefile
muncipialities <- sf::st_read(
"data/20250424/20250424_communesgemeente-belgium/georef-belgium-municipality-millesime.shp",
options = "ENCODING=UTF-8"
)
# 2. Create a geospatial data.frame called heracleum starting from heracleuam_df.
heracleum <- heracleum_df %>%
sf::st_as_sf(
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326
)
# 3. Which layers does the geospatial file 20250424_protected_areas_BE.gpkg contain? What differences can you see between the first and the other layers?
# find layers in geospatial file 20250424_protected_areas_BE.gpkg
sf::st_layers(
"data/20250424/20250424_protected_areas_BE.gpkg"
)
# 4. Read the layer NaturaSite_polygon: call it pa.
pa <- sf::st_read(
"data/20250424/20250424_protected_areas_BE.gpkg",
layer = "NaturaSite_polygon"
)
# 5. Read the other layers as well. Call them pa_bioregion and pa_habitats. Are these spatial data.frames?
pa_bioregion <- sf::st_read(
"data/20250424/20250424_protected_areas_BE.gpkg",
layer = "BIOREGION"
)
pa_habitats <- sf::st_read(
"data/20250424/20250424_protected_areas_BE.gpkg",
layer = "HABITATS"
)
# 6. How to retrieve information about the CRS of pa, heracleum and municipalities?
sf::st_crs(pa)
sf::st_crs(heracleum)
sf::st_crs(muncipialities)
# 7. Do pa and heracleum have the same CRS? Do heracleum and municipalities have the same CRS?
sf::st_crs(pa) == sf::st_crs(heracleum)
sf::st_crs(heracleum) == sf::st_crs(muncipialities)
# 8. Extract the protected areas of type (SITETYPE) A as pa_a. Hint: do it as you would do with a standard data.frame. The motto of the sf package is "Spatial data, simplified" for a reason!
pa_a <- pa %>%
dplyr::filter(SITETYPE == "A")
```
### Nele's solution
```r
# 1.1: Import shapefile and name it municipalities
municipalities <- read_sf("./data/20250424/20250424_communesgemeente-belgium", layer = "georef-belgium-municipality-millesime")
# 1.2: Create a geospatial data.frame called heracleum starting from heracleuam_df.
#Use the columns decimalLongitude and decimalLatitude to specify the coordinates.
#Note that GBIF data are stored using WGS 84 (CRS = 4326).
heracleum_df <- read_tsv("data/20250424/20250424_heracleum_BE.tsv")
heracleum <- st_as_sf(heracleuam_df, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
# 1.3 Which layers does the geospatial file 20250424_protected_areas_BE.gpkg contain?
#What differences can you see between the first and the other layers?
# "data/20250424/20250424_protected_areas_BE.gpkg"
st_layers("data/20250424/20250424_protected_areas_BE.gpkg")
# the first layer is a multi polygon
# 1.4 Read the layer NaturaSite_polygon: call it pa.
pa <- st_read("data/20250424/20250424_protected_areas_BE.gpkg", layer = "NaturaSite_polygon")
# 1.5 Read the other layers as well. Call them pa_bioregion and pa_habitats. Are these spatial data.frames?
pa_bioregion <- st_read("data/20250424/20250424_protected_areas_BE.gpkg", https://hackmd.io/s-Obj5yOTxult4YYO3XCxw?both#layer = "BIOREGION")
class(pa_bioregion)
pa_habitats <- st_read("data/20250424/20250424_protected_areas_BE.gpkg", layer = "HABITATS")
class(pa_habitats)
# 1.6 How to retrieve information about the CRS of pa, heracleum and municipalities?
st_crs(pa)
st_crs(heracleum)
st_crs(municipalities)
# 1.7 Do pa and heracleum have the same CRS? Do heracleum and municipalities have the same CRS?
st_crs(pa) == st_crs(heracleum)
#no
st_crs(municipalities) == st_crs(heracleum)
#yes
# 1.8 Extract the protected areas of type (SITETYPE) A as pa_a. Hint: do it as you would do
# with a standard data.frame. The motto of the sf package is "Spatial data, simplified" for a reason!
glimpse(pa)
pa_a <- filter(pa, SITETYPE == "A")
```
### Pieter's solution
```r
# Challenge 1 -------------------------------------------------------------
# Import the shapefile with the municipalities of Belgium: call it
# municipalities.
municipalities <-
st_read(file.path(
"data",
"20250424",
"georef-belgium-municipality-millesime.shp"
))
# Create a geospatial data.frame called heracleum starting from heracleuam_df.
# Use the columns decimalLongitude and decimalLatitude to specify the
# coordinates. Note that GBIF data are stored using WGS 84 (CRS = 4326).
heracleum <-
heracleum_df |>
st_as_sf(
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326
)
# Which layers does the geospatial file 20250424_protected_areas_BE.gpkg
# contain? What differences can you see between the first and the other layers?
st_layers(file.path("data","20250424","20250424_protected_areas_BE.gpkg"))
# Read the layer NaturaSite_polygon: call it pa.
pa <-
st_read(file.path("data", "20250424", "20250424_protected_areas_BE.gpkg"),
layer = "NaturaSite_polygon"
)
# Read the other layers as well. Call them pa_bioregion and pa_habitats. Are
# these spatial data.frames?
pa_bioregion <-
st_read(file.path("data", "20250424", "20250424_protected_areas_BE.gpkg"),
layer = "BIOREGION"
)
pa_habitats <-
st_read(file.path("data", "20250424", "20250424_protected_areas_BE.gpkg"),
layer = "HABITATS"
)
# How to retrieve information about the CRS of pa, heracleum and municipalities?
lapply(list(pa, heracleum, municipalities), st_crs)
# Do pa and heracleum have the same CRS? Do heracleum and municipalities have
# the same CRS?
st_crs(heracleum) == st_crs(pa)
waldo::compare(st_crs(heracleum), st_crs(pa))
# not exactly the same object, different WKT
waldo::compare(st_crs(heracleum), st_crs(municipalities))
## Same identifier, same CRS: we can check in different ways:
st_crs(heracleum) == st_crs(municipalities)
st_crs(heracleum)$epsg == st_crs(municipalities)$epsg
st_crs(municipalities)$epsg == 4326 # from loading heracleum
st_crs(heracleum)$proj4string == st_crs(municipalities)$proj4string
#' Get the EPSG (ID) code of a sf object
#'
#' @param sf An sf object.
#'
#' @return The EPSG (ID) code of the sf object.
get_epsg <- function(sf) {
assertthat::assert_that("sf" %in% class(sf))
st_as_text(st_crs(sf), projjson = TRUE) |>
jsonlite::fromJSON() |>
purrr::chuck("id", "code")
}
get_epsg(heracleum) == get_epsg(municipalities)
# Extract the protected areas of type (SITETYPE) A as pa_a. Hint: do it as you
# would do with a standard data.frame. The motto of the sf package is "Spatial
# data, simplified" for a reason!
pa_a <- dplyr::filter(pa, SITETYPE == "A")
```
### Jorre
```r
# Import the shapefile with the municipalities of Belgium: call it
# municipalities.
municipalities <- st_read(
'data/20250424/20250424_communesgemeente-belgium/georef-belgium-municipality-millesime.shp'
)
# Create a geospatial data.frame called heracleum starting from heracleuam_df.
# Use the columns decimalLongitude and decimalLatitude to specify the
# coordinates. Note that GBIF data are stored using WGS 84 (CRS = 4326).
heracleum <- st_as_sf(
heracleum_df,
coords = c("decimalLongitude", "decimalLatitude"),
crs = 4326
)
# Which layers does the geospatial file 20250424_protected_areas_BE.gpkg
# contain? What differences can you see between the first and the other layers?
st_layers('data/20250424/20250424_protected_areas_BE.gpkg')
# Available layers:
# layer_name geometry_type features fields crs_name
# 1 NaturaSite_polygon Multi Polygon 310 5 ETRS89-extended / LAEA Europe
# 2 BIOREGION NA 317 3 <NA>
# 3 HABITATS NA 3345 16 <NA
# Read the layer NaturaSite_polygon: call it pa.
read_layer <- function(layer) {
st_read('data/20250424/20250424_protected_areas_BE.gpkg',layer)
}
pa <- read_layer('NaturaSite_polygon')
# Read the other layers as well. Call them pa_bioregion and pa_habitats. Are
# these spatial data.frames?
pa_bioregion <- read_layer('BIOREGION')
pa_habitats <- read_layer('HABITATS')
# How to retrieve information about the CRS of pa, heracleum and municipalities?
st_crs(pa)
st_crs(heracleum)
st_crs(municipalities)
# Do pa and heracleum have the same CRS? Do heracleum and municipalities have
# the same CRS?
st_crs(pa) == st_crs(heracleum)
st_crs(heracleum) == st_crs(municipalities)
# Extract the protected areas of type (SITETYPE) A as pa_a. Hint: do it as you
# would do with a standard data.frame. The motto of the sf package is "Spatial
# data, simplified" for a reason!
pa_a <- pa |> filter(SITETYPE=='A')
```
```
```
### Sebastiaan
```r
library(sf)
library(tidyverse)
library(mapview) # Optional
heracleum_df <- readr::read_tsv(
"data/20250424/20250424_heracleum_BE.tsv"
)
# Import the shapefile with the municipalities of Belgium: call it municipalities.
municipalities <- st_read(dsn = "data/20250424/20250424_communesgemeente-belgium/georef-belgium-municipality-millesime.shp")
#View(municipalities)
# Create a geospatial data.frame called heracleum starting from heracleuam_df. Use the columns decimalLongitude and decimalLatitude to specify the coordinates. Note that GBIF data are stored using WGS 84 (CRS = 4326).
#plot(heracleum_df)
class(heracleum_df)
heracleum <- st_as_sf(heracleum_df,coords = c("decimalLongitude","decimalLatitude"))
st_crs(heracleum) <- st_crs("EPSG:4326")
# Better: heracleum <- heracleum_df %>% sf::st_as_sf(coords= c("decimalLongitude","decimalLatitude"), crs = 4326)
class(heracleum)
# Which layers does the geospatial file 20250424_protected_areas_BE.gpkg contain? What differences can you see between the first and the other layers?
st_layers("data/20250424/20250424_protected_areas_BE.gpkg")
# Read the layer NaturaSite_polygon: call it pa.
# Read the other layers as well. Call them pa_bioregion and pa_habitats. Are these spatial data.frames?
pa <- st_read(dsn="data/20250424/20250424_protected_areas_BE.gpkg",layer="NaturaSite_polygon")
pa_habitats <- st_read(dsn="data/20250424/20250424_protected_areas_BE.gpkg",layer="HABITATS")
pa_bioregion <- st_read(dsn="data/20250424/20250424_protected_areas_BE.gpkg",layer="BIOREGION")
# How to retrieve information about the CRS of pa, heracleum and municipalities?
st_crs(pa)
st_crs(heracleum)
st_crs(municipalities)
# Do pa and heracleum have the same CRS? Do heracleum and municipalities have the same CRS?
if (st_crs(pa) == st_crs(heracleum)){
print("True")
} else {
print("False")
}
if (st_crs(pa) == st_crs(municipalities)){
print("True")
} else {
print("False")
}
# Extract the protected areas of type (SITETYPE) A as pa_a. Hint: do it as you would do with a standard data.frame. The motto of the sf package is "Spatial data, simplified" for a reason!
pa_a <- pa |> filter(SITETYPE=='A')
```
## Challenge 2
### Jorre
```r
# Create pa_a_wgs84 by transforming pa_a to WGS 84 (EPSG code: 4326).
pa_a_wgs84 <- st_transform(pa_a, crs = 4326)
# Write a geopackage called pa_a_heracleum_4326.gpkg with two layers:
# NaturaSite_A, containing the protected areas of type A (pa_a_wgs84) and
# heracleum_occs, containing heracleum. Hint: check
# https://r-spatial.github.io/sf/articles/sf2.html#using-st_write.
file <- 'data/20250424/pa_a_heracleum_4326.gpkg'
st_write(pa_a_wgs84,file,layer="NaturaSite_A" )
st_write(heracleum ,file,layer="heracleum_occs",append=TRUE)
# 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
# heracleum. Call this sf object heracleum_circles.
heracleum_circles <-
st_buffer(heracleum,dist=heracleum$coordinateUncertaintyInMeters)
# Find the centroids of the municipalities. Call it municipalities_centroids.
municipalities_centroids <- st_centroid(municipalities)
```
### Nele's solution
```r
# 2.1 Create pa_a_wgs84 by transforming pa_a to WGS 84 (EPSG code: 4326).
pa_a_wgs84 <- st_transform(pa_a, crs = 4326)
st_crs(pa_a_wgs84)$epsg
# 2.2 Write a geopackage called pa_a_heracleum_4326.gpkg with two layers: NaturaSite_A,
# containing the protected areas of type A (pa_a_wgs84) and heracleum_occs, containing
# heracleum. Hint: check https://r-spatial.github.io/sf/articles/sf2.html#using-st_write.
st_write(pa_a_wgs84, "data/20250424/pa_a_heracleum_4326.gpkg", layer = "NaturaSite_A")
st_write(heracleum, "data/20250424/pa_a_heracleum_4326.gpkg", layer = "heracleum_occs")
# 2.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 heracleum. Call
# this sf object heracleum_circles.
heracleum_circles <- st_buffer(heracleum, dist = heracleum$coordinateUncertaintyInMeters)
mapview(heracleum_circles)
# 2.4 Find the centroids of the municipalities. Call it municipalities_centroids.
municipalities_centroids <- st_centroid(municipalities)
mapview(municipalities_centroids)
```
### Lawrence
```r
# 1
pa_a_wgs84 <- st_transform(pa_a, crs = "epsg:4326")
st_crs(pa_a_wgs84) #"EPSG",4326
# 2
st_write(pa_a_wgs84, dsn = "pa_a_heracleum_4326.gpkg", layer = "NaturaSite_A.gpkg", driver = "GPKG")
st_write(heracleum, dsn = "pa_a_heracleum_4326.gpkg", layer = "heracleum_occs.gpkg", driver = "GPKG")
# 3
heracleum_circles <- st_buffer(heracleum, dist = heracleum$coordinateUncertaintyInMeters)
mapview(heracleum_circles)
# 4
municipalities_centroids <- st_centroid(municipalities)
mapview(municipalities_centroids)
```
### Pieter
```r
# Create pa_a_wgs84 by transforming pa_a to WGS 84 (EPSG code: 4326).
pa_a_wgs84 <- st_transform(pa_a, crs = st_crs(4326))
# Write a geopackage called pa_a_heracleum_4326.gpkg with two layers:
# NaturaSite_A, containing the protected areas of type A (pa_a_wgs84) and
# heracleum_occs, containing heracleum. Hint: check
# https://r-spatial.github.io/sf/articles/sf2.html#using-st_write.
st_write(pa_a_wgs84, dsn = "pa_a_heracleum_4326.gpkg", layer = "NaturaSite_A")
st_write(heracleum, dsn = "pa_a_heracleum_4326.gpkg", layer = "heracleum_occs")
# 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
# heracleum. Call this sf object heracleum_circles.
heracleum_circles <-
st_buffer(heracleum, dist = heracleum$coordinateUncertaintyInMeters)
mapview(heracleum_circles)
# Find the centroids of the municipalities. Call it municipalities_centroids.
municipalities_centroids <- st_centroid(municipalities)
mapview(municipalities_centroids)
```
### Falk's attempt
```r
# read the crs ID
get_crs_id <- function (data) stringr::str_extract(
gsub(" ", "\t", sf::st_crs(data)$wkt),
"\n\tID\\[\"?(EPSG|[a-zA-Z]+)\",([0-9]+)\\]",
group = 2)
print(get_crs_id(pa))
# Create pa_a_wgs84 by transforming pa_a to WGS 84 (EPSG code: 4326)
pa_a_wgs84 <- sf::st_transform(pa_a, 4326)
# Write a geopackage called pa_a_heracleum_4326.gpkg with two layers: NaturaSite_A, containing the protected areas of type A (pa_a_wgs84) and heracleum_occs, containing heracleum. Hint: check https://r-spatial.github.io/sf/articles/sf2.html#using-st_write.
gpkg_datafile <- "./codingclub/data/20250424_pa_a_heracleum_4326.gpkg"
# NaturaSite_A
sf::st_write(pa_a_wgs84,
dsn = gpkg_datafile,
layer = "NaturaSite_A",
append = FALSE, quiet = TRUE
)
# heracleum_occs
sf::st_write(heracleum,
dsn = gpkg_datafile,
layer = "heracleum_occs",
append = FALSE, quiet = TRUE
)
# 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 heracleum. Call this sf object heracleum_circles
heracleum_proj <- sf::st_transform(heracleum, 31370) # Projecting first will noticably increase computation speed!
heracleum_circles <- sf::st_buffer(
heracleum_proj,
heracleum_proj$coordinateUncertaintyInMeters,
nQuadSegs = 1 # not necessary
)
mapview(heracleum_circles)
# Find the centroids of the municipalities. Call it municipalities_centroids.
mapview(sf::st_centroid(municipalities))
```
## Challenge 3
### Pieter
```r
# Challenge 3 -------------------------------------------------------------
# For each occurrence, find the nearest municipality as centroid.
nearest_centroids <- st_nearest_feature(heracleum, municipalities_centroids)
# For each occurrence, find the nearest municipality as polygon. Do you find
# some differences with the previous results?
nearest_polygons <- st_nearest_feature(heracleum, municipalities)
all(nearest_centroids == nearest_polygons)
waldo::compare(nearest_centroids, nearest_polygons)
## Make a little table
tibble::tibble(
nearest_centr = nearest_centroids,
nearest_poly = nearest_polygons,
occurrenceID = heracleum$occurrenceID
)
# Protected areas extend often over several municipalities. For each protected
# area of type A, find (the index of) the municipalities it belongs to.
pa_a_municipalities <-
st_intersects(pa_a, st_transform(municipalities, crs = st_crs(pa_a)))
pa_a_municipalities
# Create a new sf object called pa_durme_kruibeke with the part of the protected
# area "Durme en Middenloop van de Schelde" (SITECODE: "BE2301235") falling
# within the municipality of Kruibeke (mun_name_up: "KRUIBEKE").
municipalities_3035 <- st_transform(municipalities, crs = st_crs(3035))
pa_durme_kruibeke <-
st_intersection(
dplyr::filter(pa_a, SITECODE == "BE2301235"),
dplyr::filter(municipalities_3035, mun_name_up == "KRUIBEKE")
)
mapview(pa_durme_kruibeke)
# We would expect that the distance between pa_durme_kruibeke and the
# municipality of Kruibeke is 0. How to check it? Can you calculate the distance
# between the municipalities and the protected areas of type A? Calculate the
# distance using Lambert72 (EPSG: 31370) projection: do you get the same
# distances?
kruibeke <- dplyr::filter(municipalities_3035, mun_name_up == "KRUIBEKE")
st_distance(pa_durme_kruibeke, kruibeke)
dists_3035 <- st_distance(pa_a, municipalities_3035)
#' Transform sf object to Belge Lambert 72 (EPSG: 31370)
#'
#' @param sf An sf object.
#' @param target_epsg The target EPSG code to transform the sf object to. By
#' default 31370 (Belge Lambert 72).
#'
#' @return The transformed sf object.
to_bl72 <- function(sf, target_epsg = 31370){
assertthat::assert_that("sf" %in% class(sf))
st_transform(sf, crs = target_epsg)
}
dists_bl72 <- st_distance(to_bl72(pa_a), to_bl72(municipalities_3035))
## There are small differences
waldo::compare(dists_3035, dists_bl72)
## But not if we allow 1cm tolerance
waldo::compare(dists_3035, dists_bl72, tolerance = 0.01)
# For each protected area of type A, get (the index of) the occurrences as
# circles that intersect within the protected area. How to get only the
# occurrences that are totally contained in the protected area?
st_intersects(
st_transform(pa_a, crs = st_crs(heracleum_circles)),
heracleum_circles
)
st_covered_by(
heracleum_circles,
st_transform(pa_a, crs = st_crs(heracleum_circles))
)
# Sometimes you need to grid your polygons. Examples: you need to do a transect
# survey with a standardized research effort. Create a grid with 5kmx5km cells.
# 🐝🐝🐝 #
st_make_grid(
to_bl72(kruibeke),
cellsize = units::set_units(5, "km"),
what = "polygons",
square = FALSE,
flat_topped = TRUE) %>%
st_sf() %>%
mapview()
```
### Jorre
```r
# For each occurrence, find the nearest municipality as centroid.
tmp_nearest <- st_nearest_feature(heracleum,municipalities_centroids)
tmp_muni <- municipalities |>
st_drop_geometry(municipalities) |>
magrittr::extract(tmp_nearest,) |>
select(mun_name_nl_centroid=mun_name_nl)
heracleum_nearest_muni <- cbind(heracleum,tmp_muni)
rm(list=ls(pattern="^tmp"))
# For each occurrence, find the nearest municipality as polygon. Do you find
# some differences with the previous results?
tmp_nearest <- st_nearest_feature(heracleum,municipalities)
tmp_muni <- municipalities |>
st_drop_geometry(municipalities) |>
magrittr::extract(tmp_nearest,) |>
select(mun_name_nl_polygon=mun_name_nl)
heracleum_nearest_muni <- cbind(heracleum_nearest_muni,tmp_muni)
rm(list=ls(pattern="^tmp"))
heracleum_nearest_muni |>
filter(mun_name_nl_centroid!=mun_name_nl_polygon) |>
View()
# Protected areas extend often over several municipalities. For each protected
# area of type A, find (the index of) the municipalities it belongs to.
tmp_muni <- municipalities |>
st_drop_geometry() |>
mutate(muninumber=row_number(),mun_name_nl,.keep='none')
pa_a_intersections <- pa_a_wgs84 %>%
mutate(muninumber=st_intersects(., municipalities)) |>
unnest(muninumber) |>
left_join(tmp_muni,join_by(muninumber)) |>
select(-muninumber)
rm(list=ls(pattern="^tmp"))
# Create a new sf object called pa_durme_kruibeke with the part of the protected
# area "Durme en Middenloop van de Schelde" (SITECODE: "BE2301235") falling
# within the municipality of Kruibeke (mun_name_up: "KRUIBEKE").
kruibeke <- filter(municipalities,mun_name_up == "KRUIBEKE")
pa_durme_kruibeke <- st_intersection(
filter(pa_a_wgs84,SITECODE == "BE2301235"),
kruibeke
)
# We would expect that the distance between pa_durme_kruibeke and the
# municipality of Kruibeke is 0. How to check it? Can you calculate the distance
# between the municipalities and the protected areas of type A? Calculate the
# distance using Lambert72 (EPSG: 31370) projection: do you get the same
# distances?
st_distance(pa_durme_kruibeke,kruibeke)
st_distance(
st_transform(pa_durme_kruibeke,crs=31370),
st_transform(kruibeke,crs=31370)
)
# For each protected area of type A, get (the index of) the occurrences as
# circles that intersect within the protected area. How to get only the
# occurrences that are totally contained in the protected area?
heracleum_intersect <- heracleum_circles |>
st_within(pa_a_wgs84) |>
enframe(value='pa_number') |>
mutate(gbifID=heracleum_circles$gbifID) |>
unnest(pa_number) |>
mutate(sitecode=st_drop_geometry(pa_a)[pa_number,'SITENAME'])
# Sometimes you need to grid your polygons. Examples: you need to do a transect
# survey with a standardized research effort. Create a grid with 5kmx5km cells.
```
Google Drive
Gist
Clipboard
Sign in with Wallet
