Geoparquet in R – Marc Weber

Geoparquet extends the parquet file format and as nicely described in this blog post by Kyle Barron and provides a powerful new way to store and share geospatial data in a cloud-optimized format. I’ve been using it for more and more of my spatial data, and below is just a quick example using Overture Maps buildings data and doing spatial and attribute filtering to subset the data prior to reading in.

First open a connection to a cloud-hosted GeoParquet file

For this example we use Overture Maps buildings (public S3 bucket). We’ll open a connection (but we are not actually reading it in yet)

library(arrow)
library(dplyr)
library(sf)
library(sfarrow)
library(jsonlite)

latest <- fromJSON("https://stac.overturemaps.org/catalog.json")$latest

ds_path <- s3_bucket(
  paste0("overturemaps-us-west-2/release/", latest, "/theme=buildings/type=building/"),
  region = "us-west-2",
  anonymous = TRUE
)

ds <- open_dataset(ds_path, format = "parquet")
buildings_ds <- open_dataset(ds_path, format = "parquet")
# Inspect available columns to confirm tile partitioning
print(buildings_ds$schema$names)  # look for "z", "x", "y" as partition columns

 [1] "id"                     "names"                  "sources"               
 [4] "level"                  "height"                 "min_height"            
 [7] "is_underground"         "num_floors"             "num_floors_underground"
[10] "min_floor"              "subtype"                "class"                 
[13] "facade_color"           "facade_material"        "roof_material"         
[16] "roof_shape"             "roof_direction"         "roof_orientation"      
[19] "roof_color"             "roof_height"            "geometry"              
[22] "has_parts"              "version"                "bbox"

bbox_field <- buildings_ds$schema$GetFieldByName("bbox")

Look at a slice to see structure of data

collect(
  buildings_ds %>%
    select(id, bbox) %>%
    slice_head(n = 5)
)

# A tibble: 5 × 2
  id                                   bbox$xmin $xmax $ymin $ymax
  <chr>                                    <dbl> <dbl> <dbl> <dbl>
1 76c4a544-9ad3-4da0-be7d-22892dcbeb58     -180. -180. -78.0 -78.0
2 73d1649d-7c29-4307-9563-96989c159e02     -180. -180. -19.0 -19.0
3 4c1671f8-9fa9-4a44-b240-62febf4f1b34     -180. -180. -19.0 -19.0
4 c4f97f7d-49bd-4e20-96cc-1f4e99ee8f84     -180. -180. -19.0 -19.0
5 63e18425-986f-4402-9a5b-1ef3c9aaaa68     -180. -180. -19.0 -19.0

Materialize a bbox as columns for inspection

# Get the Field and its type
bbox_field <- buildings_ds$schema$GetFieldByName("bbox")
print(bbox_field)

Field
bbox: struct<xmin: double, xmax: double, ymin: double, ymax: double>

print(bbox_field$type)

StructType
struct<xmin: double, xmax: double, ymin: double, ymax: double>

Prune the data

We define your query bbox in EPSG:4326 lon/lat

xmin <- -123.32; ymin <- 44.52
xmax <- -123.20; ymax <- 44.63
bbox_sfc <- st_as_sfc(st_bbox(c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax), crs = 4326))


ds_pruned <- buildings_ds %>%
  filter(
    bbox$xmin <= xmax,      # feature xmin <= query xmax
    bbox$xmax >= xmin,      # feature xmax >= query xmin
    bbox$ymin <= ymax,      # feature ymin <= query ymax
    bbox$ymax >= ymin,      # feature ymax >= query ymin
    !is.na(height), height >= 4,
    !is.na(num_floors), num_floors >= 2
  ) %>%
  # Keep simple scalar columns plus geometry (avoid nested list/struct columns like 'names', 'sources', 'has_parts')
  select(id, geometry, height, num_floors, subtype, class)

Materialize and convert to sf

# Materialize and convert to sf
tbl <- collect(ds_pruned)
geom_sfc <- st_as_sfc(tbl$geometry, crs = 4326)
buildings_sf <- st_sf(tbl[, setdiff(names(tbl), "geometry")], geometry = geom_sfc)

# Optional: exact bbox cut in sf for precise geometry clipping
bbox_sfc <- st_as_sfc(st_bbox(c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax), crs = 4326))
buildings_bbox <- st_filter(buildings_sf, bbox_sfc)

mapview::mapview(buildings_bbox)