Skip to content

Commit

Permalink
A CoordTrait experiment (#103)
Browse files Browse the repository at this point in the history
Essentially:

- Traitify coordinate access
- Add dim() as required impl for `CoordinateSet`

In detail:

- A CoordTrait experiment
- Simplified, and added dimension+measure
- Trait impl of ISO-19111 coordinate tuple
- Angular conversions in CoordinateTuple
- Extend and use CoordinateTuple trait. dim in CoordinateSet
- Prefix unchecked to postfix
- xy/set_xy for CoordinateSet. Tests via tmerc/merc
  • Loading branch information
busstoptaktik authored Apr 5, 2024
1 parent 5ed8773 commit 5582bac
Show file tree
Hide file tree
Showing 27 changed files with 704 additions and 612 deletions.
1 change: 1 addition & 0 deletions Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ default-run = "kp"
[dependencies]
# Library functionality
uuid = { version = "1", features = ["v4"] }
num-traits = { version = "0.2" }

# Command line program helpers
clap = { version = "4.3.11", features = ["derive"], optional = true }
Expand Down
7 changes: 3 additions & 4 deletions examples/00-transformations.rs
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,7 @@ fn main() -> anyhow::Result<()> {
// But since a coordinate tuple is really just an array of double
// precision numbers, you may also generate it directly using plain
// Rust syntax. Note that Coor2D, like f64, provides the to_radians
// method. So compared to cph_raw above, we can use a slightly more
// compact notation.
// method, but it operates in place, so we need two steps in this case.
let cph_direct = Coor2D([12., 55.]).to_radians();
// The three versions of Copenhagen coordinates should be identical.
assert_eq!(cph, cph_raw);
Expand Down Expand Up @@ -79,8 +78,8 @@ fn main() -> anyhow::Result<()> {
// Note the use of `to_geo`, which transforms lon/lat in radians
// to lat/lon in degrees. It is defined for Coor2D as well as for
// arrays, vectors and slices of Coor2D
for coord in data.to_geo() {
println!(" {:?}", coord);
for coord in data {
println!(" {:?}", coord.to_geo());
}

// To geo again, but using slices - in two different ways
Expand Down
4 changes: 2 additions & 2 deletions examples/01-geometric_geodesy.rs
Original file line number Diff line number Diff line change
Expand Up @@ -46,8 +46,8 @@ fn main() -> Result<(), Error> {
// surface of the ellipsoid. Let's compute the distance and
// azimuth between the approximate locations of the airports
// of Copenhagen (CPH) and Paris (CDG).
let CPH = Coor4D::geo(55., 12., 0., 0.);
let CDG = Coor4D::geo(49., 2., 0., 0.);
let CPH = Coor2D::geo(55., 12.);
let CDG = Coor2D::geo(49., 2.);

// By historical convention the "from A to B" situation is considered
// the inverse sense of the geodesic problem - hence `geodesic_inv`:
Expand Down
3 changes: 1 addition & 2 deletions examples/02-user_defined_macros.rs
Original file line number Diff line number Diff line change
Expand Up @@ -46,10 +46,9 @@ fn main() -> anyhow::Result<()> {

// data.to_geo() transforms all elements in data from the internal GIS
// format (lon/lat in radians) to lat/lon in degrees.
data.to_geo();
println!("Back to ed50:");
for coord in data {
println!(" {:?}", coord);
println!(" {:?}", coord.to_geo());
}

Ok(())
Expand Down
3 changes: 3 additions & 0 deletions examples/06-user_defined_coordinate_types_and_containers.rs
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,9 @@ impl CoordinateSet for AbscissaCollection {
fn len(&self) -> usize {
4
}
fn dim(&self) -> usize {
1
}
}

fn main() -> Result<(), anyhow::Error> {
Expand Down
3 changes: 1 addition & 2 deletions examples/07-examples_from_ruminations.rs
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,8 @@ fn rumination_000() -> Result<(), anyhow::Error> {

// [6] And go back, i.e. utm -> geo
ctx.apply(utm32, Inv, &mut data)?;
data.to_geo();
for coord in data {
println!("{:?}", coord);
println!("{:?}", coord.to_geo());
}

Ok(())
Expand Down
5 changes: 3 additions & 2 deletions examples/08-user_defined_operators_using_proj.rs
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,7 @@ pub fn proj_constructor(parameters: &RawParameters, _ctx: &dyn Context) -> Resul
fn main() -> anyhow::Result<()> {
let mut prv = geodesy::Minimal::new();
prv.register_op("proj", OpConstructor(proj_constructor));
let e = Ellipsoid::default();

// Check that we can access the `proj` binary - if not, just ignore
if Command::new("proj").stderr(Stdio::piped()).spawn().is_err() {
Expand All @@ -195,7 +196,7 @@ fn main() -> anyhow::Result<()> {
println!("projected: {:?}", geo[0]);

ctx.apply(op, Inv, &mut geo)?;
assert!(rtp[0].default_ellps_dist(&geo[0]) < 1e-5);
assert!(e.distance(&rtp[0], &geo[0]) < 1e-5);
println!("roundtrip: {:?}", geo[0].to_degrees());

// Inverted invocation - note "proj inv ..."
Expand All @@ -208,7 +209,7 @@ fn main() -> anyhow::Result<()> {

// Now, we get the inverse utm projection when calling the operator in the Fwd direction
ctx.apply(op, Fwd, &mut utm)?;
assert!(geo[0].default_ellps_dist(&utm[0]) < 1e-5);
assert!(e.distance(&utm[0], &geo[0]) < 1e-5);
// ...and roundtrip back to utm
ctx.apply(op, Inv, &mut utm)?;
assert!(rtp[0].hypot2(&utm[0]) < 1e-5);
Expand Down
16 changes: 8 additions & 8 deletions src/context/minimal.rs
Original file line number Diff line number Diff line change
Expand Up @@ -143,16 +143,16 @@ mod tests {
assert_eq!(steps[2], "addone inv");

let mut data = some_basic_coor2dinates();
assert_eq!(data[0][0], 55.);
assert_eq!(data[1][0], 59.);
assert_eq!(data[0].x(), 55.);
assert_eq!(data[1].x(), 59.);

assert_eq!(2, ctx.apply(op, Fwd, &mut data)?);
assert_eq!(data[0][0], 56.);
assert_eq!(data[1][0], 60.);
assert_eq!(data[0].x(), 56.);
assert_eq!(data[1].x(), 60.);

ctx.apply(op, Inv, &mut data)?;
assert_eq!(data[0][0], 55.);
assert_eq!(data[1][0], 59.);
assert_eq!(data[0].x(), 55.);
assert_eq!(data[1].x(), 59.);

let params = ctx.params(op, 1)?;
let ellps = params.ellps(0);
Expand All @@ -168,8 +168,8 @@ mod tests {
let op = ctx.op("geo:in | utm zone=32 | neu:out")?;

let mut data = some_basic_coor2dinates();
assert_eq!(data[0][0], 55.);
assert_eq!(data[1][0], 59.);
assert_eq!(data[0].x(), 55.);
assert_eq!(data[1].x(), 59.);

ctx.apply(op, Fwd, &mut data)?;
let expected = [6098907.825005002, 691875.6321396609];
Expand Down
24 changes: 12 additions & 12 deletions src/context/plain.rs
Original file line number Diff line number Diff line change
Expand Up @@ -312,16 +312,16 @@ mod tests {

// ...and it works as expected?
let mut data = some_basic_coor2dinates();
assert_eq!(data[0][0], 55.);
assert_eq!(data[1][0], 59.);
assert_eq!(data[0].x(), 55.);
assert_eq!(data[1].x(), 59.);

ctx.apply(op, Fwd, &mut data)?;
assert_eq!(data[0][0], 56.);
assert_eq!(data[1][0], 60.);
assert_eq!(data[0].x(), 56.);
assert_eq!(data[1].x(), 60.);

ctx.apply(op, Inv, &mut data)?;
assert_eq!(data[0][0], 55.);
assert_eq!(data[1][0], 59.);
assert_eq!(data[0].x(), 55.);
assert_eq!(data[1].x(), 59.);

// Now test that the look-up functionality works in general

Expand All @@ -344,8 +344,8 @@ mod tests {
let mut data = some_basic_coor2dinates();

ctx.apply(op, Fwd, &mut data)?;
assert_eq!(data[0][0], 57.);
assert_eq!(data[1][0], 61.);
assert_eq!(data[0].x(), 57.);
assert_eq!(data[1].x(), 61.);

// 3 Console tests from stupid.md
let op = ctx.op("stupid:bad");
Expand All @@ -354,14 +354,14 @@ mod tests {
let op = ctx.op("stupid:addthree")?;
let mut data = some_basic_coor2dinates();
ctx.apply(op, Fwd, &mut data)?;
assert_eq!(data[0][0], 58.);
assert_eq!(data[1][0], 62.);
assert_eq!(data[0].x(), 58.);
assert_eq!(data[1].x(), 62.);

let op = ctx.op("stupid:addthree_one_by_one")?;
let mut data = some_basic_coor2dinates();
ctx.apply(op, Fwd, &mut data)?;
assert_eq!(data[0][0], 58.);
assert_eq!(data[1][0], 62.);
assert_eq!(data[0].x(), 58.);
assert_eq!(data[1].x(), 62.);

// Make sure we can access "sigil-less runtime defined resources"
ctx.register_resource("foo", "bar");
Expand Down
96 changes: 9 additions & 87 deletions src/coordinate/coor2d.rs
Original file line number Diff line number Diff line change
@@ -1,12 +1,11 @@
use super::*;
use crate::math::angular;
use std::ops::{Index, IndexMut};

/// Generic 2D Coordinate tuple, with no fixed interpretation of the elements
#[derive(Debug, Default, PartialEq, Copy, Clone)]
pub struct Coor2D(pub [f64; 2]);

// ----- O P E R A T O R T R A I T S -------------------------------------------------
use std::ops::{Index, IndexMut};

impl Index<usize> for Coor2D {
type Output = f64;
Expand All @@ -21,34 +20,6 @@ impl IndexMut<usize> for Coor2D {
}
}

// ----- A N G U L A R U N I T S -------------------------------------------

impl AngularUnits for Coor2D {
/// Transform the elements of a `Coor2D` from degrees to radians
#[must_use]
fn to_radians(self) -> Self {
Coor2D([self[0].to_radians(), self[1].to_radians()])
}

/// Transform the elements of a `Coor2D` from radians to degrees
#[must_use]
fn to_degrees(self) -> Self {
Coor2D([self[0].to_degrees(), self[1].to_degrees()])
}

/// Transform the elements of a `Coor2D` from radians to seconds of arc.
#[must_use]
fn to_arcsec(self) -> Self {
Coor2D([self[0].to_degrees() * 3600., self[1].to_degrees() * 3600.])
}

/// Transform the internal lon/lat-in-radians to lat/lon-in-degrees
#[must_use]
fn to_geo(self) -> Self {
Coor2D([self[1].to_degrees(), self[0].to_degrees()])
}
}

// ----- C O N S T R U C T O R S ---------------------------------------------

/// Constructors
Expand Down Expand Up @@ -127,62 +98,13 @@ impl Coor2D {
/// Multiply by a scalar
#[must_use]
pub fn scale(&self, factor: f64) -> Coor2D {
Coor2D([self[0] * factor, self[1] * factor])
Coor2D([self.x() * factor, self.y() * factor])
}

/// Scalar product
#[must_use]
pub fn dot(&self, other: Coor2D) -> f64 {
self[0] * other[0] + self[1] * other[1]
}
}

// ----- D I S T A N C E S ---------------------------------------------------

impl Coor2D {
/// Euclidean distance between two points in the 2D plane.
///
/// Primarily used to compute the distance between two projected points
/// in their projected plane. Typically, this distance will differ from
/// the actual distance in the real world.
///
/// # See also:
///
/// [`distance`](crate::ellipsoid::Ellipsoid::distance)
///
/// # Examples
///
/// ```
/// use geodesy::prelude::*;
/// let t = 1000 as f64;
/// let p0 = Coor2D::origin();
/// let p1 = Coor2D::raw(t, t);
/// assert_eq!(p0.hypot2(&p1), t.hypot(t));
/// ```
#[must_use]
pub fn hypot2(&self, other: &Self) -> f64 {
(self[0] - other[0]).hypot(self[1] - other[1])
}

/// The Geodesic distance on the default ellipsoid. Mostly a shortcut
/// for test authoring
pub fn default_ellps_dist(&self, other: &Self) -> f64 {
Ellipsoid::default().distance(
&Coor4D([self[0], self[1], 0., 0.]),
&Coor4D([other[0], other[1], 0., 0.]),
)
}
}

impl From<Coor2D> for Coor4D {
fn from(c: Coor2D) -> Self {
Coor4D([c[0], c[1], 0.0, 0.0])
}
}

impl From<Coor4D> for Coor2D {
fn from(xyzt: Coor4D) -> Self {
Coor2D([xyzt[0], xyzt[1]])
self.x() * other.x() + self.y() * other.y()
}
}

Expand All @@ -193,28 +115,28 @@ mod tests {
use super::*;
#[test]
fn distances() {
let e = Ellipsoid::default();
let lat = angular::dms_to_dd(55, 30, 36.);
let lon = angular::dms_to_dd(12, 45, 36.);
let dms = Coor2D::geo(lat, lon);
let geo = Coor2D::geo(55.51, 12.76);
assert!(geo.default_ellps_dist(&dms) < 1e-10);
assert!(e.distance(&geo, &dms) < 1e-10);
}

#[test]
fn coor2d() {
let c = Coor2D::raw(12., 55.).to_radians();
let d = Coor2D::gis(12., 55.);
assert_eq!(c, d);
assert_eq!(d[0], 12f64.to_radians());
let e = d.to_degrees();
assert_eq!(e[0], c.to_degrees()[0]);
assert_eq!(d.x(), 12f64.to_radians());
assert_eq!(d.x().to_degrees(), c.x().to_degrees());
}

#[test]
fn array() {
let b = Coor2D::raw(7., 8.);
let c = [b[0], b[1], f64::NAN, f64::NAN];
assert_eq!(b[0], c[0]);
let c = [b.x(), b.y(), f64::NAN, f64::NAN];
assert_eq!(b.x(), c[0]);
}

#[test]
Expand Down
Loading

0 comments on commit 5582bac

Please sign in to comment.