HideBa opened a new issue, #10822:
URL: https://github.com/apache/arrow-rs/issues/10822

   ### 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<T = f64>) -> Result<i32, 
ArrowError> {
   ```
   
   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<u8>, 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<u8> {
       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<u8>]) -> Vec<u8> {
       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<dyn std::error::Error>> {
       // 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<u8> = 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.
   


-- 
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.

To unsubscribe, e-mail: [email protected]

For queries about this service, please contact Infrastructure at:
[email protected]

Reply via email to