Reminders

1. Did we confirm the room reservation on the roomie?
2. Did we start the recording?

The sf package

The sf homepage is the best place to start your "spatial" voyage.

The sf package

sf stays for: simple features, a scary way to express "geographic features", where feature is a scary term to refer to an abstraction of real world phenomena with geometrical properties such as lat/lon (or x-y).

I use sp. Why do you bother me with sf?

Because Thierry says so :-)

I use sp. Why do you bother me with sf?

sp is getting old: some of its core dependencies are getting archived as obsolete.

sf is the present and the future, bro!

This sp-sf migration table will help you refactoring your code.

Install the package:

install.packages("sf")

Load the package:

library(sf)

Introduction: sf package

Spatial coordinates

Two types of spatial coordinates:

• projected (or Cartesian) coordinates: refer to points on a flat space. The coordinates are the distance along the x and y axis
• unprojected or geographic coordinates: refer to angles (latitude, longitude) pointing to locations on a sphere (or ellipsoid)

For the math lovers: the flat space is also referred to as R², the sphere as S².

No geospatial data without specifying the Coordinate Reference System (CRS) you work with!

And what is a projection? A way to visualize a surface of a sphere on a flat space.

Projections

A lot of projections exist. Each one has its own purpose. Check this projection reference page by Bill Rankin, 2006 containing "(almost) all the projections available in ArcGIS".

Example: Lambert Azimuthal Equal Area, used by European Environment Agency (EEA) to provide standard grids of Europe. Parameters: latitude of origin 52° N, longitude of origin 10° E

How to get started?

Check the Each session setup to get started.

First time coding club?

Check the First time setup section to setup.

Download data and code

You can download the material of today:

• automatically via inborutils::setup_codingclub_session()*

• manually** from GitHub folders coding-club/data/20230627 and coding-club/src/20230627

* Note: you can use the date in "YYYYMMDD" format to download the coding club material of a specific day, e.g. run setup_codingclub_session("20220428") to download the coding club material of April, 28 2022. If date is omitted, i.e. setup_codingclub_session(), the date of today is used. For all options, check the tutorial online.
** Note: check the getting started instructions on how to download a single file

Data and scripts description

Data

• 20230627_lepidoptera_iNaturalist_2023.csv: observations of Lepidoptera taken in 2023 via iNaturalist. Derived from: GBIF.org (22 June 2023) GBIF Occurrence Download https://doi.org/10.15468/dl.4z38cf
• 20230627_protected_areas.gpkg: geopackage containing the Natura2000 protected areas located in Flanders
• 20230627_provinces_be.rds: R object containing the Belgian provinces

Scripts

• 20230627_challenges.R: R script to start with.

Challenge 1

Using the cheatsheet and starting from the code provided in 20230627_challenges.R:

1. Create a geospatial data.frame called lepidoptera starting from lepidoptera_df. Use the columns decimalLongitude and decimalLatitude. Note that GBIF data are stored using WGS 84 (CRS = 4326).
2. Which layers does the geospatial file 20230627_protected_areas.gpkg contain?
3. Import the layer ps_hbtrl from 20230627_protected_areas.gpkg: call it prot_areas
4. How to retrieve information about the Coordinate Reference System (CRS) of prot_areas?
5. Do prot_areas and lepidoptera have the same CRS?
6. Read the Belgian provinces rds file as provinces_be (the code is given!). What is the class of this variable? From which package does this class come from? How to transform it to a sf object?
7. Extract the Flemish provinces as provinces_fl. Hint: do it as you would do it in a standard data.frame. The motto of the sf package is "Spatial data, simplified" for a reason!

Intermezzo: CRS & EPSG codes

A common way to specify the CRS is by providing the EPSG (numeric) code. EPSG stands for European Petroleum Survey Group and is an organization that maintains a geodetic parameter database with standard codes, the EPSG codes, for coordinate systems, datums, spheroids, units and such alike. There are a lot of EPSG codes! Full list: https://spatialreference.org/

Challenge 2

1. Transform both prot_areas and lepidoptera to the European Terrestrial Reference System 1989 (EPSG: 3035), the coordinate reference system used at EU level. Call the transformed sf objects prot_areas_3035 and lepidoptera_3035 respectively
2. Write the transformed data as a geopackage file called prot_areas_and_lepidoptera_3035.gpkg with two layers: the first called prot_areas, containing the protected areas, the second layer, lepidoptera_obs, containing the observations of lepidoptera
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 lepidoptera_3035. Call this sf object lepidoptera_3035_circles.

Challenge 3

Using data saved in CRS 3035:

1. Which observations in lepidoptera_3035 (points), are contained in which protected area?
2. But we should maybe better check which observations, as circles (!), intersect each protected area. How to do it?
3. So, how many observations as circles intersect the Sonian Forest ("Zoniënwoud")?
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?
5. How to get the portion of protected areas located in province Antwerp?
6. How to get the union of all the protected areas? How such union is expressed in sf?

The package of the month. Damiano's choice

I discovered Floris is the maintainer of qgisprocess, a package to provide an R interface to the popular and open source desktop geographic information system (GIS) program QGIS. Something I think most of you will find very helpful.

And it has a cheat sheet too!

Floris wrote also some slides presenting the package.

Bonus challenge

1. How to get the observations, as circles, totally covered in protected areas? Hint: it could be not present on the cheat sheet
2. Not a sf question, but still nice to solve: how to add to prot_areas_3035 a column called n_lepidoptera_pts with the number of observations as points in each protected area? To have an idea about the effects of coordinate uncertainty on number of observations, add also the number of observations as circles intersecting protected areas and the ones totally contained as n_lepidoptera_intersect and n_lepidoptera_all_in
3. Let's use qgisprocess! To run some models, you would like to create 50 random points in each protected area to simulate some set of observations. How to do it using qgisprocess?

Announcement

A Journey through Arrow in R

ROpenSci organizes a community call at 28 June (= tomorrow) at 6pm about the arrow R package. arrow exposes an interface to the Arrow C++ library, which is one of the libraries implementing the Apache Arrow cross-language development platform, extremely useful for in-memory and larger-than-memory data.

Everyone is welcome. No RSVP needed. Check details on event page.

Resources

• Challenges solutions are available
• Edited video recording available at the INBO coding club vimeo channel: https://vimeo.com/840403950
• Geocomputation with R
• New book of Edzer Pebesma Spatial Data Science - With Applications in R
• Spatial Data Science
• INBO tutorial on how to use open vector file formats such as geopackage in R
• sf cheatsheet even if it is far from being the best cheatsheet ever
• All articles on sf package documentation webpage are great material and very well written
• qgisprocess package homepage
• Floris' slides presenting qgisprocess package
• What does your favorite map projection say about you? Extremely funny!
• arrow R package
• ROpenSci event about arrow on Jun 28
• Keynote "Spherical geometry" of Edzer Pebesma at the OpenGeoHub summer school 2021: slides
• Keynote "R spatial" of Edzer Pebesma at the useR! 2021 conference: video and slides. Very interesting to learn the principles and basics of spatial data in R
• everything what Edzer Pebesma says or write, actually

Summer break

July = summer = no INBO coding club

library(cowsay)
cowsay::say(what = "Enjoy the summeR!", by = "egret", what_color = c("red", "orange", "green", "blue", "purple"))

Voting session for the next coding club will be launched soon. Mail will follow.