library(dplyr)
library(sf)
#> Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(terra)
#> There appears to be a problem with the PROJ installation
#> terra 1.8.93
library(gdalraster)
#> GDAL 3.8.5 (released 2024-04-02), GEOS 3.13.0, PROJ 9.5.1
library(stars)
earthdatalogin::gdal_cloud_config()
idx_url <- "https://nrs.objectstore.gov.bc.ca/sihgte/FAIB_PLFP_Index/FAIB_ALS_PROGRESS_bcgs_20k_grid.zip"
## Read the index file, use /vsizip//vsicurl/ to read directly from the zip file on the web
idx <- st_read(paste0("/vsizip//vsicurl/", idx_url))
#> Reading layer `FAIB_ALS_PROGRESS_bcgs_20k_grid' from data source
#> `/vsizip//vsicurl/https://nrs.objectstore.gov.bc.ca/sihgte/FAIB_PLFP_Index/FAIB_ALS_PROGRESS_bcgs_20k_grid.zip'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 232 features and 17 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 774674.7 ymin: 356869.5 xmax: 1375855 ymax: 1367392
#> Projected CRS: NAD83(CSRS) / BC Albers
## Filter the index for the desired tile and get the URL of the data file
data_url <- filter(idx, MAP_TILE == "092C066") |>
pull(objstore)
## List the files in the zip archive using gdalraster. Again, use /vsizip//vsicurl/ to read directly from the web
files_in_zip <- gdalraster::vsi_read_dir(paste0("/vsizip//vsicurl/", data_url))
files_in_zip
#> [1] "bc_092C066_xli_bcalb_20231002_20231030_bem_chm_202509290323.tif"
#> [2] "bc_092C066_xli_bcalb_20231002_20231030_csm_202509290324.tif"
#> [3] "bc_092C066_xli_bcalb_20231002_20231030_csm_202509290324_names.csv"
#> [4] "bc_092C066_xli_bcalb_20231002_20231030_treetops_202509290324.cpg"
#> [5] "bc_092C066_xli_bcalb_20231002_20231030_treetops_202509290324.dbf"
#> [6] "bc_092C066_xli_bcalb_20231002_20231030_treetops_202509290324.prj"
#> [7] "bc_092C066_xli_bcalb_20231002_20231030_treetops_202509290324.shp"
#> [8] "bc_092C066_xli_bcalb_20231002_20231030_treetops_202509290324.shx"
#> [9] "check_ground_092c066.qmd"
#> [10] "check_ground_092c066.sqlite"
## Read the BEM/CHM .tif file from the zip archive using stars.
## Use proxy = TRUE to read it as a proxy object, to avoid loading the entire raster into memory at once
bem_chm <- stars::read_stars(
paste0("/vsizip//vsicurl/", data_url, "/", files_in_zip[1]),
proxy = TRUE
)
plot(bem_chm)
#> downsample set to 33
## Or with terra
bem_chm_terra <- terra::rast(paste0(
"/vsizip//vsicurl/",
data_url,
"/",
files_in_zip[1]
))
bem_chm_terra
#> class : SpatRaster
#> size : 7218, 14934, 2 (nrow, ncol, nlyr)
#> resolution : 1, 1 (x, y)
#> extent : 1073734, 1088668, 402854, 410072 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(CSRS) / BC Albers
#> source : bc_092C066_xli_bcalb_20231002_20231030_bem_chm_202509290323.tif
#> names : BEM, CHM
#> min values : -2.39, 0.00
#> max values : 375.73, 81.33
plot(bem_chm_terra$BEM)
Created on 2026-02-06 with reprex v2.1.1