From 6e9c7c74eb4a53fbfb54f072dcdd3fa5ee5c278e Mon Sep 17 00:00:00 2001 From: "Brendan C. Ward" Date: Thu, 12 May 2022 16:40:19 -0700 Subject: [PATCH 1/4] Add CoordTransform::transform_bounds() --- src/spatial_ref/srs.rs | 57 +++++++++++++++++++++++++++++++++++- src/spatial_ref/tests.rs | 62 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 118 insertions(+), 1 deletion(-) diff --git a/src/spatial_ref/srs.rs b/src/spatial_ref/srs.rs index b4368dd38..7eca66db7 100644 --- a/src/spatial_ref/srs.rs +++ b/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 diff --git a/src/spatial_ref/tests.rs b/src/spatial_ref/tests.rs index 03e723731..495e0cdc1 100644 --- 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(); From 9061914a31d1447fa1ea36030179d2e36ee74e4a Mon Sep 17 00:00:00 2001 From: "Brendan C. Ward" Date: Thu, 12 May 2022 16:45:34 -0700 Subject: [PATCH 2/4] Add changes --- CHANGES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index d4bdd844d..28bd36bec 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -23,6 +23,10 @@ - +- Add `gdal::srs::CoordTransform::transform_bounds` as wrapper for `OCTTransformBounds` for GDAL 3.4 + + - + ## 0.12 - Bump Rust edition to 2021 From 79464ab02961e21ad8779825b8a24796f8e5e03b Mon Sep 17 00:00:00 2001 From: "Brendan C. Ward" Date: Sat, 21 May 2022 07:16:58 -0700 Subject: [PATCH 3/4] use match for error, try without casts to c_double --- src/spatial_ref/srs.rs | 41 +++++++++++++++++++---------------------- 1 file changed, 19 insertions(+), 22 deletions(-) diff --git a/src/spatial_ref/srs.rs b/src/spatial_ref/srs.rs index 7eca66db7..b7836d4f3 100644 --- a/src/spatial_ref/srs.rs +++ b/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_double, c_int}; +use libc::{c_char, c_int}; use std::ffi::{CStr, CString}; use std::ptr::{self, null_mut}; use std::str::FromStr; @@ -55,37 +55,34 @@ impl CoordTransform { 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, + 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 { - 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); + 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) }; - Err(GdalError::InvalidCoordinateRange { + 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. From 3ed047442608e29da355e5c4bbfd1c223de35edc Mon Sep 17 00:00:00 2001 From: "Brendan C. Ward" Date: Mon, 23 May 2022 09:17:22 -0700 Subject: [PATCH 4/4] Fix formatting issues --- src/spatial_ref/srs.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/spatial_ref/srs.rs b/src/spatial_ref/srs.rs index b7836d4f3..c21f92c94 100644 --- a/src/spatial_ref/srs.rs +++ b/src/spatial_ref/srs.rs @@ -71,15 +71,15 @@ impl CoordTransform { let msg = match _last_cpl_err(CPLErr::CE_Failure) { GdalError::CplError { msg, .. } => match msg.is_empty() { false => Some(msg), - _ => None + _ => None, }, - err => return Err(err) + 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])