Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add ability to convert 3d coordinates #176

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Add ability to convert 3d coordinates #176

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8da156e
Add ability to convert 3d coordinates. Fixes #150
PaulWagener Aug 6, 2023
ee68bf8
Fixes CI
PaulWagener Aug 13, 2023
01ae18e
Merge remote-tracking branch 'upstream/main' into add_z
PaulWagener Dec 19, 2024
8411d39
Update comments
PaulWagener Dec 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 8 additions & 2 deletions src/geo_types.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,10 @@ impl<T: crate::proj::CoordinateType> crate::Coord<T> for geo_types::Coord<T> {
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 }
}
}
Expand All @@ -50,7 +53,10 @@ impl<T: crate::proj::CoordinateType> crate::Coord<T> for geo_types::Point<T> {
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)
}
}
Expand Down
5 changes: 4 additions & 1 deletion src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 }
//! }
//! }
Expand Down
81 changes: 63 additions & 18 deletions src/proj.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,8 @@ fn result_from_create<T>(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)
Expand All @@ -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<T: CoordinateType> Coord<T> for (T, T) {
Expand All @@ -77,11 +78,29 @@ impl<T: CoordinateType> Coord<T> 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<T: CoordinateType> Coord<T> 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]
Expand Down Expand Up @@ -741,33 +760,37 @@ 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 {
proj_errno_reset(self.c_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)?))
Expand Down Expand Up @@ -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
Expand All @@ -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)?))
Expand Down Expand Up @@ -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,
},
})
Expand All @@ -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)?,
)
}
}
Expand Down Expand Up @@ -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 }
}
}
Expand Down Expand Up @@ -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) {
Expand Down