busstoptaktik/geodesy

Rust geodesy

Geodesy

Abstract

Rust Geodesy is - unsurprisingly - a geodesy library written in the Rust programming language.

Rust Geodesy provides a number of features to support a number of objectives.

The most important features are

• a set of more than 20 geodetic transformation primitives
• a set of more than 40 primitives for operations on the ellipsoid
• a means for composing these primitives into more complex operations.

The most important objectives are

• to support new, and hopefully better, abstractions,
• to use these abstractions to build better, simpler, and more tractable, geospatial standards, transformations, and software.

If any of this resonates with you, read on after this minimal usage example...

Usage

Initialize a new project, using Geodesy:

$ cargo new foo
     Created binary (application) `foo` package

$ cd foo
$ cargo add geodesy

Then copy this to the foo/src/main.rs file: A minimal example, computing the UTM coordinates of some Scandinavian capitals

use geodesy::prelude::*;

fn main() -> Result<(), Box<Error>> {
    let mut context = Minimal::new();
    let utm33 = context.op("utm zone=33")?;

    let cph = Coor2D::geo(55., 12.); // Copenhagen
    let sth = Coor2D::geo(59., 18.); // Stockholm
    let mut data = [cph, sth];

    context.apply(utm33, Fwd, &mut data)?;
    println!("{:?}", data);
    Ok(())
}

and try it out:

$ cargo r
    Finished dev [unoptimized + debuginfo] target(s) in 0.11s
     Running `C:\FLOW\AD\RG\foo\target\debug\foo.exe`

[Coor2D([308124.36786033923, 6098907.825005002]), Coor2D([672319.9640879404, 6543920.334127973])]

Concrete

Rust Geodesy (RG), is a platform for experiments with geodetic software, transformations, and standards. RG vaguely resembles the PROJ transformation system, and was built in part on the basis of experiments with alternative data flow models for PROJ. The fundamental transformation functionality of RG is fairly complete (i.e. on par with the datum shift/reference frame transformation capability of PROJ), while the number of projections supported is a far cry from PROJ's enormous gamut. It does, however, support a suite of the most important ones:

• Transverse Mercator
• Universal Transverse Mercator (UTM)
• Web Mercator
• Mercator
• Oblique Mercator
• Lambert Conformal Conic
• Lambert Azimuthal Equal Area

But fundamentally, RG is born as a geodesy, rather than cartography library. And while PROJ benefits from four decades of reality hardening, RG, being a platform for experiments, does not have operational robustness as a main focus. Hence, viewing RG as another PROJ, or PROJ RiiR, will lead to bad disappointment. At best, you may catch a weak mirage of a potential shape of jazz to come for the PROJ internal dataflow.

That said, being written in Rust, with all the memory safety guarantees Rust provides, RG by design avoids a number of pitfalls that are explicitly worked around in the PROJ code base. So the miniscule size of RG compared to PROJ is not just a matter of functional pruning. It is also a matter of development using a tool wonderfully suited for the task at hand.

Also, having the advantage of learning from PROJ experience, both from a user's and a developer's perspective, RG is designed to be significantly more extensible than PROJ. So perhaps for a number of applications, and despite its limitations, RG may be sufficient, and perhaps even useful.

Aims

Dataflow experimentation is just one aspect of RG. Overall, the aims are (at least) fourfold:

1. Support experiments for evolution of geodetic standards.
2. Support development of geodetic transformations.
3. Hence, provide easy access to a number of basic geodetic operations, not limited to coordinate operations.
4. Support experiments with data flow and alternative abstractions. Mostly as a tool for aims (1, 2, 3)

All four aims are guided by a wish to amend explicitly identified shortcomings in the existing geodetic system landscape.

Documentation

The documentation is currently limited, but take a look at:

• The coordinate operator documentation
• The description of kp, the Rust Geodesy coordinate processing program
• This essayistic rumination, outlining the overall philosophy and architecture of Rust Geodesy, and this related comparison between PROJ and RG
• The API documentation at Docs.rs
• The examples
• The tests embedded in the source code
• This rather concrete and this more philosophical description of the main discrepancy between geodesy and geomatics, RG tries to elucidate and amend.

License

Rust Geodesy: Copyright 2020, 2021, 2022, 2023, 2024 by Thomas Knudsen [email protected] and contributors.

Licensed under either of

• Apache License, Version 2.0 (LICENSE-APACHE or here)
• MIT license (LICENSE-MIT or here)

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

geodesy's Issues

Geodesic center lines

See: A MANUAL ON
TECHNICAL ASPECTS OF THE
UNITED NATIONS CONVENTION ON
THE LAW OF THE SEA, Special Publication No. 51
4th Edition - March 2006
Published by the
International Hydrographic Bureau
MONACO

Documentation for the `deformation` operator

The documentation for the deformation operator misses entirely from Rumination 002

gridshift operator

Using Gas

Support some physical geodesy and generic geophysical operations

Explore some roads toward a more smooth handling of parametric and gravity related height systems in ISO 19111/OGC Topic 2, by providing access to some geophysical context.

As a start:

• Normal gravity models, for conversion between gravity and gravity anomalies
  • https://encyclopedia.pub/entry/33797
• Conversion between the "tide free", "mean tide", and "zero tide" systems
  • ICEsat geodesy guide
• Normal atmosphere and ocean models for approximate parametric/metric conversions
  • https://en.wikipedia.org/wiki/Atmospheric_pressure
• Deflections of the vertical (ξ, η) from the astronomical (Φ, Λ) and geodetic (φ, λ) latitude and longitude https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5017334/ the local gradient of a geoid model

Move the `proj` operator to examples

The proj operator is somewhat clunky and should not be part of the library. Move it to examples, so anyone who may need the functionality may see how to obtain it and run-time register it for transparent interoperation with the built-in operators.

Reduce parallel processing boilerplate

Model parallel processing after this example from the Rayon docs

// https://docs.rs/rayon/1.1.0/rayon/slice/trait.ParallelSliceMut.html
use rayon::prelude::*;
let mut array = [1, 2, 3, 4, 5];
array.par_chunks_mut(2)
     .for_each(|slice| slice.reverse());
assert_eq!(array, [2, 1, 4, 3, 5]);

Support oblique mercator

(As per request from @pka, on GeoRust Discord, for use in tile-grid):

@pka:

Non-mercator transformations are based on Proj. But would be great, if we could cover selected projections with Geodesy!

@busstoptaktik:

That would be awesome, although the selection of projections is very limited compared to PROJ's huge gamut

@pka:

I'm mainly interested in ETRS89 and https://spatialreference.org/ref/epsg/2056/ - would these be possible?

@busstoptaktik:

epsg:4258 to epsg:2056 is a two step process (cf. projinfo -s epsg:4258 -t epsg:2056), involving first a datum shift from ETRS89 to CH1903+, then a projection to "Hotine oblique mercator variant B".

The datum shift is registered in EPSG as a 3 parameter Helmert shift, which is supported by Geodesy, while the second step, the projection to Oblique Mercator is not (yet).

It is, however, documented in "EPSG Guidance Note 7 part 2", page 64ff, and should not be too hard to implement in Rust - I would happily accept a pull request 🙂

Support NMEA pseudodecimal representation of DMS as used in navigation

The NMEA protocol uses the (D)DDMM.mmm format, which is quite useful for use in situations like:

echo 5542.25 1234.50 | kp "NMEA | whatever | geo inv"

A corresponding DMS format (D)DDMMSS.sss is also useful

echo 554215.1234 123430.5678 | kp "DMS | whatever | geo inv"

Architecture and philosophy documentation

Provide a brief overview of the architecture and the thoughts behind

Better documentation of the context providers

Minimal and Plain both need documentation. Perhaps taking off from some documentation of the Context trait itself.

Note that the Maximal example in the tests-folder comprises a fully self contained demonstration of how to implement a user defined context provider.

Consider switching to minimal-yaml

cf. https://docs.rs/crate/minimal-yaml/0.1.5

Rename nmea Op

NMEA is a bad name for something that really is ISO 6709.

Should RG accept ellipse params as well as an ellps=<name>?

While testing various output from projinfo I've come across an example where the ellipse is described in terms of rf and a parameters instead of ellps=<name>. For example a transformation between EPSG:2001 and EPSG:3857 looks like:

+proj=pipeline
  +step +inv +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995 +x_0=400000 +y_0=0 +a=6378249.145 +rf=293.465
  +step +inv +proj=longlat +a=6378249.145 +rf=293.465 # Note: I've implemented longlat as a noop in geodesy-wasm.
  +step +proj=push +v_3
  +step +proj=cart +a=6378249.145 +rf=293.465
  +step +proj=helmert +x=-255 +y=-15 +z=71
  +step +inv +proj=cart +ellps=WGS84
  +step +proj=pop +v_3
  +step +proj=webmerc +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84

# With the equivalent RG pipeline being:

  | tmerc inv lat_0=0 lon_0=-62 k=0.9995 x_0=400000 y_0=0  a=6378249.145 rf=293.465
  | noop inv a=6378249.145 rf=293.465
  | push v_3
  | cart a=6378249.145 rf=293.465
  | helmert x=-255 y=-15 z=71
  | cart inv ellps=WGS84
  | pop v_3
  | webmerc lat_0=0 lon_0=0 x_0=0 y_0=0 ellps=WGS84

When converting the same coordinate with PROJ and RG I get differences in the order of 10's of meters. The a=6378249.145 rf=293.465 parameters describe the clrk80 ellipse, so if I update both pipelines to use ellps=clrk80 instead of using a and rf the results are the same.

The question is, should RG be able to handle ellipse input via the rf and a parameters? If yes, I'm happy to attempt a PR with some guidance on where to do this and what else needs to be covered.

Additional `kp` functionality and flags

With the reincarnation of the Plain context provider and its support for handling ancillary data, kp has become rather capable recently, but there are still some obvious holes in the functionality.

Obvious because in an abstract geodetic tool landscape, we may consider kp to sit as the third vertex in a triangle, where the two other vertices comprise the PROJ command line tools proj and cct. Or perhaps even the third vertex in a tetrahedron, with the PROJ geod program sitting in the fourth out-of-plane vertex. So we should expect kp to cover more or less everything these two (three) do. And more or less none of what their cousin cs2cs does, since that is entirely out of scope for Rust Geodesy.

Hence kp should gain at least some of these items:

• Handle binary input and output: Conversion from text to floating-point binary is a very slow proces, and as kp is intended a.o. as a filter in a data processing pipeline, feeding binary stuff from and to neighbours in the pipeline should speed things somewhat up
• Have better facilities for selecting 2D and 3D input and output
• Provision for including the input record (either in its entirety or through a tag column) into the output record for information and/or for synchronisation
• Prolong the road already taken with the implementation of the geodesic and the curvature operators: More access to the ellipsoidal primitives!
• Show scale factors (i.e. the proj's -S option, and the proj_factors entry of the PROJ 4D-API.
• ...and probably a bunch more

Grids...

gridshift/vgridshift

Consider "Keep a Changelog"

The CHANGES.md file is a mostly-manual regurgitation of elements from the git log. Keep a changelog proposes an opinionated recipe for optimizing the end-user value of this manual process. Kind of like cargo fmt for changelog entries

Consider adding a feature for logging

Currently when coordinate transformation fails log warnings are written for every such point. While is is very useful for console applications and debugging, it can be quite annoying for different kinds of applications. For example, I use your library to render maps in different projections. And then I allow the user to query the objects by the mouse pointer. And when the user moves the mouse outside the projections scope, I get 1 million log warnings:

[2024-01-16T19:25:43Z WARN  geodesy::inner_op::laea] LAEA: (17613175.607089374, 7290876.352930736) outside domain
[2024-01-16T19:25:43Z WARN  geodesy::inner_op::laea] LAEA: (17585327.302897472, 7250347.963004524) outside domain
[2024-01-16T19:25:43Z WARN  geodesy::inner_op::laea] LAEA: (17582293.47894309, 7228565.651940692) outside domain
...

It would be nice if I could opt out from logging, since I handle all transformation failures by my code.

Molodensky with correction

A method for roundtrip-correction of the Molodensky transformation was described in:

A. C. Ruffhead, 2016,
The SMITSWAM method of datum transformations consisting of Standard Molodensky in two stages with applied misclosures,
Survey Review, 48:350, 376-384,
DOI: 10.1080/00396265.2016.1191748

It is simple to implement, and could be fun to try out

Remove file I/O from the Minimal-context provider

... for better WASM compatibility. Perhaps base another context (Web?) on Reqwest or ureq

Support @null and @optional grid inputs

Based on discussion here: #65 (comment)

CoordinateTuple should be a trait, rather than a type

Problem

Currently, Rust Geodesy forces input data to be in the form of a slice of CoordinateTuples, which really is forcing the coordinate representation conventions of the library on to the application code calling it.

In PROJ, I eliminated this unfortunate characteristic in 2017, by introducing the proj_trans_generic(...) function in the 4D API (cf OSGeo/PROJ#530).

With Rust having much better support for generic programming than C, one would think that it would be possible to do something better in the case of Rust Geodesy. I believe the demo code below shows that this is actually the case.

Solution

Essentially, we do four things:

1. Introduce the internal coordinate representation Coord, formerly known as CoordinateTuple:

   pub struct Coord(pub [f64; 4]);

2. Introduce the minimal CoordinateTuple trait, describing anything that can be both constructed from a Coord, turned into a Coord, and copied:

   pub trait CoordinateTuple: Copy + From<Coord> + Into<Coord> {}
   impl<T: Copy + From<Coord> + Into<Coord>> CoordinateTuple for T {}

3. Introduce the CoordinateSet trait, describing a set of any kind of CoordinateTuples. Here CoordinateSet refers to the ISO19100 concept of coordinatesets, i.e. a collection of CoordinateTuples, all of the same base type.

   pub trait CoordinateSet {
       fn len(&self) -> usize;
       fn into_coord(&self, index: usize) -> Coord;
       fn from_coord(&mut self, index: usize, coord: Coord);
   }

4. Switch the internals to work on generic CoordinateSets, rather than slices of CoordinateTuples

Conclusion

With this done, users of Rust geodesy may implement the two tiny traits above, for their particular data type, and RG will automagically support that type.

The code below includes two examples:

1. A generic plain coordinate set type, demoed as container for a 2-D single precision coordinate type.
2. A generic "Point Feature" coordinate set type, combining coordinates and attributes.

I believe the demo suggests that this will make Rust Geodesy much more generally useful.

Demo

/// Demo implementation of `CoordinateTuple` as trait. Makes it possible to
/// support transformations *directly on generic user-defined data types*,
/// as in the PROJ function `proj_trans_generic(...)`.
///
/// Compared to `proj_trans_generic` (introduced by yours truly on 2017-07-07
/// in https://github.com/OSGeo/PROJ/pull/530), `CoordinateTuple` as trait
/// is cleaner, clearer, and much safer.
///
/// Thomas Knudsen, 2022-01-02

// ----------------------------------------------------------------------------

/// `Coord` is the generic 4-D, RG-internal, double precision, coordinate
/// representation formerly known as `CoordinateTuple`
#[derive(Clone, Copy, Debug, Default)]
pub struct Coord(pub [f64; 4]);

/// `coord[0]` is a bit more readable than `coord.0[0]`, so we implement
/// the indexing traits for `Coord`
impl std::ops::Index<usize> for Coord {
    type Output = f64;
    fn index(&self, i: usize) -> &Self::Output {
        &self.0[i]
    }
}

impl std::ops::IndexMut<usize> for Coord {
    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
        &mut self.0[i]
    }
}

// ----------------------------------------------------------------------------

/// A `CoordinateTuple` in the new trait sense, is anything that can be both
/// 1: copied, 2: constructed from a `Coord`, and 3: turned into a `Coord`
pub trait CoordinateTuple: Copy + From<Coord> + Into<Coord> {}
impl<T: Copy + From<Coord> + Into<Coord>> CoordinateTuple for T {}

// ----------------------------------------------------------------------------

/// `Xy` is an example of a user-defined coordinate representation.
/// `Xy` implements the CoordinateTuple trait by deriving `Copy`, and
/// providing `Into<Coord>` and `From<Coord>`.
/// `Xy` is 2-D, single precision, as may make sense for e.g. small
/// scale map data, where it will provide ample resolution, while
/// weighing only 1/4 of a `Coord`
#[derive(Clone, Copy, Debug, Default)]
pub struct Xy(pub [f32; 2]);

impl Into<Coord> for Xy {
    fn into(self) -> Coord {
        Coord([self.0[0] as f64, self.0[1] as f64, 0., f64::NAN])
    }
}

impl From<Coord> for Xy {
    fn from(coord: Coord) -> Xy {
        Xy([coord.0[0] as f32, coord.0[1] as f32])
    }
}

// ----------------------------------------------------------------------------

/// In the ISO-19100 standards series sense, a `CoordinateSet` is a set of
/// `CoordinateTuple`s, all of the same type.
/// The trait implementation is surprisingly simple: We just need to be able
/// to iterate over the contents, to extract each element as a `Coord`, and
/// to re-introduce the contents of the transformed `Coord` into the original
/// data element
pub trait CoordinateSet {
    fn len(&self) -> usize;
    fn into_coord(&self, index: usize) -> Coord;
    fn from_coord(&mut self, index: usize, coord: Coord);
}


// ----------------------------------------------------------------------------

/// In the simplest case, our CoordinateSet consists of a slice of raw
/// `CoordinateTuple`s
#[derive(Debug, Default)]
struct PlainCoordinateSet<'a, C: CoordinateTuple> (
    &'a mut [C]
);

/// Since the data type `T` implements the `CoordinateTuple` trait,
/// it is extremely simple to implement the `CoordinateSet` trait,
/// as we already have the necessary `From<Coord>` and `Into<Coord>`
/// implementations.
///
/// Note that due to the from/into duality, `From<Coord>` is not
/// explicitly used: The compiler auto-derives `Into<T> for Coord`
/// when it sees `From<Coord> for T`
impl<'a, T: CoordinateTuple> CoordinateSet for PlainCoordinateSet<'a, T> {
    fn len(&self) -> usize {
        self.0.len()
    }

    fn into_coord(&self, index: usize) -> Coord {
        self.0[index].into()
    }

    fn from_coord(&mut self, index: usize, coord: Coord) {
        self.0[index] = coord.into()
    }
}

// ----------------------------------------------------------------------------

/// This is a more realistic case of implementing `CoordinateSet` for a rather
/// generic class of potential user defined types, composed from a payload and
/// a coordinate tuple.

/// The element `PointFeature`
pub struct PointFeature<T: CoordinateTuple, P> {
    pub coordinatetuple: T,
    pub payload: P
}

/// ...a set of `PointFeature`s
struct PointFeatureSet<'a, T: CoordinateTuple, P> (
    &'a mut [PointFeature<T, P>]
);

/// And the implementation of the `CoordinateSet` trait for generic `PointFeatureSet`s
impl<'a, T: CoordinateTuple, P> CoordinateSet for PointFeatureSet<'a, T, P> {
    fn len(&self) -> usize {
        self.0.len()
    }

    fn into_coord(&self, index: usize) -> Coord {
        self.0[index].coordinatetuple.into()
    }

    fn from_coord(&mut self, index: usize, coord: Coord) {
        self.0[index].coordinatetuple = coord.into()
    }
}

// ----------------------------------------------------------------------------

/// Code taking a generic `CoordinateTuple` as arg
fn print_coordinatetuple<T: CoordinateTuple>(ct: T) {
    let coord: Coord = ct.into();
    println!("{:#?}", coord);
}

/// Dummy transformation, here using a vtable
fn transformation_of_a_coordinate_set(cs: &mut dyn CoordinateSet) {
    for i in 0..cs.len() {
        let mut co = cs.into_coord(i);
        co[0] += 2.;
        cs.from_coord(i, co);
    }
}

fn main() {
    // A vector of `CoordinateTuple`s
    let mut xxxx = vec![Xy([1.,2.]), Xy([3.,4.])];
    print_coordinatetuple(xxxx[1]);

    // Now turn the vector into a `PlainCoordinateSet` (or perhaps any other
    // type, T, implementing the `CoordinateSet` trait), in order to be able
    // to transform it.
    //
    // Note that the life time of the `PlainCoordinateSet` struct need not be
    // any longer than the run time of the transformation, so we generate it
    // inline with the function call.
    transformation_of_a_coordinate_set(&mut PlainCoordinateSet(&mut xxxx));
    dbg!(&xxxx);

    // But of course we may also generate it prior to the call
    let mut xyxy = PlainCoordinateSet(&mut xxxx);
    transformation_of_a_coordinate_set(&mut xyxy);
    dbg!(&xyxy);
}

Better documentation of `kp`

The documentation of kp in Rumination 003 lags heavily - kp is way more useful now than at time of documentation

LAEA projection does not preserve latitude on inverse transformation

This fails with longitude being equal before and after transformations, but latitude == 0.

#[test]
fn lambert_projection() {
    let mut ctx = Minimal::new();
    let op = ctx.op("laea lon_0=10 lat_0=52 x_0=4321000 y_0=3210000").unwrap();
    let mut data = [Coor2D::geo(52.0, 10.0)];
    let clone = data.clone();
    ctx.apply(op, Fwd, &mut data).unwrap();
    ctx.apply(op, Inv, &mut data).unwrap();
    assert_eq!(data, clone);
}

Consider alternative parameter names for the Helmert operator

The Helmert operator can take up to 16 parameters in the dynamic case. This makes expressions involving Helmert steps quite hard to read. And as the parameters typically come in families of 3, we might as well introduce an optional set of parameter names, each taking a group of 3 numbers.

For instance, we may replace the group

x=0.03054 y=0.04606 z=-0.07944
rx=0.00141958 ry=0.00015132 rz=0.00150337
s=0.003002

by the more easily eaten

translation = 0.03054, 0.04606, -0.07944
rotation = 0.00141958, 0.00015132, 0.00150337
scale = 0.003002

perhaps adding homologous names like velocity and angular_velocity for the dx, dy, dz and drx, dry, drz pairs.

Or what about translation_trend and rotation_trend?

For PROJ compatibility, the existing parameters should be kept as alternatives.

Tasks

• Implement functionality (@nrhill1 in #91)
• Update tests (@nrhill1 in #91)
• Update documentation

Support translation from PROJ strings to Geodesy format

Over at geodesy-wasm, @Rennzie identifies the need for a

Proj string to geodesy_rs parser so that proj strings can be used with geodesy_wasm

Given that the Geodesy crate contains all necessary material for doing this, we really should support it directly in the library. Whether it should be directly supported as input format is less evident.

Documentation for KP

Implementing Swiss Oblique Mercator

For some reason I've not been able to tag you successfully in my PR over at geodesy-wasm. I've implemented the somerc operator and wondered if you'd rather I PR it into RG?

Rennzie/geodesy-wasm#27

Constructing `OpHandle` for EPSG:4326 results in error

EPSG:4326 from crs-definitions:

Def { code: 4326, proj4: "+proj=longlat +datum=WGS84 +no_defs", wkt: "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]" }

The relevant part being the proj4 string: "+proj=longlat +datum=WGS84 +no_defs".

When using parse_proj on that string, this returns "longlat datum=WGS84 no_defs".

Attempting to construct an OpHandle from that via geodesy::Plain::new().op(...).unwrap() results in this error:

Result::unwrap()` on an `Err` value: NotFound("longlat", ": longlat datum=WGS84 no_defs")

UTM: support "south" aspect

y_0 = 10000000 yada yada...

Documentation for operators

Operators are below the API surface, hence do not (and should not) get handled by the docs built for docs.rs

In due course, they need their own book, but for now, it suffices with a RG Comment

Parallel

utilize the minions in the context-definition

Implement improved meridian distance series from Charles Karney

In PROJ PR 3516, @cffk demonstrated new series expansions for the meridian distance. Rust geodesy currently uses an older implementation by B.R.Bowring. As we have already implemented most of the machinery needed for Karney's work in constants.rs, and in latitudes.rs, we really should implement Karney's work here as well

`CoordinateSet` trait: `get_coord/set_coord` should be `get/set`

According to the Rust API guidelines, when only one thing could be reasonably gotten/sat, the corresponding methods should be named get and set respectively.

Hence the CoordinateSet methods get_coord(i) and set_coord(i, value) should be deprecated, replaced by get(i) and set(i, value), and eventually removed.

Imperfect handling of flag parameters for macros

Demonstrated in the test macro_expansion_with_defaults_provided(): Providing inv=true works,

geodesy/src/op/mod.rs

Line 359 in 018fc20

let op = prv.op("addone|helmert:one x=2 inv=true")?;

while just plain inv doesn't

Improve error messaging

Use thiserror/anyhow

Clean up YAML-handling

Likely, the pre-processing into PROJ-style tokens is superfluous (it was probably the result of early struggling with the crude YAML-Rust documentation, while still learning the basics of Rust)

External definitions in zipfile

a zipfile rather than a tree of read only pipeline-definitions?
cf.
https://www.reddit.com/r/rust/comments/d5i74n/accessing_zipped_file/
https://doc.rust-lang.org/std/str/fn.from_utf8.html

Bug: tmerc operation returning incorrect results

I've been testing things with a EPSG:27700 (BNG) to EPSG:3857 (webmerc) pipeline and the round trip result for the tmerc operation is incorrect.

# Full pipeline roundtrip
echo 544748. 258372. | cargo run --features binary --bin kp -- "enu:in | tmerc inv lat_0=49 lon_0=-2 k_0=0.9996012717 x_0=400000 y_0=-100000 ellps=airy | webmerc lat_0=0 lon_0=0 x_0=0 y_0=0 ellps=WGS84 | enu:out" --roundtrip

-190686. 5461518. 0.. NaN

# tmerc only roundtrip
echo 544748. 258372.| cargo run --features binary --bin kp -- "enu:in | tmerc inv lat_0=49 lon_0=-2 k_0=0.9996012717 x_0=400000 y_0=-100000 ellps=airy | gis:out" --roundtrip

-190686. 5461518. 0. NaN

# tmerc Cart -> Geo
echo 544748. 258372. | cargo run --features binary --bin kp -- "enu:in | tmerc inv lat_0=49 lon_0=-2 k_0=0.9996012717 x_0=400000 y_0=-100000 ellps=airy | gis:out

0. 52. 0. NaN  # What I expect

# tmerc Geo -> Cart - Incorrect result
echo 0.118 52. | cargo run --features binary --bin kp -- "gis:in | tmerc  lat_0=49 lon_0=-2 k_0=0.9996012717 x_0=400000 y_0=-100000 ellps=airy | enu:out"

354062. 5719890. 0. NaN # NOT what I expect. Should be 544748. 258372.

I've tested this with kp and with geodesy-wasm which implements the Minimal context internally.

Possibly useful pointers:

• The following pipeline produces the correct result in the forward direction but incorrect in the inverse (same control coordinate as above):
  • tmerc inv lat_0=49 lon_0=-2 k_0=0.9996012717 x_0=400000 y_0=-100000 ellps=airy | webmerc lat_0=0 lon_0=0 x_0=0 y_0=0 ellps=WGS84
• Evaluating the individual parts:
  • tmerc and webmerc correctly transform the equivalent coordinate to WGS84[gis] (0.118 52.).
  • tmerc does NOT transform that WGS84 coordinate to cartesian coords correctly while webmerc does.

The last point I think suggests the tmerc FWD operation is producing incorrect results?

---

Tested on latest main as of 18/08/2023.
I've truncated all results to integers for readability

Alternative interpolation algorithms

For grid interpolation, this article and its code companion shows Constrained Bicubic, which is "bicubic without overshoots"

// Constrained bicubic interpolation
template<typename T>
void cbi (const Matrix<T>& in, Matrix<T>& out)
{
  //scaling factors
  double xr = (double)out.cols () / (in.cols () - 1);
  double yr = (double)out.rows () / (in.rows () - 1);

  assert (xr >= 1 && yr >= 1);

  //weighting function
  auto w = [](double x, double y) -> double {
    return x * x * y * y * (9 - 6 * x - 6 * y + 4 * x * y);
  };
  
  for (int i = 0; i < out.rows (); i++)
  {
    int ii = (int)floor (i / xr);
    for (int j = 0; j < out.cols (); j++)
    {
      int jj = (int)floor (j / yr);
      T v00 = in (ii, jj), v01 = in (ii, jj + 1),
        v10 = in (ii + 1, jj), v11 = in (ii + 1, jj + 1);
      double fi = i / xr - ii, fj = j / yr - jj;
      out (i, j) = w(1 - fi, 1 - fj) * v00 + w(1 - fi, fj) * v01 +
        w(fi, 1 - fj) * v10 + w(fi, fj) * v11;
    }
  }
}

Improve test support

A CoordinateTuple::dms(latd, latm, lats, lond, lonm, lons, h, t) would vastly simplify harvesting of test coordinates from published literature.

A point1.rough_distance_2D(point2), giving spherical distance between geographical coordinates on a 6400 km sphere, and a point1.rough_distance_3D(point2) giving "hypot3D supposing the ellipsoid is a 6400 km sphere" would be very useful in tests, providing upper bounds for deviations between computed and expected coordinates.

Flesh out kp

Compact cmd line notation, actual functionality, help text

Support Jacobian/Factors in kp

With #54 merged (5eab7f2), it is now possible to support geometrical factors in kp.

It should work much like the -S and -SV options in proj:

$ echo 12 55 | proj -S +proj=utm +zone=32
691875.63       6098907.83      <1.00005 1.00005 1.0001 1.20736e-06 1.00005 1.00005>

$ echo 12 55 | proj -SV +proj=utm +zone=32
Longitude: 12dE [ 12 ]
Latitude:  55dN [ 55 ]
Easting (x):   691875.63
Northing (y):  6098907.83
Meridian scale (h) : 1.00005168  ( 0.005168 % error )
Parallel scale (k) : 1.00005168  ( 0.005168 % error )
Areal scale (s):     1.00010336  ( 0.01034 % error )
Angular distortion (w): 0.000
Meridian/Parallel angle: 90.00000
Convergence : 2d27'29.52" [ 2.45819987 ]
Max-min (Tissot axis a-b) scale error: 1.00005 1.00005

(also see #46)

context api

clean up,

parallel + vector operations

Intersection of geodesics

Consider algorithms

Make better use of `split_into_steps()`

Move split_into_steps() from inner_op/pipeline.rs to the Ancillary functions section in op/mod.rs, and use it as primary building block for the other functions in that section.

Support Geodetic TIFF grids

As a question. Should we consider supporting TIFF grids alongside Gravsoft and Ntv2?

• https://proj.org/en/9.3/specifications/geodetictiffgrids.html

Logging

Use Rust's excellent logging facilities