From 616c95460482d696a1c571c1e5a3ac7b9d63a662 Mon Sep 17 00:00:00 2001 From: "Simeon H.K. Fitch" Date: Fri, 5 Jan 2024 15:51:16 -0500 Subject: [PATCH 1/2] Implemented raster::warp::reproject_into and raster::warp::create_and_reproject. --- CHANGES.md | 6 + fixtures/labels-nad.tif | Bin 0 -> 32199 bytes src/errors.rs | 2 + src/lib.rs | 1 + src/raster/buffer.rs | 1 + src/raster/{processing => }/dem/aspect.rs | 4 +- .../{processing => }/dem/color_relief.rs | 4 +- src/raster/{processing => }/dem/hillshade.rs | 6 +- src/raster/{processing => }/dem/mod.rs | 14 +- src/raster/{processing => }/dem/options.rs | 0 src/raster/{processing => }/dem/roughness.rs | 6 +- src/raster/{processing => }/dem/slope.rs | 6 +- src/raster/{processing => }/dem/tpi.rs | 6 +- src/raster/{processing => }/dem/tri.rs | 4 +- src/raster/mod.rs | 17 +- src/raster/processing/mod.rs | 3 - src/raster/warp.rs | 27 -- src/raster/warp/mod.rs | 254 +++++++++++++ src/raster/warp/reproject_options.rs | 228 +++++++++++ src/raster/warp/resample.rs | 72 ++++ src/raster/warp/warp_options.rs | 357 ++++++++++++++++++ src/xml.rs | 89 +++++ 22 files changed, 1047 insertions(+), 60 deletions(-) create mode 100644 fixtures/labels-nad.tif rename src/raster/{processing => }/dem/aspect.rs (97%) rename src/raster/{processing => }/dem/color_relief.rs (98%) rename src/raster/{processing => }/dem/hillshade.rs (97%) rename src/raster/{processing => }/dem/mod.rs (97%) rename src/raster/{processing => }/dem/options.rs (100%) rename src/raster/{processing => }/dem/roughness.rs (92%) rename src/raster/{processing => }/dem/slope.rs (96%) rename src/raster/{processing => }/dem/tpi.rs (92%) rename src/raster/{processing => }/dem/tri.rs (96%) delete mode 100644 src/raster/processing/mod.rs delete mode 100644 src/raster/warp.rs create mode 100644 src/raster/warp/mod.rs create mode 100644 src/raster/warp/reproject_options.rs create mode 100644 src/raster/warp/resample.rs create mode 100644 src/raster/warp/warp_options.rs create mode 100644 src/xml.rs diff --git a/CHANGES.md b/CHANGES.md index 6f4541a67..9581e1644 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -21,6 +21,12 @@ - Added `DriverManager::get_output_driver_for_dataset_name` and `DriverManager::get_output_drivers_for_dataset_name` for the ability to auto detect compatible `Driver`(s) for writing data. - + +- **Breaking**: Moved `raster::processing::dem` to `raster::dem` +- **Breaking**: Moved `raster::reproject` to `raster::warp::reproject_into` and added more options +- Added `raster::warp::create_and_reproject_image`. + + - - Added `Feature::unset_field` diff --git a/fixtures/labels-nad.tif b/fixtures/labels-nad.tif new file mode 100644 index 0000000000000000000000000000000000000000..b9be2c0baf7e0c01200d3ff976c39fdc2e819c06 GIT binary patch literal 32199 zcmeI0y=xRf7{=ed3o*tZDyW5uXDkGph%q2|#vF-gG}@enr6;EIg2Yr};TK$&LInST zAR>rKV__lw0aht3Ec`?Rwqj#^=dyb{dpA2f!)D{`<+*UT@6P+oyU+8R5Ryzz9;H4a z+C}6Q$>VC7-=1GT<+|Y4LRIu@q3f0Eon8<1@Y*ka@0Qm0`1KK9-^;boeB<|iNmG^T z?GFd2$o&iafCrBTsF&;OT)%1*sl@eTt`{2vRN#7<>!n74ws8H4>w60Yy2$lke&6R> zJnQ(#2!Hb0LGkrv-`5J!57vG^MD)2vR6aLRnYuJRSG!T0tBv-TO6AG9>3V6VR-ZmS zI6hUWR?k)kOSfzFTQ_;(*x=}Jc}OTm`^!VotwbMxue|?o{=}K!+I+C)C zR-PPKm=*KJo?W?Gxm=kzIeeE^@d5!5009sH0T2KI5C8!X009sH0T2KI5C8!X0DBFS(AJ2QAyV4TaYrACA>Fv+(9(&NAzE6_F_C+oi7Kro8l!cW zcC5eOYkRHJSnZfr2C2|)Flx3132u${rWHv+f->dEz`o*|a%K&*ibw_^9F}Z` zKpXNKUHUWTqfiY|oRC_m+Wi04Lfai4!d6&yh@4vpy{=V~!S9GE(C*iGSGKgnOqwFm+RH?3t5JH^FWNq$Lm%YHB$eK~vfE+MdLfd? zdbfUq(sIJ0exgZNN@OCp(Ij|Rwr%FeXtPVCbeyyp z+O5HAf?Kn;pUF&Yt}?K%keJWcDVvdWf>b?)?f9isPm(rsB` manages cell values in in raster I/O operations. /// /// It conceptually represents a 2-D array backed by a `Vec` with row-major organization diff --git a/src/raster/processing/dem/aspect.rs b/src/raster/dem/aspect.rs similarity index 97% rename from src/raster/processing/dem/aspect.rs rename to src/raster/dem/aspect.rs index 181bdfe03..99d6bb59f 100644 --- a/src/raster/processing/dem/aspect.rs +++ b/src/raster/dem/aspect.rs @@ -3,7 +3,7 @@ use std::num::NonZeroUsize; use super::options::{common_dem_options, CommonOptions}; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::DemSlopeAlg; +use crate::raster::dem::DemSlopeAlg; /// Configuration options for [`aspect()`][super::aspect()]. #[derive(Debug, Clone, Default)] @@ -71,7 +71,7 @@ mod tests { use crate::assert_near; use crate::cpl::CslStringList; use crate::errors::Result; - use crate::raster::processing::dem::aspect; + use crate::raster::dem::aspect; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/processing/dem/color_relief.rs b/src/raster/dem/color_relief.rs similarity index 98% rename from src/raster/processing/dem/color_relief.rs rename to src/raster/dem/color_relief.rs index 904646439..3b25fdf30 100644 --- a/src/raster/processing/dem/color_relief.rs +++ b/src/raster/dem/color_relief.rs @@ -3,7 +3,7 @@ use std::path::{Path, PathBuf}; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::options::common_dem_options; +use crate::raster::dem::options::common_dem_options; use super::options::CommonOptions; @@ -125,7 +125,7 @@ mod tests { use crate::assert_near; use crate::cpl::CslStringList; use crate::errors::Result; - use crate::raster::processing::dem::color_relief; + use crate::raster::dem::color_relief; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/processing/dem/hillshade.rs b/src/raster/dem/hillshade.rs similarity index 97% rename from src/raster/processing/dem/hillshade.rs rename to src/raster/dem/hillshade.rs index e4ebd9b6c..064c961b1 100644 --- a/src/raster/processing/dem/hillshade.rs +++ b/src/raster/dem/hillshade.rs @@ -2,8 +2,8 @@ use std::num::NonZeroUsize; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::options::common_dem_options; -use crate::raster::processing::dem::DemSlopeAlg; +use crate::raster::dem::options::common_dem_options; +use crate::raster::dem::DemSlopeAlg; use super::options::CommonOptions; @@ -160,7 +160,7 @@ mod tests { use crate::assert_near; use crate::cpl::CslStringList; use crate::errors::Result; - use crate::raster::processing::dem::hillshade; + use crate::raster::dem::hillshade; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/processing/dem/mod.rs b/src/raster/dem/mod.rs similarity index 97% rename from src/raster/processing/dem/mod.rs rename to src/raster/dem/mod.rs index 7557227ac..407de98a7 100644 --- a/src/raster/processing/dem/mod.rs +++ b/src/raster/dem/mod.rs @@ -75,7 +75,7 @@ mod tri; /// use gdal::Dataset; /// # fn main() -> gdal::errors::Result<()> { /// use std::path::Path; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; /// let mut opts = AspectOptions::new(); /// opts @@ -135,7 +135,7 @@ pub fn aspect>( /// use gdal::Dataset; /// # fn main() -> gdal::errors::Result<()> { /// use std::path::Path; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; /// let mut opts = ColorReliefOptions::new("fixtures/color-relief.clr"); /// opts.with_alpha(true); @@ -189,7 +189,7 @@ pub fn color_relief>( /// /// ```rust, no_run /// use gdal::Dataset; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// # fn main() -> gdal::errors::Result<()> { /// use std::path::Path; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; @@ -247,7 +247,7 @@ pub fn hillshade>( /// # fn main() -> gdal::errors::Result<()> { /// use gdal::Dataset; /// use std::path::Path; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; /// let roughness_ds = roughness(&ds, Path::new("target/dem-hills-roughness.tiff"), &RoughnessOptions::default())?; /// let stats = roughness_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); @@ -301,7 +301,7 @@ pub fn roughness>( /// # fn main() -> gdal::errors::Result<()> { /// use std::path::Path; /// use gdal::Dataset; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; /// let mut opts = SlopeOptions::new(); /// opts @@ -356,7 +356,7 @@ pub fn slope>( /// # fn main() -> gdal::errors::Result<()> { /// use std::path::Path; /// use gdal::Dataset; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; /// let tpi_ds = topographic_position_index(&ds, Path::new("target/dem-hills-tpi.tiff"), &TpiOptions::default())?; /// let stats = tpi_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); @@ -404,7 +404,7 @@ pub fn topographic_position_index>( /// # fn main() -> gdal::errors::Result<()> { /// use std::path::Path; /// use gdal::Dataset; -/// use gdal::raster::processing::dem::*; +/// use gdal::raster::dem::*; /// let ds = Dataset::open("fixtures/dem-hills.tiff")?; /// let mut opts = TriOptions::new(); /// opts.with_algorithm(DemTriAlg::Wilson); diff --git a/src/raster/processing/dem/options.rs b/src/raster/dem/options.rs similarity index 100% rename from src/raster/processing/dem/options.rs rename to src/raster/dem/options.rs diff --git a/src/raster/processing/dem/roughness.rs b/src/raster/dem/roughness.rs similarity index 92% rename from src/raster/processing/dem/roughness.rs rename to src/raster/dem/roughness.rs index 063c6fb3d..00e83b76e 100644 --- a/src/raster/processing/dem/roughness.rs +++ b/src/raster/dem/roughness.rs @@ -2,7 +2,7 @@ use std::num::NonZeroUsize; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::options::common_dem_options; +use crate::raster::dem::options::common_dem_options; use super::options::CommonOptions; @@ -35,8 +35,8 @@ mod tests { use crate::assert_near; use crate::cpl::CslStringList; use crate::errors::Result; - use crate::raster::processing::dem::roughness; - use crate::raster::processing::dem::roughness::RoughnessOptions; + use crate::raster::dem::roughness; + use crate::raster::dem::roughness::RoughnessOptions; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/processing/dem/slope.rs b/src/raster/dem/slope.rs similarity index 96% rename from src/raster/processing/dem/slope.rs rename to src/raster/dem/slope.rs index 568e7e0c8..25637db4e 100644 --- a/src/raster/processing/dem/slope.rs +++ b/src/raster/dem/slope.rs @@ -2,8 +2,8 @@ use std::num::NonZeroUsize; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::options::common_dem_options; -use crate::raster::processing::dem::DemSlopeAlg; +use crate::raster::dem::options::common_dem_options; +use crate::raster::dem::DemSlopeAlg; use super::options::CommonOptions; @@ -85,7 +85,7 @@ impl SlopeOptions { mod tests { use crate::cpl::CslStringList; use crate::errors::Result; - use crate::raster::processing::dem::slope; + use crate::raster::dem::slope; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/processing/dem/tpi.rs b/src/raster/dem/tpi.rs similarity index 92% rename from src/raster/processing/dem/tpi.rs rename to src/raster/dem/tpi.rs index 504e78016..ef6496579 100644 --- a/src/raster/processing/dem/tpi.rs +++ b/src/raster/dem/tpi.rs @@ -2,7 +2,7 @@ use std::num::NonZeroUsize; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::options::common_dem_options; +use crate::raster::dem::options::common_dem_options; use super::options::CommonOptions; @@ -34,8 +34,8 @@ mod tests { use crate::assert_near; use crate::cpl::CslStringList; use crate::errors::Result; - use crate::raster::processing::dem::topographic_position_index; - use crate::raster::processing::dem::tpi::TpiOptions; + use crate::raster::dem::topographic_position_index; + use crate::raster::dem::tpi::TpiOptions; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/processing/dem/tri.rs b/src/raster/dem/tri.rs similarity index 96% rename from src/raster/processing/dem/tri.rs rename to src/raster/dem/tri.rs index cff6024b9..0adf635a5 100644 --- a/src/raster/processing/dem/tri.rs +++ b/src/raster/dem/tri.rs @@ -2,7 +2,7 @@ use std::num::NonZeroUsize; use crate::cpl::CslStringList; use crate::errors; -use crate::raster::processing::dem::options::common_dem_options; +use crate::raster::dem::options::common_dem_options; use super::options::CommonOptions; @@ -80,7 +80,7 @@ impl DemTriAlg { mod tests { use crate::assert_near; use crate::errors::Result; - use crate::raster::processing::dem::terrain_ruggedness_index; + use crate::raster::dem::terrain_ruggedness_index; use crate::raster::StatisticsAll; use crate::test_utils::{fixture, InMemoryFixture}; use crate::Dataset; diff --git a/src/raster/mod.rs b/src/raster/mod.rs index 1fd9947f0..80ab06ebd 100644 --- a/src/raster/mod.rs +++ b/src/raster/mod.rs @@ -75,7 +75,6 @@ //! ``` pub use buffer::{Buffer, ByteBuffer}; -pub use create_options::RasterCreationOptions; #[cfg(all(major_ge_3, minor_ge_1))] pub use mdarray::{ Attribute, Dimension, ExtendedDataType, ExtendedDataTypeClass, Group, MDArray, MdStatisticsAll, @@ -86,16 +85,24 @@ pub use rasterband::{ }; pub use rasterize::{rasterize, BurnSource, MergeAlgorithm, OptimizeMode, RasterizeOptions}; pub use types::{AdjustedValue, GdalDataType, GdalType}; -pub use warp::reproject; mod buffer; -mod create_options; +pub mod dem; #[cfg(all(major_ge_3, minor_ge_1))] mod mdarray; -pub mod processing; mod rasterband; mod rasterize; #[cfg(test)] mod tests; mod types; -mod warp; +pub mod warp; + +/// Key/value pair for passing driver-specific creation options to +/// [`Driver::create_with_band_type_wth_options`](crate::Driver::create_with_band_type_with_options`). +/// +/// See `papszOptions` in [GDAL's `Create(...)` API documentation](https://gdal.org/api/gdaldriver_cpp.html#_CPPv4N10GDALDriver6CreateEPKciii12GDALDataType12CSLConstList). +#[derive(Debug)] +pub struct RasterCreationOption<'a> { + pub key: &'a str, + pub value: &'a str, +} diff --git a/src/raster/processing/mod.rs b/src/raster/processing/mod.rs deleted file mode 100644 index 6ec54fb82..000000000 --- a/src/raster/processing/mod.rs +++ /dev/null @@ -1,3 +0,0 @@ -//! GDAL processing routines. - -pub mod dem; diff --git a/src/raster/warp.rs b/src/raster/warp.rs deleted file mode 100644 index e803fb360..000000000 --- a/src/raster/warp.rs +++ /dev/null @@ -1,27 +0,0 @@ -use crate::dataset::Dataset; -use crate::utils::_last_cpl_err; -use gdal_sys::{self, CPLErr, GDALResampleAlg}; -use std::ptr::{null, null_mut}; - -use crate::errors::*; - -pub fn reproject(src: &Dataset, dst: &Dataset) -> Result<()> { - let rv = unsafe { - gdal_sys::GDALReprojectImage( - src.c_dataset(), - null(), - dst.c_dataset(), - null(), - GDALResampleAlg::GRA_Bilinear, - 0.0, - 0.0, - None, - null_mut(), - null_mut(), - ) - }; - if rv != CPLErr::CE_None { - return Err(_last_cpl_err(rv)); - } - Ok(()) -} diff --git a/src/raster/warp/mod.rs b/src/raster/warp/mod.rs new file mode 100644 index 000000000..cce886ce0 --- /dev/null +++ b/src/raster/warp/mod.rs @@ -0,0 +1,254 @@ +//! GDAL Warp API Bindings +//! +//! See also: +//! * [Warper C++ API](https://gdal.org/api/gdalwarp_cpp.html) +//! * [Warp API Tutorial](https://gdal.org/tutorials/warp_tut.html) +//! * [`gdalwarp` Program](https://gdal.org/programs/gdalwarp.html#gdalwarp) + +mod reproject_options; +mod resample; +mod warp_options; + +use gdal_sys::CPLErr; +pub use reproject_options::*; +pub use resample::*; +use std::ffi::CString; +use std::path::Path; +use std::ptr; +pub use warp_options::*; + +use crate::dataset::Dataset; +use crate::DriverManager; + +use crate::errors::*; +use crate::spatial_ref::SpatialRef; +use crate::utils::{_last_cpl_err, _path_to_c_string}; + +/// Reproject raster dataset into the given [`SpatialRef`] and save result to `dst_file`. +pub fn create_and_reproject>( + ds: &Dataset, + dst_file: P, + dst_srs: &SpatialRef, + options: &ReprojectOptions, +) -> Result<()> { + let dest_file = dst_file.as_ref(); + fn reproject( + src: &Dataset, + dst_file: &Path, + dst_srs: &SpatialRef, + options: &ReprojectOptions, + ) -> Result<()> { + let dest = _path_to_c_string(dst_file)?; + // Format the destination projection. + let dst_wkt = CString::new(dst_srs.to_wkt()?)?; + // Format the source projection, if specified. + let src_wkt = options + .src_spatial_ref() + .map(|s| s.to_wkt()) + .transpose()? + .map(CString::new) + .transpose()?; + let src_wkt_ptr = src_wkt.map(|s| s.as_ptr()).unwrap_or(ptr::null()); + + let driver = options + .output_format() + .as_ref() + .map(|f| DriverManager::get_driver_by_name(f)) + .transpose()? + .unwrap_or(src.driver()); + + let mut warp_options = options.clone_and_init_warp_options(src.raster_count())?; + + let rv = unsafe { + // See: https://github.com/OSGeo/gdal/blob/7b6c3fe71d61699abe66ea372bcd110701e38ff3/alg/gdalwarper.cpp#L235 + gdal_sys::GDALCreateAndReprojectImage( + src.c_dataset(), + src_wkt_ptr, + dest.as_ptr(), + dst_wkt.as_ptr(), + driver.c_driver(), + ptr::null_mut(), // create options + warp_options.resampling_alg().to_gdal(), + warp_options.memory_limit() as f64, + options.max_error().unwrap_or(0.0), + None, // progress fn + ptr::null_mut(), // progress arg + warp_options.as_ptr_mut(), + ) + }; + + if rv != CPLErr::CE_None { + return Err(_last_cpl_err(rv)); + } + + // See https://lists.osgeo.org/pipermail/gdal-dev/2023-November/057887.html for + // why this is required. To get around it We should rewrite this function to use the + // lower-level `GDALWarp` API. + if options.dst_nodata().is_some() { + let ds = Dataset::open(dst_file)?; + for b in 1..=ds.raster_count() { + let mut rb = ds.rasterband(b)?; + rb.set_no_data_value(options.dst_nodata())?; + } + } + + Ok(()) + } + reproject(ds, dest_file, dst_srs, options) +} + +/// Reproject one dataset into another dataset. +/// +/// Assumes destination dataset is properly sized and setup with a [`SpatialRef`], +/// [`GeoTransform`][crate::GeoTransform], [`RasterBand`][crate::raster::RasterBand], etc. +/// +/// See [`create_and_reproject`] for a more flexible alternative. +pub fn reproject_into( + src: &Dataset, + dst: &mut Dataset, + options: &ReprojectIntoOptions, +) -> Result<()> { + // Format the source projection, if specified. + let src_wkt = options + .src_spatial_ref() + .map(|s| s.to_wkt()) + .transpose()? + .map(CString::new) + .transpose()?; + let src_wkt_ptr = src_wkt.map(|s| s.as_ptr()).unwrap_or(ptr::null()); + + // Format the destination projection, if specified. + let dst_wkt = options + .src_spatial_ref() + .map(|s| s.to_wkt()) + .transpose()? + .map(CString::new) + .transpose()?; + let dst_wkt_ptr = dst_wkt.map(|s| s.as_ptr()).unwrap_or(ptr::null()); + + // GDALCreateAndReprojectImage requires a mutable pointer to + // an GDALWarpOptions instance. We could either propagate mutability up the call chain + // or clone the given options. Given the user may want to reuse settings for consistent + // application across multiple files and may find mutation unexpected, we clone make a clone. + let mut warp_options = options.clone_and_init_warp_options(dst.raster_count())?; + + let rv = unsafe { + gdal_sys::GDALReprojectImage( + src.c_dataset(), + src_wkt_ptr, + dst.c_dataset(), + dst_wkt_ptr, + warp_options.resampling_alg().to_gdal(), + warp_options.memory_limit() as f64, + options.max_error().unwrap_or(0.0), + None, // progress fn + ptr::null_mut(), // progress arg + warp_options.as_ptr_mut(), + ) + }; + if rv != CPLErr::CE_None { + return Err(_last_cpl_err(rv)); + } + + // See https://lists.osgeo.org/pipermail/gdal-dev/2023-November/057887.html for + // why this is required. To get around it We should rewrite this function to use the + // lower-level `GDALWarp` API. + if options.dst_nodata().is_some() { + for b in 1..=dst.raster_count() { + let mut rb = dst.rasterband(b)?; + rb.set_no_data_value(options.dst_nodata())?; + } + } + + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::errors::Result; + use crate::raster::GdalDataType; + use crate::spatial_ref::SpatialRef; + use crate::test_utils::{fixture, InMemoryFixture, TempFixture}; + use crate::{assert_near, Dataset}; + + // TODO: For some unknown reason this test fails on GDAL < 3.4 + #[cfg(any(all(major_ge_3, minor_ge_4), major_ge_4))] + #[test] + fn test_create_reproject() -> Result<()> { + use std::path::Path; + let dst_srs = SpatialRef::from_epsg(4269)?; + let source = TempFixture::fixture("labels.tif"); + + let dest = Path::new("target").join("labels-proj.tif"); + let ds = Dataset::open(&source)?; + + let mut opts = ReprojectOptions::default(); + opts.with_output_format("GTiff") + .with_dst_nodata(255.0) + .warp_options_mut() + .with_initial_value(InitValue::NoData) + .with_resampling_alg(WarpResampleAlg::NearestNeighbour); + + create_and_reproject(&ds, &dest, &dst_srs, &opts)?; + + let result = Dataset::open(dest)?; + let rb = result.rasterband(1)?; + let result_stats = rb.get_statistics(true, false)?.unwrap(); + + // Expected raster created with: + // gdalwarp -overwrite -t_srs EPSG:4269 -dstnodata 255 -r near -of GTiff fixtures/labels.tif fixtures/labels-nad.tif + let expected = Dataset::open(fixture("labels-nad.tif"))?; + let erb = expected.rasterband(1)?; + assert_eq!(erb.no_data_value(), Some(255.0)); + assert_eq!(erb.band_type(), GdalDataType::UInt8); + + let expected_stats = erb.get_statistics(true, false)?.unwrap(); + assert_near!(StatisticsAll, result_stats, expected_stats, epsilon = 1e-2); + + Ok(()) + } + + #[test] + fn test_reproject_into() -> Result<()> { + let source = TempFixture::fixture("labels.tif"); + let source_ds = Dataset::open(&source)?; + + let drv = DriverManager::get_driver_by_name("GTiff")?; + let outfile = InMemoryFixture::new("foo.tif"); + let dst_srs = SpatialRef::from_epsg(4269)?; + let mut dest_ds = drv.create_with_band_type::(outfile.path(), 210, 151, 1)?; + dest_ds.set_spatial_ref(&dst_srs)?; + dest_ds.set_geo_transform(&[ + -78.66496151541256, + 0.0003095182591293914, + 0.0, + 38.41639646432918, + 0.0, + -0.0003095182591293914, + ])?; + + let mut opts = ReprojectIntoOptions::default(); + opts.with_dst_nodata(255.0) + .warp_options_mut() + .with_initial_value(InitValue::NoData) + .with_resampling_alg(WarpResampleAlg::NearestNeighbour); + + reproject_into(&source_ds, &mut dest_ds, &opts)?; + + let rb = dest_ds.rasterband(1)?; + let result_stats = rb.get_statistics(true, false)?.unwrap(); + + // Expected raster created with: + // gdalwarp -overwrite -t_srs EPSG:4269 -dstnodata 255 -r near -of GTiff fixtures/labels.tif fixtures/labels-nad.tif + let expected = Dataset::open(fixture("labels-nad.tif"))?; + let erb = expected.rasterband(1)?; + assert_eq!(erb.no_data_value(), Some(255.0)); + assert_eq!(erb.band_type(), GdalDataType::UInt8); + + let expected_stats = erb.get_statistics(true, false)?.unwrap(); + assert_near!(StatisticsAll, result_stats, expected_stats, epsilon = 1e-2); + + Ok(()) + } +} diff --git a/src/raster/warp/reproject_options.rs b/src/raster/warp/reproject_options.rs new file mode 100644 index 000000000..885a5c4a6 --- /dev/null +++ b/src/raster/warp/reproject_options.rs @@ -0,0 +1,228 @@ +use crate::errors::Result; +use crate::raster::warp::GdalWarpOptions; +use crate::spatial_ref::SpatialRef; + +/// Injects methods associated with specifying warp no-data values. +macro_rules! nodata_accessors { + () => { + /// Specify the source no-data value. + /// + /// Overrides any no-data value specified in the source dataset. + pub fn with_src_nodata(&mut self, nodata_value: f64) -> &mut Self { + self.src_nodata = Some(nodata_value); + self + } + + /// Get the specified source no-data value, if any. + pub fn src_nodata(&self) -> Option { + self.src_nodata + } + + /// Specify the destination no-data value. + pub fn with_dst_nodata(&mut self, nodata_value: f64) -> &mut Self { + self.dst_nodata = Some(nodata_value); + self + } + + /// Get the specified destination no-data value, if any. + pub fn dst_nodata(&self) -> Option { + self.dst_nodata + } + }; +} + +/// Injects methods around [`GdalWarpOptions`]. +macro_rules! warp_options_accessors { + () => { + /// Set the general Warp options. + pub fn with_warp_options(&mut self, warp_options: GdalWarpOptions) -> &mut Self { + self.warp_options = warp_options; + self + } + + /// Fetch an immutable reference to the general Warp options. + pub fn warp_options(&self) -> &GdalWarpOptions { + &self.warp_options + } + + /// Fetch a mutable reference to the general Warp options. + pub fn warp_options_mut(&mut self) -> &mut GdalWarpOptions { + &mut self.warp_options + } + + /// Clone `warp_options` and initialize any required sub-structures. + /// + /// We clone because `GDALCreateAndReprojectImage` and siblings require a mutable pointer to + /// an GDALWarpOptions instance. We could either propagate mutability up the call chain + /// or clone the given options. Given the user may want to reuse settings for consistent + /// application across multiple files and may find mutation unexpected, we clone make a clone. + pub(crate) fn clone_and_init_warp_options( + &self, + band_count: usize, + ) -> Result { + let mut warp_options = self.warp_options().clone(); + + // If nodata values are specified, we need to initialize some state in + // `GdalWarpOptions` first. + if self.src_nodata().is_some() || self.dst_nodata().is_some() { + warp_options.with_band_count(band_count); + } + + if let Some(src_nodata) = self.src_nodata() { + warp_options.apply_src_nodata(src_nodata)?; + } + + if let Some(dst_nodata) = self.dst_nodata() { + warp_options.apply_dst_nodata(dst_nodata)?; + } + + if warp_options.working_datatype().is_none() { + warp_options.with_auto_working_datatype(); + } + + Ok(warp_options) + } + }; +} + +macro_rules! src_sr_accessors { + () => { + /// Set the source spatial reference system. + /// + /// If not specified here, the source [`SpatialRef`] is read from the source dataset. + /// + /// If specified here, any [`SpatialRef`] in the source dataset is overridden. + pub fn with_src_spatial_ref(&mut self, srs: SpatialRef) -> &mut Self { + self.src_srs = Some(srs); + self + } + + /// Fetch the source spatial reference system, if set. + pub fn src_spatial_ref(&self) -> Option<&SpatialRef> { + self.src_srs.as_ref() + } + }; +} + +/// Settings for [`create_and_reproject`][super::create_and_reproject]. +#[derive(Debug, Clone, Default)] +pub struct ReprojectOptions { + warp_options: GdalWarpOptions, + max_error: Option, + src_srs: Option, + src_nodata: Option, + dst_nodata: Option, + output_format: Option, +} + +impl ReprojectOptions { + pub fn new() -> Self { + Default::default() + } + + /// Set the maximum error. + /// + /// Measured in input pixels, it is the allowed in approximating + /// transformations. + /// + /// `0.0` indicates exact calculations. + pub fn with_max_error(&mut self, max_error: f64) -> &mut Self { + self.max_error = Some(max_error); + self + } + + /// Fetch the specified maximum error. + /// + /// Returns `None` if unset. + pub fn max_error(&self) -> Option { + self.max_error + } + + /// Explicitly specify output raster format. + /// + /// This is equivalent to the `-of ` CLI flag accepted by many GDAL tools. + /// + /// The value of `format` must be the identifier of a driver supported by the runtime + /// environment's GDAL library (e.g. `COG`, `JPEG`, `VRT`, etc.). A list of these identifiers + /// is available from `gdalinfo --formats`: + /// + /// ```text + /// ❯ gdalinfo --formats + /// Supported Formats: + /// VRT -raster,multidimensional raster- (rw+v): Virtual Raster + /// DERIVED -raster- (ro): Derived datasets using VRT pixel functions + /// GTiff -raster- (rw+vs): GeoTIFF + /// COG -raster- (wv): Cloud optimized GeoTIFF generator + /// NITF -raster- (rw+vs): National Imagery Transmission Format + /// ... + /// ``` + /// If no output format is specified, then the driver from the source dataset is used. + /// + pub fn with_output_format(&mut self, format: &str) -> &mut Self { + self.output_format = Some(format.to_owned()); + self + } + + /// Fetch the specified output format driver identifier, if any. + pub fn output_format(&self) -> Option { + self.output_format.clone() + } + + src_sr_accessors!(); + nodata_accessors!(); + warp_options_accessors!(); +} + +/// Settings for [`reproject_into`][super::reproject_into]. +#[derive(Debug, Clone, Default)] +pub struct ReprojectIntoOptions { + warp_options: GdalWarpOptions, + max_error: Option, + src_srs: Option, + dst_srs: Option, + src_nodata: Option, + dst_nodata: Option, +} + +impl ReprojectIntoOptions { + pub fn new() -> Self { + Default::default() + } + + /// Set the maximum error. + /// + /// Measured in input pixels, it is the allowed in approximating + /// transformations. + /// + /// `0.0` indicates exact calculations. + pub fn with_max_error(&mut self, max_error: f64) -> &mut Self { + self.max_error = Some(max_error); + self + } + + /// Fetch the specified maximum error. + /// + /// Returns `None` if unset. + pub fn max_error(&self) -> Option { + self.max_error + } + + /// Set the destination spatial reference system. + /// + /// If not specified here, the source [`SpatialRef`] is read from the destination dataset. + /// + /// If specified here, any [`SpatialRef`] in the destination dataset is overridden. + pub fn with_dst_spatial_ref(&mut self, srs: SpatialRef) -> &mut Self { + self.dst_srs = Some(srs); + self + } + + /// Fetch the destination spatial reference system, if set. + pub fn dst_spatial_ref(&self) -> Option<&SpatialRef> { + self.dst_srs.as_ref() + } + + src_sr_accessors!(); + nodata_accessors!(); + warp_options_accessors!(); +} diff --git a/src/raster/warp/resample.rs b/src/raster/warp/resample.rs new file mode 100644 index 000000000..df46df919 --- /dev/null +++ b/src/raster/warp/resample.rs @@ -0,0 +1,72 @@ +use crate::errors::{GdalError, Result}; +use gdal_sys::GDALResampleAlg; + +/// GDAL Warp Resampling Algorithm +#[derive(Debug, Copy, Clone, Default)] +#[repr(u32)] +#[allow(clippy::upper_case_acronyms)] +pub enum WarpResampleAlg { + /// Nearest neighbour (select on one input pixel) + #[default] + NearestNeighbour = GDALResampleAlg::GRA_NearestNeighbour, + /// Bilinear (2x2 kernel) + Bilinear = GDALResampleAlg::GRA_Bilinear, + /// Cubic Convolution Approximation (4x4 kernel) + Cubic = GDALResampleAlg::GRA_Cubic, + /// Cubic B-Spline Approximation (4x4 kernel) + CubicSpline = GDALResampleAlg::GRA_CubicSpline, + /// Lanczos windowed sinc interpolation (6x6 kernel) + Lanczos = GDALResampleAlg::GRA_Lanczos, + /// Average (computes the weighted average of all non-NODATA contributing pixels) + Average = GDALResampleAlg::GRA_Average, + /// Mode (selects the value which appears most often of all the sampled points) + Mode = GDALResampleAlg::GRA_Mode, + /// Max (selects maximum of all non-NODATA contributing pixels) + Max = GDALResampleAlg::GRA_Max, + /// Min (selects minimum of all non-NODATA contributing pixels) + Min = GDALResampleAlg::GRA_Min, + /// Med (selects median of all non-NODATA contributing pixels) + Med = GDALResampleAlg::GRA_Med, + /// Q1 (selects first quartile of all non-NODATA contributing pixels) + Q1 = GDALResampleAlg::GRA_Q1, + /// Q3 (selects third quartile of all non-NODATA contributing pixels) + Q3 = GDALResampleAlg::GRA_Q3, + #[cfg(any(all(major_ge_3, minor_ge_1), major_ge_4))] + /// Sum (weighed sum of all non-NODATA contributing pixels). + Sum = GDALResampleAlg::GRA_Sum, + #[cfg(any(all(major_ge_3, minor_ge_3), major_ge_4))] + /// RMS (weighted root mean square (quadratic mean) of all non-NODATA contributing pixels) + RMS = GDALResampleAlg::GRA_RMS, +} + +impl WarpResampleAlg { + pub fn to_gdal(self) -> GDALResampleAlg::Type { + self as GDALResampleAlg::Type + } + + pub fn from_gdal(alg: GDALResampleAlg::Type) -> Result { + Ok(match alg { + GDALResampleAlg::GRA_NearestNeighbour => WarpResampleAlg::NearestNeighbour, + GDALResampleAlg::GRA_Bilinear => WarpResampleAlg::Bilinear, + GDALResampleAlg::GRA_Cubic => WarpResampleAlg::Cubic, + GDALResampleAlg::GRA_CubicSpline => WarpResampleAlg::CubicSpline, + GDALResampleAlg::GRA_Lanczos => WarpResampleAlg::Lanczos, + GDALResampleAlg::GRA_Average => WarpResampleAlg::Average, + GDALResampleAlg::GRA_Mode => WarpResampleAlg::Mode, + GDALResampleAlg::GRA_Max => WarpResampleAlg::Max, + GDALResampleAlg::GRA_Min => WarpResampleAlg::Min, + GDALResampleAlg::GRA_Med => WarpResampleAlg::Med, + GDALResampleAlg::GRA_Q1 => WarpResampleAlg::Q1, + GDALResampleAlg::GRA_Q3 => WarpResampleAlg::Q3, + #[cfg(any(all(major_ge_3, minor_ge_1), major_ge_4))] + GDALResampleAlg::GRA_Sum => WarpResampleAlg::Sum, + #[cfg(any(all(major_ge_3, minor_ge_3), major_ge_4))] + GDALResampleAlg::GRA_RMS => WarpResampleAlg::RMS, + o => { + return Err(GdalError::BadArgument(format!( + "Ordinal {o} does not map to a supported WarpResampleAlg" + ))) + } + }) + } +} diff --git a/src/raster/warp/warp_options.rs b/src/raster/warp/warp_options.rs new file mode 100644 index 000000000..b70a6b909 --- /dev/null +++ b/src/raster/warp/warp_options.rs @@ -0,0 +1,357 @@ +use std::fmt::{Debug, Display, Formatter}; +use std::ptr::NonNull; + +use crate::cpl::CslStringList; +use crate::errors::{GdalError, Result}; +use crate::raster::warp::resample::WarpResampleAlg; +use crate::raster::GdalDataType; +use crate::utils::_last_null_pointer_err; +use crate::xml::GdalXmlNode; +use gdal_sys::{ + GDALCloneWarpOptions, GDALCreateWarpOptions, GDALDeserializeWarpOptions, + GDALDestroyWarpOptions, GDALSerializeWarpOptions, GDALWarpInitDefaultBandMapping, + GDALWarpInitDstNoDataReal, GDALWarpInitSrcNoDataReal, GDALWarpOptions, + GDALWarpResolveWorkingDataType, +}; +use libc::c_char; + +/// Container for options provided to GDAL Warp routines. +/// +/// See: [`GDALWarpOptions`](https://gdal.org/api/gdalwarp_cpp.html#_CPPv415GDALWarpOptions) +/// for details. +pub struct GdalWarpOptions(NonNull); + +impl GdalWarpOptions { + pub fn new() -> Self { + unsafe { Self::from_ptr(GDALCreateWarpOptions()) } + } + + /// Create Self from a raw pointer. + /// + /// # Safety + /// Caller is responsible for ensuring `ptr` is not null, and + /// ownership of `ptr` is properly transferred + pub unsafe fn from_ptr(ptr: *mut GDALWarpOptions) -> Self { + Self(NonNull::new_unchecked(ptr)) + } + + /// Specify the resampling algorithm to use in Warp operation. + pub fn with_resampling_alg(&mut self, alg: WarpResampleAlg) -> &mut Self { + unsafe { (*self.as_ptr_mut()).eResampleAlg = alg.to_gdal() }; + self + } + + /// Get the resampling algorithm to be used in Warp operation. + pub fn resampling_alg(&self) -> WarpResampleAlg { + // `unwrap` below is ok because `with_resampling_alg` is the only way it got set, + // aside from the GDAL default, which is `GRA_NearestNeighbour`. + WarpResampleAlg::from_gdal(unsafe { (*self.as_ptr()).eResampleAlg }).unwrap_or_default() + } + + /// Set the datatype used during processing. + /// + /// If unset, the algorithm picks the datatype. + pub fn with_working_datatype(&mut self, dt: GdalDataType) -> &mut Self { + unsafe { (*self.as_ptr_mut()).eWorkingDataType = dt.gdal_ordinal() }; + self + } + + /// Fetch working datatype, if specified. + pub fn working_datatype(&self) -> Option { + let c_dt = unsafe { (*self.as_ptr()).eWorkingDataType }; + let dt: GdalDataType = c_dt.try_into().ok()?; + + // Default is `Unknown`, so we consider that "unspecified". + if dt == GdalDataType::Unknown { + None + } else { + Some(dt) + } + } + + /// This option is a convenience method for setting `INIT_DEST` setting via [`extra_options_mut`]. + /// + /// It forces the destination image to be initialized to the indicated value (for all bands), + /// or indicates that it should be initialized to the band's no-data value. + /// + /// If this value isn't set the destination image will be read and overlaid. + /// + /// See: [`GDALWarpOptions::papszWarpOptions](https://gdal.org/api/gdalwarp_cpp.html#_CPPv4N15GDALWarpOptions16papszWarpOptionsE) + pub fn with_initial_value(&mut self, init: InitValue) -> &mut Self { + self.extra_options_mut() + .set_name_value("INIT_DEST", &init.to_string()) + .expect("Failed to set INIT_DEST"); // We can unwrap because we know the strings are valid. + self + } + + /// Fetch the initial value setting, if any. + /// + /// See [`with_initial_value`][Self::with_initial_value] for details. + pub fn initial_value(&self) -> Option { + let init = self.extra_options().fetch_name_value("INIT_DEST"); + + match init.as_deref() { + Some("NO_DATA") => Some(InitValue::NoData), + Some(s) => s.parse::().ok().map(InitValue::Value), + None => None, + } + } + + /// If the working data type is unknown, this method will determine a valid working + /// data type to support the data in the src and dest data sets and any noData values. + pub fn with_auto_working_datatype(&mut self) -> &mut Self { + unsafe { GDALWarpResolveWorkingDataType(self.as_ptr_mut()) }; + self + } + + /// Memory limit in in bytes, + /// + /// Use `0` to specify GDAL default. + pub fn with_memory_limit(&mut self, limit_bytes: usize) -> &mut Self { + unsafe { (*self.as_ptr_mut()).dfWarpMemoryLimit = limit_bytes as f64 }; + self + } + + /// Fetch the memory limit setting in bytes. + /// + /// Zero means use GDAL default. + pub fn memory_limit(&self) -> usize { + unsafe { (*self.as_ptr()).dfWarpMemoryLimit as usize } + } + + /// Number of bands to process + /// + /// `0` selects all bands. + pub fn with_band_count(&mut self, num_bands: usize) -> &mut Self { + unsafe { GDALWarpInitDefaultBandMapping(self.as_ptr_mut(), num_bands as libc::c_int) }; + self + } + + /// Get the specified number of bands to process + /// + /// `0` indicates all bands. + pub fn band_count(&mut self) -> usize { + let cnt = unsafe { (*self.as_ptr()).nBandCount }; + cnt as usize + } + + /// Sets the source Dataset no-data value. Internal use only. + /// + /// This method exists to facilitate delaying the specification of a homogeneous no-data value + /// until the number of bands is known (at the point of warp call), via `with_band_count`, + /// which initializes the band mapping, as required by GDALWarp. + /// + /// Therefore, the caller of this method is responsible for ensuring that a non-zero number of bands has been + /// first set via [`with_band_count`][Self::with_band_count]. If not, this method will return + /// `Err(GdalError::UnexpectedLogicError(...))` + pub(super) fn apply_src_nodata(&mut self, no_data_value: f64) -> Result<&mut Self> { + if self.band_count() == 0 { + return Err(GdalError::UnexpectedLogicError( + "Specification of source no-data value prior to initializing band mapping via `with_band_count`".into()) + ); + } + + // GDALWarpOptions destructor frees this. See: + // https://github.com/OSGeo/gdal/blob/a9635785a2db8f575328326f2b1833e743ec8828/alg/gdalwarper.cpp#L1293 + unsafe { GDALWarpInitSrcNoDataReal(self.as_ptr_mut(), no_data_value) }; + + Ok(self) + } + + /// Sets the destination Dataset no-data value. Internal use only. + /// + /// See [`Self::apply_src_nodata`] for additional details. + pub(super) fn apply_dst_nodata(&mut self, no_data_value: f64) -> Result<&mut Self> { + if self.band_count() == 0 { + return Err(GdalError::UnexpectedLogicError( + "Specification of destination no-data value prior to initializing band mapping via `with_band_count`".into()) + ); + } + + // The GDALWarpOptions destructor frees this. See: + // https://github.com/OSGeo/gdal/blob/a9635785a2db8f575328326f2b1833e743ec8828/alg/gdalwarper.cpp#L1295 + unsafe { GDALWarpInitDstNoDataReal(self.as_ptr_mut(), no_data_value) }; + + Ok(self) + } + + /// Get any extra options attached to the Warp options. + pub fn extra_options(&self) -> &CslStringList { + let opts_array: &*mut *mut c_char = unsafe { &(*self.as_ptr()).papszWarpOptions }; + // Proof that GDALWarpOptions owns the CslStringList, and we just need to wrap it: + // https://github.com/OSGeo/gdal/blob/a9635785a2db8f575328326f2b1833e743ec8828/alg/gdalwarper.cpp#L1290 + // `CslStringList` is `transparent` with a single field, so this should be ok. + unsafe { std::mem::transmute(opts_array) } + } + + /// Get a mutable reference to extra options attached to the Warp options. + pub fn extra_options_mut(&mut self) -> &mut CslStringList { + let opts_array: &*mut *mut c_char = unsafe { &(*self.as_ptr()).papszWarpOptions }; + // See `extra_options` for rationale on transmute. + let csl: *mut CslStringList = opts_array as *const *mut *mut i8 as *mut CslStringList; + // `unwrap` should be ok because `opts_array` points to an offset against `self`, and + // we can assume `self` is not null. + unsafe { csl.as_mut().unwrap() } + } + + /// Serialize settings to GDAL XML. + pub fn to_xml(&self) -> Result { + let c_xml = unsafe { GDALSerializeWarpOptions(self.as_ptr()) }; + Ok(unsafe { GdalXmlNode::from_ptr(c_xml) }) + } + + /// Deserialize options from GDAL XML + pub fn from_xml(xml: &GdalXmlNode) -> Result { + let c_opts = unsafe { GDALDeserializeWarpOptions(xml.as_ptr_mut()) }; + if c_opts.is_null() { + Err(_last_null_pointer_err("GDALDeserializeWarpOptions")) + } else { + Ok(unsafe { Self::from_ptr(c_opts) }) + } + } + + /// Get a immutable pointer to C API options. + pub fn as_ptr(&self) -> *const GDALWarpOptions { + self.0.as_ptr() + } + + /// Get a mutable pointer to C API options. + pub fn as_ptr_mut(&mut self) -> *mut GDALWarpOptions { + self.0.as_ptr() + } +} + +impl Clone for GdalWarpOptions { + fn clone(&self) -> Self { + unsafe { Self::from_ptr(GDALCloneWarpOptions(self.as_ptr())) } + } +} + +impl Drop for GdalWarpOptions { + fn drop(&mut self) { + unsafe { GDALDestroyWarpOptions(self.as_ptr_mut()) } + } +} + +impl Default for GdalWarpOptions { + fn default() -> Self { + Self::new() + } +} + +impl Debug for GdalWarpOptions { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + let xml = self.to_xml().map_err(|_| std::fmt::Error)?; + Debug::fmt(&xml, f) + } +} + +/// Specifies the initial value cells in the destination dataset during a warp operation. +/// +/// See [`GdalWarpOptions::with_initial_value`]. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum InitValue { + NoData, + Value(f64), +} + +impl From for InitValue { + fn from(v: f64) -> Self { + InitValue::Value(v) + } +} + +impl Display for InitValue { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let str = match self { + InitValue::NoData => "NO_DATA".to_string(), + InitValue::Value(v) => v.to_string(), + }; + write!(f, "{}", str) + } +} + +#[cfg(test)] +mod tests { + use crate::errors::Result; + use crate::raster::warp::*; + use crate::raster::GdalDataType; + + fn ops_str(ops: &GdalWarpOptions) -> String { + let xml = ops.to_xml().unwrap(); + xml.to_string() + } + + #[test] + fn defaults() { + let opts = GdalWarpOptions::default(); + assert!(ops_str(&opts).contains("NearestNeighbour")); + } + + #[test] + fn with_settings() -> Result<()> { + let mut opts = GdalWarpOptions::default(); + assert_eq!(opts.memory_limit(), 0); + opts.with_memory_limit(1 << 16) + .with_working_datatype(GdalDataType::UInt16) + .with_band_count(2) + .with_resampling_alg(WarpResampleAlg::Cubic); + + let extra = opts.extra_options_mut(); + extra.set_name_value("NUM_THREADS", "4")?; + extra.set_name_value("SOURCE_EXTRA", "2")?; + + assert_eq!(opts.memory_limit(), 1 << 16); + + let expected = r#" + 65536 + Cubic + UInt16 + + + + + + +"#; + assert_eq!(ops_str(&opts), expected); + + Ok(()) + } + + #[test] + fn band_count() -> Result<()> { + let mut opts = GdalWarpOptions::default(); + assert_eq!(opts.band_count(), 0); + opts.with_band_count(3); + assert_eq!(opts.band_count(), 3); + + assert!(!ops_str(&opts).contains("")); + + opts.apply_src_nodata(255.0)?; + + assert!(ops_str(&opts).contains("255")); + + opts.apply_dst_nodata(0.0)?; + + assert!(ops_str(&opts).contains("0")); + + Ok(()) + } + + #[test] + fn init_value() -> Result<()> { + let rendered = |v: &str| format!(""); + + let mut opts = GdalWarpOptions::default(); + assert_eq!(opts.initial_value(), None); + + opts.with_initial_value(InitValue::NoData); + assert!(ops_str(&opts).contains(&rendered("NO_DATA"))); + + opts.with_initial_value(InitValue::Value(255.0)); + assert!(ops_str(&opts).contains(&rendered("255"))); + + Ok(()) + } +} diff --git a/src/xml.rs b/src/xml.rs new file mode 100644 index 000000000..ade2af645 --- /dev/null +++ b/src/xml.rs @@ -0,0 +1,89 @@ +use std::ffi::CString; +use std::fmt::{Debug, Display, Formatter}; +use std::ptr::NonNull; +use std::str::FromStr; + +use gdal_sys::{ + CPLCloneXMLTree, CPLDestroyXMLNode, CPLErr, CPLParseXMLString, CPLSerializeXMLTree, CPLXMLNode, +}; + +use crate::utils::_last_cpl_err; + +/// An XML node, as captured from GDAL serialization APIs. +pub struct GdalXmlNode(NonNull); + +impl GdalXmlNode { + /// Create a Self from a raw pointer. + /// + /// # Safety + /// Caller is responsible for ensuring `ptr` is not null, and + /// ownership of `ptr` is properly transferred. + pub unsafe fn from_ptr(ptr: *mut CPLXMLNode) -> GdalXmlNode { + Self(NonNull::new_unchecked(ptr)) + } + + pub fn as_ptr(&self) -> *const CPLXMLNode { + self.0.as_ptr() + } + + pub fn as_ptr_mut(&self) -> *mut CPLXMLNode { + self.0.as_ptr() + } +} + +impl Clone for GdalXmlNode { + fn clone(&self) -> Self { + unsafe { GdalXmlNode::from_ptr(CPLCloneXMLTree(self.as_ptr())) } + } +} + +impl Drop for GdalXmlNode { + fn drop(&mut self) { + unsafe { CPLDestroyXMLNode(self.0.as_mut()) }; + } +} + +impl FromStr for GdalXmlNode { + type Err = crate::errors::GdalError; + + fn from_str(s: &str) -> Result { + let s = CString::new(s)?; + let c_xml = unsafe { CPLParseXMLString(s.as_ptr()) }; + if c_xml.is_null() { + Err(_last_cpl_err(CPLErr::CE_Failure)) + } else { + Ok(unsafe { GdalXmlNode::from_ptr(c_xml) }) + } + } +} + +impl Display for GdalXmlNode { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + let s = unsafe { CString::from_raw(CPLSerializeXMLTree(self.as_ptr_mut())) }; + f.write_str(s.to_string_lossy().trim_end()) + } +} + +impl Debug for GdalXmlNode { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + let mut s = self.to_string(); + if !f.alternate() { + // Flatten the display string to fit on one line. + s = s.replace(['\n', '\r', ' ', '\t'], ""); + } + f.write_str(&s) + } +} + +#[cfg(test)] +mod tests { + use crate::xml::GdalXmlNode; + + #[test] + fn serde() { + let src = r#"yeet"#; + let xml: GdalXmlNode = src.parse().unwrap(); + let s = xml.to_string(); + assert_eq!(src, s.trim()); + } +} From 74edd4851a41bee204a926ae910b89b3af93ce53 Mon Sep 17 00:00:00 2001 From: "Simeon H.K. Fitch" Date: Wed, 31 Jan 2024 15:47:48 -0500 Subject: [PATCH 2/2] Scaffolding for GDALWarp binding. --- src/raster/mod.rs | 12 +---- src/raster/warp/mod.rs | 96 +++++++++++++++++++++++++++++++-- src/raster/warp/warp_options.rs | 48 +++++++++++++++-- 3 files changed, 140 insertions(+), 16 deletions(-) diff --git a/src/raster/mod.rs b/src/raster/mod.rs index 80ab06ebd..24ee43e3f 100644 --- a/src/raster/mod.rs +++ b/src/raster/mod.rs @@ -75,6 +75,7 @@ //! ``` pub use buffer::{Buffer, ByteBuffer}; +pub use create_options::RasterCreationOptions; #[cfg(all(major_ge_3, minor_ge_1))] pub use mdarray::{ Attribute, Dimension, ExtendedDataType, ExtendedDataTypeClass, Group, MDArray, MdStatisticsAll, @@ -87,6 +88,7 @@ pub use rasterize::{rasterize, BurnSource, MergeAlgorithm, OptimizeMode, Rasteri pub use types::{AdjustedValue, GdalDataType, GdalType}; mod buffer; +mod create_options; pub mod dem; #[cfg(all(major_ge_3, minor_ge_1))] mod mdarray; @@ -96,13 +98,3 @@ mod rasterize; mod tests; mod types; pub mod warp; - -/// Key/value pair for passing driver-specific creation options to -/// [`Driver::create_with_band_type_wth_options`](crate::Driver::create_with_band_type_with_options`). -/// -/// See `papszOptions` in [GDAL's `Create(...)` API documentation](https://gdal.org/api/gdaldriver_cpp.html#_CPPv4N10GDALDriver6CreateEPKciii12GDALDataType12CSLConstList). -#[derive(Debug)] -pub struct RasterCreationOption<'a> { - pub key: &'a str, - pub value: &'a str, -} diff --git a/src/raster/warp/mod.rs b/src/raster/warp/mod.rs index cce886ce0..a560b77be 100644 --- a/src/raster/warp/mod.rs +++ b/src/raster/warp/mod.rs @@ -9,11 +9,11 @@ mod reproject_options; mod resample; mod warp_options; -use gdal_sys::CPLErr; +use gdal_sys::{CPLErr, GDALDatasetH, GDALWarp}; pub use reproject_options::*; pub use resample::*; use std::ffi::CString; -use std::path::Path; +use std::path::{Path, PathBuf}; use std::ptr; pub use warp_options::*; @@ -22,7 +22,7 @@ use crate::DriverManager; use crate::errors::*; use crate::spatial_ref::SpatialRef; -use crate::utils::{_last_cpl_err, _path_to_c_string}; +use crate::utils::{_last_cpl_err, _last_null_pointer_err, _path_to_c_string}; /// Reproject raster dataset into the given [`SpatialRef`] and save result to `dst_file`. pub fn create_and_reproject>( @@ -163,6 +163,79 @@ pub fn reproject_into( Ok(()) } +pub fn warp(source: &Dataset, dest: D, options: &GdalWarpOptions) -> Result +where + D: Into, +{ + warp_multiple(&[source], dest, options) +} + +pub fn warp_multiple(source: &[&Dataset], dest: D, options: &GdalWarpOptions) -> Result +where + D: Into, +{ + let app_opts = GdalWarpAppOptions::default(); + + if true { + todo!("how the hell do you go from {options:?} to GdalWarpAppOptions?"); + } + + let mut source = source.iter().map(|ds| ds.c_dataset()).collect::>(); + + let dest = dest.into(); + match dest { + WarpDestination::Dataset(ds) => { + let ds_c = unsafe { + GDALWarp( + ptr::null_mut(), + ds.c_dataset(), + source.len() as libc::c_int, + source.as_mut_ptr(), + app_opts.as_ptr(), + ptr::null_mut(), + ) + }; + if ds_c.is_null() { + Err(_last_null_pointer_err("GDALWarp")) + } else { + Ok(ds) + } + } + WarpDestination::Path(p) => { + let path = _path_to_c_string(&p)?; + let ds_c = unsafe { + GDALWarp( + path.as_ptr(), + ptr::null_mut(), + source.len() as libc::c_int, + source.as_ptr() as *mut GDALDatasetH, + app_opts.as_ptr(), + ptr::null_mut(), + ) + }; + Ok(unsafe { Dataset::from_c_dataset(ds_c) }) + } + } +} + +#[derive(Debug)] +pub enum WarpDestination { + Dataset(Dataset), + Path(PathBuf), +} + +impl From for WarpDestination { + fn from(ds: Dataset) -> Self { + WarpDestination::Dataset(ds) + } +} + +impl From for WarpDestination { + fn from(path: PathBuf) -> Self { + WarpDestination::Path(path) + } +} + #[cfg(test)] mod tests { use super::*; @@ -251,4 +324,21 @@ mod tests { Ok(()) } + + #[test] + #[ignore] + fn test_warp() -> Result<()> { + let source = TempFixture::fixture("labels.tif"); + let source_ds = Dataset::open(&source)?; + let dest = Path::new("target").join("labels-warp.tif"); + + let mut options = GdalWarpOptions::default(); + options + .with_band_count(source_ds.raster_count()) + .with_initial_value(InitValue::NoData) + .with_resampling_alg(WarpResampleAlg::NearestNeighbour); + + warp(&source_ds, dest, &options)?; + Ok(()) + } } diff --git a/src/raster/warp/warp_options.rs b/src/raster/warp/warp_options.rs index b70a6b909..fc9625bc8 100644 --- a/src/raster/warp/warp_options.rs +++ b/src/raster/warp/warp_options.rs @@ -1,4 +1,5 @@ use std::fmt::{Debug, Display, Formatter}; +use std::ptr; use std::ptr::NonNull; use crate::cpl::CslStringList; @@ -9,9 +10,9 @@ use crate::utils::_last_null_pointer_err; use crate::xml::GdalXmlNode; use gdal_sys::{ GDALCloneWarpOptions, GDALCreateWarpOptions, GDALDeserializeWarpOptions, - GDALDestroyWarpOptions, GDALSerializeWarpOptions, GDALWarpInitDefaultBandMapping, - GDALWarpInitDstNoDataReal, GDALWarpInitSrcNoDataReal, GDALWarpOptions, - GDALWarpResolveWorkingDataType, + GDALDestroyWarpOptions, GDALSerializeWarpOptions, GDALWarpAppOptions, GDALWarpAppOptionsFree, + GDALWarpAppOptionsNew, GDALWarpInitDefaultBandMapping, GDALWarpInitDstNoDataReal, + GDALWarpInitSrcNoDataReal, GDALWarpOptions, GDALWarpResolveWorkingDataType, }; use libc::c_char; @@ -246,6 +247,47 @@ impl Debug for GdalWarpOptions { } } +#[derive(Debug)] +pub struct GdalWarpAppOptions(NonNull); + +impl GdalWarpAppOptions { + /// Instatiate and empty set of warp options. + pub fn new() -> Self { + unsafe { Self::from_ptr(GDALWarpAppOptionsNew(ptr::null_mut(), ptr::null_mut())) } + } + + /// Create Self from a raw pointer. + /// + /// # Safety + /// Caller is responsible for ensuring `ptr` is not null, and + /// ownership of `ptr` is properly transferred + pub unsafe fn from_ptr(ptr: *mut GDALWarpAppOptions) -> Self { + Self(NonNull::new_unchecked(ptr)) + } + + /// Get a immutable pointer to C API options. + pub fn as_ptr(&self) -> *const GDALWarpAppOptions { + self.0.as_ptr() + } + + /// Get a mutable pointer to C API options. + pub fn as_ptr_mut(&mut self) -> *mut GDALWarpAppOptions { + self.0.as_ptr() + } +} + +impl Drop for GdalWarpAppOptions { + fn drop(&mut self) { + unsafe { GDALWarpAppOptionsFree(self.as_ptr_mut()) } + } +} + +impl Default for GdalWarpAppOptions { + fn default() -> Self { + Self::new() + } +} + /// Specifies the initial value cells in the destination dataset during a warp operation. /// /// See [`GdalWarpOptions::with_initial_value`].