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.

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

A CoordTrait experiment #103

Merged
merged 10 commits into from
Apr 5, 2024
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
Loading