[parquet-geospatial] GeospatialStatistics mislabels PolyhedralSurface WKB as GeometryCollection, and mis-bounds Triangle
- Dominant language
- Rust
- Stars
- 3.6k
- Forks
- 1.3k
- Avg merge
- 2d 16h
- Merged PRs (30d)
- 168
Description
### Describe the bug
`parquet-geospatial`'s statistics accumulator derives the `geospatial_types` code by
re-deriving it through `geo-traits` rather than reading it from the WKB header. Because
the underlying `wkb` reader classifies a geometry with `code & 0x7`, ISO WKB codes outside
the seven basic types are silently folded onto one of those seven instead of being
rejected:
| Input | `code & 0x7` | Classified as | Result |
| --- | --- | --- | --- |
| `PolyhedralSurface Z` (1015) | 7 | `GeometryCollection Z` | statistics say `[1007]`; bbox correct |
| `Triangle Z` (1017) | 1 | `Point Z` | statistics say `[1001]`; **bbox is garbage** |
| `TIN Z` (1016) | 0 | — | rejected with `WKB type code out of range` |
The `PolyhedralSurface` case writes a false `geospatial_types` into the file. The
`Triangle` case is worse: the reader consumes the ring counts as coordinate bytes, so the
bounding box written to `GeospatialStatistics` does not contain the geometry. Since the
bbox is exactly what row-group pruning consults, that silently drops matching rows.
`PolyhedralSurface` parses "successfully" and bounds correctly by coincidence — its WKB
body is layout-identical to a `GeometryCollection` of `Polygon`s, so walking it as a
collection visits the same coordinates. `Triangle` has no such coincidence.
The relevant code is [`geometry_type()` in `parquet-geospatial/src/bounding.rs`](https://github.com/apache/arrow-rs/blob/e1e71ea89b211758e20254296008f96e744ed138/parquet-geospatial/src/bounding.rs#L286-L323),
whose own doc comment already notes the alternative:
```rust
/// This can also be derived from bytes 2-5 (possibly endian-swapped according to byte 1)
/// of the input WKB buffer but is slightly clearer recomputed.
fn geometry_type(geom: &impl GeometryTrait) -> Result {
```
The `_ => Err(...)` arm in that match is unreachable for these inputs, because
`geo-traits` has already been told the geometry is a `GeometryCollection` / `Point`.
### To Reproduce
`parquet = { version = "58", features = ["geospatial"] }`, `parquet-geospatial = "58"`,
`arrow-array = "58"`, `arrow-schema = "58"`, `bytes = "1"` — resolving to `parquet` 58.4.0,
`parquet-geospatial` 58.4.0, `wkb` 0.9.2, `geo-traits` 0.3.0.
```rust
use std::sync::Arc;
use arrow_array::{ArrayRef, BinaryArray, RecordBatch};
use arrow_schema::{DataType, Field, Schema};
use parquet::arrow::ArrowWriter;
use parquet::file::reader::{FileReader, SerializedFileReader};
use parquet_geospatial::bounding::GeometryBounder;
use parquet_geospatial::interval::IntervalTrait;
use parquet_geospatial::{WkbMetadata, WkbType};
/// ISO WKB header: little-endian byte order + a 4-byte geometry type code.
fn hdr(buf: &mut Vec, code: u32) {
buf.push(0x01);
buf.extend_from_slice(&code.to_le_bytes());
}
/// One ISO `Polygon Z` (code 1003) with a single exterior ring.
fn polygon_z(ring: &[(f64, f64, f64)]) -> Vec {
let mut b = Vec::new();
hdr(&mut b, 1003);
b.extend_from_slice(&1u32.to_le_bytes()); // numRings
b.extend_from_slice(&(ring.len() as u32).to_le_bytes());
for (x, y, z) in ring {
b.extend_from_slice(&x.to_le_bytes());
b.extend_from_slice(&y.to_le_bytes());
b.extend_from_slice(&z.to_le_bytes());
}
b
}
/// One ISO `PolyhedralSurface Z` (code 1015) built from `Polygon Z` members.
fn polyhedral_surface_z(polys: &[Vec]) -> Vec {
let mut b = Vec::new();
hdr(&mut b, 1015);
b.extend_from_slice(&(polys.len() as u32).to_le_bytes());
for p in polys {
b.extend_from_slice(p);
}
b
}
fn main() -> Result<(), Box> {
// Two triangles spanning x[0,1] y[0,1] z[0,2].
let wkb = polyhedral_surface_z(&[
polygon_z(&[(0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 0.0)]),
polygon_z(&[(0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 0.0, 2.0), (0.0, 0.0, 0.0)]),
]);
println!("WKB header type code = {}", u32::from_le_bytes(wkb[1..5].try_into()?));
// ---- 1. the accumulator directly -------------------------------------
let mut bounder = GeometryBounder::empty();
bounder.update_wkb(&wkb)?;
println!(
"GeometryBounder::geometry_types() = {:?} (expected [1015])",
bounder.geometry_types()
);
println!(
"GeometryBounder bbox = x{:?} y{:?} z{:?}",
(bounder.x().lo(), bounder.x().hi()),
(bounder.y().lo(), bounder.y().hi()),
(bounder.z().lo(), bounder.z().hi())
);
// ---- 2. end to end through ArrowWriter --------------------------------
let field = Field::new("geom", DataType::Binary, true)
.with_extension_type(WkbType::new(Some(WkbMetadata::new(Some("EPSG:7415"), None))));
let schema = Arc::new(Schema::new(vec![field]));
let array: ArrayRef = Arc::new(BinaryArray::from(vec![Some(wkb.as_slice())]));
let batch = RecordBatch::try_new(schema.clone(), vec![array])?;
let mut buf: Vec = Vec::new();
let mut writer = ArrowWriter::try_new(&mut buf, schema, None)?;
writer.write(&batch)?;
writer.close()?;
let reader = SerializedFileReader::new(bytes::Bytes::from(buf))?;
let meta = reader.metadata();
let col = meta.row_group(0).column(0);
println!("Parquet logical type = {:?}", col.column_descr().logical_type_ref());
println!("GeospatialStatistics = {:?}", col.geo_statistics());
Ok(())
}
```
Actual output:
```
WKB header type code = 1015
GeometryBounder::geometry_types() = [1007] (expected [1015])
GeometryBounder bbox = x(0.0, 1.0) y(0.0, 1.0) z(0.0, 2.0)
Parquet logical type = Some(Geometry { crs: Some("\"EPSG:7415\"") })
GeospatialStatistics = Some(GeospatialStatistics { bbox: Some(BoundingBox { x_range: (0.0, 1.0), y_range: (0.0, 1.0), z_range: Some((0.0, 2.0)), m_range: None }), geospatial_types: Some([1007]) })
```
The `Triangle` / `TIN` rows in the table above come from the same harness, feeding
`GeometryBounder::update_wkb` a hand-built `Triangle Z` whose true bounds are
`x(10,11) y(20,21) z(30,30)`:
```
Triangle Z (1017) -> types [1001] bbox x(8.4879831644e-314, 8.4879831644e-314) y(10.0, 10.0) z(20.0, 20.0)
TIN Z (1016) -> Err: External error: General error: WKB type code out of range. Got: 1016
```
### Expected behavior
For `PolyhedralSurface Z`, either of these would be correct; the current behaviour is the
one clearly wrong option, because it states something false about the data:
1. **Report the header code**, i.e. `[1015]`. [`Geospatial.md`](https://github.com/apache/parquet-format/blob/master/Geospatial.md#geospatial-types)
defines `geospatial_types` as "[WKB (ISO-variant) integer codes]" and presents its table
as "the most common geospatial types", not a closed vocabulary — so 1015 appears
admissible. This is what the existing `bytes 2-5` comment suggests, and it is a local
change to `geometry_type()`.
2. **Reject the geometry**, so the caller can fall back to omitting statistics. Note that
`update_wkb`'s own doc comment already anticipates this: "clients may wish to ignore
such an error for the purposes of writing statistics."
For `Triangle Z`, only rejection or correct bounding is acceptable — the bbox currently
written is not a superset of the geometry, which breaks the pruning contract.
### Additional context
The root cause is in [`georust/wkb`](https://github.com/georust/wkb), whose
`WkbGeometryCode::get_type` masks with `code & 0x7` and so never rejects a code outside
the basic seven. That is not something `arrow-rs` can fix directly, which is why deriving
the type code from the WKB header inside `geometry_type()` looks like the right fix at
this layer. Whether `parquet-geospatial` should additionally validate the header code
before handing the buffer to `Wkb::try_new` is the part I am unsure about — the
`Triangle` bbox result suggests it might have to.
I found this writing `PolyhedralSurface Z` columns from
[CityParquet](https://github.com/cityjson/cityparquet), a Parquet encoding for 3D city
models where solid geometry is the common case. The bbox being exact means this is not
currently blocking us — we consume the bounding box and ignore `geospatial_types` — but
the file does carry a false type code, and readers that trust it will be misled.
I'm happy to open a PR for the header-derived `geometry_type()` if that is the direction
you'd prefer.
Contributor guide
Research direction
Start with geometry_type() in parquet-geospatial/src/bounding.rs and trace GeometryBounder::update_wkb, then run the provided PolyhedralSurface, Triangle, and TIN harness. Check the existing bounding and geospatial statistics tests for coverage. Done means unsupported WKB geometries cannot produce misleading type codes or bounding boxes, and the resulting statistics follow one of the stated expected behaviors.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- rust
- Domain
- data-engineering
- Issue type
- Bug
- Difficulty
- 3/5
- Estimated time
- 1-2 days
- Activity status
- Active
- Clarity
- Mostly clear
- Newbie friendliness
- 68/100