Reference Grid is required for `grid_grab`

[slacasm]

Excellent update to make this also pull Fuel Types from the national layer for a simple use-case. I would suggest an additional function be created that enables a user to build a reference grid from an input polygon. I have provided a template function here that does this, though it is weighty and was largely done with co-pilot, so caution is urged:

`suppressPackageStartupMessages({
library(sf)
library(terra)
library(stars)
})

.snap_extent <- function(ext, res, origin) {
xmin <- floor((ext[1] - origin[1]) / res[1]) * res[1] + origin[1]
ymin <- floor((ext[3] - origin[2]) / res[2]) * res[2] + origin[2]
xmax <- ceiling((ext[2] - origin[1]) / res[1]) * res[1] + origin[1]
ymax <- ceiling((ext[4] - origin[2]) / res[2]) * res[2] + origin[2]
c(xmin, xmax, ymin, ymax)
}

infer_grid_from_raster <- function(reference_grid) {
r <- if (inherits(reference_grid, "SpatRaster")) reference_grid else rast(reference_grid)
list(
crs = crs(r, proj = TRUE),
res = res(r),
origin = terra::origin(r)
)
}

infer_grid_from_wcs <- function(aoi_sf, wcs_base, coverage_id) {
aoi_ll <- st_transform(aoi_sf, 4326)
cxy <- st_coordinates(st_centroid(st_union(aoi_ll)))[1,]
dx <- 0.005; dy <- 0.005
url <- paste0(
wcs_base,
"?SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage",
"&coverageId=", utils::URLencode(coverage_id, reserved = TRUE),
"&subset=Long(", cxy[1]-dx, ",", cxy[1]+dx, ")",
"&subset=Lat(", cxy[2]-dy, ",", cxy[2]+dy, ")"
)
cov <- try(read_stars(url, proxy = FALSE), silent = TRUE)
if (inherits(cov, "try-error")) stop("WCS probe failed; check URL, version, or coverageId.")
r <- rast(cov)
list(
crs = crs(r, proj = TRUE),
res = res(r),
origin = terra::origin(r)
)
}

make_reference_grid_from_aoi <- function(
aoi_path,
out_dir,
out_stem = "aoi_grid",
buffer_m = 0,
reference_grid = NULL,
wcs_base = NULL,
wcs_coverage = NULL,
crs_fallback = "EPSG:3347",
cellsize_fallback = 100,
grid_type = c("polygons","centroids","both"),
clip_to_aoi = FALSE
) {
grid_type <- match.arg(grid_type)
if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)

aoi <- st_read(aoi_path, quiet = TRUE)
if (!inherits(aoi, "sf")) stop("AOI read failed.")
if (!any(sf::st_geometry_type(aoi) %in% c("POLYGON","MULTIPOLYGON"))) {
stop("AOI must contain polygon geometry.")
}

spec <- NULL
if (!is.null(reference_grid)) {
spec <- infer_grid_from_raster(reference_grid)
} else if (!is.null(wcs_base) && !is.null(wcs_coverage)) {
spec <- infer_grid_from_wcs(aoi, wcs_base, wcs_coverage)
} else {
spec <- list(crs = crs_fallback, res = c(cellsize_fallback, cellsize_fallback), origin = c(0,0))
}

Ensure origin is numeric length-2

if (is.null(spec$origin) || any(is.na(spec$origin))) spec$origin <- c(0,0)

aoi_tgt <- st_transform(aoi, spec$crs)
if (buffer_m != 0) {
aoi_tgt <- st_buffer(st_union(aoi_tgt), dist = buffer_m)
aoi_tgt <- st_as_sf(aoi_tgt)
} else {
aoi_tgt <- st_as_sf(st_union(aoi_tgt))
}

bb <- st_bbox(aoi_tgt)
ext0 <- c(bb["xmin"], bb["xmax"], bb["ymin"], bb["ymax"])
ext1 <- .snap_extent(ext0, res = spec$res, origin = spec$origin)

Build raster template

ext_r <- terra::ext(ext1[1], ext1[2], ext1[3], ext1[4]) # xmin, xmax, ymin, ymax
r <- terra::rast(ext = ext_r, resolution = spec$res, crs = spec$crs)
values(r) <- NA_real_

out_list <- list(template = r)

=== robust bbox polygon from raster extent ===

e <- terra::ext(r)
stopifnot(!any(is.na(c(e$xmin, e$xmax, e$ymin, e$ymax))))
coords <- matrix(
c(
e$xmin, e$ymin,
e$xmax, e$ymin,
e$xmax, e$ymax,
e$xmin, e$ymax,
e$xmin, e$ymin
),
ncol = 2, byrow = TRUE
)
bbox_sfc <- sf::st_sfc(sf::st_polygon(list(coords)), crs = sf::st_crs(spec$crs))

if (grid_type %in% c("polygons","both")) {
grd_poly <- sf::st_make_grid(
x = bbox_sfc,
cellsize = spec$res,
what = "polygons",
square = TRUE
) |>
sf::st_as_sf()
if (clip_to_aoi) grd_poly <- suppressWarnings(sf::st_intersection(grd_poly, aoi_tgt))
grd_poly <- sf::st_collection_extract(grd_poly, "POLYGON", warn = FALSE)
out_list$grid_polygons <- grd_poly
}

if (grid_type %in% c("centroids","both")) {
tmp_poly <- if (!"grid_polygons" %in% names(out_list)) {
tp <- sf::st_make_grid(
x = bbox_sfc,
cellsize = spec$res,
what = "polygons",
square = TRUE
) |>
sf::st_as_sf()
if (clip_to_aoi) tp <- suppressWarnings(sf::st_intersection(tp, aoi_tgt))
sf::st_collection_extract(tp, "POLYGON", warn = FALSE)
} else {
out_list$grid_polygons
}
cents <- sf::st_centroid(tmp_poly)
out_list$grid_centroids <- cents
}

tmpl_path <- file.path(out_dir, paste0(out_stem, "_template.tif"))
writeRaster(r, tmpl_path, overwrite = TRUE)

gpkg_path <- file.path(out_dir, paste0(out_stem, ".gpkg"))
if ("grid_polygons" %in% names(out_list)) {
st_write(out_list$grid_polygons, gpkg_path, layer = "grid_polygons", delete_layer = TRUE, quiet = TRUE)
}
if ("grid_centroids" %in% names(out_list)) {
st_write(out_list$grid_centroids, gpkg_path, layer = "grid_centroids", delete_layer = TRUE, quiet = TRUE)
}
st_write(st_as_sf(aoi_tgt), gpkg_path, layer = "aoi_used", delete_layer = TRUE, quiet = TRUE)

message("Wrote: ", tmpl_path)
if (file.exists(gpkg_path)) message("Wrote: ", gpkg_path)

invisible(out_list)
}
`