ENH: Add CoordTransform::transform_bounds()

CHANGES.md

@@ -39,6 +39,10 @@
 
   - <https://github.com/georust/gdal/pull/273>
 
+- Add `gdal::srs::CoordTransform::transform_bounds` as wrapper for  `OCTTransformBounds` for GDAL 3.4
+
+  - <https://github.com/georust/gdal/pull/272>
+
 ## 0.12
 
 - Bump Rust edition to 2021

src/spatial_ref/srs.rs

@@ -1,6 +1,6 @@
 use crate::utils::{_last_cpl_err, _last_null_pointer_err, _string};
 use gdal_sys::{self, CPLErr, OGRCoordinateTransformationH, OGRErr, OGRSpatialReferenceH};
-use libc::{c_char, c_int};
+use libc::{c_char, c_double, c_int};
 use std::ffi::{CStr, CString};
 use std::ptr::{self, null_mut};
 use std::str::FromStr;
@@ -33,6 +33,61 @@ impl CoordTransform {
         })
     }
 
+    /// Transform bounding box, densifying the edges to account for nonlinear
+    /// transformations.
+    ///
+    /// # Arguments
+    /// * bounds - array of [axis0_min, axis1_min, axis0_max, axis1_min],
+    ///            interpreted in the axis order of the source SpatialRef,
+    ///            typically [xmin, ymin, xmax, ymax]
+    /// * densify_pts - number of points per edge (recommended: 21)
+    ///
+    /// # Returns
+    /// Some([f64; 4]) with bounds in axis order of target SpatialRef
+    /// None if there is an error.
+    #[cfg(all(major_ge_3, minor_ge_4))]
+    pub fn transform_bounds(&self, bounds: &[f64; 4], densify_pts: i32) -> Result<([f64; 4])> {
+        let mut out_xmin: f64 = 0.;
+        let mut out_ymin: f64 = 0.;
+        let mut out_xmax: f64 = 0.;
+        let mut out_ymax: f64 = 0.;
+
+        let ret_val = unsafe {
+            gdal_sys::OCTTransformBounds(
+                self.inner,
+                bounds[0] as c_double,
+                bounds[1] as c_double,
+                bounds[2] as c_double,
+                bounds[3] as c_double,
+                &mut out_xmin as *mut c_double,
+                &mut out_ymin as *mut c_double,
+                &mut out_xmax as *mut c_double,
+                &mut out_ymax as *mut c_double,
+                densify_pts as c_int,
+            ) == 1
+        };
+
+        if ret_val {
+            Ok([out_xmin, out_ymin, out_xmax, out_ymax])
+        } else {
+            let err = _last_cpl_err(CPLErr::CE_Failure);
+            let msg = if let GdalError::CplError { msg, .. } = err {
+                if msg.trim().is_empty() {
+                    None
+                } else {
+                    Some(msg)
+                }
+            } else {
+                return Err(err);
+            };
+            Err(GdalError::InvalidCoordinateRange {
+                from: self.from.clone(),
+                to: self.to.clone(),
+                msg,
+            })
+        }
+    }
+
     /// Transform coordinates in place.
     ///
     /// # Arguments

src/spatial_ref/tests.rs

@@ -64,6 +64,68 @@ fn comparison() {
     assert!(spatial_ref5 == spatial_ref4);
 }
 
+#[cfg(all(major_ge_3, minor_ge_4))]
+#[test]
+fn transform_bounds() {
+    let bounds: [f64; 4] = [-180., -80., 180., 80.];
+    // bounds for y,x ordered SpatialRefs
+    let yx_bounds: [f64; 4] = [-80.0, -180.0, 80.0, 180.];
+
+    let spatial_ref1 = SpatialRef::from_definition("OGC:CRS84").unwrap();
+
+    // transforming between the same SpatialRef should return existing bounds
+    let mut transform = CoordTransform::new(&spatial_ref1, &spatial_ref1).unwrap();
+    let mut out_bounds = transform.transform_bounds(&bounds, 21).unwrap();
+    assert_almost_eq(out_bounds[0], bounds[0]);
+    assert_almost_eq(out_bounds[1], bounds[1]);
+    assert_almost_eq(out_bounds[2], bounds[2]);
+    assert_almost_eq(out_bounds[3], bounds[3]);
+
+    // EPSG:4326 is in y,x order by default; returned bounds are [ymin, xmin, ymax, xmax]
+    let mut spatial_ref2 = SpatialRef::from_epsg(4326).unwrap();
+    transform = CoordTransform::new(&spatial_ref1, &spatial_ref2).unwrap();
+    out_bounds = transform.transform_bounds(&bounds, 21).unwrap();
+    assert_almost_eq(out_bounds[0], yx_bounds[0]);
+    assert_almost_eq(out_bounds[1], yx_bounds[1]);
+    assert_almost_eq(out_bounds[2], yx_bounds[2]);
+    assert_almost_eq(out_bounds[3], yx_bounds[3]);
+
+    // if source SpatialRef is in y,x order and and target SpatialRef is in x,y order
+    // input bounds are interpreted as [ymin, xmin, ymax, xmax] and returns
+    // [xmin, ymin, xmax, ymax]
+    transform = CoordTransform::new(&spatial_ref2, &spatial_ref1).unwrap();
+    out_bounds = transform.transform_bounds(&yx_bounds, 21).unwrap();
+    assert_almost_eq(out_bounds[0], bounds[0]);
+    assert_almost_eq(out_bounds[1], bounds[1]);
+    assert_almost_eq(out_bounds[2], bounds[2]);
+    assert_almost_eq(out_bounds[3], bounds[3]);
+
+    // force EPSG:4326 into x,y order to match source SpatialRef
+    spatial_ref2
+        .set_axis_mapping_strategy(gdal_sys::OSRAxisMappingStrategy::OAMS_TRADITIONAL_GIS_ORDER);
+    transform = CoordTransform::new(&spatial_ref1, &spatial_ref2).unwrap();
+    out_bounds = transform.transform_bounds(&bounds, 21).unwrap();
+    assert_almost_eq(out_bounds[0], bounds[0]);
+    assert_almost_eq(out_bounds[1], bounds[1]);
+    assert_almost_eq(out_bounds[2], bounds[2]);
+    assert_almost_eq(out_bounds[3], bounds[3]);
+
+    spatial_ref2 = SpatialRef::from_epsg(3857).unwrap();
+    transform = CoordTransform::new(&spatial_ref1, &spatial_ref2).unwrap();
+    out_bounds = transform.transform_bounds(&bounds, 21).unwrap();
+
+    let expected_bounds: [f64; 4] = [
+        -20037508.342789244,
+        -15538711.096309224,
+        20037508.342789244,
+        15538711.09630923,
+    ];
+    assert_almost_eq(out_bounds[0], expected_bounds[0]);
+    assert_almost_eq(out_bounds[1], expected_bounds[1]);
+    assert_almost_eq(out_bounds[2], expected_bounds[2]);
+    assert_almost_eq(out_bounds[3], expected_bounds[3]);
+}
+
 #[test]
 fn transform_coordinates() {
     let spatial_ref1 = SpatialRef::from_wkt("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],UNIT[\"DMSH\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",4326]]").unwrap();