Raster to a5 grid in R

raster-to-a5.R

# ---- raster zonal stats example ----
# Extract mean elevation per A5 cell from an ALOS DEM COG tile
# using GDAL's zonal-stats algorithm via gdalraster.

# pak::pkg_install("pak::pak("belian-earth/a5R")")
library(a5R)
library(gdalraster)

# Open a cloud-optimised GeoTIFF via /vsicurl (Planetary Computer)
src <- "/vsicurl?pc_url_signing=yes&url=https://ai4edataeuwest.blob.core.windows.net/alos-dem/AW3D30_global/ALPSMLC30_N047E010_DSM.tif"
rast <- new(GDALRaster, src)
rast_wgs84 <- autoCreateWarpedVRT(
  rast,
  dst_wkt = epsg_to_wkt(4326),
  resample_alg = "NearestNeighbour"
)

# Generate A5 grid covering the raster extent
cells <- a5_grid(rast_wgs84$bbox(), resolution = 12)

# Write cell boundaries to a temp Parquet file for GDAL
zones_file <- tempfile(fileext = ".parquet")
defn <- ogr_def_layer("Polygon", srs = epsg_to_wkt(4326))
defn$cell_id <- ogr_def_field("OFTString")
lyr <- ogr_ds_create(
  "PARQUET",
  zones_file,
  layer = "zones",
  layer_defn = defn,
  return_obj = TRUE,
  srs = epsg_to_wkt(4326)
)
lyr$batchCreateFeature(data.frame(
  cell_id = as.character(cells),
  geom = as.character(a5_cell_to_boundary(cells, format = "wkt"))
))
lyr$close()

# Run zonal stats
out_file <- tempfile(fileext = ".parquet")
alg <- gdal_run(
  "raster zonal-stats",
  list(
    input = rast_wgs84,
    zones = zones_file,
    output = out_file,
    format = "PARQUET",
    stat = "mean",
    strategy = "raster",
    "include-field" = "cell_id",
    "skip-errors" = TRUE,
    overwrite = TRUE
  )
)
alg$release()


# Read results and plot
df <- arrow::read_parquet(out_file)
df$geometry <- a5_cell_to_boundary(df$cell_id)

pal <- hcl.colors(100, "Inferno")
cols <- pal[as.numeric(cut(df$mean, 100))]
mean_lat <- mean(rast_wgs84$bbox()[c(2, 4)])
wk::wk_plot(
  df$geometry,
  col = cols,
  border = NA,
  asp = 1 / cos(mean_lat * pi / 180)
)