This is an automated email from the ASF dual-hosted git repository.
kontinuation pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/sedona-db.git
The following commit(s) were added to refs/heads/main by this push:
new f110d44a feat(sedona-gdal): add raster operations and vrt support
(#698)
f110d44a is described below
commit f110d44abba11486ed63a91226db132cd971f7fe
Author: Kristin Cowalcijk <[email protected]>
AuthorDate: Wed Apr 1 11:38:49 2026 +0800
feat(sedona-gdal): add raster operations and vrt support (#698)
## Summary
- add VRT dataset support plus rasterize, rasterize-affine, and polygonize
wrappers
- layer higher-level raster operations on top of the dataset/raster/vector
wrapper stack
- remove raster wrapper re-export aliases so call sites import concrete
module paths directly
---
c/sedona-gdal/src/lib.rs | 1 +
c/sedona-gdal/src/raster.rs | 3 +
c/sedona-gdal/src/raster/polygonize.rs | 275 ++++++++++++++++++++
c/sedona-gdal/src/raster/rasterize.rs | 260 +++++++++++++++++++
c/sedona-gdal/src/raster/rasterize_affine.rs | 374 +++++++++++++++++++++++++++
c/sedona-gdal/src/vrt.rs | 321 +++++++++++++++++++++++
6 files changed, 1234 insertions(+)
diff --git a/c/sedona-gdal/src/lib.rs b/c/sedona-gdal/src/lib.rs
index c331f5df..f6138451 100644
--- a/c/sedona-gdal/src/lib.rs
+++ b/c/sedona-gdal/src/lib.rs
@@ -35,4 +35,5 @@ pub mod geo_transform;
pub mod raster;
pub mod spatial_ref;
pub mod vector;
+pub mod vrt;
pub mod vsi;
diff --git a/c/sedona-gdal/src/raster.rs b/c/sedona-gdal/src/raster.rs
index 389d9d73..83967429 100644
--- a/c/sedona-gdal/src/raster.rs
+++ b/c/sedona-gdal/src/raster.rs
@@ -15,5 +15,8 @@
// specific language governing permissions and limitations
// under the License.
+pub mod polygonize;
pub mod rasterband;
+pub mod rasterize;
+pub mod rasterize_affine;
pub mod types;
diff --git a/c/sedona-gdal/src/raster/polygonize.rs
b/c/sedona-gdal/src/raster/polygonize.rs
new file mode 100644
index 00000000..9291fd71
--- /dev/null
+++ b/c/sedona-gdal/src/raster/polygonize.rs
@@ -0,0 +1,275 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements. See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership. The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License. You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied. See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+use std::ptr;
+
+use crate::cpl::CslStringList;
+use crate::errors::Result;
+use crate::gdal_api::{call_gdal_api, GdalApi};
+use crate::gdal_dyn_bindgen::*;
+use crate::raster::rasterband::RasterBand;
+use crate::vector::layer::Layer;
+
+#[derive(Clone, Debug, Default)]
+pub struct PolygonizeOptions {
+ /// Use 8 connectedness (diagonal pixels are considered connected).
+ ///
+ /// If `false` (default), 4 connectedness is used.
+ pub eight_connected: bool,
+
+ /// Name of a dataset from which to read the geotransform.
+ ///
+ /// This is useful if the source band has no related dataset, which is
typical for mask bands.
+ ///
+ /// Corresponds to GDAL's `DATASET_FOR_GEOREF=dataset_name` option.
+ pub dataset_for_georef: Option<String>,
+
+ /// Interval in number of features at which transactions must be flushed.
+ ///
+ /// - `0` means that no transactions are opened.
+ /// - a negative value means a single transaction.
+ ///
+ /// Corresponds to GDAL's `COMMIT_INTERVAL=num` option.
+ pub commit_interval: Option<i64>,
+}
+
+impl PolygonizeOptions {
+ /// Build a GDAL option list from these polygonize options.
+ pub fn to_options_list(&self) -> Result<CslStringList> {
+ let mut options = CslStringList::new();
+
+ if self.eight_connected {
+ options.set_name_value("8CONNECTED", "8")?;
+ }
+
+ if let Some(ref ds) = self.dataset_for_georef {
+ options.set_name_value("DATASET_FOR_GEOREF", ds)?;
+ }
+
+ if let Some(interval) = self.commit_interval {
+ options.set_name_value("COMMIT_INTERVAL", &interval.to_string())?;
+ }
+
+ Ok(options)
+ }
+}
+
+/// Polygonize a raster band into a vector layer.
+/// This uses `GDALPolygonize`, which reads source pixels as integers.
+pub fn polygonize(
+ api: &'static GdalApi,
+ src_band: &RasterBand<'_>,
+ mask_band: Option<&RasterBand<'_>>,
+ out_layer: &Layer<'_>,
+ pixel_value_field: i32,
+ options: &PolygonizeOptions,
+) -> Result<()> {
+ let mask = mask_band.map_or(ptr::null_mut(), |b| b.c_rasterband());
+ let csl = options.to_options_list()?;
+
+ let rv = unsafe {
+ call_gdal_api!(
+ api,
+ GDALPolygonize,
+ src_band.c_rasterband(),
+ mask,
+ out_layer.c_layer(),
+ pixel_value_field,
+ csl.as_ptr(),
+ ptr::null_mut(), // pfnProgress
+ ptr::null_mut() // pProgressData
+ )
+ };
+ if rv != CE_None {
+ return Err(api.last_cpl_err(rv as u32));
+ }
+ Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ #[test]
+ fn test_polygonizeoptions_as_ptr() {
+ let c_options =
PolygonizeOptions::default().to_options_list().unwrap();
+ assert_eq!(c_options.fetch_name_value("8CONNECTED"), None);
+ assert_eq!(c_options.fetch_name_value("DATASET_FOR_GEOREF"), None);
+ assert_eq!(c_options.fetch_name_value("COMMIT_INTERVAL"), None);
+
+ let c_options = PolygonizeOptions {
+ eight_connected: true,
+ dataset_for_georef: Some("/vsimem/georef.tif".to_string()),
+ commit_interval: Some(12345),
+ }
+ .to_options_list()
+ .unwrap();
+ assert_eq!(c_options.fetch_name_value("8CONNECTED"), Some("8".into()));
+ assert_eq!(
+ c_options.fetch_name_value("DATASET_FOR_GEOREF"),
+ Some("/vsimem/georef.tif".into())
+ );
+ assert_eq!(
+ c_options.fetch_name_value("COMMIT_INTERVAL"),
+ Some("12345".into())
+ );
+ }
+
+ #[cfg(feature = "gdal-sys")]
+ #[test]
+ fn test_polygonize_connectivity_affects_regions() {
+ use crate::dataset::LayerOptions;
+ use crate::driver::DriverManager;
+ use crate::global::with_global_gdal_api;
+ use crate::raster::types::Buffer;
+ use crate::vector::feature::FieldDefn;
+ use crate::vsi::unlink_mem_file;
+
+ with_global_gdal_api(|api| {
+ let mem_driver = DriverManager::get_driver_by_name(api,
"MEM").unwrap();
+ let raster_ds = mem_driver.create("", 3, 3, 1).unwrap();
+ let band = raster_ds.rasterband(1).unwrap();
+
+ // 3x3 raster with diagonal 1s:
+ // 1 0 0
+ // 0 1 0
+ // 0 0 1
+ let mut data = Buffer::new((3, 3), vec![1u8, 0, 0, 0, 1, 0, 0, 0,
1]);
+ band.write((0, 0), (3, 3), &mut data).unwrap();
+
+ let gpkg_path = "/vsimem/test_polygonize_connectivity.gpkg";
+ let gpkg_driver = DriverManager::get_driver_by_name(api,
"GPKG").unwrap();
+ let vector_ds = gpkg_driver.create_vector_only(gpkg_path).unwrap();
+
+ // 4-connected output
+ let mut layer_4 = vector_ds
+ .create_layer(LayerOptions {
+ name: "four",
+ srs: None,
+ ty: OGRwkbGeometryType::wkbPolygon,
+ options: None,
+ })
+ .unwrap();
+ let field_defn = FieldDefn::new(api, "val",
OGRFieldType::OFTInteger).unwrap();
+ layer_4.create_field(&field_defn).unwrap();
+
+ polygonize(api, &band, None, &layer_4, 0,
&PolygonizeOptions::default()).unwrap();
+
+ let ones_4 = layer_4
+ .features()
+ .filter_map(|f| f.field_as_integer(0))
+ .filter(|v| *v == 1)
+ .count();
+ assert_eq!(ones_4, 3);
+
+ // 8-connected output
+ let mut layer_8 = vector_ds
+ .create_layer(LayerOptions {
+ name: "eight",
+ srs: None,
+ ty: OGRwkbGeometryType::wkbPolygon,
+ options: None,
+ })
+ .unwrap();
+ let field_defn = FieldDefn::new(api, "val",
OGRFieldType::OFTInteger).unwrap();
+ layer_8.create_field(&field_defn).unwrap();
+
+ polygonize(
+ api,
+ &band,
+ None,
+ &layer_8,
+ 0,
+ &PolygonizeOptions {
+ eight_connected: true,
+ dataset_for_georef: None,
+ commit_interval: None,
+ },
+ )
+ .unwrap();
+
+ let ones_8 = layer_8
+ .features()
+ .filter_map(|f| f.field_as_integer(0))
+ .filter(|v| *v == 1)
+ .count();
+ assert_eq!(ones_8, 1);
+
+ unlink_mem_file(api, gpkg_path).unwrap();
+ })
+ .unwrap();
+ }
+
+ #[cfg(feature = "gdal-sys")]
+ #[test]
+ fn test_polygonize_with_mask_band_restricts_output() {
+ use crate::dataset::LayerOptions;
+ use crate::driver::DriverManager;
+ use crate::global::with_global_gdal_api;
+ use crate::raster::types::Buffer;
+ use crate::vector::feature::FieldDefn;
+ use crate::vsi::unlink_mem_file;
+
+ with_global_gdal_api(|api| {
+ let mem_driver = DriverManager::get_driver_by_name(api,
"MEM").unwrap();
+ let raster_ds = mem_driver.create("", 3, 3, 2).unwrap();
+
+ let value_band = raster_ds.rasterband(1).unwrap();
+ let mask_band = raster_ds.rasterband(2).unwrap();
+
+ // Value band: all 7s.
+ let mut values = Buffer::new((3, 3), vec![7u8; 9]);
+ value_band.write((0, 0), (3, 3), &mut values).unwrap();
+
+ // Mask: only the center pixel is included.
+ let mut mask = Buffer::new((3, 3), vec![0u8, 0, 0, 0, 1, 0, 0, 0,
0]);
+ mask_band.write((0, 0), (3, 3), &mut mask).unwrap();
+
+ let gpkg_path = "/vsimem/test_polygonize_mask.gpkg";
+ let gpkg_driver = DriverManager::get_driver_by_name(api,
"GPKG").unwrap();
+ let vector_ds = gpkg_driver.create_vector_only(gpkg_path).unwrap();
+
+ let mut layer = vector_ds
+ .create_layer(LayerOptions {
+ name: "masked",
+ srs: None,
+ ty: OGRwkbGeometryType::wkbPolygon,
+ options: None,
+ })
+ .unwrap();
+ let field_defn = FieldDefn::new(api, "val",
OGRFieldType::OFTInteger).unwrap();
+ layer.create_field(&field_defn).unwrap();
+
+ polygonize(
+ api,
+ &value_band,
+ Some(&mask_band),
+ &layer,
+ 0,
+ &PolygonizeOptions::default(),
+ )
+ .unwrap();
+
+ assert_eq!(layer.feature_count(true), 1);
+ let only_val =
layer.features().next().unwrap().field_as_integer(0);
+ assert_eq!(only_val, Some(7));
+
+ unlink_mem_file(api, gpkg_path).unwrap();
+ })
+ .unwrap();
+ }
+}
diff --git a/c/sedona-gdal/src/raster/rasterize.rs
b/c/sedona-gdal/src/raster/rasterize.rs
new file mode 100644
index 00000000..a7a89a3c
--- /dev/null
+++ b/c/sedona-gdal/src/raster/rasterize.rs
@@ -0,0 +1,260 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements. See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership. The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License. You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied. See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! Ported (and contains copied code) from georust/gdal:
+//! <https://github.com/georust/gdal/blob/v0.19.0/src/raster/rasterize.rs>.
+//! Original code is licensed under MIT.
+
+use std::ptr;
+
+use crate::cpl::CslStringList;
+use crate::dataset::Dataset;
+use crate::errors::{GdalError, Result};
+use crate::gdal_api::{call_gdal_api, GdalApi};
+use crate::gdal_dyn_bindgen::*;
+use crate::vector::geometry::Geometry;
+
+/// Source of burn values.
+#[derive(Copy, Clone, Debug)]
+pub enum BurnSource {
+ /// Use whatever `burn_values` argument is supplied to `rasterize`.
+ UserSupplied,
+
+ /// Add the geometry's Z value to whatever `burn_values` argument
+ /// is supplied to `rasterize`.
+ Z,
+}
+
+/// Algorithm for merging new raster values with existing values.
+#[derive(Copy, Clone, Debug)]
+pub enum MergeAlgorithm {
+ /// Overwrite existing value (default).
+ Replace,
+ /// Add new value to existing value (useful for heatmaps).
+ Add,
+}
+
+/// Optimization mode for rasterization.
+#[derive(Copy, Clone, Debug)]
+pub enum OptimizeMode {
+ /// Let GDAL decide (default).
+ Automatic,
+ /// Force raster-based scan (iterates over pixels).
+ Raster,
+ /// Force vector-based scan (iterates over geometry edges).
+ Vector,
+}
+
+/// Options that specify how to rasterize geometries.
+#[derive(Copy, Clone, Debug)]
+pub struct RasterizeOptions {
+ /// Set to `true` to set all pixels touched by the line or polygons,
+ /// not just those whose center is within the polygon or that are
+ /// selected by Bresenham's line algorithm. Defaults to `false`.
+ pub all_touched: bool,
+
+ /// May be set to `BurnSource::Z` to use the Z values of the geometries.
+ /// `burn_value` is added to this before burning. Defaults to
+ /// `BurnSource::UserSupplied` in which case just the `burn_value` is
burned.
+ pub source: BurnSource,
+
+ /// May be `MergeAlgorithm::Replace` (default) or `MergeAlgorithm::Add`.
+ /// `Replace` overwrites existing values; `Add` adds to them.
+ pub merge_algorithm: MergeAlgorithm,
+
+ /// The height in lines of the chunk to operate on. `0` (default) lets GDAL
+ /// choose based on cache size. Not used in `OPTIM=RASTER` mode.
+ pub chunk_y_size: usize,
+
+ /// Optimization mode for rasterization.
+ pub optimize: OptimizeMode,
+}
+
+impl Default for RasterizeOptions {
+ fn default() -> Self {
+ RasterizeOptions {
+ all_touched: false,
+ source: BurnSource::UserSupplied,
+ merge_algorithm: MergeAlgorithm::Replace,
+ chunk_y_size: 0,
+ optimize: OptimizeMode::Automatic,
+ }
+ }
+}
+
+impl RasterizeOptions {
+ /// Build a GDAL option list from these rasterize options.
+ pub fn to_options_list(self) -> Result<CslStringList> {
+ let mut options = CslStringList::with_capacity(5);
+
+ options.set_name_value(
+ "ALL_TOUCHED",
+ if self.all_touched { "TRUE" } else { "FALSE" },
+ )?;
+
+ options.set_name_value(
+ "MERGE_ALG",
+ match self.merge_algorithm {
+ MergeAlgorithm::Replace => "REPLACE",
+ MergeAlgorithm::Add => "ADD",
+ },
+ )?;
+
+ options.set_name_value("CHUNKYSIZE", &self.chunk_y_size.to_string())?;
+
+ options.set_name_value(
+ "OPTIM",
+ match self.optimize {
+ OptimizeMode::Automatic => "AUTO",
+ OptimizeMode::Raster => "RASTER",
+ OptimizeMode::Vector => "VECTOR",
+ },
+ )?;
+
+ if let BurnSource::Z = self.source {
+ options.set_name_value("BURN_VALUE_FROM", "Z")?;
+ }
+
+ Ok(options)
+ }
+}
+
+/// Rasterize geometries into the selected dataset bands.
+/// Supply one burn value per geometry; values are replicated across the
target bands.
+pub fn rasterize(
+ api: &'static GdalApi,
+ dataset: &Dataset,
+ band_list: &[i32],
+ geometries: &[&Geometry],
+ burn_values: &[f64],
+ options: Option<RasterizeOptions>,
+) -> Result<()> {
+ if band_list.is_empty() {
+ return Err(GdalError::BadArgument(
+ "`band_list` must not be empty".to_string(),
+ ));
+ }
+ if burn_values.len() != geometries.len() {
+ return Err(GdalError::BadArgument(format!(
+ "burn_values length ({}) must match geometries length ({})",
+ burn_values.len(),
+ geometries.len()
+ )));
+ }
+ let raster_count = dataset.raster_count();
+ for &band in band_list {
+ let is_good = band > 0 && (band as usize) <= raster_count;
+ if !is_good {
+ return Err(GdalError::BadArgument(format!(
+ "Band index {} is out of bounds",
+ band
+ )));
+ }
+ }
+
+ let geom_handles: Vec<OGRGeometryH> = geometries.iter().map(|g|
g.c_geometry()).collect();
+
+ // Replicate each burn value across all bands, matching the GDAL C API
+ // contract that expects nGeomCount * nBandCount burn values.
+ let expanded_burn_values: Vec<f64> = burn_values
+ .iter()
+ .flat_map(|burn| std::iter::repeat_n(burn, band_list.len()))
+ .copied()
+ .collect();
+
+ let opts = options.unwrap_or_default();
+ let csl = opts.to_options_list()?;
+
+ let n_band_count: i32 = band_list.len().try_into()?;
+ let n_geom_count: i32 = geom_handles.len().try_into()?;
+
+ let rv = unsafe {
+ call_gdal_api!(
+ api,
+ GDALRasterizeGeometries,
+ dataset.c_dataset(),
+ n_band_count,
+ band_list.as_ptr(),
+ n_geom_count,
+ geom_handles.as_ptr(),
+ ptr::null_mut(), // pfnTransformer
+ ptr::null_mut(), // pTransformArg
+ expanded_burn_values.as_ptr(),
+ csl.as_ptr(),
+ ptr::null_mut(), // pfnProgress
+ ptr::null_mut() // pProgressData
+ )
+ };
+ if rv != CE_None {
+ return Err(api.last_cpl_err(rv as u32));
+ }
+ Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use crate::driver::DriverManager;
+ use crate::global::with_global_gdal_api;
+
+ #[test]
+ fn test_rasterizeoptions_as_ptr() {
+ let c_options = RasterizeOptions::default().to_options_list().unwrap();
+ assert_eq!(
+ c_options.fetch_name_value("ALL_TOUCHED"),
+ Some("FALSE".to_string())
+ );
+ assert_eq!(c_options.fetch_name_value("BURN_VALUE_FROM"), None);
+ assert_eq!(
+ c_options.fetch_name_value("MERGE_ALG"),
+ Some("REPLACE".to_string())
+ );
+ assert_eq!(
+ c_options.fetch_name_value("CHUNKYSIZE"),
+ Some("0".to_string())
+ );
+ assert_eq!(
+ c_options.fetch_name_value("OPTIM"),
+ Some("AUTO".to_string())
+ );
+ }
+
+ #[cfg(feature = "gdal-sys")]
+ #[test]
+ fn test_rasterize() {
+ with_global_gdal_api(|api| {
+ let wkt = "POLYGON ((2 2, 2 4.25, 4.25 4.25, 4.25 2, 2 2))";
+ let poly = Geometry::from_wkt(api, wkt).unwrap();
+
+ let driver = DriverManager::get_driver_by_name(api,
"MEM").unwrap();
+ let dataset = driver.create("", 5, 5, 1).unwrap();
+
+ let bands = [1];
+ let geometries = [&poly];
+ let burn_values = [1.0];
+ rasterize(api, &dataset, &bands, &geometries, &burn_values,
None).unwrap();
+
+ let rb = dataset.rasterband(1).unwrap();
+ let values = rb.read_as::<u8>((0, 0), (5, 5), (5, 5),
None).unwrap();
+ assert_eq!(
+ values.data(),
+ [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0,
0, 0, 0, 0, 0]
+ );
+ })
+ .unwrap();
+ }
+}
diff --git a/c/sedona-gdal/src/raster/rasterize_affine.rs
b/c/sedona-gdal/src/raster/rasterize_affine.rs
new file mode 100644
index 00000000..58b347b8
--- /dev/null
+++ b/c/sedona-gdal/src/raster/rasterize_affine.rs
@@ -0,0 +1,374 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements. See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership. The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License. You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied. See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! Fast affine-transformer rasterize wrapper.
+//!
+//! GDALRasterizeGeometries() will internally call
GDALCreateGenImgProjTransformer2()
+//! if pfnTransformer is NULL, even in the common case where only a
GeoTransform-based
+//! affine conversion from georeferenced coords to pixel/line is needed.
+//!
+//! This module supplies a minimal GDALTransformerFunc that applies the dataset
+//! GeoTransform (and its inverse), avoiding expensive transformer creation.
+
+use std::ffi::{c_char, c_int, c_void};
+use std::ptr;
+
+use crate::dataset::Dataset;
+use crate::errors::{GdalError, Result};
+use crate::gdal_api::{call_gdal_api, GdalApi};
+use crate::gdal_dyn_bindgen::CE_None;
+use crate::geo_transform::{GeoTransform, GeoTransformEx};
+use crate::vector::geometry::Geometry;
+
+#[repr(C)]
+struct AffineTransformArg {
+ gt: GeoTransform,
+ inv_gt: GeoTransform,
+}
+
+unsafe extern "C" fn affine_transformer(
+ p_transformer_arg: *mut c_void,
+ b_dst_to_src: c_int,
+ n_point_count: c_int,
+ x: *mut f64,
+ y: *mut f64,
+ _z: *mut f64,
+ pan_success: *mut c_int,
+) -> c_int {
+ if p_transformer_arg.is_null() || x.is_null() || y.is_null() ||
pan_success.is_null() {
+ return 0;
+ }
+ if n_point_count < 0 {
+ return 0;
+ }
+
+ // Treat transformer arg as immutable.
+ let arg = &*(p_transformer_arg as *const AffineTransformArg);
+ let affine = if b_dst_to_src == 0 {
+ &arg.inv_gt
+ } else {
+ &arg.gt
+ };
+
+ let n = n_point_count as usize;
+ for i in 0..n {
+ // SAFETY: x/y/pan_success are assumed to point to arrays of length
n_point_count.
+ let xin = unsafe { *x.add(i) };
+ let yin = unsafe { *y.add(i) };
+ let (xout, yout) = affine.apply(xin, yin);
+ unsafe {
+ *x.add(i) = xout;
+ *y.add(i) = yout;
+ *pan_success.add(i) = 1;
+ }
+ }
+
+ 1
+}
+
+/// Rasterize geometries using the dataset geotransform as the transformer.
+/// Geometry coordinates must already be in the dataset georeferenced
coordinate space.
+pub fn rasterize_affine(
+ api: &'static GdalApi,
+ dataset: &Dataset,
+ bands: &[usize],
+ geometries: &[Geometry],
+ burn_values: &[f64],
+ all_touched: bool,
+) -> Result<()> {
+ if bands.is_empty() {
+ return Err(GdalError::BadArgument(
+ "`bands` must not be empty".to_string(),
+ ));
+ }
+ if burn_values.len() != geometries.len() {
+ return Err(GdalError::BadArgument(format!(
+ "Burn values length ({}) must match geometries length ({})",
+ burn_values.len(),
+ geometries.len()
+ )));
+ }
+
+ let raster_count = dataset.raster_count();
+ for band in bands {
+ let is_good = *band > 0 && *band <= raster_count;
+ if !is_good {
+ return Err(GdalError::BadArgument(format!(
+ "Band index {} is out of bounds",
+ *band
+ )));
+ }
+ }
+
+ let bands_i32: Vec<c_int> = bands.iter().map(|&band| band as
c_int).collect();
+
+ // Keep this stack-allocated option array instead of going through
`CslStringList`
+ // so the affine fast path avoids extra allocation overhead for tiny
geometries.
+ let mut c_options: [*mut c_char; 2] = if all_touched {
+ [c"ALL_TOUCHED=TRUE".as_ptr() as *mut c_char, ptr::null_mut()]
+ } else {
+ [
+ c"ALL_TOUCHED=FALSE".as_ptr() as *mut c_char,
+ ptr::null_mut(),
+ ]
+ };
+
+ let geometries_c: Vec<_> = geometries.iter().map(|geo|
geo.c_geometry()).collect();
+ let burn_values_expanded: Vec<f64> = burn_values
+ .iter()
+ .flat_map(|burn| std::iter::repeat_n(burn, bands_i32.len()))
+ .copied()
+ .collect();
+
+ let gt = dataset.geo_transform().map_err(|_e| {
+ GdalError::BadArgument(
+ "Missing geotransform: only geotransform-based affine rasterize is
supported"
+ .to_string(),
+ )
+ })?;
+ let inv_gt = gt.invert().map_err(|_e| {
+ GdalError::BadArgument(
+ "Non-invertible geotransform: only geotransform-based affine
rasterize is supported"
+ .to_string(),
+ )
+ })?;
+ let mut arg = AffineTransformArg { gt, inv_gt };
+
+ unsafe {
+ let error = call_gdal_api!(
+ api,
+ GDALRasterizeGeometries,
+ dataset.c_dataset(),
+ bands_i32.len() as c_int,
+ bands_i32.as_ptr(),
+ geometries_c.len() as c_int,
+ geometries_c.as_ptr(),
+ // SAFETY: `affine_transformer` has the GDAL transformer callback
ABI.
+ // This binding models the raw C API slot as `void*`, and on
supported
+ // targets we rely on the platform ABI allowing this callback
pointer to
+ // pass through that raw field unchanged.
+ (affine_transformer as *const ()).cast::<c_void>() as *mut c_void,
+ (&mut arg as *mut AffineTransformArg).cast::<c_void>(),
+ burn_values_expanded.as_ptr(),
+ // SAFETY: GDAL reads option strings through this mutable `char**`
API.
+ // The pointed-to string literals remain valid for the call and
are not
+ // expected to be modified by `GDALRasterizeGeometries`.
+ c_options.as_mut_ptr(),
+ ptr::null_mut(),
+ ptr::null_mut()
+ );
+ if error != CE_None {
+ return Err(api.last_cpl_err(error as u32));
+ }
+ }
+ Ok(())
+}
+
+#[cfg(all(test, feature = "gdal-sys"))]
+mod tests {
+ use super::*;
+
+ use crate::driver::{Driver, DriverManager};
+ use crate::global::with_global_gdal_api;
+ use crate::raster::rasterize::{rasterize, RasterizeOptions};
+ use crate::raster::types::Buffer;
+
+ fn mem_driver(api: &'static GdalApi) -> Driver {
+ DriverManager::get_driver_by_name(api, "MEM").unwrap()
+ }
+
+ fn make_dataset_u8(
+ api: &'static GdalApi,
+ width: usize,
+ height: usize,
+ gt: GeoTransform,
+ ) -> Result<Dataset> {
+ let driver = mem_driver(api);
+ let ds = driver.create_with_band_type::<u8>("", width, height, 1)?;
+ ds.set_geo_transform(>)?;
+ let band = ds.rasterband(1)?;
+ let mut buf = Buffer::new((width, height), vec![0u8; width * height]);
+ band.write((0, 0), (width, height), &mut buf)?;
+ Ok(ds)
+ }
+
+ fn read_u8(ds: &Dataset, width: usize, height: usize) -> Vec<u8> {
+ let band = ds.rasterband(1).unwrap();
+ let buf = band
+ .read_as::<u8>((0, 0), (width, height), (width, height), None)
+ .unwrap();
+ buf.data().to_vec()
+ }
+
+ fn poly_from_pixel_rect(
+ api: &'static GdalApi,
+ gt: &GeoTransform,
+ x0: f64,
+ y0: f64,
+ x1: f64,
+ y1: f64,
+ ) -> Geometry {
+ let (wx0, wy0) = gt.apply(x0, y0);
+ let (wx1, wy1) = gt.apply(x1, y0);
+ let (wx2, wy2) = gt.apply(x1, y1);
+ let (wx3, wy3) = gt.apply(x0, y1);
+ let wkt =
+ format!("POLYGON (({wx0} {wy0}, {wx1} {wy1}, {wx2} {wy2}, {wx3}
{wy3}, {wx0} {wy0}))");
+ Geometry::from_wkt(api, &wkt).unwrap()
+ }
+
+ fn line_from_pixel_points(
+ api: &'static GdalApi,
+ gt: &GeoTransform,
+ pts: &[(f64, f64)],
+ ) -> Geometry {
+ assert!(pts.len() >= 2);
+ let mut s = String::from("LINESTRING (");
+ for (i, (px, py)) in pts.iter().copied().enumerate() {
+ let (wx, wy) = gt.apply(px, py);
+ if i > 0 {
+ s.push_str(", ");
+ }
+ s.push_str(&format!("{wx} {wy}"));
+ }
+ s.push(')');
+ Geometry::from_wkt(api, &s).unwrap()
+ }
+
+ #[test]
+ fn test_rasterize_affine_matches_baseline_north_up() {
+ with_global_gdal_api(|api| {
+ let (w, h) = (32usize, 24usize);
+ let gt: GeoTransform = [100.0, 2.0, 0.0, 200.0, 0.0, -2.0];
+
+ let geom_baseline = poly_from_pixel_rect(api, >, 3.2, 4.7, 20.4,
18.1);
+ let geom_affine = poly_from_pixel_rect(api, >, 3.2, 4.7, 20.4,
18.1);
+
+ let ds_baseline = make_dataset_u8(api, w, h, gt).unwrap();
+ let ds_affine = make_dataset_u8(api, w, h, gt).unwrap();
+
+ let geom_refs: Vec<&Geometry> = vec![&geom_baseline];
+ rasterize(api, &ds_baseline, &[1], &geom_refs, &[1.0],
None).unwrap();
+ rasterize_affine(api, &ds_affine, &[1], &[geom_affine], &[1.0],
false).unwrap();
+
+ assert_eq!(read_u8(&ds_affine, w, h), read_u8(&ds_baseline, w, h));
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_rasterize_affine_matches_baseline_rotated_gt_all_touched() {
+ with_global_gdal_api(|api| {
+ let (w, h) = (40usize, 28usize);
+ // Rotated/skewed GeoTransform.
+ let gt: GeoTransform = [10.0, 1.2, 0.15, 50.0, -0.1, -1.1];
+
+ let geom_baseline = poly_from_pixel_rect(api, >, 5.25, 4.5,
25.75, 20.25);
+ let geom_affine = poly_from_pixel_rect(api, >, 5.25, 4.5, 25.75,
20.25);
+
+ let ds_baseline = make_dataset_u8(api, w, h, gt).unwrap();
+ let ds_affine = make_dataset_u8(api, w, h, gt).unwrap();
+
+ let geom_refs: Vec<&Geometry> = vec![&geom_baseline];
+ rasterize(
+ api,
+ &ds_baseline,
+ &[1],
+ &geom_refs,
+ &[1.0],
+ Some(RasterizeOptions {
+ all_touched: true,
+ ..Default::default()
+ }),
+ )
+ .unwrap();
+ rasterize_affine(api, &ds_affine, &[1], &[geom_affine], &[1.0],
true).unwrap();
+
+ assert_eq!(read_u8(&ds_affine, w, h), read_u8(&ds_baseline, w, h));
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_rasterize_affine_matches_baseline_linestring() {
+ with_global_gdal_api(|api| {
+ let (w, h) = (64usize, 48usize);
+ // Rotated/skewed GeoTransform.
+ let gt: GeoTransform = [5.0, 1.0, 0.2, 100.0, -0.15, -1.05];
+
+ // A polyline with many vertices, defined in pixel/line space.
+ let mut pts: Vec<(f64, f64)> = Vec::new();
+ for i in 0..200 {
+ let t = i as f64 / 199.0;
+ let x = 2.625 + t * ((w as f64) - 5.25);
+ let y = 5.25 + (t * 6.0).sin() * 8.0 + t * ((h as f64) -
12.25);
+ pts.push((x, y));
+ }
+ let geom_baseline = line_from_pixel_points(api, >, &pts);
+ let geom_affine = line_from_pixel_points(api, >, &pts);
+
+ let ds_baseline = make_dataset_u8(api, w, h, gt).unwrap();
+ let ds_affine = make_dataset_u8(api, w, h, gt).unwrap();
+
+ let geom_refs: Vec<&Geometry> = vec![&geom_baseline];
+ rasterize(api, &ds_baseline, &[1], &geom_refs, &[1.0],
None).unwrap();
+ rasterize_affine(api, &ds_affine, &[1], &[geom_affine], &[1.0],
false).unwrap();
+
+ let got = read_u8(&ds_affine, w, h);
+ let expected = read_u8(&ds_baseline, w, h);
+ if got != expected {
+ let mut diffs = Vec::new();
+ for (i, (a, b)) in got
+ .iter()
+ .copied()
+ .zip(expected.iter().copied())
+ .enumerate()
+ {
+ if a != b {
+ let x = i % w;
+ let y = i / w;
+ diffs.push((x, y, a, b));
+ }
+ }
+ panic!(
+ "raster mismatch: {} differing pixels; first 10: {:?}",
+ diffs.len(),
+ &diffs[..diffs.len().min(10)]
+ );
+ }
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_rasterize_affine_fails_on_noninvertible_gt() {
+ with_global_gdal_api(|api| {
+ let (w, h) = (8usize, 8usize);
+ let gt: GeoTransform = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
+ let ds = make_dataset_u8(api, w, h, gt).unwrap();
+ let geom = Geometry::from_wkt(api, "POINT (0 0)").unwrap();
+ let err = rasterize_affine(api, &ds, &[1], &[geom], &[1.0],
true).unwrap_err();
+ match err {
+ GdalError::BadArgument(msg) => {
+ assert!(msg.contains("Non-invertible geotransform"));
+ }
+ other => panic!("Unexpected error: {other:?}"),
+ }
+ })
+ .unwrap();
+ }
+}
diff --git a/c/sedona-gdal/src/vrt.rs b/c/sedona-gdal/src/vrt.rs
new file mode 100644
index 00000000..4287915f
--- /dev/null
+++ b/c/sedona-gdal/src/vrt.rs
@@ -0,0 +1,321 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements. See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership. The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License. You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied. See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! GDAL VRT (Virtual Raster) API wrappers.
+
+use std::ffi::CString;
+use std::ops::{Deref, DerefMut};
+use std::ptr::{null, null_mut};
+
+use crate::cpl::CslStringList;
+use crate::dataset::Dataset;
+use crate::errors::Result;
+use crate::gdal_api::{call_gdal_api, GdalApi};
+use crate::raster::rasterband::RasterBand;
+use crate::{gdal_dyn_bindgen::*, raster::types::GdalDataType};
+
+/// Special value indicating that nodata is not set for a VRT source.
+/// Matches `VRT_NODATA_UNSET` from GDAL's `gdal_vrt.h`.
+pub const NODATA_UNSET: f64 = -1234.56;
+
+/// A VRT (Virtual Raster) dataset.
+pub struct VrtDataset {
+ dataset: Dataset,
+}
+
+impl VrtDataset {
+ /// Create an empty VRT dataset with the given raster size.
+ pub fn create(api: &'static GdalApi, x_size: usize, y_size: usize) ->
Result<Self> {
+ let x: i32 = x_size.try_into()?;
+ let y: i32 = y_size.try_into()?;
+ let c_dataset = unsafe { call_gdal_api!(api, VRTCreate, x, y) };
+
+ if c_dataset.is_null() {
+ return Err(api.last_null_pointer_err("VRTCreate"));
+ }
+
+ Ok(VrtDataset {
+ dataset: Dataset::new(api, c_dataset),
+ })
+ }
+
+ /// Return the underlying `Dataset`, transferring ownership.
+ pub fn as_dataset(self) -> Dataset {
+ let VrtDataset { dataset } = self;
+ dataset
+ }
+
+ /// Add a band to this VRT dataset.
+ /// Return the 1-based index of the new band.
+ pub fn add_band(&mut self, data_type: GdalDataType, options:
Option<&[&str]>) -> Result<usize> {
+ let csl =
CslStringList::try_from_iter(options.unwrap_or(&[]).iter().copied())?;
+
+ // Preserve null semantics: pass null when no options given.
+ let opts_ptr = if csl.is_empty() {
+ null_mut()
+ } else {
+ csl.as_ptr()
+ };
+
+ let rv = unsafe {
+ call_gdal_api!(
+ self.dataset.api(),
+ GDALAddBand,
+ self.dataset.c_dataset(),
+ data_type.to_c(),
+ opts_ptr
+ )
+ };
+
+ if rv != CE_None {
+ return Err(self.dataset.api().last_cpl_err(rv as u32));
+ }
+
+ Ok(self.raster_count())
+ }
+
+ /// Fetch a VRT band by 1-indexed band number.
+ pub fn rasterband(&self, band_index: usize) -> Result<VrtRasterBand<'_>> {
+ let band = self.dataset.rasterband(band_index)?;
+ Ok(VrtRasterBand { band })
+ }
+}
+
+impl Deref for VrtDataset {
+ type Target = Dataset;
+
+ fn deref(&self) -> &Self::Target {
+ &self.dataset
+ }
+}
+
+impl DerefMut for VrtDataset {
+ fn deref_mut(&mut self) -> &mut Self::Target {
+ &mut self.dataset
+ }
+}
+
+impl AsRef<Dataset> for VrtDataset {
+ fn as_ref(&self) -> &Dataset {
+ &self.dataset
+ }
+}
+
+/// A raster band within a VRT dataset.
+pub struct VrtRasterBand<'a> {
+ band: RasterBand<'a>,
+}
+
+impl<'a> VrtRasterBand<'a> {
+ /// Return the raw GDAL raster band handle.
+ pub fn c_rasterband(&self) -> GDALRasterBandH {
+ self.band.c_rasterband()
+ }
+
+ /// Add a simple source to this VRT band.
+ /// Map a source window to a destination window, with optional resampling
and nodata.
+ pub fn add_simple_source(
+ &self,
+ source_band: &RasterBand<'_>,
+ src_window: (i32, i32, i32, i32),
+ dst_window: (i32, i32, i32, i32),
+ resampling: Option<&str>,
+ nodata: Option<f64>,
+ ) -> Result<()> {
+ let c_resampling = resampling.map(CString::new).transpose()?;
+
+ let resampling_ptr = c_resampling.as_ref().map(|s|
s.as_ptr()).unwrap_or(null());
+
+ let nodata_value = nodata.unwrap_or(NODATA_UNSET);
+
+ let rv = unsafe {
+ call_gdal_api!(
+ self.band.api(),
+ VRTAddSimpleSource,
+ self.band.c_rasterband(),
+ source_band.c_rasterband(),
+ src_window.0,
+ src_window.1,
+ src_window.2,
+ src_window.3,
+ dst_window.0,
+ dst_window.1,
+ dst_window.2,
+ dst_window.3,
+ resampling_ptr,
+ nodata_value
+ )
+ };
+
+ if rv != CE_None {
+ return Err(self.band.api().last_cpl_err(rv as u32));
+ }
+ Ok(())
+ }
+
+ /// Set the nodata value for this VRT band.
+ pub fn set_no_data_value(&self, nodata: f64) -> Result<()> {
+ let rv = unsafe {
+ call_gdal_api!(
+ self.band.api(),
+ GDALSetRasterNoDataValue,
+ self.band.c_rasterband(),
+ nodata
+ )
+ };
+
+ if rv != CE_None {
+ return Err(self.band.api().last_cpl_err(rv as u32));
+ }
+ Ok(())
+ }
+}
+
+impl<'a> Deref for VrtRasterBand<'a> {
+ type Target = RasterBand<'a>;
+
+ fn deref(&self) -> &Self::Target {
+ &self.band
+ }
+}
+
+impl<'a> DerefMut for VrtRasterBand<'a> {
+ fn deref_mut(&mut self) -> &mut Self::Target {
+ &mut self.band
+ }
+}
+
+impl<'a> AsRef<RasterBand<'a>> for VrtRasterBand<'a> {
+ fn as_ref(&self) -> &RasterBand<'a> {
+ &self.band
+ }
+}
+
+#[cfg(all(test, feature = "gdal-sys"))]
+mod tests {
+ use crate::dataset::Dataset;
+ use crate::gdal_dyn_bindgen::GDAL_OF_READONLY;
+ use crate::global::with_global_gdal_api;
+ use crate::raster::types::GdalDataType;
+ use crate::vrt::{VrtDataset, NODATA_UNSET};
+
+ fn fixture(name: &str) -> String {
+ sedona_testing::data::test_raster(name).unwrap()
+ }
+
+ #[test]
+ fn test_vrt_create() {
+ with_global_gdal_api(|api| {
+ let vrt = VrtDataset::create(api, 100, 100).unwrap();
+ assert_eq!(vrt.raster_count(), 0);
+ assert!(!vrt.c_dataset().is_null());
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_vrt_add_band() {
+ with_global_gdal_api(|api| {
+ let mut vrt = VrtDataset::create(api, 100, 100).unwrap();
+ let band_idx = vrt.add_band(GdalDataType::Float32, None).unwrap();
+ assert_eq!(band_idx, 1);
+ assert_eq!(vrt.raster_count(), 1);
+
+ let band_idx = vrt.add_band(GdalDataType::UInt8, None).unwrap();
+ assert_eq!(band_idx, 2);
+ assert_eq!(vrt.raster_count(), 2);
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_vrt_set_geo_transform() {
+ with_global_gdal_api(|api| {
+ let vrt = VrtDataset::create(api, 100, 100).unwrap();
+ let transform = [0.0, 1.0, 0.0, 100.0, 0.0, -1.0];
+ vrt.set_geo_transform(&transform).unwrap();
+ assert_eq!(vrt.geo_transform().unwrap(), transform);
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_vrt_set_projection() {
+ with_global_gdal_api(|api| {
+ let vrt = VrtDataset::create(api, 100, 100).unwrap();
+ vrt.set_projection("EPSG:4326").unwrap();
+ assert!(vrt.projection().contains("4326"));
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_vrt_add_simple_source() {
+ with_global_gdal_api(|api| {
+ let source = Dataset::open_ex(
+ api,
+ &fixture("tinymarble.tif"),
+ GDAL_OF_READONLY,
+ None,
+ None,
+ None,
+ )
+ .unwrap();
+ let source_band_type = source.rasterband(1).unwrap().band_type();
+
+ let mut vrt = VrtDataset::create(api, 1, 1).unwrap();
+ vrt.add_band(source_band_type, None).unwrap();
+
+ let source_band = source.rasterband(1).unwrap();
+ let vrt_band = vrt.rasterband(1).unwrap();
+
+ vrt_band
+ .add_simple_source(&source_band, (0, 0, 1, 1), (0, 0, 1, 1),
None, None)
+ .unwrap();
+
+ let source_px = source_band
+ .read_as::<f64>((0, 0), (1, 1), (1, 1), None)
+ .unwrap()
+ .data()[0];
+ let vrt_px = vrt_band
+ .read_as::<f64>((0, 0), (1, 1), (1, 1), None)
+ .unwrap()
+ .data()[0];
+
+ assert_eq!(vrt_px, source_px);
+ })
+ .unwrap();
+ }
+
+ #[test]
+ fn test_vrt_nodata_unset() {
+ assert_eq!(NODATA_UNSET, -1234.56);
+ }
+
+ #[test]
+ #[allow(clippy::float_cmp)]
+ fn test_vrt_set_no_data_value() {
+ with_global_gdal_api(|api| {
+ let mut vrt = VrtDataset::create(api, 1, 1).unwrap();
+ vrt.add_band(GdalDataType::UInt8, None).unwrap();
+ let band = vrt.rasterband(1).unwrap();
+ band.set_no_data_value(-9999.0).unwrap();
+ assert_eq!(band.no_data_value(), Some(-9999.0));
+ })
+ .unwrap();
+ }
+}
