library(freestiler)
library(mapgl)
library(duckdb)
library(duckspatial)
con <- dbConnect(duckdb::duckdb(), dbdir = "mydb.duckdb")
dbExecute(con, "INSTALL spatial; LOAD spatial; INSTALL httpfs; LOAD httpfs;")

## [1] 0

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] 0

muni_qc <- ddbs_read_table(con, 'muni_qc', quiet = FALSE)

## ✔ view muni_qc successfully imported.

Create a PMTiles file from the muni_qc view using the freestiler package. The freestile function takes the input data, the output file name, the layer name, and the minimum and maximum zoom levels for the tiles.

freestile(muni_qc, "pmtiles/muni_qc.pmtiles", layer_name = "municipalities", min_zoom = 4,  max_zoom = 12)

## Creating MVT tiles (zoom 4-12) for 1285 features across 1 layer...

##   Parsed 1285 features across 1 layer in 0.1s
##   Zoom  4/12:      4 tiles ...
##                 4 encoded (0.2s)
##   Zoom  5/12:     12 tiles ...
##                11 encoded (0.1s)
##   Zoom  6/12:     29 tiles ...
##                26 encoded (0.1s)
##   Zoom  7/12:     83 tiles ...
##                76 encoded (0.2s)
##   Zoom  8/12:    279 tiles ...
##               242 encoded (0.6s)
##   Zoom  9/12:   1000 tiles ...
##               855 encoded (1.6s)
##   Zoom 10/12:   3821 tiles ...
##              3152 encoded (6.2s)
##   Zoom 11/12:  14970 tiles ...
##             12016 encoded (29.1s)
##   Zoom 12/12:  59228 tiles ...
##             46880 encoded (126.3s)
##   Total: 63262 tiles in 164.5s
##   Writing PMTiles archive (63262 tiles) ...
##   PMTiles write: 0.4s

## Created /home/glaroc/Repos/QCBS-Repos/space-clouds-workshop/pmtiles/muni_qc.pmtiles (24.7 MB)

## View with: view_tiles("muni_qc.pmtiles")

dbExecute(con, "CREATE OR REPLACE VIEW atlas_inat AS SELECT * FROM atlas WHERE dataset_name = 'iNaturalist Research-grade Observations'")

## [1] 0

atlas_inat <- ddbs_read_table(con, 'atlas_inat', quiet = FALSE)

## ✔ view atlas_inat successfully imported.

freestile(atlas_inat, "pmtiles/atlas_inat.pmtiles", layer_name = "iNat", min_zoom = 3,  max_zoom = 12, 
        drop_rate = 2.5, cluster_distance = 1, base_zoom = 11)

## Creating MVT tiles (zoom 3-12) for 1013971 features across 1 layer...

##   Parsed 1013971 features across 1 layer in 4.4s
##   Zoom  3/12:      1 tiles ...
##                 1 encoded (0.2s)
##   Zoom  4/12:      4 tiles ...
##                 4 encoded (0.3s)
##   Zoom  5/12:     10 tiles ...
##                10 encoded (0.3s)
##   Zoom  6/12:     24 tiles ...
##                24 encoded (0.3s)
##   Zoom  7/12:     69 tiles ...
##                69 encoded (0.3s)
##   Zoom  8/12:    190 tiles ...
##               189 encoded (0.4s)
##   Zoom  9/12:    520 tiles ...
##               505 encoded (0.5s)
##   Zoom 10/12:   1304 tiles ...
##              1275 encoded (0.4s)
##   Zoom 11/12:   3163 tiles ...
##              3106 encoded (0.6s)
##   Zoom 12/12:   7448 tiles ...
##              7308 encoded (2.6s)
##   Total: 12491 tiles in 5.9s
##   Writing PMTiles archive (12491 tiles) ...
##   PMTiles write: 4.8s

## Created /home/glaroc/Repos/QCBS-Repos/space-clouds-workshop/pmtiles/atlas_inat.pmtiles (548.9 MB)

## View with: view_tiles("atlas_inat.pmtiles")

categories <- unique(atlas_inat$kingdom)

random_colors <- paste0("#", sapply(1:length(categories), function(i) {
  paste0(sprintf("%02X", sample(0:255, 3)), collapse = "")
}))

Display this on a map using the mapgl package. The serve_tiles function is used to serve the PMTiles file locally.

library(mapgl)

serve_tiles("pmtiles", port = 8080)

## Serving tiles at http://localhost:8080

##   Directory: pmtiles

## Use stop_server() to stop

maplibre(hash = TRUE, 
  style = carto_style("voyager"),
  center = c(-72, 55),
  zoom = 4) |>
  add_pmtiles_source(
    id = "muni",
    url = "http://localhost:8080/muni_qc.pmtiles",
    promote_id = "fid"
  ) |>
  add_pmtiles_source(
    id = "atlas_inat",
    url = "http://localhost:8080/atlas_inat.pmtiles",
  ) |>
  add_fill_layer(
    id = "muni_layer",
    source = "muni",
    source_layer = "municipalities",
    fill_color = "#ececec",
    fill_opacity = 0.5,
    fill_outline_color = "#333333",
    hover_options = list(
      fill_color = "#ffffcc",
      fill_opacity = 0.9
    )
  ) |> 
  add_circle_layer(
    id = "atlas_inat_layer",
    source = "atlas_inat",
    source_layer = "iNat",
    circle_radius = 3,
    circle_opacity = 0.7,
    circle_color = match_expr(
      column = "kingdom",
      values = categories,
      stops  = random_colors
    ),
  ) |> 
  add_categorical_legend(
    legend_title = "kingdom",
    values       = categories,
    colors       = random_colors
  )

View a COG using titiler from Biodiversité Québec. We need to specify the rescale and colormap_name parameters in the URL to visualize the data properly.

maplibre(hash = TRUE, 
  style = carto_style("dark-matter"),
  center = c(-72, 55),
  zoom = 4) |>
  add_raster_source(
        id = "ndvi",
        tileSize = 256,
        tiles = c("https://tiler.biodiversite-quebec.ca/cog/tiles/{z}/{x}/{y}?url=https://object-arbutus.cloud.computecanada.ca/bq-io/io/NDVImax/NDVImax_Quebec_2022.tif&rescale=0,1&colormap_name=speed")
  ) |> 
  add_raster_layer(
    id = "ndvi_layer",
    source = "ndvi",
    raster_opacity = 0.7
  )