From 8da156e57721994705a9c40d86183c178de4851e Mon Sep 17 00:00:00 2001 From: Paul Wagener Date: Sun, 6 Aug 2023 15:36:24 +0200 Subject: [PATCH 1/3] Add ability to convert 3d coordinates. Fixes #150 --- src/geo_types.rs | 10 +++++-- src/lib.rs | 5 +++- src/proj.rs | 75 ++++++++++++++++++++++++++++++++++++++---------- 3 files changed, 72 insertions(+), 18 deletions(-) diff --git a/src/geo_types.rs b/src/geo_types.rs index 3c7d2ddc..6a42486e 100644 --- a/src/geo_types.rs +++ b/src/geo_types.rs @@ -23,7 +23,10 @@ impl crate::Coord for geo_types::Coord { fn y(&self) -> T { self.y } - fn from_xy(x: T, y: T) -> Self { + fn z(&self) -> T { + T::zero() + } + fn from_xyz(x: T, y: T, _z: T) -> Self { coord! { x: x, y: y } } } @@ -50,7 +53,10 @@ impl crate::Coord for geo_types::Point { fn y(&self) -> T { geo_types::Point::y(*self) } - fn from_xy(x: T, y: T) -> Self { + fn z(&self) -> T { + T::zero() + } + fn from_xyz(x: T, y: T, _z: T) -> Self { Self::new(x, y) } } diff --git a/src/lib.rs b/src/lib.rs index 03f604e6..8830eb50 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -158,7 +158,10 @@ //! fn y(&self) -> f64 { //! self.lat //! } -//! fn from_xy(x: f64, y: f64) -> Self { +//! fn z(&self) -> f64 { +//! 0.0 +//! } +//! fn from_xyz(x: f64, y: f64, _z: f64) -> Self { //! Self { lon: x, lat: y } //! } //! } diff --git a/src/proj.rs b/src/proj.rs index ce37e4a6..455f4f62 100644 --- a/src/proj.rs +++ b/src/proj.rs @@ -67,7 +67,8 @@ where { fn x(&self) -> T; fn y(&self) -> T; - fn from_xy(x: T, y: T) -> Self; + fn z(&self) -> T; + fn from_xyz(x: T, y: T, z: T) -> Self; } impl Coord for (T, T) { @@ -77,11 +78,29 @@ impl Coord for (T, T) { fn y(&self) -> T { self.1 } - fn from_xy(x: T, y: T) -> Self { + fn z(&self) -> T { + T::zero() + } + fn from_xyz(x: T, y: T, _z: T) -> Self { (x, y) } } +impl Coord for (T, T, T) { + fn x(&self) -> T { + self.0 + } + fn y(&self) -> T { + self.1 + } + fn z(&self) -> T { + self.2 + } + fn from_xyz(x: T, y: T, z: T) -> Self { + (x, y, z) + } +} + /// Errors originating in PROJ which can occur during projection and conversion #[derive(Error, Debug)] pub enum ProjError { @@ -733,8 +752,10 @@ impl Proj { }; let c_x: c_double = point.x().to_f64().ok_or(ProjError::FloatConversion)?; let c_y: c_double = point.y().to_f64().ok_or(ProjError::FloatConversion)?; + let c_z: c_double = point.z().to_f64().ok_or(ProjError::FloatConversion)?; let new_x; let new_y; + let new_z; let err; // Input coords are defined in terms of lambda & phi, using the PJ_LP struct. // This signals that we wish to project geodetic coordinates. @@ -744,7 +765,7 @@ impl Proj { let coords = PJ_LPZT { lam: c_x, phi: c_y, - z: 0.0, + z: c_z, t: f64::INFINITY, }; unsafe { @@ -752,14 +773,16 @@ impl Proj { // PJ_DIRECTION_* determines a forward or inverse projection let trans = proj_trans(self.c_proj, inv, PJ_COORD { lpzt: coords }); // output of coordinates uses the PJ_XY struct - new_x = trans.xy.x; - new_y = trans.xy.y; + new_x = trans.xyz.x; + new_y = trans.xyz.y; + new_z = trans.xyz.z; err = proj_errno(self.c_proj); } if err == 0 { - Ok(Coord::from_xy( + Ok(Coord::from_xyz( F::from(new_x).ok_or(ProjError::FloatConversion)?, F::from(new_y).ok_or(ProjError::FloatConversion)?, + F::from(new_z).ok_or(ProjError::FloatConversion)?, )) } else { Err(ProjError::Projection(error_message(err)?)) @@ -807,8 +830,10 @@ impl Proj { { let c_x: c_double = point.x().to_f64().ok_or(ProjError::FloatConversion)?; let c_y: c_double = point.y().to_f64().ok_or(ProjError::FloatConversion)?; + let c_z: c_double = point.z().to_f64().ok_or(ProjError::FloatConversion)?; let new_x; let new_y; + let new_z; let err; // This doesn't seem strictly correct, but if we set PJ_XY or PJ_LP here, the @@ -817,20 +842,22 @@ impl Proj { let xyzt = PJ_XYZT { x: c_x, y: c_y, - z: 0.0, + z: c_z, t: f64::INFINITY, }; unsafe { proj_errno_reset(self.c_proj); let trans = proj_trans(self.c_proj, PJ_DIRECTION_PJ_FWD, PJ_COORD { xyzt }); - new_x = trans.xy.x; - new_y = trans.xy.y; + new_x = trans.xyz.x; + new_y = trans.xyz.y; + new_z = trans.xyz.z; err = proj_errno(self.c_proj); } if err == 0 { - Ok(C::from_xy( + Ok(C::from_xyz( F::from(new_x).ok_or(ProjError::FloatConversion)?, F::from(new_y).ok_or(ProjError::FloatConversion)?, + F::from(new_z).ok_or(ProjError::FloatConversion)?, )) } else { Err(ProjError::Conversion(error_message(err)?)) @@ -1022,11 +1049,12 @@ impl Proj { .map(|point| { let c_x: c_double = point.x().to_f64().ok_or(ProjError::FloatConversion)?; let c_y: c_double = point.y().to_f64().ok_or(ProjError::FloatConversion)?; + let c_z: c_double = point.z().to_f64().ok_or(ProjError::FloatConversion)?; Ok(PJ_COORD { xyzt: PJ_XYZT { x: c_x, y: c_y, - z: 0.0, + z: c_z, t: f64::INFINITY, }, }) @@ -1053,9 +1081,10 @@ impl Proj { // feels a bit clunky, but we're guaranteed that pj and points have the same length unsafe { for (i, coord) in pj.iter().enumerate() { - points[i] = Coord::from_xy( - F::from(coord.xy.x).ok_or(ProjError::FloatConversion)?, - F::from(coord.xy.y).ok_or(ProjError::FloatConversion)?, + points[i] = Coord::from_xyz( + F::from(coord.xyz.x).ok_or(ProjError::FloatConversion)?, + F::from(coord.xyz.y).ok_or(ProjError::FloatConversion)?, + F::from(coord.xyz.z).ok_or(ProjError::FloatConversion)?, ) } } @@ -1179,7 +1208,11 @@ mod test { self.y } - fn from_xy(x: f64, y: f64) -> Self { + fn z(&self) -> f64 { + 0.0 + } + + fn from_xyz(x: f64, y: f64, _z: f64) -> Self { MyPoint { x, y } } } @@ -1285,6 +1318,18 @@ mod test { assert_relative_eq!(t.y(), 1141263.0111604782); } + #[test] + fn test_3d_crs() { + let from = "EPSG:4978"; + let to = "EPSG:7415"; + let proj = Proj::new_known_crs(from, to, None).unwrap(); + let t = proj.convert((3968419.4, 379068.0, 4962142.2)).unwrap(); + + assert_relative_eq!(t.x(), 159814.79923647654); + assert_relative_eq!(t.y(), 380108.089631567); + assert_relative_eq!(t.z(), -27.239349158361257); + } + #[test] fn test_from_crs_nul_error() { match Proj::new_known_crs("\0", "EPSG:4326", None) { From ee68bf828360fa64102bf0c0bbc67cef38e80669 Mon Sep 17 00:00:00 2001 From: Paul Wagener Date: Sun, 13 Aug 2023 09:44:26 +0200 Subject: [PATCH 2/3] Fixes CI --- src/proj.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/proj.rs b/src/proj.rs index 455f4f62..f0afc644 100644 --- a/src/proj.rs +++ b/src/proj.rs @@ -1325,9 +1325,9 @@ mod test { let proj = Proj::new_known_crs(from, to, None).unwrap(); let t = proj.convert((3968419.4, 379068.0, 4962142.2)).unwrap(); - assert_relative_eq!(t.x(), 159814.79923647654); - assert_relative_eq!(t.y(), 380108.089631567); - assert_relative_eq!(t.z(), -27.239349158361257); + assert_relative_eq!(t.x(), 159783.64355817687, epsilon = 1e-5); + assert_relative_eq!(t.y(), 380007.7906745551, epsilon = 1e-5); + assert_relative_eq!(t.z(), 16.769676147028804, epsilon = 1e-5); } #[test] From 8411d39160739087c8d1e2c46f051e9f56674522 Mon Sep 17 00:00:00 2001 From: Paul Wagener Date: Thu, 19 Dec 2024 23:46:39 +0100 Subject: [PATCH 3/3] Update comments --- src/proj.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/proj.rs b/src/proj.rs index bff55a04..6fbeed73 100644 --- a/src/proj.rs +++ b/src/proj.rs @@ -56,8 +56,8 @@ fn result_from_create(context: *mut PJ_CONTEXT, ptr: *mut T) -> Result<*mut T /// A point in two dimensional space. The primary unit of input/output for proj. /// -/// By default, any numeric `(x, y)` tuple implements `Coord`, but you can conform your type to -/// `Coord` to pass it directly into proj. +/// By default, any numeric `(x, y)` or `(x, y, z)` tuple implements `Coord`, but you can conform +/// your type to `Coord` to pass it directly into proj. /// /// See the [`geo-types` feature](#feature-flags) for interop with the [`geo-types` /// crate](https://docs.rs/crate/geo-types) @@ -769,7 +769,7 @@ impl Proj { // This signals that we wish to project geodetic coordinates. // For conversion (i.e. between projected coordinates) you should use // PJ_XY {x: , y: } - // We also initialize z and t in case libproj tries to read them. + // We also initialize t in case libproj tries to read it. let coords = PJ_LPZT { lam: c_x, phi: c_y,