Accessing the STAC Catalogue

1 What is a STAC Catalog?

From Bon in a Box’s “Working with STAC” documentation, which you can see here:

STAC stands for SpatioTemporal Asset Catalog. It is a common language to describe geospatial information, so it can be more easily worked with, indexed, and discovered. Spatial information is stored in a geoJSON format, which can be easily adapted to different domains. The goal of STAC is to enable a global index of all imagery derived data products with an easily implementable standard for organization to expose their data in a persistent and reliable way. STAC objects are hosted in geoJSON format as HTML pages.

2 Biodiversité Québec’s STAC Catalog

Biodiversité Québec hosts many layers for biodiversity maps and models for Canada and Québec. Most of the layers we will use for challenges will be from this STAC Catalog.

You can view the layers at the following link:

Optionally, you can also browse the GEO BON STAC Catalog for global layers and access them using the documentation here.

iNaturalist density maps

A series of density maps of iNaturalist observations is hosted on the Biodiversité Québec STAC. This means you can access density maps and work with them as a proxy, or download it locally to use as the basis for a Blitz the Gap challenge. There are two resolutions available for each taxonomic group: 100m, and 1km.

First, load the R packages you will need to use to access the density map layers:

library(gdalcubes)
library(rstac)
library(knitr)
library(stars)
library(terra)
library(ggplot2)
library(mapview)

Set up the connection to the STAC catalog:

# set up
stac_obj <- stac("https://io.biodiversite-quebec.ca/stac/")

# get iNaturalist heatmaps for Canada
it_obj <- stac_obj |>
  stac_search(collections = "inat_canada_heatmaps") |>
  post_request() |> items_fetch()
it_obj

Get a table showing which layers are available in the catalog:

# See layers in the object
df <- data.frame(description=character())
for (f in it_obj[['features']]){
  df <- rbind(df,data.frame(description=f$properties$description)) 
}
df$feature_number = 1:nrow(df)

description	feature_number
Density of observations in iNaturalist Canada for species of taxonomic group Reptilia at 1km resolution	1
Density of observations in iNaturalist Canada for species of taxonomic group Reptilia at 100m resolution	2
Density of observations in iNaturalist Canada for species of taxonomic group Protozoa at 1km resolution	3
Density of observations in iNaturalist Canada for species of taxonomic group Protozoa at 100m resolution	4
Density of observations in iNaturalist Canada for species of taxonomic group Plantae at 1km resolution	5
Density of observations in iNaturalist Canada for species of taxonomic group Plantae at 100m resolution	6
Density of observations in iNaturalist Canada for species of taxonomic group Mollusca at 1km resolution	7
Density of observations in iNaturalist Canada for species of taxonomic group Mollusca at 100m resolution	8
Density of observations in iNaturalist Canada for species of taxonomic group Mammalia at 1km resolution	9
Density of observations in iNaturalist Canada for species of taxonomic group Mammalia at 100m resolution	10
Density of observations in iNaturalist Canada for species of taxonomic group Insecta at 1km resolution	11
Density of observations in iNaturalist Canada for species of taxonomic group Insecta at 100m resolution	12
Density of observations in iNaturalist Canada for species of taxonomic group Fungi at 1km resolution	13
Density of observations in iNaturalist Canada for species of taxonomic group Fungi at 100m resolution	14
Density of observations in iNaturalist Canada for species of taxonomic group Chromista at 1km resolution	15
Density of observations in iNaturalist Canada for species of taxonomic group Chromista at 100m resolution	16
Density of observations in iNaturalist Canada for species of taxonomic group Aves at 1km resolution	17
Density of observations in iNaturalist Canada for species of taxonomic group Aves at 100m resolution	18
Density of observations in iNaturalist Canada for species of taxonomic group Arachnida at 1km resolution	19
Density of observations in iNaturalist Canada for species of taxonomic group Arachnida at 100m resolution	20
Density of observations in iNaturalist Canada for species of taxonomic group Animalia at 1km resolution	21
Density of observations in iNaturalist Canada for species of taxonomic group Animalia at 100m resolution	22
Density of observations in iNaturalist Canada for species of taxonomic group Amphibia at 1km resolution	23
Density of observations in iNaturalist Canada for species of taxonomic group Amphibia at 100m resolution	24
Density of observations in iNaturalist Canada for species of taxonomic group All at 1km resolution	25
Density of observations in iNaturalist Canada for species of taxonomic group All at 100m resolution	26
Density of observations in iNaturalist Canada for species of taxonomic group Actinopterygii at 1km resolution	27
Density of observations in iNaturalist Canada for species of taxonomic group Actinopterygii at 100m resolution	28

df

Let’s access one of the layers to map it, as an example. We’ll map the density of iNaturalist observations across all taxa by selecting one of the features we listed in the table above.

How do we select one of the layers in the STAC Catalog? You can index them by checking the feature_number column above.

Here, we want “Density of observations in iNaturalist Canada for species of taxonomic group All at 1km resolution”, which is the 25th layer:

# 37 is the number to change if you want another layer from the STAC
inat <- read_stars(paste0('/vsicurl/', it_obj[['features']][[25]]$assets[[1]]$href), 
                   proxy = TRUE) # using it as a proxy (not a local download,
                                 # to save time and space!)

Now, let’s map it!

First, load a polygon of Canada that you can download from the data/raw folder on the Sharepoint (or the 00_rawdata folder on the GitHub repository). Make sure to update the file path to match the folder in which you placed the file when you downloaded it!

# read the Canada polygon
canada = sf::read_sf("data/raw/base-layers/canada-polygon/canada.outline.shp")

Prepare the data for plotting, and use ggplot2 to map it!

# crop the raster to canada (to cut off the ocean, USA, etc.)
inat = sf::st_crop(inat, canada)

# (optional) transform data for easier plotting
# this will square root each raster cell's value. this could also be log10, etc.
inat_sqrt = inat |> sqrt()

# map!
ggplot() +
  geom_stars(data = inat_sqrt, 
             # plotting at a coarser resolution to make this faster
             downsample = 10, 
             na.action = na.omit) +
  scale_fill_viridis_c(option = "turbo", 
                       end = .5,
                       na.value = "transparent") +
  theme_void() +
  theme(legend.position = "top") +
  labs(fill = "Density") +
  coord_equal() 

To make an interactive map that we can zoom into, we can use the package mapview:

First, let’s some mapview options to apply to all the maps we will make here:

mapviewOptions(basemaps = c("OpenStreetMap"),
               na.color = "transparent")

Now, here’s the map!

pal = viridis::turbo(5)
mapview(inat_sqrt, 
        col.regions = pal, # assign color palette
        layer.name = "Density (sqrt)") 

Number of pixels is above 5e+05.Only about 5e+05 pixels will be shown.
You can increase the value of `maxpixels` to 24627120 to avoid this.