Support for Web Clients?

[mmgeorge]

Hi! We've been working on adding display optimization that builds on top of GeoParquet, with the goal of enabling web clients to read and visualize GeoParquet files directly by streaming via range requests. I was wondering if there was any interest in standardizing any of this into the GeoParquet spec directly.

The approach we're using I think is actually pretty straightforward:

• For points we simply cluster the dataset using Z-order.
• For non-points, we cluster using an XZ-order code derived from feature extents, using the approach outlined in https://doi.org/10.1007/3-540-48482-5_7. Then we generate several multiscale columns that include quantized, delta-encoded geometries to handle different LODs. We write these multiscale columns, and the XZ-key, as an additional structured column.

Because the XZ-code is an actual column, clients can use statistics to query the XZ-code ranges relevant for a given spatial query (like an on-screen display tile). Provided the file was generated with a Page Index (we currently require this), we are seeing pretty decent performance.

This is the full Overture building footprints dataset, (2.7 billion features spread across 16 files, each about ~25 GB):

footprints2.mp4

Overture Building Footprints

And here's a single file containing 2010 US Census Blocks (11 million features in a single 10 GB file):

census2.mp4

Census Blocks

Right now we are storing this provisional "Spatial Optimization" metadata as a key-value pair in addition to the GeoParquet metadata (which covers the full resolution geometry).

I noticed there's some other exploratory work going on in this area (#279). Where this approach potentially differs from that is that while clustering spatially in this way does lead to better spatial query performance by colocating nearby features, it does fundamentally come at the expense of other variables you might want to cluster on. Storing additional LOD columns also means increasing the size of your dataset, so there are some real tradeoffs here. That said, if your data is spatial and you are using spatial filtering or performing spatial analysis, I imagine standard query engine predicate skipping should be able to benefit from this as well, provided we also write per-feature bounds that they could take advantage of (right now we are not doing this).

I was following some of the Iceberg / native geospatial types discussions, and I think clustering by XZ-code is something that was being explored there as well? So it seems like there might be some standardization possible here. I'm also planning to look at spatio-temporal clustering, something like a Z(XZ,T) for discrete, or XZ(x0, x1, y0, y1, t0, t1) for time intervals.

[Maxxen]

Hello! Im curious how you generate the "quantized delta-encoded geometries"? Is this just regular quantized WKB with DELTA_(LENGTH)_BYTE_ARRAY or something bespoke?

[mmgeorge]

@Maxxen, currently it's a custom binary format using PBF, based on Esri's Feature proto spec:

message Geometry {
  repeated uint32 lengths = 2 [packed = true]; // coordinate structure in lengths
  repeated sint64 coords = 3 [packed = true]; // delta-encoded integer values
}

Geometries are flattened, with each length in the lengths array pointing to the start of a ring. So a polygon like this:

Polygon {
  rings: [ 
    [ [256, 256], [258, 256], [258, 258], [256, 258], [256, 256] ], 
    [ [56, 56], [58, 56], [58, 58], [56, 58], [56, 56] ] 
}

Becomes:

geometry: {
  lengths: [5, 5]
  coords: [256, 256, 2, 0, 0, 2, -2, 0, 0, -2, 56, 56, 2, 0, 0, 2, -2, 0, 0, -2]
}

Since sint64 in PBF is varint encoding, the smaller values for the deltas can use fewer bytes. Some more information here: https://protobuf.dev/programming-guides/encoding/. For the actual quantization part we select several LODs and compute the transform needed to snap vertices to pixels at the LOD. Haven't yet implemented 3D support yet, which would use a slightly different approach since we also need to consider Z values.

I was wondering if we could take a similar approach to the native encoding in 1.1, but using DELTA_BINARY_PACKED. You would need to store the x's & y's in separate columns, and then either denote rings using repetition levels, or using a third column of offsets. The problem with breaking up the geometry into different columns though is that you end up needing to pull in multiple different pages (of different columns) to rehydrate the geometry for a given set of features. Page sizes are determined in bytes, and a page is the smallest granularity that we can download. So separate X & Y columns means that each page covers twice the number of features than if they were combined. You also end up having redundant repetition and definition levels, since those are duplicated for each column. It could still be worth exploring though.

[mmgeorge]

Opened a draft pull request here: #294

[mmgeorge]

Split into two proposals, one for LOD and another for spatial sorting