Read GeoTIFF and TIFF through the nd pipeline

CHANGELOG.md

@@ -45,6 +45,19 @@ tag. Releases before 2.0.0 are recorded in the
 
 ### Changed
 
+- **GeoTIFF and TIFF read through the nd pipeline**, as netCDF, HDF5 and Zarr already did. A raster
+  is a grid over the axes `y` (image rows) and `x` (image columns), and the reader now hands the
+  plan its columns un-broadcast — the bands on `y × x`, `geo.lat` on `y`, `geo.lon` on `x`, and the
+  TIFF tags as scalars — for a decode and a broadcast above the scan. The pixel cross-product is
+  therefore built once, at the top of the plan, instead of per column in the reader, and the two nd
+  optimizer rules now reach a raster: a `WHERE` on a coordinate selects the grid before it is
+  materialized, and an element-wise `SELECT` expression on a coordinate is evaluated on that axis
+  alone rather than on every pixel. `read_tiff` also takes the optional `dimensions` argument the
+  other nd readers take: `read_tiff('raster.tif', ['y'])` returns one row per image row. The same
+  list is accepted as `OPTIONS (read_dimensions 'y')` on `CREATE EXTERNAL TABLE … STORED AS TIFF`.
+  As for every nd format, the projected columns define the grid, so selecting only 1-d coordinate
+  columns returns their own axes rather than the full pixel grid. See
+  [GeoTIFF](docs/docs/2.0.0-rc2/formats/geotiff.md).
 - **Minimum supported Rust is 1.94**, up from 1.91. `iceberg` and `iceberg-datafusion` 0.10 — the
   only release line built against the DataFusion 53 and Arrow 58 this workspace unifies on —
   declare `rust-version = "1.94"`, so the workspace floor follows. Beacon's own code uses no

beacon-db/beacon-core/tests/read_functions.rs

@@ -34,6 +34,8 @@ fn seed(rt: TestRuntime) -> TestRuntime {
         .expect("copy netcdf fixture as hdf5");
     std::fs::copy(nested_hdf5_fixture(), rt.datasets_dir().join("nested.h5"))
         .expect("copy nested hdf5 fixture");
+    std::fs::copy(tiff_fixture(), rt.datasets_dir().join("raster.tif"))
+        .expect("copy geotiff fixture");
     rt
 }
 
@@ -62,6 +64,19 @@ fn nested_hdf5_fixture() -> std::path::PathBuf {
         .join("beacon-file-formats/beacon-arrow-hdf5/test_files/nested-groups.h5")
 }
 
+/// A stripped single-band GeoTIFF, shipped with the TIFF reader: 1287 x 380
+/// float32 pixels on the axes `x` and `y`.
+fn tiff_fixture() -> std::path::PathBuf {
+    std::path::Path::new(env!("CARGO_MANIFEST_DIR"))
+        .parent()
+        .unwrap()
+        .join("beacon-file-formats/beacon-arrow-tiff/test-files/test.tif")
+}
+
+/// The GeoTIFF fixture's grid: `y` (image rows) x `x` (image columns).
+const TIFF_HEIGHT: i64 = 380;
+const TIFF_WIDTH: i64 = 1287;
+
 #[tokio::test(flavor = "multi_thread")]
 async fn read_csv_scans_filters_and_projects() {
     let rt = seeded("read-csv").await;
@@ -197,6 +212,68 @@ async fn read_hdf5_schema_lists_the_nested_columns() {
     );
 }
 
+/// A GeoTIFF rides the same nd pipeline as netCDF, HDF5 and zarr: the raster is
+/// a `y` x `x` grid, and the 1-d coordinate axes broadcast over it.
+#[tokio::test(flavor = "multi_thread")]
+async fn read_tiff_scans_the_fixture_as_a_grid() {
+    let rt = seeded("read-tiff").await;
+
+    assert_eq!(
+        scalar_i64(&rt.sql("SELECT count(*) FROM read_tiff('raster.tif')").await),
+        TIFF_HEIGHT * TIFF_WIDTH,
+        "the raster is one row per pixel"
+    );
+
+    // `geo.lat` lives on `y` alone, so it is a broadcast column of the grid the
+    // full-rank band establishes. The band is co-selected because, as for every
+    // nd format, the projected columns are what define the grid.
+    let broadcast = rt
+        .sql(
+            r#"SELECT count("geo.lat")          AS lat_rows,
+                      count(DISTINCT "geo.lat") AS lat_values,
+                      count("band.0")           AS band_values
+               FROM read_tiff('raster.tif')"#,
+        )
+        .await;
+    let column = |name: &str| {
+        let index = broadcast[0].schema().index_of(name).expect(name);
+        arrow::array::as_primitive_array::<arrow::datatypes::Int64Type>(broadcast[0].column(index))
+            .value(0)
+    };
+    assert_eq!(column("lat_rows"), TIFF_HEIGHT * TIFF_WIDTH);
+    assert_eq!(column("lat_values"), TIFF_HEIGHT);
+    // The band's nodata pixels come back as nulls, so it counts fewer.
+    assert!(column("band_values") > 0);
+    assert!(column("band_values") < column("lat_rows"));
+}
+
+/// The optional second argument sets the grid for a raster too: `['y']` keeps
+/// the latitude axis and drops the band, which needs both axes.
+#[tokio::test(flavor = "multi_thread")]
+async fn read_tiff_takes_a_dimensions_argument() {
+    let rt = seeded("read-tiff-dimensions").await;
+
+    let narrowed = rt
+        .sql("SELECT * FROM read_tiff(['raster.tif'], ['y'])")
+        .await;
+    let columns: Vec<String> = narrowed[0]
+        .schema()
+        .fields()
+        .iter()
+        .map(|f| f.name().clone())
+        .collect();
+    assert!(columns.contains(&"geo.lat".to_string()), "{columns:?}");
+    assert!(
+        !columns.contains(&"band.0".to_string()),
+        "a 2-d band does not fit a 1-d grid: {columns:?}"
+    );
+    assert_eq!(
+        total_rows(&narrowed) as i64,
+        TIFF_HEIGHT,
+        "one row for each image row"
+    );
+}
+
 #[tokio::test(flavor = "multi_thread")]
 async fn read_of_a_missing_file_is_an_error_not_an_empty_result() {
     let rt = seeded("read-missing").await;

beacon-db/beacon-datafusion-ext/src/fast_object/mod.rs

@@ -38,10 +38,10 @@
 //! # The format still plans its own scan
 //!
 //! `create_physical_plan` is what turns the config into a plan, so every format
-//! keeps the shape it wants: netCDF and HDF5 stack decode and broadcast nodes
-//! over their scan, and Zarr and Atlas expand a store directory into partitions
-//! and reduce it to the marker at its root. Nothing here knows about any of
-//! that.
+//! keeps the shape it wants: netCDF, HDF5, Zarr and GeoTIFF stack decode and
+//! broadcast nodes over their scan, and Zarr and Atlas expand a store directory
+//! into partitions and reduce it to the marker at its root. Nothing here knows
+//! about any of that.
 //!
 //! # What `EXPLAIN` shows
 //!

beacon-db/beacon-datafusion-ext/src/fast_object/table.rs

@@ -281,8 +281,8 @@ impl FastObjectTable {
 /// blocking a later repartition or limit pushdown.
 ///
 /// The scan is found by descending the single-child chain, because an nd format
-/// returns a stack: netCDF and HDF5 hand back decode and broadcast nodes above
-/// their scan.
+/// returns a stack: netCDF, HDF5, Zarr and GeoTIFF hand back decode and
+/// broadcast nodes above their scan.
 fn record_counters(plan: &Arc<dyn ExecutionPlan>, considered: usize, dropped: usize) {
     let mut node: &dyn ExecutionPlan = plan.as_ref();
     loop {