Skip to content

Instantly share code, notes, and snippets.

@h-a-graham
Last active December 9, 2025 15:03
Show Gist options
  • Select an option

  • Save h-a-graham/f1bddd6f8fc192d14d7d91cea9acffa8 to your computer and use it in GitHub Desktop.

Select an option

Save h-a-graham/f1bddd6f8fc192d14d7d91cea9acffa8 to your computer and use it in GitHub Desktop.
Generate an a5 grid ysing R and duckdb
#' Generate an A5 grid over a specified bounding box or geometry
#' @param x Either a numeric vector of length 4 specifying a bounding box
#' (xmin, ymin, xmax, ymax), or a WKT string or WKB blob specifying a geometry
#' to cover with the A5 grid.
#' @param resolution Integer from 0 (lowest resolution) to 30
#' (highest resolution)
#' @param output Character string specifying the output Parquet file path
#' to save the generated A5 grid.
#' @return Character string specifying the output Parquet file path
a5_grid <- function(x, resolution, output) {
assertthat::assert_that(
inherits(resolution, "numeric")
)
assertthat::assert_that(resolution >= 0 & resolution <= 30)
assertthat::assert_that((resolution %% 1) == 0)
assertthat::assert_that(
inherits(x, "numeric") || inherits(x, "character") || inherits(x, "blob")
)
if (inherits(x, "numeric")) {
assertthat::assert_that(length(x) == 4)
}
assertthat::assert_that(inherits(output, "character"))
assertthat::assert_that(fs::path_ext(output) == "parquet")
if (inherits(x, "character")) {
bbox <- gdalraster::bbox_from_wkt(x)
} else if (inherits(x, "blob")) {
bbox <- gdalraster::g_wk2wk(x) |> gdalraster::bbox_from_wkt()
} else {
bbox <- x
}
# Build WHERE clause for spatial intersection if geometry provided
if (!inherits(x, "numeric")) {
if (inherits(x, "blob")) {
x <- gdalraster::g_wk2wk(x)
}
where_clause <- glue::glue(
"WHERE ST_Intersects(geometry, ST_GeomFromText('{x}'))"
)
} else {
where_clause <- ""
}
con <- DBI::dbConnect(duckdb::duckdb(), dbdir = ":memory:")
on.exit(DBI::dbDisconnect(con, shutdown = TRUE))
# Install and load extensions
DBI::dbExecute(
con,
"INSTALL a5 FROM community;
LOAD a5;
INSTALL spatial;
LOAD spatial;
"
)
sql <- glue::glue(
"
COPY (
WITH params AS (
SELECT CASE
WHEN {resolution} = 0 THEN 35.0
WHEN {resolution} = 1 THEN 18.0
WHEN {resolution} = 2 THEN 10.0
WHEN {resolution} = 3 THEN 5.0
ELSE 5.0 * POWER(0.5, {resolution} - 3)
END AS step,
{bbox[1]} AS xmin,
{bbox[3]} AS xmax,
{bbox[2]} AS ymin,
{bbox[4]} AS ymax
),
xs AS (
SELECT xmin + (i * step) AS x
FROM params,
unnest(range(0, CAST(CEIL((xmax - xmin) / step) + 1 AS INTEGER))) AS t(i)
),
ys AS (
SELECT ymin + (i * step) AS y
FROM params,
unnest(range(0, CAST(CEIL((ymax - ymin) / step) + 1 AS INTEGER))) AS t(i)
),
cells AS (
SELECT DISTINCT a5_lonlat_to_cell(x, y, {resolution}) AS cell
FROM xs CROSS JOIN ys
)
SELECT cell,
ST_MakePolygon(
ST_MakeLine(
list_transform(
a5_cell_to_boundary(cell),
pt -> ST_Point(pt[1], pt[2])
)
)
) AS geometry
FROM cells
{where_clause}
) TO '{output}' WITH (FORMAT 'parquet', COMPRESSION 'zstd');
"
)
DBI::dbExecute(con, sql)
return(output)
}
a5_grid(c(-179, -89, 179, 89), 2, tempfile(fileext = ".parquet")) |>
sf::read_sf() |>
plot()
@h-a-graham
Copy link
Author

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment