Add ability to convert 3d coordinates

[PaulWagener]

Fixes #150

This PR adds the ability to convert 3d coordinates.

The Coord trait is changed to add a Z coordinate so at least a minor version bump is warranted with this change.

[urschrei]

Hi @PaulWagener, thanks for the PR!

As ever with 2D -> 3D, there's likely to be some discussion…

At minimum, I don't think we can remove the xy functions in one fell swoop as that would break a great deal of downstream code. I think marking them as deprecated in favour of the new ones is better, here.

As for Coord, we're looking at the same problem. In a sense, this is our fault for using an overly general name, but changing the signature will break a lot of code. How do people feel about a new Coord3 (or maybe Coord4 to include M values) struct and deprecating Coord.

Even more radically: we could adopt Rust-Geodesy's coordinate struct…

[PaulWagener]

As long as 3D coordinates will be supported I will be happy :)

[oreilles]

What more is needed for this feature to be added ? I'd be happy to help.

[lx-88]

I just wanted to bring up that this issue hurt me on a recent project. It would be great to see a solution to 3d transformations, or at least allowing transformation into ECEF, into the project and released.

[michaelkirk]

@lx-88 - does using this branch work for you?

[michaelkirk]

How do people feel about a new Coord3 (or maybe Coord4 to include M values) struct and deprecating Coord.

A new struct or a new trait?

I thought the whole "bring your own struct" thing was kind of nice, so users don't have to first convert their data to an intermediate proj-crate format.

Just brainstorming... if avoiding api-breakage is a concern, could we introduce the new 3d method with a default impl?

 trait Coord {
   ...
+  #[deprecated(note = "implement from_xyz instead, this method will be removed in a future release")]
   fn from_xy(x: f64, y: f64, z:f64) -> Self;

   // at some point we'd remove the default impl, which would be a breaking release - presumably in the same release as we remove `from_xy`
+  fn from_xyz(x: f64, y: f64, z:f64) -> Self {
+     Self::from_xy(x, y)
+  }
 }

[nokonoko1203]

Hello!
I would be very happy if this crate could be made to correspond to three-dimensional coordinates!

I know that it is very difficult to do this, but what is the reason that this PR has not been merged?
Is there anything I can do to help?

[balazsdukai]

In the meantime, while this PR is being worked on, here is a fork that reprojects 3d coordinates . Here is an example how I use it.
I did the minimum necessary for my use case, but it might be useful for others.

[nokonoko1203]

This is great! I'll use it as a reference!
I was just thinking of using it for a similar project.
(PointCloud to 3D Tiles v1.1: https://github.com/MIERUNE/point-tiler)

[nokonoko1203]

incorporated the contents of this PR into my fork.
I used it to convert from the projection coordinate system used in Japan, and it works very smoothly!

https://github.com/nokonoko1203/proj/blob/f643a9a197921fc7f31571972f522bbd5149a20e/src/proj.rs#L1330

#[test]
fn test_3d_crs() {
	let from = "EPSG:6677";
	let to = "EPSG:6697";
	let proj = Proj::new_known_crs(from, to, None).unwrap();
	let (x, y, z) = (-5998.998, -35838.918, 3.901);
	let (lng, lat, height) = proj.convert((x, y, z)).unwrap();
	assert_relative_eq!(lng, 139.76706139226548, epsilon = 1e-3);
	assert_relative_eq!(lat, 35.67694831658619, epsilon = 1e-3);
	assert_relative_eq!(height, 3.901, epsilon = 1e-3);
	
	let from = "EPSG:6697";
	let to = "EPSG:4979";
	let proj = Proj::new_known_crs(from, to, None).unwrap();
	let (lng, lat, height) = proj.convert((lng, lat, height)).unwrap();
	assert_relative_eq!(lng, 139.76706139226548, epsilon = 1e-3);
	assert_relative_eq!(lat, 35.67694831658619, epsilon = 1e-3);
	assert_relative_eq!(height, 40.53393140934767, epsilon = 1e-3);
	
	let from = "EPSG:4979";
	let to = "EPSG:4978";
	let proj = Proj::new_known_crs(from, to, None).unwrap();
	let (x, y, z) = proj.convert((lng, lat, height)).unwrap();
	assert_relative_eq!(x, -3959898.9249523925, epsilon = 1e-3);
	assert_relative_eq!(y, 3350278.1607454875, epsilon = 1e-3);
	assert_relative_eq!(z, 3699157.243055472, epsilon = 1e-3);
}

[arpadav]

I just encountered this lack-of-feature, and implemented the same solution, so I am glad this is getting attention.

However, from an API standpoint, sometimes I just want to change my CRS for degrees of an ellipsoid with no 3d component, but sometimes I require 3d for operations like ECEF

So should users be required to add 0 z's for 2d coordinates in their from_xyz implementation, or should there be some sort of CoordXy/Coord2 and CoordXyz/Coord3 traits?

I'm honestly fine with either solution, but am just curious on others thoughts

[susu]

I've stumbled upon this missing functionality, trying to convert from ECEF (EPSG:4978) to EPSG:4326. Tried nokonoko1203's fork which works well (rebased to main). What leftover tasks remaining here to merge this?

[nokonoko1203]

What is the work required to merge this PR?
I'm very happy if this feature is merged!

[busstoptaktik]

(Comming back to look at this after a l-o-o-n-g time)

In 2023, @urschrei said:

Even more radically: we could adopt Rust-Geodesy's coordinate struct…

While much later @michaelkirk added:

I thought the whole "bring your own struct" thing was kind of nice, so users don't have to first convert their data to an intermediate proj-crate format.

I just want to add, that in the meantime, Rust Geodesy has gone the "bring your own struct" way, so implementing the CoordinateSet trait for your bag of data structures implementing CoordinateTuple, brings Rust Geodesy directly into the guts of your own data structures. As said here:

/// `CoordinateTuple` is the ISO-19111 atomic spatial/spatiotemporal
/// referencing element. So loosely speaking, a CoordinateSet is a
/// collection of CoordinateTuples.
///
/// Note that (despite the formal name) the underlying data structure
/// need not be a Rust tuple: It can be any item, for which it makes
///  sense to implement the CoordinateTuple trait.
///
/// The CoordinateTuple trait provides a number of convenience accessors
/// for accessing single coordinate elements or tuples of subsets.
/// These accessors are pragmatically named (x, y, xy, etc.). While these
/// names may be geodetically naïve, they are suggestive, practical, and
/// aligns well with the internal coordinate order convention of most
/// Geodesy operators.
///
/// All accessors have default implementations, except the 3 methods
/// [`nth_unchecked()`](Self::nth_unchecked()),
/// [`set_nth_unchecked()`](Self::set_nth_unchecked) and
/// [`dim()`](Self::dim()),
/// which must be provided by the implementer.
///
/// When accessing dimensions outside of the domain of the CoordinateTuple,
/// [NaN](f64::NAN) will be returned.

Returning NaNs for out-of-domain ensures that computations unintentionally needing access to e.g. height or time, will ensure NaN propagating in cases of data sets not supporting them.

For convenience, Rust Geodesy provides ready-to-eat implementations for these cases:

pub struct Coor2D(pub [f64; 2]);
pub struct Coor3D(pub [f64; 3]);
pub struct Coor4D(pub [f64; 4]);
pub struct Coor32(pub [f32; 2]);

and CoordinateSet implementations for Arrays, slices and Vecs of data structures implementing CoordinateTuple.

It would be nice (but non-trivial, due to integer indexing, and need for length being the total number of coordinates, not of geometries) to implement those for GeometryCollection, for in-place transformation of geometry coordinates.