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 bf852fbb..905e8ed1 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 a24e64b3..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) @@ -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)] #[non_exhaustive] @@ -741,18 +760,20 @@ 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. // 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, - z: 0.0, + z: c_z, t: f64::INFINITY, }; unsafe { @@ -760,14 +781,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)?)) @@ -815,8 +838,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 @@ -825,20 +850,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)?)) @@ -1030,11 +1057,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, }, }) @@ -1061,9 +1089,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)?, ) } } @@ -1187,7 +1216,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 } } } @@ -1293,6 +1326,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(), 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] fn test_from_crs_nul_error() { match Proj::new_known_crs("\0", "EPSG:4326", None) {