Script examples using DuckDB and DuckSpatial to read Parquet files and perform spatial operations
Guillaume Larocque
2026-05-06
library(dplyr)
library(duckdb)
library(duckspatial)
library(duckdbfs)
library(sf)Create a DuckDB database connection. The database will be stored in
the file mydb.duckdb in the current working directory. If
the file does not exist, it will be created.
con <- dbConnect(duckdb::duckdb(), dbdir = "mydb.duckdb")
dbExecute(con, "INSTALL spatial; LOAD spatial; INSTALL httpfs; LOAD httpfs;")## [1] 0Create a view named atlas that reads the Parquet file
from the specified URL. The read_parquet function is used
to read the Parquet file and create a view that can be queried using
SQL. A view is a virtual table that is defined by a SQL query. It does
not store data itself, but rather provides a way to access data from one
or more tables or other views. In this case, the view atlas
will allow us to query the data in the Parquet file as if it were a
regular table in the database.
dbExecute(con, "CREATE OR REPLACE VIEW atlas AS SELECT * FROM read_parquet('https://object-arbutus.cloud.computecanada.ca/bq-io/atlas/parquet/atlas_public_2026-05-04.parquet')")## [1] 0Describe the structure of the atlas view.
atl <- dbGetQuery(con, "DESCRIBE atlas;")
atl## column_name
## 1 code_flags
## 2 valid_scientific_name
## 3 observation_value
## 4 observation_type
## 5 effort_sampling_value
## 6 effort_sampling_method
## 7 effort_sampling_unit
## 8 longitude
## 9 latitude
## 10 coordinate_uncertainty
## 11 coordinate_uncertainty_unit
## 12 year_obs
## 13 month_obs
## 14 day_obs
## 15 time_obs
## 16 vernacular_fr
## 17 vernacular_en
## 18 group_en
## 19 group_fr
## 20 observed_rank
## 21 kingdom
## 22 phylum
## 23 class
## 24 order
## 25 family
## 26 genus
## 27 dataset_name
## 28 dataset_creator
## 29 dataset_publisher
## 30 dataset_url_information
## 31 dataset_url_download
## 32 dataset_doi
## 33 source_parent_event_id
## 34 source_event_id
## 35 source_occurrence_id
## 36 license
## 37 within_quebec
## 38 geom
## 39 geom_bbox
## column_type null key default
## 1 VARCHAR[] YES <NA> <NA>
## 2 VARCHAR YES <NA> <NA>
## 3 DOUBLE YES <NA> <NA>
## 4 VARCHAR YES <NA> <NA>
## 5 DOUBLE YES <NA> <NA>
## 6 VARCHAR YES <NA> <NA>
## 7 VARCHAR YES <NA> <NA>
## 8 DOUBLE YES <NA> <NA>
## 9 DOUBLE YES <NA> <NA>
## 10 VARCHAR YES <NA> <NA>
## 11 VARCHAR YES <NA> <NA>
## 12 INTEGER YES <NA> <NA>
## 13 INTEGER YES <NA> <NA>
## 14 INTEGER YES <NA> <NA>
## 15 VARCHAR YES <NA> <NA>
## 16 VARCHAR YES <NA> <NA>
## 17 VARCHAR YES <NA> <NA>
## 18 VARCHAR YES <NA> <NA>
## 19 VARCHAR YES <NA> <NA>
## 20 VARCHAR YES <NA> <NA>
## 21 VARCHAR YES <NA> <NA>
## 22 VARCHAR YES <NA> <NA>
## 23 VARCHAR YES <NA> <NA>
## 24 VARCHAR YES <NA> <NA>
## 25 VARCHAR YES <NA> <NA>
## 26 VARCHAR YES <NA> <NA>
## 27 VARCHAR YES <NA> <NA>
## 28 VARCHAR YES <NA> <NA>
## 29 VARCHAR YES <NA> <NA>
## 30 VARCHAR[] YES <NA> <NA>
## 31 VARCHAR[] YES <NA> <NA>
## 32 VARCHAR YES <NA> <NA>
## 33 VARCHAR YES <NA> <NA>
## 34 VARCHAR YES <NA> <NA>
## 35 VARCHAR YES <NA> <NA>
## 36 VARCHAR YES <NA> <NA>
## 37 BOOLEAN YES <NA> <NA>
## 38 GEOMETRY('OGC:CRS84') YES <NA> <NA>
## 39 STRUCT(xmin FLOAT, ymin FLOAT, xmax FLOAT, ymax FLOAT) YES <NA> <NA>
## extra
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## 7 <NA>
## 8 <NA>
## 9 <NA>
## 10 <NA>
## 11 <NA>
## 12 <NA>
## 13 <NA>
## 14 <NA>
## 15 <NA>
## 16 <NA>
## 17 <NA>
## 18 <NA>
## 19 <NA>
## 20 <NA>
## 21 <NA>
## 22 <NA>
## 23 <NA>
## 24 <NA>
## 25 <NA>
## 26 <NA>
## 27 <NA>
## 28 <NA>
## 29 <NA>
## 30 <NA>
## 31 <NA>
## 32 <NA>
## 33 <NA>
## 34 <NA>
## 35 <NA>
## 36 <NA>
## 37 <NA>
## 38 <NA>
## 39 <NA>Connect using duckdbfs allowing dplyr syntax
library(duckdbfs)
atlas_rem <- open_dataset("https://object-arbutus.cloud.computecanada.ca/bq-io/atlas/parquet/atlas_public_2026-05-04.parquet", tblname = 'atlas')See the available datasets in the atlas by grouping the records by
dataset_name and counting the number of records for each
dataset. The results are ordered in descending order of count.
datasets <- atlas_rem |> group_by(dataset_name) |> summarize(cnt=count()) |> arrange(desc(cnt))
datasets## # Source: SQL [?? x 2]
## # Database: DuckDB 1.5.2 [glaroc@Linux 6.17.0-20-generic:R 4.6.0/:memory:]
## # Ordered by: desc(cnt)
## dataset_name cnt
## <chr> <dbl>
## 1 EOD – eBird Observation Dataset 41235823
## 2 iNaturalist Research-grade Observations 1013971
## 3 Données de biodiversité de l'Observatoire global du St-Laurent 829995
## 4 Points d'observation écologique 808022
## 5 Placettes-échantillons temporaires du quatrième inventaire 526616
## 6 Placettes-échantillons temporaires du deuxième inventaire 372418
## 7 Placettes-échantillons temporaires du troisième inventaire 349620
## 8 Placettes-échantillons permanentes 204311
## 9 Placettes-échantillons temporaires du cinquième inventaire 181020
## 10 Herbier du Québec (QUE) - Collection de plantes vasculaires 121263
## # ℹ more rowsObtain Eastern bluebird (Sialia sialis) records from the atlas by filtering the dataset for records with the valid scientific name ‘Sialia sialis’ and collecting the results into a data frame.
sialia_sialis <- atlas_rem |> filter(valid_scientific_name == 'Sialia sialis') |> select(valid_scientific_name, latitude, longitude, dataset_name, year_obs, day_obs) |> collect()
plot(sialia_sialis$longitude, sialia_sialis$latitude, cex=0.5)
GBIF parquet dataset
Connect to the GBIF dataset. Note that the GBIF dataset is not a GeoParquet file yet. https://registry.opendata.aws/gbif/
dbExecute(con, "SET s3_region='us-east-1';SET s3_use_ssl=false;")## [1] 0dbExecute(con, "CREATE OR REPLACE VIEW gbif AS SELECT * FROM read_parquet('s3://gbif-open-data-us-east-1/occurrence/2026-04-01/occurrence.parquet/*')")## [1] 0gbif <- dbGetQuery(con, "DESCRIBE gbif;")Query GBIF dataset to get the number of records per kingdom.
gbif_kingdom <- dbGetQuery(con, "SELECT kingdom, COUNT(*) AS cnt FROM gbif GROUP BY kingdom ORDER BY cnt DESC;")Spatial operations in DuckDB
dbExecute(con, "CREATE OR REPLACE VIEW muni_qc AS SELECT * FROM 'https://data.fieldmaps.io/edge-matched/humanitarian/intl/adm3_polygons.parquet' WHERE adm1_name = 'Quebec'")## [1] 0dbExecute(con, "CREATE OR REPLACE VIEW longueuil AS SELECT * FROM muni_qc WHERE adm3_name = 'Longueuil'")## [1] 0Get bluebird records in the municipality of Longueuil, Quebec, and plot their locations. Not that the resulting output is not an sf object, but a regular data frame.
bluebird_mtl <- dbGetQuery(con,"SELECT atlas.* FROM atlas, longueuil WHERE ST_Contains(longueuil.geometry, atlas.geom) AND atlas.valid_scientific_name = 'Sialia sialis'")
plot(bluebird_mtl$longitude, bluebird_mtl$latitude, cex=0.5)
Spatial operations with duckspatial
If we want to get an sf object, we can use the
duckspatial package to read the view as an sf object.
longueuil_sf <- ddbs_read_table(con, 'longueuil', quiet = FALSE)## ✔ view longueuil successfully imported.plot(longueuil_sf["fid"])
Extract the bluebird records in Atlas as an sf object and plot them.
bluebird_atlas <- ddbs_read_table(con, 'atlas', clauses = "WHERE valid_scientific_name = 'Sialia sialis'", quiet = FALSE)## ✔ view atlas successfully imported.plot(bluebird_atlas["dataset_name"], cex = 0.5, key.pos = 4, key.width = lcm(5), key.length = 0.5)
Use duckspatial to get all the observations from iNaturalist in Atlas from Longueuil.
atlas_longueuil_inat <- ddbs_join('atlas', 'longueuil', join = "within", conn=con) |>
filter(dataset_name == 'iNaturalist Research-grade Observations') |>
ddbs_collect()
plot(longueuil_sf["fid"], color = "lightblue", main = "iNaturalist Observations in Longueuil", key.pos = 4, key.width = lcm(5), key.length = 0.5, reset = FALSE)## Warning in title(...): "color" is not a graphical parameter
## Warning in title(...): "color" is not a graphical parameterplot(atlas_longueuil_inat['valid_scientific_name'], cex=0.5, add=TRUE)