diff --git a/c/sedona-gdal/src/dataset.rs b/c/sedona-gdal/src/dataset.rs new file mode 100644 index 000000000..a7119a6b4 --- /dev/null +++ b/c/sedona-gdal/src/dataset.rs @@ -0,0 +1,458 @@ +// 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: +//! . +//! Original code is licensed under MIT. + +use std::ffi::{CStr, CString}; +use std::ptr; + +use crate::cpl::CslStringList; +use crate::driver::Driver; +use crate::errors::Result; +use crate::gdal_api::{call_gdal_api, GdalApi}; +use crate::gdal_dyn_bindgen::*; +use crate::raster::rasterband::RasterBand; +use crate::raster::types::{DatasetOptions, GdalDataType as RustGdalDataType}; +use crate::spatial_ref::SpatialRef; +use crate::vector::layer::Layer; + +/// A GDAL dataset. +pub struct Dataset { + api: &'static GdalApi, + c_dataset: GDALDatasetH, + owned: bool, +} + +// SAFETY: `Dataset` carries an opaque GDAL dataset handle plus an ownership flag. +// Moving the wrapper across threads only transfers ownership of that handle; it does +// not permit concurrent shared access. The handle is closed at most once on drop when +// `owned` is true, so `Send` is sound while `Sync` remains intentionally unimplemented. +unsafe impl Send for Dataset {} + +impl Drop for Dataset { + fn drop(&mut self) { + if self.owned && !self.c_dataset.is_null() { + unsafe { call_gdal_api!(self.api, GDALClose, self.c_dataset) }; + } + } +} + +impl Dataset { + /// Open a dataset with extended options. + pub fn open_ex( + api: &'static GdalApi, + path: &str, + open_flags: GDALOpenFlags, + allowed_drivers: Option<&[&str]>, + open_options: Option<&[&str]>, + sibling_files: Option<&[&str]>, + ) -> Result { + let c_path = CString::new(path)?; + + // Build CslStringLists from Option<&[&str]>. + // None → null pointer (use GDAL default). + // Some(&[]) → pointer to [null] (explicitly empty list). + let drivers_csl = allowed_drivers + .map(|v| CslStringList::try_from_iter(v.iter().copied())) + .transpose()?; + let options_csl = open_options + .map(|v| CslStringList::try_from_iter(v.iter().copied())) + .transpose()?; + let siblings_csl = sibling_files + .map(|v| CslStringList::try_from_iter(v.iter().copied())) + .transpose()?; + + let c_dataset = unsafe { + call_gdal_api!( + api, + GDALOpenEx, + c_path.as_ptr(), + open_flags, + drivers_csl + .as_ref() + .map_or(ptr::null(), |csl| csl.as_ptr() as *const *const _), + options_csl + .as_ref() + .map_or(ptr::null(), |csl| csl.as_ptr() as *const *const _), + siblings_csl + .as_ref() + .map_or(ptr::null(), |csl| csl.as_ptr() as *const *const _) + ) + }; + + if c_dataset.is_null() { + return Err(api.last_cpl_err(CE_Failure as u32)); + } + + Ok(Self { + api, + c_dataset, + owned: true, + }) + } + + /// Create a new owned Dataset from a C handle. + pub(crate) fn new_owned(api: &'static GdalApi, c_dataset: GDALDatasetH) -> Self { + Self { + api, + c_dataset, + owned: true, + } + } + + /// Wrap an existing C dataset handle (non-owning). + /// + /// # Safety + /// + /// The caller must ensure the handle is valid and outlives this `Dataset`. + pub unsafe fn from_c_dataset(api: &'static GdalApi, c_dataset: GDALDatasetH) -> Self { + Self { + api, + c_dataset, + owned: false, + } + } + + /// Return the raw C dataset handle. + pub fn c_dataset(&self) -> GDALDatasetH { + self.c_dataset + } + + /// Return raster size as (x_size, y_size). + pub fn raster_size(&self) -> (usize, usize) { + let x = unsafe { call_gdal_api!(self.api, GDALGetRasterXSize, self.c_dataset) }; + let y = unsafe { call_gdal_api!(self.api, GDALGetRasterYSize, self.c_dataset) }; + (x as usize, y as usize) + } + + /// Return the number of raster bands. + pub fn raster_count(&self) -> usize { + unsafe { call_gdal_api!(self.api, GDALGetRasterCount, self.c_dataset) as usize } + } + + /// Fetch a raster band by 1-indexed band number. + /// Band numbers start at 1, as in GDAL. + pub fn rasterband(&self, band_index: usize) -> Result> { + let band_index_i32 = i32::try_from(band_index)?; + let c_band = + unsafe { call_gdal_api!(self.api, GDALGetRasterBand, self.c_dataset, band_index_i32) }; + if c_band.is_null() { + return Err(self.api.last_null_pointer_err("GDALGetRasterBand")); + } + Ok(RasterBand::new(self.api, c_band, self)) + } + + /// Fetch the dataset geotransform coefficients. + /// Return an error if no geotransform is available. + pub fn geo_transform(&self) -> Result<[f64; 6]> { + let mut gt = [0.0f64; 6]; + let rv = unsafe { + call_gdal_api!( + self.api, + GDALGetGeoTransform, + self.c_dataset, + gt.as_mut_ptr() + ) + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(gt) + } + + /// Set the geo-transform. + pub fn set_geo_transform(&self, gt: &[f64; 6]) -> Result<()> { + let rv = unsafe { + call_gdal_api!( + self.api, + GDALSetGeoTransform, + self.c_dataset, + gt.as_ptr() as *mut f64 + ) + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Fetch the projection definition string for this dataset. + /// Return an empty string if no projection is available. + pub fn projection(&self) -> String { + unsafe { + let ptr = call_gdal_api!(self.api, GDALGetProjectionRef, self.c_dataset); + if ptr.is_null() { + String::new() + } else { + CStr::from_ptr(ptr).to_string_lossy().into_owned() + } + } + } + + /// Set the projection definition string for this dataset. + pub fn set_projection(&self, projection: &str) -> Result<()> { + let c_projection = CString::new(projection)?; + let rv = unsafe { + call_gdal_api!( + self.api, + GDALSetProjection, + self.c_dataset, + c_projection.as_ptr() + ) + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Fetch the spatial reference for this dataset. + /// GDAL returns a borrowed handle; this method clones it. + pub fn spatial_ref(&self) -> Result { + let c_srs = unsafe { call_gdal_api!(self.api, GDALGetSpatialRef, self.c_dataset) }; + if c_srs.is_null() { + return Err(self.api.last_null_pointer_err("GDALGetSpatialRef")); + } + // GDALGetSpatialRef returns a borrowed reference — clone it via OSRClone. + unsafe { SpatialRef::from_c_srs_clone(self.api, c_srs) } + } + + /// Set the spatial reference. + pub fn set_spatial_ref(&self, srs: &SpatialRef) -> Result<()> { + let rv = + unsafe { call_gdal_api!(self.api, GDALSetSpatialRef, self.c_dataset, srs.c_srs()) }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Create a copy of this dataset to a new file using the given driver. + pub fn create_copy( + &self, + driver: &Driver, + filename: &str, + options: &[&str], + ) -> Result { + let c_filename = CString::new(filename)?; + let csl = CslStringList::try_from_iter(options.iter().copied())?; + + let c_ds = unsafe { + call_gdal_api!( + self.api, + GDALCreateCopy, + driver.c_driver(), + c_filename.as_ptr(), + self.c_dataset, + 0, // bStrict + csl.as_ptr(), + ptr::null_mut(), + ptr::null_mut() + ) + }; + if c_ds.is_null() { + return Err(self.api.last_cpl_err(CE_Failure as u32)); + } + Ok(Dataset { + api: self.api, + c_dataset: c_ds, + owned: true, + }) + } + + /// Create a new vector layer. + pub fn create_layer(&self, options: LayerOptions<'_>) -> Result> { + let c_name = CString::new(options.name)?; + let c_srs = options.srs.map_or(ptr::null_mut(), |s| s.c_srs()); + + let csl = CslStringList::try_from_iter(options.options.unwrap_or(&[]).iter().copied())?; + + let c_layer = unsafe { + call_gdal_api!( + self.api, + GDALDatasetCreateLayer, + self.c_dataset, + c_name.as_ptr(), + c_srs, + options.ty, + csl.as_ptr() + ) + }; + if c_layer.is_null() { + return Err(self.api.last_null_pointer_err("GDALDatasetCreateLayer")); + } + Ok(Layer::new(self.api, c_layer, self)) + } + + /// Get the GDAL API reference. + pub fn api(&self) -> &'static GdalApi { + self.api + } + + /// Open a dataset using a `DatasetOptions` struct (georust-compatible convenience). + pub fn open_ex_with_options( + api: &'static GdalApi, + path: &str, + options: DatasetOptions<'_>, + ) -> Result { + Self::open_ex( + api, + path, + options.open_flags, + options.allowed_drivers, + options.open_options, + options.sibling_files, + ) + } + + /// Add a band backed by an existing memory buffer. + /// Pass `DATAPOINTER`, `PIXELOFFSET`, and `LINEOFFSET` to `GDALAddBand`. + /// + /// # Safety + /// + /// `data_ptr` must point to valid band data that outlives this dataset. + pub unsafe fn add_band_with_data( + &self, + data_type: RustGdalDataType, + data_ptr: *const u8, + pixel_offset: Option, + line_offset: Option, + ) -> Result<()> { + let data_pointer = format!("DATAPOINTER={data_ptr:p}"); + + let mut options = CslStringList::with_capacity(3); + options.add_string(&data_pointer)?; + + if let Some(pixel) = pixel_offset { + options.set_name_value("PIXELOFFSET", &pixel.to_string())?; + } + + if let Some(line) = line_offset { + options.set_name_value("LINEOFFSET", &line.to_string())?; + } + + let err = call_gdal_api!( + self.api, + GDALAddBand, + self.c_dataset, + data_type.to_c(), + options.as_ptr() + ); + if err != CE_None { + return Err(self.api.last_cpl_err(err as u32)); + } + Ok(()) + } + + /// Mark this dataset as owning its handle (for `Drop`). + pub fn set_owned(&mut self, owned: bool) { + self.owned = owned; + } +} + +/// Options for creating a vector layer. +pub struct LayerOptions<'a> { + pub name: &'a str, + pub srs: Option<&'a SpatialRef>, + pub ty: OGRwkbGeometryType, + /// Additional driver-specific options, in the form `"name=value"`. + pub options: Option<&'a [&'a str]>, +} + +impl Default for LayerOptions<'_> { + fn default() -> Self { + Self { + name: "", + srs: None, + ty: OGRwkbGeometryType::wkbUnknown, + options: None, + } + } +} + +#[cfg(all(test, feature = "gdal-sys"))] +mod tests { + use crate::driver::DriverManager; + use crate::global::with_global_gdal_api; + + #[test] + fn test_geo_transform_roundtrip() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create("", 256, 256, 1).unwrap(); + + let gt = [0.0, 1.0, 0.0, 0.0, 0.0, -1.0]; + ds.set_geo_transform(>).unwrap(); + let got = ds.geo_transform().unwrap(); + assert_eq!(gt, got); + }) + .unwrap(); + } + + #[test] + fn test_geo_transform_unset() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create("", 256, 256, 1).unwrap(); + + // MEM driver without an explicit set_geo_transform returns an error + assert!(ds.geo_transform().is_err()); + }) + .unwrap(); + } + + #[test] + fn test_set_projection_roundtrip() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create("", 256, 256, 1).unwrap(); + + let wkt = r#"GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]"#; + ds.set_projection(wkt).unwrap(); + let got = ds.projection(); + // The returned WKT may be reformatted by GDAL, so just check it contains WGS 84 + assert!(got.contains("WGS 84"), "Expected WGS 84 in: {got}"); + }) + .unwrap(); + } + + #[test] + fn test_dataset_raster_count() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + + let ds1 = driver.create("", 64, 64, 1).unwrap(); + assert_eq!(ds1.raster_count(), 1); + + let ds3 = driver.create("", 64, 64, 3).unwrap(); + assert_eq!(ds3.raster_count(), 3); + }) + .unwrap(); + } + + #[test] + fn test_dataset_raster_size() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create("", 123, 456, 1).unwrap(); + assert_eq!(ds.raster_size(), (123, 456)); + }) + .unwrap(); + } +} diff --git a/c/sedona-gdal/src/driver.rs b/c/sedona-gdal/src/driver.rs new file mode 100644 index 000000000..c9ed92c2d --- /dev/null +++ b/c/sedona-gdal/src/driver.rs @@ -0,0 +1,242 @@ +// 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: +//! . +//! Original code is licensed under MIT. + +use std::ffi::CString; +use std::ptr; + +use crate::dataset::Dataset; +use crate::errors::Result; +use crate::gdal_api::{call_gdal_api, GdalApi}; +use crate::gdal_dyn_bindgen::*; +use crate::raster::types::GdalDataType as RustGdalDataType; +use crate::raster::types::GdalType; + +/// A GDAL driver. +pub struct Driver { + api: &'static GdalApi, + c_driver: GDALDriverH, +} + +impl Driver { + /// Wrap an existing C driver handle. + /// + /// # Safety + /// + /// The caller must ensure the handle is valid. + pub unsafe fn from_c_driver(api: &'static GdalApi, c_driver: GDALDriverH) -> Self { + Self { api, c_driver } + } + + /// Return the raw C driver handle. + pub fn c_driver(&self) -> GDALDriverH { + self.c_driver + } + + /// Create a new raster dataset (with u8 band type). + pub fn create( + &self, + filename: &str, + size_x: usize, + size_y: usize, + bands: usize, + ) -> Result { + self.create_with_band_type::(filename, size_x, size_y, bands) + } + + /// Create a new raster dataset with a specific band type. + pub fn create_with_band_type( + &self, + filename: &str, + size_x: usize, + size_y: usize, + bands: usize, + ) -> Result { + let c_filename = CString::new(filename)?; + let x: i32 = size_x.try_into()?; + let y: i32 = size_y.try_into()?; + let b: i32 = bands.try_into()?; + let c_ds = unsafe { + call_gdal_api!( + self.api, + GDALCreate, + self.c_driver, + c_filename.as_ptr(), + x, + y, + b, + T::gdal_ordinal(), + ptr::null_mut() + ) + }; + if c_ds.is_null() { + return Err(self.api.last_cpl_err(CE_Failure as u32)); + } + Ok(Dataset::new_owned(self.api, c_ds)) + } + + /// Create a new raster dataset with a runtime data type. + /// + /// Unlike [`create_with_band_type`](Self::create_with_band_type), this accepts a + /// [`GdalDataType`](RustGdalDataType) enum value instead of a compile-time generic, + /// which is useful when the data type is only known at runtime. + pub fn create_with_data_type( + &self, + filename: &str, + size_x: usize, + size_y: usize, + bands: usize, + data_type: RustGdalDataType, + ) -> Result { + let c_filename = CString::new(filename)?; + let x: i32 = size_x.try_into()?; + let y: i32 = size_y.try_into()?; + let b: i32 = bands.try_into()?; + let c_ds = unsafe { + call_gdal_api!( + self.api, + GDALCreate, + self.c_driver, + c_filename.as_ptr(), + x, + y, + b, + data_type.to_c(), + ptr::null_mut() + ) + }; + if c_ds.is_null() { + return Err(self.api.last_cpl_err(CE_Failure as u32)); + } + Ok(Dataset::new_owned(self.api, c_ds)) + } + + /// Create a new dataset (vector-only, no raster bands). + pub fn create_vector_only(&self, filename: &str) -> Result { + let c_filename = CString::new(filename)?; + let c_ds = unsafe { + call_gdal_api!( + self.api, + GDALCreate, + self.c_driver, + c_filename.as_ptr(), + 1, + 1, + 0, + GDALDataType::GDT_Unknown, + ptr::null_mut() + ) + }; + if c_ds.is_null() { + return Err(self.api.last_cpl_err(CE_Failure as u32)); + } + Ok(Dataset::new_owned(self.api, c_ds)) + } +} + +/// Driver manager for looking up drivers by name. +pub struct DriverManager; + +impl DriverManager { + pub fn get_driver_by_name(api: &'static GdalApi, name: &str) -> Result { + let c_name = CString::new(name)?; + let c_driver = unsafe { call_gdal_api!(api, GDALGetDriverByName, c_name.as_ptr()) }; + if c_driver.is_null() { + return Err(api.last_null_pointer_err("GDALGetDriverByName")); + } + Ok(Driver { api, c_driver }) + } +} + +#[cfg(all(test, feature = "gdal-sys"))] +mod tests { + use crate::driver::DriverManager; + use crate::errors::GdalError; + use crate::global::with_global_gdal_api; + use crate::raster::types::GdalDataType; + + #[test] + fn test_get_driver_by_name() { + with_global_gdal_api(|api| { + let gtiff = DriverManager::get_driver_by_name(api, "GTiff").unwrap(); + assert!(!gtiff.c_driver().is_null()); + let mem = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + assert!(!mem.c_driver().is_null()); + }) + .unwrap(); + } + + #[test] + fn test_get_driver_by_name_invalid() { + with_global_gdal_api(|api| { + let err = DriverManager::get_driver_by_name(api, "NO_SUCH_DRIVER"); + assert!(matches!(err, Err(GdalError::NullPointer { .. }))); + }) + .unwrap(); + } + + #[test] + fn test_driver_create() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create("", 32, 16, 2).unwrap(); + assert_eq!(ds.raster_size(), (32, 16)); + assert_eq!(ds.raster_count(), 2); + }) + .unwrap(); + } + + #[test] + fn test_driver_create_with_band_type() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create_with_band_type::("", 10, 20, 1).unwrap(); + assert_eq!(ds.raster_count(), 1); + let ds = driver.create_with_band_type::("", 10, 20, 2).unwrap(); + assert_eq!(ds.raster_count(), 2); + let ds = driver.create_with_band_type::("", 10, 20, 3).unwrap(); + assert_eq!(ds.raster_count(), 3); + }) + .unwrap(); + } + + #[test] + fn test_driver_create_with_data_type() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver + .create_with_data_type("", 8, 8, 1, GdalDataType::UInt16) + .unwrap(); + assert_eq!(ds.raster_count(), 1); + }) + .unwrap(); + } + + #[test] + fn test_driver_create_vector_only() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let ds = driver.create_vector_only("").unwrap(); + assert_eq!(ds.raster_count(), 0); + assert_eq!(ds.raster_size(), (1, 1)); + }) + .unwrap(); + } +} diff --git a/c/sedona-gdal/src/errors.rs b/c/sedona-gdal/src/errors.rs index f344667ec..3ff5d03db 100644 --- a/c/sedona-gdal/src/errors.rs +++ b/c/sedona-gdal/src/errors.rs @@ -63,6 +63,9 @@ pub enum GdalError { #[error(transparent)] IntConversionError(#[from] TryFromIntError), + + #[error("Buffer length {0} does not match raster size {1:?}")] + BufferSizeMismatch(usize, (usize, usize)), } pub type Result = std::result::Result; diff --git a/c/sedona-gdal/src/gdal.rs b/c/sedona-gdal/src/gdal.rs new file mode 100644 index 000000000..e1afbd1f5 --- /dev/null +++ b/c/sedona-gdal/src/gdal.rs @@ -0,0 +1,237 @@ +// 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. + +//! High-level convenience wrapper around [`GdalApi`]. +//! +//! [`Gdal`] bundles a `&'static GdalApi` reference and exposes ergonomic +//! methods that delegate to the lower-level constructors and free functions +//! scattered across the crate, eliminating the need to pass `api` explicitly +//! at every call site. + +use crate::config; +use crate::dataset::Dataset; +use crate::driver::{Driver, DriverManager}; +use crate::errors::Result; +use crate::gdal_api::GdalApi; +use crate::gdal_dyn_bindgen::{GDALOpenFlags, OGRFieldType}; +use crate::mem::create_mem_dataset; +use crate::raster::polygonize::{polygonize, PolygonizeOptions}; +use crate::raster::rasterband::RasterBand; +use crate::raster::rasterize::{rasterize, RasterizeOptions}; +use crate::raster::rasterize_affine::rasterize_affine; +use crate::raster::types::DatasetOptions; +use crate::raster::types::GdalDataType; +use crate::spatial_ref::SpatialRef; +use crate::vector::feature::FieldDefn; +use crate::vector::geometry::Geometry; +use crate::vector::layer::Layer; +use crate::vrt::VrtDataset; +use crate::vsi; + +/// High-level convenience wrapper around [`GdalApi`]. +/// +/// Stores a `&'static GdalApi` reference and provides ergonomic methods that +/// delegate to the various constructors and free functions in the crate. +pub struct Gdal { + api: &'static GdalApi, +} + +impl Gdal { + /// Create a `Gdal` wrapper for the given API reference. + pub(crate) fn new(api: &'static GdalApi) -> Self { + Self { api } + } + + // -- Info ---------------------------------------------------------------- + + /// Return the name of the loaded GDAL library. + pub fn name(&self) -> &str { + self.api.name() + } + + /// Fetch GDAL version information for a standard request key. + pub fn version_info(&self, request: &str) -> String { + self.api.version_info(request) + } + + // -- Config -------------------------------------------------------------- + + /// Set a thread-local GDAL configuration option. + pub fn set_thread_local_config_option(&self, key: &str, value: &str) -> Result<()> { + config::set_thread_local_config_option(self.api, key, value) + } + + // -- Driver -------------------------------------------------------------- + + /// Fetch a GDAL driver by its short name. + pub fn get_driver_by_name(&self, name: &str) -> Result { + DriverManager::get_driver_by_name(self.api, name) + } + + // -- Dataset ------------------------------------------------------------- + + /// Open a dataset with extended options. + pub fn open_ex( + &self, + path: &str, + open_flags: GDALOpenFlags, + allowed_drivers: Option<&[&str]>, + open_options: Option<&[&str]>, + sibling_files: Option<&[&str]>, + ) -> Result { + Dataset::open_ex( + self.api, + path, + open_flags, + allowed_drivers, + open_options, + sibling_files, + ) + } + + /// Open a dataset using a [`DatasetOptions`] struct. + pub fn open_ex_with_options(&self, path: &str, options: DatasetOptions<'_>) -> Result { + Dataset::open_ex_with_options(self.api, path, options) + } + + // -- Spatial Reference --------------------------------------------------- + + /// Create a spatial reference from a WKT string. + pub fn spatial_ref_from_wkt(&self, wkt: &str) -> Result { + SpatialRef::from_wkt(self.api, wkt) + } + + // -- VRT ----------------------------------------------------------------- + + /// Create an empty VRT dataset with the given raster size. + pub fn create_vrt(&self, x_size: usize, y_size: usize) -> Result { + VrtDataset::create(self.api, x_size, y_size) + } + + // -- Geometry ------------------------------------------------------------ + + /// Create a geometry from WKB bytes. + pub fn geometry_from_wkb(&self, wkb: &[u8]) -> Result { + Geometry::from_wkb(self.api, wkb) + } + + /// Create a geometry from a WKT string. + pub fn geometry_from_wkt(&self, wkt: &str) -> Result { + Geometry::from_wkt(self.api, wkt) + } + + // -- Vector -------------------------------------------------------------- + + /// Create an OGR field definition. + pub fn create_field_defn(&self, name: &str, field_type: OGRFieldType) -> Result { + FieldDefn::new(self.api, name, field_type) + } + + // -- VSI (Virtual File System) ------------------------------------------- + + /// Create a VSI in-memory file from the given bytes. + pub fn create_mem_file(&self, file_name: &str, data: &[u8]) -> Result<()> { + vsi::create_mem_file(self.api, file_name, data) + } + + /// Delete a VSI in-memory file. + pub fn unlink_mem_file(&self, file_name: &str) -> Result<()> { + vsi::unlink_mem_file(self.api, file_name) + } + + /// Copy the bytes of a VSI in-memory file, taking ownership of the GDAL buffer. + pub fn get_vsi_mem_file_bytes_owned(&self, file_name: &str) -> Result> { + vsi::get_vsi_mem_file_bytes_owned(self.api, file_name) + } + + // -- Raster operations --------------------------------------------------- + + /// Create a bare in-memory MEM dataset with GDAL-owned bands. + /// For a higher-level builder with external bands and metadata, use `MemDatasetBuilder`. + pub fn create_mem_dataset( + &self, + width: usize, + height: usize, + n_owned_bands: usize, + owned_bands_data_type: GdalDataType, + ) -> Result { + create_mem_dataset( + self.api, + width, + height, + n_owned_bands, + owned_bands_data_type, + ) + } + + /// Rasterize geometries using the dataset geotransform as the transformer. + pub fn rasterize_affine( + &self, + dataset: &Dataset, + bands: &[usize], + geometries: &[Geometry], + burn_values: &[f64], + all_touched: bool, + ) -> Result<()> { + rasterize_affine( + self.api, + dataset, + bands, + geometries, + burn_values, + all_touched, + ) + } + + /// Rasterize geometries into the selected dataset bands. + pub fn rasterize( + &self, + dataset: &Dataset, + band_list: &[i32], + geometries: &[&Geometry], + burn_values: &[f64], + options: Option, + ) -> Result<()> { + rasterize( + self.api, + dataset, + band_list, + geometries, + burn_values, + options, + ) + } + + /// Polygonize a raster band into a vector layer. + pub fn polygonize( + &self, + src_band: &RasterBand<'_>, + mask_band: Option<&RasterBand<'_>>, + out_layer: &Layer<'_>, + pixel_value_field: i32, + options: &PolygonizeOptions, + ) -> Result<()> { + polygonize( + self.api, + src_band, + mask_band, + out_layer, + pixel_value_field, + options, + ) + } +} diff --git a/c/sedona-gdal/src/gdal_api.rs b/c/sedona-gdal/src/gdal_api.rs index c4baa6b19..cef854d4c 100644 --- a/c/sedona-gdal/src/gdal_api.rs +++ b/c/sedona-gdal/src/gdal_api.rs @@ -33,6 +33,7 @@ use crate::gdal_dyn_bindgen::SedonaGdalApi; /// initialization of [`GdalApi`] via [`GdalApi::try_from_shared_library`] or /// [`GdalApi::try_from_current_process`], and you cannot obtain a `&GdalApi` /// without successful initialization. +#[macro_export] macro_rules! call_gdal_api { ($api:expr, $func:ident $(, $arg:expr)*) => { if let Some(func) = $api.inner.$func { diff --git a/c/sedona-gdal/src/global.rs b/c/sedona-gdal/src/global.rs index 9365ef89a..5b9581163 100644 --- a/c/sedona-gdal/src/global.rs +++ b/c/sedona-gdal/src/global.rs @@ -16,6 +16,7 @@ // under the License. use crate::errors::GdalInitLibraryError; +use crate::gdal::Gdal; use crate::gdal_api::GdalApi; use std::path::PathBuf; use std::sync::{Mutex, OnceLock}; @@ -203,6 +204,16 @@ where Ok(func(api)) } +/// Execute a closure with the process-global high-level [`Gdal`] handle. +/// The global API is initialized lazily on first use and then reused. +pub fn with_global_gdal(func: F) -> Result +where + F: FnOnce(&Gdal) -> R, +{ + let api = get_global_gdal_api()?; + Ok(func(&Gdal::new(api))) +} + /// Verify that the GDAL library meets the minimum version requirement. /// /// We use `GDALVersionInfo("VERSION_NUM")` instead of `GDALCheckVersion` because diff --git a/c/sedona-gdal/src/lib.rs b/c/sedona-gdal/src/lib.rs index 0646a241a..6f8f70275 100644 --- a/c/sedona-gdal/src/lib.rs +++ b/c/sedona-gdal/src/lib.rs @@ -23,14 +23,19 @@ pub mod gdal_dyn_bindgen; pub mod errors; // --- Core API --- +pub mod gdal; pub mod gdal_api; pub mod global; // --- High-level wrappers --- pub mod config; pub mod cpl; +pub mod dataset; +pub mod driver; pub mod geo_transform; +pub mod mem; pub mod raster; pub mod spatial_ref; pub mod vector; +pub mod vrt; pub mod vsi; diff --git a/c/sedona-gdal/src/mem.rs b/c/sedona-gdal/src/mem.rs new file mode 100644 index 000000000..9c2701697 --- /dev/null +++ b/c/sedona-gdal/src/mem.rs @@ -0,0 +1,460 @@ +// 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. + +//! High-level builder for creating in-memory (MEM) GDAL datasets. +//! +//! [`MemDatasetBuilder`] provides a fluent, type-safe API for constructing GDAL MEM +//! datasets with zero-copy band attachment, optional geo-transform, projection, and +//! per-band nodata values. +//! +//! # Example +//! +//! ```rust,ignore +//! use sedona_gdal::global::with_global_gdal; +//! use sedona_gdal::mem::{MemDatasetBuilder, Nodata}; +//! use sedona_gdal::GdalDataType; +//! +//! with_global_gdal(|gdal| { +//! let data: Vec = vec![0u8; 256 * 256]; +//! let dataset = unsafe { +//! MemDatasetBuilder::new(256, 256) +//! .add_band(GdalDataType::UInt8, data.as_ptr()) +//! .geo_transform([0.0, 1.0, 0.0, 0.0, 0.0, -1.0]) +//! .projection("EPSG:4326") +//! .build(gdal) +//! .unwrap() +//! }; +//! assert_eq!(dataset.raster_count(), 1); +//! }).unwrap(); +//! ``` + +use crate::dataset::Dataset; +use crate::errors::Result; +use crate::gdal::Gdal; +use crate::gdal_api::{call_gdal_api, GdalApi}; +use crate::gdal_dyn_bindgen::CE_Failure; +use crate::raster::types::GdalDataType; + +/// Nodata value for a raster band. +/// +/// GDAL has three separate APIs for setting nodata depending on the band data type: +/// - [`f64`] for most types (UInt8 through Float64, excluding Int64/UInt64) +/// - [`i64`] for Int64 bands +/// - [`u64`] for UInt64 bands +/// +/// This enum encapsulates all three variants so callers don't need to match on +/// the band type when setting nodata. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum Nodata { + F64(f64), + I64(i64), + U64(u64), +} + +/// A band specification for [`MemDatasetBuilder`]. +struct MemBand { + data_type: GdalDataType, + data_ptr: *const u8, + pixel_offset: Option, + line_offset: Option, + nodata: Option, +} + +/// A builder for constructing in-memory (MEM) GDAL datasets. +/// +/// This creates datasets using `MEMDataset::Create` (bypassing GDAL's open-dataset-list +/// mutex for better concurrency) and attaches bands via `GDALAddBand` with `DATAPOINTER` +/// options for zero-copy operation. +/// +/// # Safety +/// +/// All `add_band*` methods are `unsafe` because the caller must ensure that the +/// provided data pointers remain valid for the lifetime of the built [`Dataset`]. +pub struct MemDatasetBuilder { + width: usize, + height: usize, + n_owned_bands: usize, + owned_bands_data_type: Option, + bands: Vec, + geo_transform: Option<[f64; 6]>, + projection: Option, +} + +impl MemDatasetBuilder { + /// Create a builder for a MEM dataset with the given dimensions. + pub fn new(width: usize, height: usize) -> Self { + Self { + width, + height, + n_owned_bands: 0, + owned_bands_data_type: None, + bands: Vec::new(), + geo_transform: None, + projection: None, + } + } + + /// Create a builder for a MEM dataset with GDAL-owned bands. + pub fn new_with_owned_bands( + width: usize, + height: usize, + n_owned_bands: usize, + owned_bands_data_type: GdalDataType, + ) -> Self { + Self { + width, + height, + n_owned_bands, + owned_bands_data_type: Some(owned_bands_data_type), + bands: Vec::new(), + geo_transform: None, + projection: None, + } + } + + /// Create a MEM dataset with GDAL-owned bands. + /// This is a safe shortcut for `new_with_owned_bands(...).build(gdal)`. + pub fn create( + gdal: &Gdal, + width: usize, + height: usize, + n_owned_bands: usize, + owned_bands_data_type: GdalDataType, + ) -> Result { + // SAFETY: `new_with_owned_bands` creates a builder with zero external bands, + // so no data pointers need to outlive the dataset. + unsafe { + Self::new_with_owned_bands(width, height, n_owned_bands, owned_bands_data_type) + .build(gdal) + } + } + + /// Add a zero-copy band from a raw data pointer. + /// Use default contiguous row-major pixel and line offsets. + /// + /// # Safety + /// + /// The caller must ensure `data_ptr` points to a valid buffer of at least + /// `height * width * data_type.byte_size()` bytes, and that the buffer + /// outlives the built [`Dataset`]. + pub unsafe fn add_band(self, data_type: GdalDataType, data_ptr: *const u8) -> Self { + self.add_band_with_options(data_type, data_ptr, None, None, None) + } + + /// Add a zero-copy band with custom offsets and optional nodata. + /// + /// # Safety + /// + /// The caller must ensure `data_ptr` points to a valid buffer of sufficient size + /// for the given dimensions and offsets, and that the buffer outlives the built + /// [`Dataset`]. + pub unsafe fn add_band_with_options( + mut self, + data_type: GdalDataType, + data_ptr: *const u8, + pixel_offset: Option, + line_offset: Option, + nodata: Option, + ) -> Self { + self.bands.push(MemBand { + data_type, + data_ptr, + pixel_offset, + line_offset, + nodata, + }); + self + } + + /// Set the dataset geotransform coefficients. + pub fn geo_transform(mut self, gt: [f64; 6]) -> Self { + self.geo_transform = Some(gt); + self + } + + /// Set the dataset projection definition string. + pub fn projection(mut self, wkt: impl Into) -> Self { + self.projection = Some(wkt.into()); + self + } + + /// Build the MEM dataset and attach the configured bands and metadata. + /// + /// # Safety + /// + /// This method is unsafe because the built dataset references memory provided via + /// the `add_band*` methods. The caller must ensure all data pointers remain valid + /// for the lifetime of the returned [`Dataset`]. + pub unsafe fn build(self, gdal: &Gdal) -> Result { + let dataset = gdal.create_mem_dataset( + self.width, + self.height, + self.n_owned_bands, + self.owned_bands_data_type.unwrap_or(GdalDataType::UInt8), + )?; + + // Attach bands (zero-copy via DATAPOINTER). + for band_spec in &self.bands { + dataset.add_band_with_data( + band_spec.data_type, + band_spec.data_ptr, + band_spec.pixel_offset, + band_spec.line_offset, + )?; + } + + // Set geo-transform. + if let Some(gt) = &self.geo_transform { + dataset.set_geo_transform(gt)?; + } + + // Set projection/CRS. + if let Some(proj) = &self.projection { + dataset.set_projection(proj)?; + } + + // Set per-band nodata values. + for (i, band_spec) in self.bands.iter().enumerate() { + if let Some(nodata) = &band_spec.nodata { + let raster_band = dataset.rasterband(i + 1 + self.n_owned_bands)?; + match nodata { + Nodata::F64(v) => raster_band.set_no_data_value(Some(*v))?, + Nodata::I64(v) => raster_band.set_no_data_value_i64(Some(*v))?, + Nodata::U64(v) => raster_band.set_no_data_value_u64(Some(*v))?, + } + } + } + + Ok(dataset) + } +} + +/// Create a bare in-memory MEM dataset with GDAL-owned bands. +/// For a higher-level builder with external bands and metadata, use `MemDatasetBuilder`. +pub(crate) fn create_mem_dataset( + api: &'static GdalApi, + width: usize, + height: usize, + n_owned_bands: usize, + owned_bands_data_type: GdalDataType, +) -> Result { + let empty_filename = c""; + let c_data_type = owned_bands_data_type.to_c(); + let handle = unsafe { + call_gdal_api!( + api, + MEMDatasetCreate, + empty_filename.as_ptr(), + width.try_into()?, + height.try_into()?, + n_owned_bands.try_into()?, + c_data_type, + std::ptr::null_mut() + ) + }; + + if handle.is_null() { + return Err(api.last_cpl_err(CE_Failure as u32)); + } + Ok(Dataset::new_owned(api, handle)) +} + +#[cfg(all(test, feature = "gdal-sys"))] +mod tests { + use crate::global::with_global_gdal; + use crate::mem::{MemDatasetBuilder, Nodata}; + use crate::raster::types::GdalDataType; + + #[test] + fn test_mem_builder_single_band() { + with_global_gdal(|gdal| { + let data = vec![42u8; 64 * 64]; + let dataset = unsafe { + MemDatasetBuilder::new(64, 64) + .add_band(GdalDataType::UInt8, data.as_ptr()) + .build(gdal) + .unwrap() + }; + assert_eq!(dataset.raster_size(), (64, 64)); + assert_eq!(dataset.raster_count(), 1); + }) + .unwrap(); + } + + #[test] + fn test_mem_builder_multi_band() { + with_global_gdal(|gdal| { + let band1 = vec![1u16; 32 * 32]; + let band2 = vec![2u16; 32 * 32]; + let band3 = vec![3u16; 32 * 32]; + let dataset = unsafe { + MemDatasetBuilder::new(32, 32) + .add_band(GdalDataType::UInt16, band1.as_ptr() as *const u8) + .add_band(GdalDataType::UInt16, band2.as_ptr() as *const u8) + .add_band(GdalDataType::UInt16, band3.as_ptr() as *const u8) + .build(gdal) + .unwrap() + }; + assert_eq!(dataset.raster_count(), 3); + }) + .unwrap(); + } + + #[test] + fn test_mem_builder_with_geo_transform() { + with_global_gdal(|gdal| { + let data = vec![0f32; 10 * 10]; + let gt = [100.0, 0.5, 0.0, 200.0, 0.0, -0.5]; + let dataset = unsafe { + MemDatasetBuilder::new(10, 10) + .add_band(GdalDataType::Float32, data.as_ptr() as *const u8) + .geo_transform(gt) + .build(gdal) + .unwrap() + }; + let got = dataset.geo_transform().unwrap(); + assert_eq!(gt, got); + }) + .unwrap(); + } + + #[test] + fn test_mem_builder_with_projection() { + with_global_gdal(|gdal| { + let data = [0u8; 8 * 8]; + let dataset = unsafe { + MemDatasetBuilder::new(8, 8) + .add_band(GdalDataType::UInt8, data.as_ptr()) + .projection(r#"GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]"#) + .build(gdal) + .unwrap() + }; + let proj = dataset.projection(); + assert!(proj.contains("WGS 84"), "Expected WGS 84 in: {proj}"); + }) + .unwrap(); + } + + #[test] + fn test_mem_builder_with_nodata() { + with_global_gdal(|gdal| { + let data = [0f64; 4 * 4]; + let dataset = unsafe { + MemDatasetBuilder::new(4, 4) + .add_band_with_options( + GdalDataType::Float64, + data.as_ptr() as *const u8, + None, + None, + Some(Nodata::F64(-9999.0)), + ) + .build(gdal) + .unwrap() + }; + let band = dataset.rasterband(1).unwrap(); + let nodata = band.no_data_value(); + assert_eq!(nodata, Some(-9999.0)); + }) + .unwrap(); + } + + #[test] + fn test_mem_builder_zero_bands() { + with_global_gdal(|gdal| { + let dataset = unsafe { MemDatasetBuilder::new(16, 16).build(gdal).unwrap() }; + assert_eq!(dataset.raster_count(), 0); + assert_eq!(dataset.raster_size(), (16, 16)); + }) + .unwrap(); + } + + #[test] + fn test_mem_builder_mixed_band_types() { + with_global_gdal(|gdal| { + let band_u8 = [0u8; 8 * 8]; + let band_f64 = vec![0f64; 8 * 8]; + let dataset = unsafe { + MemDatasetBuilder::new(8, 8) + .add_band(GdalDataType::UInt8, band_u8.as_ptr()) + .add_band(GdalDataType::Float64, band_f64.as_ptr() as *const u8) + .build(gdal) + .unwrap() + }; + assert_eq!(dataset.raster_count(), 2); + }) + .unwrap(); + } + + #[test] + pub fn test_mem_builder_with_owned_bands() { + with_global_gdal(|gdal| { + let dataset = unsafe { + MemDatasetBuilder::new_with_owned_bands(16, 16, 2, GdalDataType::UInt16) + .build(gdal) + .unwrap() + }; + assert_eq!(dataset.raster_count(), 2); + assert_eq!( + dataset.rasterband(1).unwrap().band_type(), + GdalDataType::UInt16 + ); + assert_eq!( + dataset.rasterband(2).unwrap().band_type(), + GdalDataType::UInt16 + ); + + let dataset = MemDatasetBuilder::create(gdal, 10, 8, 1, GdalDataType::Float32).unwrap(); + assert_eq!(dataset.raster_count(), 1); + assert_eq!( + dataset.rasterband(1).unwrap().band_type(), + GdalDataType::Float32 + ); + }) + .unwrap(); + } + + #[test] + pub fn test_mem_builder_mixed_owned_and_external_bands() { + with_global_gdal(|gdal| { + let external_band = [0u8; 8 * 8]; + let dataset = unsafe { + MemDatasetBuilder::new_with_owned_bands(8, 8, 1, GdalDataType::Float32) + .add_band_with_options( + GdalDataType::UInt8, + external_band.as_ptr(), + None, + None, + Some(Nodata::U64(255)), + ) + .build(gdal) + .unwrap() + }; + assert_eq!(dataset.raster_count(), 2); + assert_eq!( + dataset.rasterband(1).unwrap().band_type(), + GdalDataType::Float32 + ); + assert_eq!( + dataset.rasterband(2).unwrap().band_type(), + GdalDataType::UInt8 + ); + let nodata = dataset.rasterband(2).unwrap().no_data_value(); + assert_eq!(nodata, Some(255.0)); + }) + .unwrap(); + } +} diff --git a/c/sedona-gdal/src/raster.rs b/c/sedona-gdal/src/raster.rs index 1ddc9b2ed..839674298 100644 --- a/c/sedona-gdal/src/raster.rs +++ b/c/sedona-gdal/src/raster.rs @@ -15,4 +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 000000000..9291fd71a --- /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, + + /// 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, +} + +impl PolygonizeOptions { + /// Build a GDAL option list from these polygonize options. + pub fn to_options_list(&self) -> Result { + 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/rasterband.rs b/c/sedona-gdal/src/raster/rasterband.rs new file mode 100644 index 000000000..9e5376692 --- /dev/null +++ b/c/sedona-gdal/src/raster/rasterband.rs @@ -0,0 +1,372 @@ +// 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: +//! . +//! Original code is licensed under MIT. + +use std::marker::PhantomData; + +use crate::dataset::Dataset; +use crate::errors::{GdalError, Result}; +use crate::gdal_api::{call_gdal_api, GdalApi}; +use crate::raster::types::{Buffer, GdalType, ResampleAlg}; +use crate::{gdal_dyn_bindgen::*, raster::types::GdalDataType}; + +/// A raster band of a dataset. +pub struct RasterBand<'a> { + api: &'static GdalApi, + c_rasterband: GDALRasterBandH, + _dataset: PhantomData<&'a Dataset>, +} + +impl<'a> RasterBand<'a> { + pub(crate) fn new( + api: &'static GdalApi, + c_rasterband: GDALRasterBandH, + _dataset: &'a Dataset, + ) -> Self { + Self { + api, + c_rasterband, + _dataset: PhantomData, + } + } + + /// Return the raw C raster band handle. + pub fn c_rasterband(&self) -> GDALRasterBandH { + self.c_rasterband + } + + /// Read a window of this band into a typed buffer. + /// If `e_resample_alg` is `None`, use nearest-neighbour resampling. + pub fn read_as( + &self, + window: (isize, isize), + window_size: (usize, usize), + size: (usize, usize), + e_resample_alg: Option, + ) -> Result> { + let len = size.0 * size.1; + // Safety: all GdalType implementations are numeric primitives (u8, i8, u16, ..., f64), + // for which zeroed memory is a valid bit pattern. + let mut data: Vec = vec![unsafe { std::mem::zeroed() }; len]; + + let resample_alg = e_resample_alg.unwrap_or(ResampleAlg::NearestNeighbour); + let mut extra_arg = GDALRasterIOExtraArg { + eResampleAlg: resample_alg.to_gdal(), + ..GDALRasterIOExtraArg::default() + }; + + let rv = unsafe { + call_gdal_api!( + self.api, + GDALRasterIOEx, + self.c_rasterband, + GF_Read, + i32::try_from(window.0)?, + i32::try_from(window.1)?, + i32::try_from(window_size.0)?, + i32::try_from(window_size.1)?, + data.as_mut_ptr() as *mut std::ffi::c_void, + i32::try_from(size.0)?, + i32::try_from(size.1)?, + T::gdal_ordinal(), + 0, // nPixelSpace (auto) + 0, // nLineSpace (auto) + &mut extra_arg + ) + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + + Ok(Buffer::new(size, data)) + } + + /// Write a buffer to this raster band. + pub fn write( + &self, + window: (isize, isize), + window_size: (usize, usize), + buffer: &mut Buffer, + ) -> Result<()> { + let expected_len = buffer.shape.0 * buffer.shape.1; + if buffer.data.len() != expected_len { + return Err(GdalError::BufferSizeMismatch( + buffer.data.len(), + buffer.shape, + )); + } + let rv = unsafe { + call_gdal_api!( + self.api, + GDALRasterIO, + self.c_rasterband, + GF_Write, + i32::try_from(window.0)?, + i32::try_from(window.1)?, + i32::try_from(window_size.0)?, + i32::try_from(window_size.1)?, + buffer.data.as_mut_ptr() as *mut std::ffi::c_void, + i32::try_from(buffer.shape.0)?, + i32::try_from(buffer.shape.1)?, + T::gdal_ordinal(), + 0, // nPixelSpace (auto) + 0 // nLineSpace (auto) + ) + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Fetch this band's data type. + pub fn band_type(&self) -> GdalDataType { + GdalDataType::from_c(self.c_band_type()).unwrap_or(GdalDataType::Unknown) + } + + /// Fetch this band's raw GDAL data type. + pub fn c_band_type(&self) -> GDALDataType { + unsafe { call_gdal_api!(self.api, GDALGetRasterDataType, self.c_rasterband) } + } + + /// Fetch band size as `(x_size, y_size)`. + pub fn size(&self) -> (usize, usize) { + let x = unsafe { call_gdal_api!(self.api, GDALGetRasterBandXSize, self.c_rasterband) }; + let y = unsafe { call_gdal_api!(self.api, GDALGetRasterBandYSize, self.c_rasterband) }; + (x as usize, y as usize) + } + + /// Fetch the natural block size as `(x_size, y_size)`. + pub fn block_size(&self) -> (usize, usize) { + let mut x: i32 = 0; + let mut y: i32 = 0; + unsafe { + call_gdal_api!( + self.api, + GDALGetBlockSize, + self.c_rasterband, + &mut x, + &mut y + ) + }; + (x as usize, y as usize) + } + + /// Fetch the band's nodata value. + /// Return `None` if no nodata value is set. + pub fn no_data_value(&self) -> Option { + let mut success: i32 = 0; + let value = unsafe { + call_gdal_api!( + self.api, + GDALGetRasterNoDataValue, + self.c_rasterband, + &mut success + ) + }; + if success != 0 { + Some(value) + } else { + None + } + } + + /// Set the band's nodata value. + /// Pass `None` to clear any existing nodata value. + pub fn set_no_data_value(&self, value: Option) -> Result<()> { + let rv = if let Some(val) = value { + unsafe { call_gdal_api!(self.api, GDALSetRasterNoDataValue, self.c_rasterband, val) } + } else { + unsafe { call_gdal_api!(self.api, GDALDeleteRasterNoDataValue, self.c_rasterband) } + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Set the band's nodata value as `u64`. + /// Pass `None` to clear any existing nodata value. + pub fn set_no_data_value_u64(&self, value: Option) -> Result<()> { + let rv = if let Some(val) = value { + unsafe { + call_gdal_api!( + self.api, + GDALSetRasterNoDataValueAsUInt64, + self.c_rasterband, + val + ) + } + } else { + unsafe { call_gdal_api!(self.api, GDALDeleteRasterNoDataValue, self.c_rasterband) } + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Set the band's nodata value as `i64`. + /// Pass `None` to clear any existing nodata value. + pub fn set_no_data_value_i64(&self, value: Option) -> Result<()> { + let rv = if let Some(val) = value { + unsafe { + call_gdal_api!( + self.api, + GDALSetRasterNoDataValueAsInt64, + self.c_rasterband, + val + ) + } + } else { + unsafe { call_gdal_api!(self.api, GDALDeleteRasterNoDataValue, self.c_rasterband) } + }; + if rv != CE_None { + return Err(self.api.last_cpl_err(rv as u32)); + } + Ok(()) + } + + /// Get the GDAL API reference. + pub fn api(&self) -> &'static GdalApi { + self.api + } +} + +/// Return the actual block size for a block index. +/// Clamp edge blocks to the raster extent. +pub fn actual_block_size( + band: &RasterBand<'_>, + block_index: (usize, usize), +) -> Result<(usize, usize)> { + let (block_x, block_y) = band.block_size(); + let (raster_x, raster_y) = band.size(); + let x_off = block_index.0 * block_x; + let y_off = block_index.1 * block_y; + if x_off >= raster_x || y_off >= raster_y { + return Err(GdalError::BadArgument(format!( + "block index ({}, {}) is out of bounds for raster size ({}, {})", + block_index.0, block_index.1, raster_x, raster_y + ))); + } + let actual_x = if x_off + block_x > raster_x { + raster_x - x_off + } else { + block_x + }; + let actual_y = if y_off + block_y > raster_y { + raster_y - y_off + } else { + block_y + }; + Ok((actual_x, actual_y)) +} + +#[cfg(all(test, feature = "gdal-sys"))] +mod tests { + use crate::dataset::Dataset; + use crate::driver::DriverManager; + use crate::gdal_dyn_bindgen::*; + use crate::global::with_global_gdal_api; + use crate::raster::types::ResampleAlg; + + fn fixture(name: &str) -> String { + sedona_testing::data::test_raster(name).unwrap() + } + + #[test] + fn test_read_raster() { + with_global_gdal_api(|api| { + let path = fixture("tinymarble.tif"); + let dataset = Dataset::open_ex(api, &path, GDAL_OF_READONLY, None, None, None).unwrap(); + let rb = dataset.rasterband(1).unwrap(); + let rv = rb.read_as::((20, 30), (2, 3), (2, 3), None).unwrap(); + assert_eq!(rv.shape, (2, 3)); + assert_eq!(rv.data(), [7, 7, 7, 10, 8, 12]); + }) + .unwrap(); + } + + #[test] + fn test_read_raster_with_default_resample() { + with_global_gdal_api(|api| { + let path = fixture("tinymarble.tif"); + let dataset = Dataset::open_ex(api, &path, GDAL_OF_READONLY, None, None, None).unwrap(); + let rb = dataset.rasterband(1).unwrap(); + let rv = rb.read_as::((20, 30), (4, 4), (2, 2), None).unwrap(); + assert_eq!(rv.shape, (2, 2)); + // Default is NearestNeighbour; exact values are GDAL-version-dependent + // when downsampling from 4x4 to 2x2. Just verify shape and non-emptiness. + assert_eq!(rv.data().len(), 4); + }) + .unwrap(); + } + + #[test] + fn test_read_raster_with_average_resample() { + with_global_gdal_api(|api| { + let path = fixture("tinymarble.tif"); + let dataset = Dataset::open_ex(api, &path, GDAL_OF_READONLY, None, None, None).unwrap(); + let rb = dataset.rasterband(1).unwrap(); + let rv = rb + .read_as::((20, 30), (4, 4), (2, 2), Some(ResampleAlg::Average)) + .unwrap(); + assert_eq!(rv.shape, (2, 2)); + // Average resampling; exact values are GDAL-version-dependent, so just + // verify that the downsampled result has the expected shape and length. + assert_eq!(rv.data().len(), 4); + }) + .unwrap(); + } + + #[test] + fn test_get_no_data_value() { + with_global_gdal_api(|api| { + // tinymarble.tif has no nodata + let path = fixture("tinymarble.tif"); + let dataset = Dataset::open_ex(api, &path, GDAL_OF_READONLY, None, None, None).unwrap(); + let rb = dataset.rasterband(1).unwrap(); + assert!(rb.no_data_value().is_none()); + + // labels.tif has nodata=255 + let path = fixture("labels.tif"); + let dataset = Dataset::open_ex(api, &path, GDAL_OF_READONLY, None, None, None).unwrap(); + let rb = dataset.rasterband(1).unwrap(); + assert_eq!(rb.no_data_value(), Some(255.0)); + }) + .unwrap(); + } + + #[test] + #[allow(clippy::float_cmp)] + fn test_set_no_data_value() { + with_global_gdal_api(|api| { + let driver = DriverManager::get_driver_by_name(api, "MEM").unwrap(); + let dataset = driver.create("", 20, 10, 1).unwrap(); + let rasterband = dataset.rasterband(1).unwrap(); + assert_eq!(rasterband.no_data_value(), None); + assert!(rasterband.set_no_data_value(Some(1.23)).is_ok()); + assert_eq!(rasterband.no_data_value(), Some(1.23)); + assert!(rasterband.set_no_data_value(None).is_ok()); + assert_eq!(rasterband.no_data_value(), None); + }) + .unwrap(); + } +} diff --git a/c/sedona-gdal/src/raster/rasterize.rs b/c/sedona-gdal/src/raster/rasterize.rs new file mode 100644 index 000000000..a7a89a3c8 --- /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: +//! . +//! 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 { + 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, +) -> 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 = 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 = 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::((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 000000000..e4358f455 --- /dev/null +++ b/c/sedona-gdal/src/raster/rasterize_affine.rs @@ -0,0 +1,362 @@ +// 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_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_Failure, 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 = bands.iter().map(|&band| band as c_int).collect(); + + let c_options = if all_touched { + [c"ALL_TOUCHED=TRUE".as_ptr(), ptr::null_mut()] + } else { + [c"ALL_TOUCHED=FALSE".as_ptr(), ptr::null_mut()] + }; + + let geometries_c: Vec<_> = geometries.iter().map(|geo| geo.c_geometry()).collect(); + let burn_values_expanded: Vec = 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(), + (affine_transformer as *const ()).cast::() as *mut c_void, + (&mut arg as *mut AffineTransformArg).cast::(), + burn_values_expanded.as_ptr(), + c_options.as_ptr() as *mut *mut i8, + ptr::null_mut(), + ptr::null_mut() + ); + if error != CE_None { + return Err(api.last_cpl_err(CE_Failure 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 { + let driver = mem_driver(api); + let ds = driver.create_with_band_type::("", 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 { + let band = ds.rasterband(1).unwrap(); + let buf = band + .read_as::((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/vector.rs b/c/sedona-gdal/src/vector.rs index 10c038edc..52ff4e1bf 100644 --- a/c/sedona-gdal/src/vector.rs +++ b/c/sedona-gdal/src/vector.rs @@ -15,4 +15,6 @@ // specific language governing permissions and limitations // under the License. +pub mod feature; pub mod geometry; +pub mod layer; diff --git a/c/sedona-gdal/src/vector/feature.rs b/c/sedona-gdal/src/vector/feature.rs new file mode 100644 index 000000000..a80648b17 --- /dev/null +++ b/c/sedona-gdal/src/vector/feature.rs @@ -0,0 +1,213 @@ +// 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: +//! . +//! Original code is licensed under MIT. + +use std::ffi::CString; +use std::marker::PhantomData; + +use crate::errors::{GdalError, Result}; +use crate::gdal_api::{call_gdal_api, GdalApi}; +use crate::gdal_dyn_bindgen::*; +use crate::vector::geometry::Envelope; + +/// An OGR feature. +pub struct Feature<'a> { + api: &'static GdalApi, + c_feature: OGRFeatureH, + _lifetime: PhantomData<&'a ()>, +} + +impl Drop for Feature<'_> { + fn drop(&mut self) { + if !self.c_feature.is_null() { + unsafe { call_gdal_api!(self.api, OGR_F_Destroy, self.c_feature) }; + } + } +} + +impl<'a> Feature<'a> { + pub(crate) fn new(api: &'static GdalApi, c_feature: OGRFeatureH) -> Self { + Self { + api, + c_feature, + _lifetime: PhantomData, + } + } + + /// Fetch the feature geometry. + /// The returned geometry is borrowed; return `None` if no geometry is set. + pub fn geometry(&self) -> Option> { + let c_geom = unsafe { call_gdal_api!(self.api, OGR_F_GetGeometryRef, self.c_feature) }; + if c_geom.is_null() { + None + } else { + Some(BorrowedGeometry { + api: self.api, + c_geom, + _lifetime: PhantomData, + }) + } + } + + /// Fetch the index of a field by name. + /// Return an error if the field is not found. + pub fn field_index(&self, name: &str) -> Result { + let c_name = CString::new(name)?; + let idx = unsafe { + call_gdal_api!( + self.api, + OGR_F_GetFieldIndex, + self.c_feature, + c_name.as_ptr() + ) + }; + if idx < 0 { + return Err(GdalError::BadArgument(format!("field '{name}' not found"))); + } + Ok(idx) + } + + /// Fetch a field value as `f64`. + pub fn field_as_double(&self, field_index: i32) -> f64 { + unsafe { + call_gdal_api!( + self.api, + OGR_F_GetFieldAsDouble, + self.c_feature, + field_index + ) + } + } + + /// Fetch a field value as `i32`. + /// Return `None` if the field is unset or null. + pub fn field_as_integer(&self, field_index: i32) -> Option { + let is_set = unsafe { + call_gdal_api!( + self.api, + OGR_F_IsFieldSetAndNotNull, + self.c_feature, + field_index + ) + }; + if is_set != 0 { + Some(unsafe { + call_gdal_api!( + self.api, + OGR_F_GetFieldAsInteger, + self.c_feature, + field_index + ) + }) + } else { + None + } + } +} + +/// A geometry borrowed from a feature (not owned — will NOT be destroyed). +pub struct BorrowedGeometry<'a> { + api: &'static GdalApi, + c_geom: OGRGeometryH, + _lifetime: PhantomData<&'a ()>, +} + +impl<'a> BorrowedGeometry<'a> { + /// Return the raw C geometry handle. + pub fn c_geometry(&self) -> OGRGeometryH { + self.c_geom + } + + /// Export to ISO WKB. + pub fn wkb(&self) -> Result> { + let size = unsafe { call_gdal_api!(self.api, OGR_G_WkbSize, self.c_geom) }; + if size < 0 { + return Err(GdalError::BadArgument(format!( + "OGR_G_WkbSize returned negative size: {size}" + ))); + } + let mut buf = vec![0u8; size as usize]; + let rv = unsafe { + call_gdal_api!( + self.api, + OGR_G_ExportToIsoWkb, + self.c_geom, + wkbNDR, + buf.as_mut_ptr() + ) + }; + if rv != OGRERR_NONE { + return Err(GdalError::OgrError { + err: rv, + method_name: "OGR_G_ExportToIsoWkb", + }); + } + Ok(buf) + } + + /// Fetch the 2D envelope of this geometry. + pub fn envelope(&self) -> Envelope { + let mut env = OGREnvelope { + MinX: 0.0, + MaxX: 0.0, + MinY: 0.0, + MaxY: 0.0, + }; + unsafe { call_gdal_api!(self.api, OGR_G_GetEnvelope, self.c_geom, &mut env) }; + Envelope { + MinX: env.MinX, + MaxX: env.MaxX, + MinY: env.MinY, + MaxY: env.MaxY, + } + } +} + +/// An OGR field definition. +pub struct FieldDefn { + api: &'static GdalApi, + c_field_defn: OGRFieldDefnH, +} + +impl Drop for FieldDefn { + fn drop(&mut self) { + if !self.c_field_defn.is_null() { + unsafe { call_gdal_api!(self.api, OGR_Fld_Destroy, self.c_field_defn) }; + } + } +} + +impl FieldDefn { + /// Create a new field definition. + pub fn new(api: &'static GdalApi, name: &str, field_type: OGRFieldType) -> Result { + let c_name = CString::new(name)?; + let c_field_defn = + unsafe { call_gdal_api!(api, OGR_Fld_Create, c_name.as_ptr(), field_type) }; + if c_field_defn.is_null() { + return Err(api.last_null_pointer_err("OGR_Fld_Create")); + } + Ok(Self { api, c_field_defn }) + } + + /// Return the raw C handle. + pub fn c_field_defn(&self) -> OGRFieldDefnH { + self.c_field_defn + } +} diff --git a/c/sedona-gdal/src/vector/layer.rs b/c/sedona-gdal/src/vector/layer.rs new file mode 100644 index 000000000..8111f231f --- /dev/null +++ b/c/sedona-gdal/src/vector/layer.rs @@ -0,0 +1,208 @@ +// 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: +//! . +//! Original code is licensed under MIT. + +use std::marker::PhantomData; + +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::feature::{Feature, FieldDefn}; + +/// An OGR layer (borrowed from a Dataset). +pub struct Layer<'a> { + api: &'static GdalApi, + c_layer: OGRLayerH, + _dataset: PhantomData<&'a Dataset>, +} + +impl<'a> Layer<'a> { + pub(crate) fn new(api: &'static GdalApi, c_layer: OGRLayerH, _dataset: &'a Dataset) -> Self { + Self { + api, + c_layer, + _dataset: PhantomData, + } + } + + /// Return the raw C layer handle. + pub fn c_layer(&self) -> OGRLayerH { + self.c_layer + } + + /// Reset reading to the first feature. + pub fn reset_reading(&self) { + unsafe { call_gdal_api!(self.api, OGR_L_ResetReading, self.c_layer) }; + } + + /// Fetch the next feature from the current read cursor. + /// Return `None` when no more features are available. + pub fn next_feature(&self) -> Option> { + let c_feature = unsafe { call_gdal_api!(self.api, OGR_L_GetNextFeature, self.c_layer) }; + if c_feature.is_null() { + None + } else { + Some(Feature::new(self.api, c_feature)) + } + } + + /// Create a field in this layer. + /// Allow the driver to approximate the definition if needed. + pub fn create_field(&self, field_defn: &FieldDefn) -> Result<()> { + let rv = unsafe { + call_gdal_api!( + self.api, + OGR_L_CreateField, + self.c_layer, + field_defn.c_field_defn(), + 1 // bApproxOK + ) + }; + if rv != OGRERR_NONE { + return Err(GdalError::OgrError { + err: rv, + method_name: "OGR_L_CreateField", + }); + } + Ok(()) + } + + /// Fetch the feature count for this layer. + /// Return `-1` if the count is unknown and `force` is `false`. + pub fn feature_count(&self, force: bool) -> i64 { + unsafe { + call_gdal_api!( + self.api, + OGR_L_GetFeatureCount, + self.c_layer, + if force { 1 } else { 0 } + ) + } + } + + /// Iterate over features from the start of the layer. + /// This resets the layer read cursor before iteration. + pub fn features(&mut self) -> FeatureIterator<'_> { + self.reset_reading(); + FeatureIterator { layer: self } + } +} + +/// Iterator over features in a layer. +pub struct FeatureIterator<'a> { + layer: &'a Layer<'a>, +} + +impl<'a> Iterator for FeatureIterator<'a> { + type Item = Feature<'a>; + + fn next(&mut self) -> Option { + self.layer.next_feature() + } +} + +#[cfg(all(test, feature = "gdal-sys"))] +mod tests { + use gdal_sys::{ + GDALDatasetGetLayer, GDALDatasetGetLayerCount, OGR_F_Create, OGR_F_SetGeometry, + OGR_L_CreateFeature, OGR_L_GetLayerDefn, + }; + + use super::Layer; + use crate::dataset::{Dataset, LayerOptions}; + use crate::driver::DriverManager; + use crate::gdal_dyn_bindgen::OGRwkbGeometryType; + use crate::global::with_global_gdal_api; + use crate::vector::geometry::Geometry; + use crate::vsi::unlink_mem_file; + + #[test] + fn test_layer_iteration_and_reset() { + with_global_gdal_api(|api| { + let path = "/vsimem/test_layer_iteration.gpkg"; + let driver = DriverManager::get_driver_by_name(api, "GPKG").unwrap(); + let dataset = driver.create_vector_only(path).unwrap(); + + let layer = dataset + .create_layer(LayerOptions { + name: "features", + srs: None, + ty: OGRwkbGeometryType::wkbPoint, + options: None, + }) + .unwrap(); + + let layer_defn = unsafe { OGR_L_GetLayerDefn(layer.c_layer()) }; + assert!(!layer_defn.is_null()); + + for x in [1.0_f64, 2.0, 3.0] { + let feature = unsafe { OGR_F_Create(layer_defn) }; + assert!(!feature.is_null()); + + let geometry = Geometry::from_wkt(api, &format!("POINT ({x} 0)")).unwrap(); + let set_geometry_err = unsafe { OGR_F_SetGeometry(feature, geometry.c_geometry()) }; + assert_eq!(set_geometry_err, gdal_sys::OGRErr::OGRERR_NONE); + + let create_feature_err = unsafe { OGR_L_CreateFeature(layer.c_layer(), feature) }; + assert_eq!(create_feature_err, gdal_sys::OGRErr::OGRERR_NONE); + + unsafe { gdal_sys::OGR_F_Destroy(feature) }; + } + + let write_count = unsafe { GDALDatasetGetLayerCount(dataset.c_dataset()) }; + assert_eq!(write_count, 1); + + let read_dataset = Dataset::open_ex( + api, + path, + crate::gdal_dyn_bindgen::GDAL_OF_VECTOR | crate::gdal_dyn_bindgen::GDAL_OF_READONLY, + None, + None, + None, + ) + .unwrap(); + + let read_count = unsafe { GDALDatasetGetLayerCount(read_dataset.c_dataset()) }; + assert_eq!(read_count, 1); + + let c_layer = unsafe { GDALDatasetGetLayer(read_dataset.c_dataset(), 0) }; + assert!(!c_layer.is_null()); + let mut read_layer = Layer::new(api, c_layer, &read_dataset); + + assert_eq!(read_layer.feature_count(true), 3); + + let mut iter = read_layer.features(); + assert!(iter.next().is_some()); + assert!(iter.next().is_some()); + assert!(iter.next().is_some()); + assert!(iter.next().is_none()); + + read_layer.reset_reading(); + assert!(read_layer.next_feature().is_some()); + + assert_eq!(read_layer.features().count(), 3); + assert_eq!(read_layer.features().count(), 3); + + unlink_mem_file(api, path).unwrap(); + }) + .unwrap(); + } +} diff --git a/c/sedona-gdal/src/vrt.rs b/c/sedona-gdal/src/vrt.rs new file mode 100644 index 000000000..200c8fe50 --- /dev/null +++ b/c/sedona-gdal/src/vrt.rs @@ -0,0 +1,326 @@ +// 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_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, +} + +unsafe impl Send for VrtDataset {} + +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 { + 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_owned(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 { + 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> { + 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 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<'a>, + src_window: (i32, i32, i32, i32), + dst_window: (i32, i32, i32, i32), + resampling: Option<&str>, + nodata: Option, + ) -> Result<()> { + let c_resampling = resampling.map(CString::new).transpose()?; + + let resampling_ptr = c_resampling + .as_ref() + .map(|s| s.as_ptr()) + .unwrap_or(null_mut()); + + 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> 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::((0, 0), (1, 1), (1, 1), None) + .unwrap() + .data()[0]; + let vrt_px = vrt_band + .read_as::((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(); + } +} diff --git a/c/sedona-gdal/src/vsi.rs b/c/sedona-gdal/src/vsi.rs index 80022c7c9..20a60e11c 100644 --- a/c/sedona-gdal/src/vsi.rs +++ b/c/sedona-gdal/src/vsi.rs @@ -276,7 +276,7 @@ mod tests { let buffer = get_vsi_mem_file_buffer_owned(api, file_name).unwrap(); assert!(buffer.is_empty()); - assert_eq!(buffer.as_ref(), &[]); + assert_eq!(buffer.as_ref(), b""); }) .unwrap(); }