Merge #272

272: ENH: Add CoordTransform::transform_bounds() r=lnicola a=brendan-ward - [x] I agree to follow the project's [code of conduct](https://github.com/georust/gdal/blob/master/CODE_OF_CONDUCT.md). - [x] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users. --- This wraps `OCTTransformBounds` to provide `CoordTransform::transform_bounds()` for GDAL >= 3.4 Output values were verified manually using [rasterio](https://github.com/rasterio/rasterio). Co-authored-by: Brendan C. Ward <[email protected]>
georust · May 24, 2022 · fb6e260 · fb6e260
2 parents 16230bd + 3ed0474
commit fb6e260
Show file tree

Hide file tree

Showing 3 changed files with 118 additions and 0 deletions.
diff --git a/CHANGES.md b/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

diff --git a/src/spatial_ref/srs.rs b/src/spatial_ref/srs.rs
@@ -33,6 +33,58 @@ 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],
+                bounds[1],
+                bounds[2],
+                bounds[3],
+                &mut out_xmin,
+                &mut out_ymin,
+                &mut out_xmax,
+                &mut out_ymax,
+                densify_pts as c_int,
+            ) == 1
+        };
+
+        if !ret_val {
+            let msg = match _last_cpl_err(CPLErr::CE_Failure) {
+                GdalError::CplError { msg, .. } => match msg.is_empty() {
+                    false => Some(msg),
+                    _ => None,
+                },
+                err => return Err(err),
+            };
+            return Err(GdalError::InvalidCoordinateRange {
+                from: self.from.clone(),
+                to: self.to.clone(),
+                msg,
+            });
+        }
+
+        Ok([out_xmin, out_ymin, out_xmax, out_ymax])
+    }
+
     /// Transform coordinates in place.
     ///
     /// # Arguments

diff --git a/src/spatial_ref/tests.rs b/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();