From dc5bcd6d9a77edfd5797cc93af23ab5ff785cd3a Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Fri, 11 Aug 2023 15:44:47 +0200 Subject: [PATCH] Simplify ParsedParameters (#58) Remove the special treatment of x_0.., y_0.., lat_0.., lon_0.. in ParsedParameters, so all parameters are treated identically. The accessors x, y, lat, lon and ellps are still available. This may make the code slightly slower, but also more accessible --- src/context/minimal.rs | 20 +++--- src/inner_op/btmerc.rs | 38 ++++++----- src/inner_op/cart.rs | 14 ++-- src/inner_op/curvature.rs | 7 +- src/inner_op/geodesic.rs | 4 +- src/inner_op/laea.rs | 29 ++++---- src/inner_op/latitude.rs | 6 +- src/inner_op/lcc.rs | 54 ++++++++------- src/inner_op/merc.rs | 29 ++++---- src/inner_op/molodensky.rs | 6 +- src/inner_op/omerc.rs | 35 +++++++--- src/inner_op/pipeline.rs | 2 +- src/inner_op/tmerc.rs | 109 ++++++++++-------------------- src/inner_op/webmerc.rs | 10 +-- src/op/mod.rs | 14 +++- src/op/parameter.rs | 2 +- src/op/parsed_parameters.rs | 128 +++++++++++++++--------------------- 17 files changed, 248 insertions(+), 259 deletions(-) diff --git a/src/context/minimal.rs b/src/context/minimal.rs index 194be4f..9b6d485 100644 --- a/src/context/minimal.rs +++ b/src/context/minimal.rs @@ -136,11 +136,17 @@ mod tests { ctx.register_resource("stupid:way", "addone | addone | addone inv"); let op = ctx.op("stupid:way")?; + let steps = ctx.steps(op)?; + assert_eq!(steps.len(), 3); + assert_eq!(steps[0], "addone"); + assert_eq!(steps[1], "addone"); + assert_eq!(steps[2], "addone inv"); + let mut data = some_basic_coor2dinates(); assert_eq!(data[0][0], 55.); assert_eq!(data[1][0], 59.); - ctx.apply(op, Fwd, &mut data)?; + assert_eq!(2, ctx.apply(op, Fwd, &mut data)?); assert_eq!(data[0][0], 56.); assert_eq!(data[1][0], 60.); @@ -148,12 +154,6 @@ mod tests { assert_eq!(data[0][0], 55.); assert_eq!(data[1][0], 59.); - let steps = ctx.steps(op)?; - assert_eq!(steps.len(), 3); - assert_eq!(steps[0], "addone"); - assert_eq!(steps[1], "addone"); - assert_eq!(steps[2], "addone inv"); - let params = ctx.params(op, 1)?; let ellps = params.ellps(0); assert_eq!(ellps.semimajor_axis(), 6378137.); @@ -195,9 +195,9 @@ mod tests { // while ellps_? defaults to GRS80 - so they are there even though we havent // set them let params = ctx.params(op, 1)?; - let ellps = params.ellps[0]; + let ellps = params.ellps(0); assert_eq!(ellps.semimajor_axis(), 6378137.); - assert_eq!(0., params.lat[0]); + assert_eq!(0., params.real("lat_0")?); // The zone id is found among the natural numbers (which here includes 0) let zone = params.natural("zone")?; @@ -225,7 +225,7 @@ mod tests { let op = ctx.op("utm zone=32")?; let steps = ctx.steps(op)?; assert!(steps.len() == 1); - let ellps = ctx.params(op, 0)?.ellps[0]; + let ellps = ctx.params(op, 0)?.ellps(0); let jac = Jacobian::new( &ctx, op, diff --git a/src/inner_op/btmerc.rs b/src/inner_op/btmerc.rs index 0d4ca47..98882b0 100644 --- a/src/inner_op/btmerc.rs +++ b/src/inner_op/btmerc.rs @@ -5,13 +5,13 @@ use crate::operator_authoring::*; // Forward transverse mercator, following Bowring (1989) fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let eps = ellps.second_eccentricity_squared(); - let lat_0 = op.params.lat[0]; - let lon_0 = op.params.lon[0]; - let x_0 = op.params.x[0]; - let y_0 = op.params.y[0]; - let k_0 = op.params.k[0]; + let lat_0 = op.params.lat(0).to_radians(); + let lon_0 = op.params.lon(0).to_radians(); + let x_0 = op.params.x(0); + let y_0 = op.params.y(0); + let k_0 = op.params.k(0); let mut successes = 0_usize; let n = operands.len(); @@ -53,13 +53,13 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // Inverse transverse mercator, following Bowring (1989) fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let eps = ellps.second_eccentricity_squared(); - let lat_0 = op.params.lat[0]; - let lon_0 = op.params.lon[0]; - let x_0 = op.params.x[0]; - let y_0 = op.params.y[0]; - let k_0 = op.params.k[0]; + let lat_0 = op.params.lat(0).to_radians(); + let lon_0 = op.params.lon(0).to_radians(); + let x_0 = op.params.x(0); + let y_0 = op.params.y(0); + let k_0 = op.params.k(0); let mut successes = 0_usize; let n = operands.len(); @@ -134,22 +134,24 @@ pub fn utm(parameters: &RawParameters, _ctx: &dyn Context) -> Result } // The scaling factor is 0.9996 by definition of UTM - params.k[0] = 0.9996; + params.real.insert("k_0", 0.9996); // The center meridian is determined by the zone - params.lon[0] = (-183. + 6. * zone as f64).to_radians(); + params + .real + .insert("lon_0", -183. + 6. * zone as f64); // The base parallel is by definition the equator - params.lat[0] = 0.0; + params.real.insert("lat_0", 0.); // The false easting is 500000 m by definition of UTM - params.x[0] = 500000.0; + params.real.insert("x_0", 500_000.); // The false northing is 0 m by definition of UTM - params.y[0] = 0.0; + params.real.insert("y_0", 0.); // or 10_000_000 m if using the southern aspect if params.boolean("south") { - params.y[0] = 10_000_000.0; + params.real.insert("y_0", 10_000_000.0); } let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv))); diff --git a/src/inner_op/cart.rs b/src/inner_op/cart.rs index b9fe5ed..bf421fe 100644 --- a/src/inner_op/cart.rs +++ b/src/inner_op/cart.rs @@ -6,9 +6,10 @@ use crate::operator_authoring::*; fn cart_fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let n = operands.len(); let mut successes = 0; + let ellps = op.params.ellps(0); for i in 0..n { let mut coord = operands.get_coord(i); - coord = op.params.ellps[0].cartesian(&coord); + coord = ellps.cartesian(&coord); if !coord.0.iter().any(|c| c.is_nan()) { successes += 1; } @@ -20,20 +21,21 @@ fn cart_fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> us // ----- I N V E R S E -------------------------------------------------------------- fn cart_inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { + let ellps = op.params.ellps(0); // eccentricity squared, Fukushima's E, Claessens' c3 = 1-c2` - let es = op.params.ellps[0].eccentricity_squared(); + let es = ellps.eccentricity_squared(); // semiminor axis - let b = op.params.ellps[0].semiminor_axis(); + let b = ellps.semiminor_axis(); // semimajor axis - let a = op.params.ellps[0].semimajor_axis(); + let a = ellps.semimajor_axis(); // reciproque of a - let ra = 1. / op.params.ellps[0].semimajor_axis(); + let ra = 1. / ellps.semimajor_axis(); // aspect ratio, b/a: Fukushima's ec, Claessens' c4 let ar = b * ra; // 1.5 times the fourth power of the eccentricity let ce4 = 1.5 * es * es; // if we're closer than this to the Z axis, we force latitude to one of the poles - let cutoff = op.params.ellps[0].semimajor_axis() * 1e-16; + let cutoff = ellps.semimajor_axis() * 1e-16; let n = operands.len(); let mut successes = 0; diff --git a/src/inner_op/curvature.rs b/src/inner_op/curvature.rs index a87443f..691165e 100644 --- a/src/inner_op/curvature.rs +++ b/src/inner_op/curvature.rs @@ -6,7 +6,7 @@ use crate::operator_authoring::*; fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let n = operands.len(); let sliced = 0..n; - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let prime = op.params.boolean("prime"); let meridional = op.params.boolean("meridian"); @@ -116,6 +116,11 @@ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { )); } + // Check that we have a proper ellipsoid + if op.params.text("ellps").is_ok() { + let _ = Ellipsoid::named(op.params.text("ellps")?.as_str())?; + } + Ok(op) } diff --git a/src/inner_op/geodesic.rs b/src/inner_op/geodesic.rs index d218b60..3e0e9d6 100644 --- a/src/inner_op/geodesic.rs +++ b/src/inner_op/geodesic.rs @@ -4,7 +4,7 @@ use crate::operator_authoring::*; // ----- F O R W A R D ----------------------------------------------------------------- fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let n = operands.len(); let sliced = 0..n; @@ -35,7 +35,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // ----- I N V E R S E ----------------------------------------------------------------- fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let reversible = op.params.boolean("reversible"); let n = operands.len(); diff --git a/src/inner_op/laea.rs b/src/inner_op/laea.rs index c530896..47f7cf9 100644 --- a/src/inner_op/laea.rs +++ b/src/inner_op/laea.rs @@ -20,9 +20,9 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let north_polar = op.params.boolean("north_polar"); let south_polar = op.params.boolean("south_polar"); - let lon_0 = op.params.lon(0); - let x_0 = op.params.x(0); - let y_0 = op.params.y(0); + let lon_0 = op.params.real("lon_0").unwrap_or(0.).to_radians(); + let x_0 = op.params.real("x_0").unwrap_or(0.); + let y_0 = op.params.real("y_0").unwrap_or(0.); let ellps = op.params.ellps(0); let e = ellps.eccentricity(); let a = ellps.semimajor_axis(); @@ -95,10 +95,10 @@ fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let north_polar = op.params.boolean("north_polar"); let south_polar = op.params.boolean("south_polar"); - let lon_0 = op.params.lon(0); - let lat_0 = op.params.lat(0); - let x_0 = op.params.x(0); - let y_0 = op.params.y(0); + let lon_0 = op.params.real("lon_0").unwrap_or(0.).to_radians(); + let lat_0 = op.params.real("lat_0").unwrap_or(0.).to_radians(); + let x_0 = op.params.real("x_0").unwrap_or(0.); + let y_0 = op.params.real("y_0").unwrap_or(0.); let ellps = op.params.ellps(0); let a = ellps.semimajor_axis(); @@ -191,7 +191,8 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result let def = ¶meters.definition; let mut params = ParsedParameters::new(parameters, &GAMUT)?; - let lat_0 = params.lat[0]; + let lat_0 = params.real("lat_0").unwrap_or(0.).to_radians(); + if lat_0.is_nan() { warn!("LAEA: Bad central latitude!"); return Err(Error::BadParam("lat_0".to_string(), def.clone())); @@ -216,7 +217,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result // --- Precompute some latitude invariant factors --- - let ellps = params.ellps[0]; + let ellps = params.ellps(0); let a = ellps.semimajor_axis(); let es = ellps.eccentricity_squared(); let e = es.sqrt(); @@ -266,6 +267,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result #[cfg(test)] mod tests { use super::*; + use float_eq::assert_float_eq; #[test] fn laea_oblique() -> Result<(), Error> { @@ -275,11 +277,16 @@ mod tests { let op = ctx.op("laea ellps=GRS80 lat_0=52 lon_0=10 x_0=4321000 y_0=3210000")?; // The test point from IOGP - let p = Coor4D::geo(50.0, 5.0, 0.0, 0.0); - let mut operands = [p]; + let p = Coor2D::geo(50.0, 5.0); + let geo = [p]; + let p = Coor2D::raw(3962799.45, 2999718.85); + let projected = [p]; + + let mut operands = geo.clone(); // Forward ctx.apply(op, Fwd, &mut operands)?; + assert_float_eq!(operands[0].0, projected[0].0, abs_all <= 0.01); assert!((operands[0][0] - 3962799.45).abs() < 0.01); assert!((operands[0][1] - 2999718.85).abs() < 0.01); ctx.apply(op, Inv, &mut operands)?; diff --git a/src/inner_op/latitude.rs b/src/inner_op/latitude.rs index 8f324e3..16f007c 100644 --- a/src/inner_op/latitude.rs +++ b/src/inner_op/latitude.rs @@ -7,7 +7,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let mut successes = 0_usize; let n = operands.len(); let sliced = 0..n; - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); if op.params.boolean("geocentric") { for i in sliced { @@ -66,7 +66,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let mut successes = 0_usize; let n = operands.len(); - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); if op.params.boolean("geocentric") { for i in 0..n { @@ -137,7 +137,7 @@ pub const GAMUT: [OpParameter; 8] = [ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { let mut op = Op::plain(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT, ctx)?; - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let mut number_of_flags = 0_usize; if op.params.boolean("geocentric") { diff --git a/src/inner_op/lcc.rs b/src/inner_op/lcc.rs index 5e4e49b..1eed9df 100644 --- a/src/inner_op/lcc.rs +++ b/src/inner_op/lcc.rs @@ -9,12 +9,13 @@ const EPS10: f64 = 1e-10; // Forward Lambert conformal conic, following the PROJ implementation, // cf. https://proj.org/operations/projections/lcc.html fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let a = op.params.ellps[0].semimajor_axis(); - let e = op.params.ellps[0].eccentricity(); - let lon_0 = op.params.lon[0]; - let k_0 = op.params.k[0]; - let x_0 = op.params.x[0]; - let y_0 = op.params.y[0]; + let ellps = op.params.ellps(0); + let a = ellps.semimajor_axis(); + let e = ellps.eccentricity(); + let lon_0 = op.params.lon(0); + let k_0 = op.params.k(0); + let x_0 = op.params.x(0); + let y_0 = op.params.y(0); let Ok(n) = op.params.real("n") else { return 0 }; let Ok(c) = op.params.real("c") else { return 0 }; let Ok(rho0) = op.params.real("rho0") else { return 0 }; @@ -48,12 +49,13 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // ----- I N V E R S E ----------------------------------------------------------------- fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let a = op.params.ellps[0].semimajor_axis(); - let e = op.params.ellps[0].eccentricity(); - let lon_0 = op.params.lon[0]; - let k_0 = op.params.k[0]; - let x_0 = op.params.x[0]; - let y_0 = op.params.y[0]; + let ellps = op.params.ellps(0); + let a = ellps.semimajor_axis(); + let e = ellps.eccentricity(); + let lon_0 = op.params.lon(0); + let k_0 = op.params.k(0); + let x_0 = op.params.x(0); + let y_0 = op.params.y(0); let Ok(n) = op.params.real("n") else { return 0 }; let Ok(c) = op.params.real("c") else { return 0 }; let Ok(rho0) = op.params.real("rho0") else { return 0 }; @@ -119,18 +121,25 @@ pub const GAMUT: [OpParameter; 9] = [ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result { let def = ¶meters.definition; let mut params = ParsedParameters::new(parameters, &GAMUT)?; - if params.lat[2].is_nan() { - params.lat[2] = params.lat[1]; + if !params.real.contains_key("lat_2") { + params.real.insert("lat_2", params.lat(1)); } - let phi1 = params.lat[1]; - let mut phi2 = params.lat[2]; + let phi1 = params.lat(1).to_radians(); + let mut phi2 = params.lat(2).to_radians(); if phi2.is_nan() { phi2 = phi1; } - params.lat[2] = phi2; - - let mut lat_0 = params.lat[0]; + params + .real + .insert("lon_0", params.real["lon_0"].to_radians()); + params + .real + .insert("lat_0", params.real["lat_0"].to_radians()); + params.real.insert("lat_1", phi1); + params.real.insert("lat_2", phi2); + + let mut lat_0 = params.lat(0); if lat_0.is_nan() { lat_0 = 0.; if (phi1 - phi2).abs() < EPS10 { @@ -140,8 +149,9 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result let sc = phi1.sin_cos(); let mut n = sc.0; - let e = params.ellps[0].eccentricity(); - let es = params.ellps[0].eccentricity_squared(); + let ellps = params.ellps(0); + let e = ellps.eccentricity(); + let es = ellps.eccentricity_squared(); if (phi1 + phi2).abs() < EPS10 { return Err(Error::General( @@ -189,7 +199,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result params.real.insert("c", c); params.real.insert("n", n); params.real.insert("rho0", rho0); - params.lat[0] = lat_0; + params.real.insert("lat_0", lat_0); let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv))); let steps = Vec::::new(); diff --git a/src/inner_op/merc.rs b/src/inner_op/merc.rs index d069ede..69e1372 100644 --- a/src/inner_op/merc.rs +++ b/src/inner_op/merc.rs @@ -4,13 +4,13 @@ use crate::operator_authoring::*; // ----- F O R W A R D ----------------------------------------------------------------- fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let a = ellps.semimajor_axis(); - let k_0 = op.params.k[0]; - let x_0 = op.params.x[0]; - let y_0 = op.params.y[0]; - let lat_0 = op.params.lat[0]; - let lon_0 = op.params.lon[0]; + let k_0 = op.params.k(0); + let x_0 = op.params.x(0); + let y_0 = op.params.y(0); + let lat_0 = op.params.lat(0); + let lon_0 = op.params.lon(0); let mut successes = 0_usize; let length = operands.len(); @@ -33,13 +33,13 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // ----- I N V E R S E ----------------------------------------------------------------- fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let a = ellps.semimajor_axis(); - let k_0 = op.params.k[0]; - let x_0 = op.params.x[0]; - let y_0 = op.params.y[0]; - let lat_0 = op.params.lat[0]; - let lon_0 = op.params.lon[0]; + let k_0 = op.params.k(0); + let x_0 = op.params.x(0); + let y_0 = op.params.y(0); + let lat_0 = op.params.lat(0); + let lon_0 = op.params.lon(0); let mut successes = 0_usize; let length = operands.len(); @@ -80,7 +80,7 @@ pub const GAMUT: [OpParameter; 8] = [ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result { let def = ¶meters.definition; let mut params = ParsedParameters::new(parameters, &GAMUT)?; - let ellps = params.ellps[0]; + let ellps = params.ellps(0); let lat_ts = params.real("lat_ts")?; if lat_ts.abs() > 90. { @@ -92,7 +92,8 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result // lat_ts trumps k_0 if lat_ts != 0.0 { let sc = lat_ts.to_radians().sin_cos(); - params.k[0] = sc.1 / (1. - ellps.eccentricity_squared() * sc.0 * sc.0).sqrt() + let k_0 = sc.1 / (1. - ellps.eccentricity_squared() * sc.0 * sc.0).sqrt(); + params.real.insert("k_0", k_0); } let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv))); diff --git a/src/inner_op/molodensky.rs b/src/inner_op/molodensky.rs index 35af18f..b6c8272 100644 --- a/src/inner_op/molodensky.rs +++ b/src/inner_op/molodensky.rs @@ -25,7 +25,7 @@ fn common( operands: &mut dyn CoordinateSet, direction: Direction, ) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let a = ellps.semimajor_axis(); let f = ellps.flattening(); let es = ellps.eccentricity_squared(); @@ -100,8 +100,8 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result let def = ¶meters.definition; let mut params = ParsedParameters::new(parameters, &GAMUT)?; - let ellps_0 = params.ellps[0]; - let ellps_1 = params.ellps[1]; + let ellps_0 = params.ellps(0); + let ellps_1 = params.ellps(1); // We may use `ellps, da, df`, to parameterize the op, but `ellps_0, ellps_1` // is a more likely set of parameters to come across in real life. diff --git a/src/inner_op/omerc.rs b/src/inner_op/omerc.rs index 0355dee..1929ce8 100644 --- a/src/inner_op/omerc.rs +++ b/src/inner_op/omerc.rs @@ -9,14 +9,14 @@ use std::f64::consts::FRAC_PI_4; #[allow(non_snake_case)] fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let es = ellps.eccentricity_squared(); let e = es.sqrt(); - let kc = op.params.k[0]; + let kc = op.params.k(0); - let FE = op.params.x[0]; - let FN = op.params.y[0]; + let FE = op.params.x(0); + let FN = op.params.y(0); let Ec = FE; let Nc = FN; @@ -92,6 +92,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { let u = A * (S * c0 + V * s0).atan2(cblon) / B; coord[0] = v * cc + u * sc + FE; coord[1] = u * cc - v * sc + FN; + successes += 1; continue; } @@ -106,6 +107,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { }; coord[0] = v * cc + u * sc + Ec; coord[1] = u * cc - v * sc + Nc; + successes += 1; continue; } @@ -125,14 +127,14 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { #[allow(non_snake_case)] fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let es = ellps.eccentricity_squared(); let e = es.sqrt(); - let kc = op.params.k[0]; + let kc = op.params.k(0); - let FE = op.params.x[0]; - let FN = op.params.y[0]; + let FE = op.params.x(0); + let FN = op.params.y(0); let latc = op.params.real["latc"].to_radians(); let lonc = op.params.real["lonc"].to_radians(); @@ -276,6 +278,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result #[cfg(test)] mod tests { use super::*; + use float_eq::assert_float_eq; #[test] fn omerc() -> Result<(), Error> { @@ -291,15 +294,27 @@ mod tests { // Validation value from EPSG let geo = [Coor2D::geo(5.3872535833, 115.8055054444)]; - let projected = [Coor2D::raw(679245.7281740266, 596562.7774687681)]; // Forward let mut operands = geo.clone(); + assert_eq!(1, ctx.apply(op, Fwd, &mut operands)?); + for i in 0..operands.len() { + assert_float_eq!(operands[i].0, projected[i].0, abs_all <= 1e-9); + } + + // Roundtrip + assert_eq!(1, ctx.apply(op, Inv, &mut operands)?); + for i in 0..operands.len() { + assert_float_eq!(operands[i].0, geo[i].0, abs_all <= 1e-9); + } + + // Forward + let mut operands = geo.clone(); + ctx.apply(op, Fwd, &mut operands)?; for i in 0..operands.len() { - dbg!(operands[i]); assert!(operands[i].hypot2(&projected[i]) < 1e-9); } diff --git a/src/inner_op/pipeline.rs b/src/inner_op/pipeline.rs index 8f3746e..2b8558b 100644 --- a/src/inner_op/pipeline.rs +++ b/src/inner_op/pipeline.rs @@ -194,7 +194,7 @@ mod tests { let op = ctx.op("addone|addone|addone")?; let mut data = some_basic_coor2dinates(); - ctx.apply(op, Fwd, &mut data)?; + assert_eq!(2, ctx.apply(op, Fwd, &mut data)?); assert_eq!(data[0][0], 58.); assert_eq!(data[1][0], 62.); diff --git a/src/inner_op/tmerc.rs b/src/inner_op/tmerc.rs index f660f4a..71543b9 100644 --- a/src/inner_op/tmerc.rs +++ b/src/inner_op/tmerc.rs @@ -7,10 +7,10 @@ use crate::operator_authoring::*; // Forward transverse mercator, following Engsager & Poder(2007) fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // Make all precomputed parameters directly accessible - let ellps = op.params.ellps[0]; - let lat_0 = op.params.lat[0]; - let lon_0 = op.params.lon[0]; - let x_0 = op.params.x[0]; + let ellps = op.params.ellps(0); + let lat_0 = op.params.lat(0).to_radians(); + let lon_0 = op.params.lon(0).to_radians(); + let x_0 = op.params.x(0); let Some(conformal) = op.params.fourier_coefficients.get("conformal") else { warn!("Missing Fourier coefficients for conformal mapping!"); return 0; @@ -94,12 +94,12 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // ----- I N V E R S E ----------------------------------------------------------------- -// Inverse Transverse Mercator, following Engsager & Poder (2007) (currently Bowring stands in!) +// Inverse Transverse Mercator, following Engsager & Poder (2007) fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // Make all precomputed parameters directly accessible - let ellps = op.params.ellps[0]; - let lon_0 = op.params.lon[0]; - let x_0 = op.params.x[0]; + let ellps = op.params.ellps(0); + let lon_0 = op.params.lon(0).to_radians(); + let x_0 = op.params.x(0); let Some(conformal) = op.params.fourier_coefficients.get("conformal") else { warn!("Missing Fourier coefficients for conformal mapping!"); return 0; @@ -201,22 +201,24 @@ pub fn utm(parameters: &RawParameters, _ctx: &dyn Context) -> Result } // The scaling factor is 0.9996 by definition of UTM - params.k[0] = 0.9996; + params.real.insert("k_0", 0.9996); // The center meridian is determined by the zone - params.lon[0] = (-183. + 6. * zone as f64).to_radians(); + params + .real + .insert("lon_0", -183. + 6. * zone as f64); // The base parallel is by definition the equator - params.lat[0] = 0.0; + params.real.insert("lat_0", 0.); // The false easting is 500000 m by definition of UTM - params.x[0] = 500_000.0; + params.real.insert("x_0", 500_000.); // The false northing is 0 m by definition of UTM - params.y[0] = 0.0; + params.real.insert("y_0", 0.); // or 10_000_000 m if using the southern aspect if params.boolean("south") { - params.y[0] = 10_000_000.0; + params.real.insert("y_0", 10_000_000.0); } let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv))); @@ -263,13 +265,13 @@ const TRANSVERSE_MERCATOR: PolynomialCoefficients = PolynomialCoefficients { // Pre-compute some of the computationally heavy prerequisites, // to get better amortization over the full operator lifetime. fn precompute(op: &mut Op) { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let n = ellps.third_flattening(); - let lat_0 = op.params.lat[0]; - let y_0 = op.params.y[0]; + let lat_0 = op.params.lat(0).to_radians(); + let y_0 = op.params.y(0); // The scaled spherical Earth radius - Qn in Engsager's implementation - let qs = op.params.k[0] * ellps.semimajor_axis() * ellps.normalized_meridian_arc_unit(); + let qs = op.params.k(0) * ellps.semimajor_axis() * ellps.normalized_meridian_arc_unit(); op.params.real.insert("scaled_radius", qs); // The Fourier series for the conformal latitude @@ -303,51 +305,10 @@ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { #[cfg(test)] mod tests { use super::*; + use float_eq::assert_float_eq; #[test] fn tmerc() -> Result<(), Error> { - // Validation values from PROJ: - // echo 12 55 0 0 | cct -d18 +proj=utm +zone=32 | clip - #[rustfmt::skip] - let geo = [ - Coor4D::geo( 55., 12., 0., 0.), - Coor4D::geo(-55., 12., 0., 0.), - Coor4D::geo( 55., -6., 0., 0.), - Coor4D::geo(-55., -6., 0., 0.) - ]; - - #[rustfmt::skip] - let projected = [ - Coor4D::raw( 691_875.632_139_661, 6_098_907.825_005_012, 0., 0.), - Coor4D::raw( 691_875.632_139_661,-6_098_907.825_005_012, 0., 0.), - Coor4D::raw(-455_673.814_189_040, 6_198_246.671_090_279, 0., 0.), - Coor4D::raw(-455_673.814_189_040,-6_198_246.671_090_279, 0., 0.) - ]; - - let mut ctx = Minimal::default(); - let definition = "tmerc k_0=0.9996 lon_0=9 x_0=500000"; - let op = ctx.op(definition)?; - - let mut operands = geo.clone(); - ctx.apply(op, Fwd, &mut operands)?; - - for i in 0..operands.len() { - dbg!(operands[i]); - dbg!(projected[i]); - assert!(operands[i].hypot2(&projected[i]) < 1e-6); - } - - ctx.apply(op, Inv, &mut operands)?; - for i in 0..operands.len() { - assert!(operands[i].hypot2(&geo[i]) < 5e-6); - } - - Ok(()) - } - - // Same as above, but using the 2D Coor2D data structure - #[test] - fn coor2d() -> Result<(), Error> { // Validation values from PROJ: // echo 12 55 0 0 | cct -d18 +proj=utm +zone=32 | clip #[rustfmt::skip] @@ -374,9 +335,7 @@ mod tests { ctx.apply(op, Fwd, &mut operands)?; for i in 0..operands.len() { - dbg!(operands[i]); - dbg!(projected[i]); - assert!(operands[i].hypot2(&projected[i]) < 1e-6); + assert_float_eq!(operands[i].0, projected[i].0, abs_all <= 1e-8); } ctx.apply(op, Inv, &mut operands)?; @@ -397,29 +356,29 @@ mod tests { // echo 12 55 0 0 | cct -d18 +proj=utm +zone=32 | clip #[rustfmt::skip] let geo = [ - Coor4D::geo( 55., 12., 0., 0.), - Coor4D::geo(-55., 12., 0., 0.), - Coor4D::geo( 55., -6., 0., 0.), - Coor4D::geo(-55., -6., 0., 0.) + Coor2D::geo( 55., 12.), + Coor2D::geo(-55., 12.), + Coor2D::geo( 55., -6.), + Coor2D::geo(-55., -6.) ]; #[rustfmt::skip] let projected = [ - Coor4D::raw( 691_875.632_139_661, 6_098_907.825_005_012, 0., 0.), - Coor4D::raw( 691_875.632_139_661,-6_098_907.825_005_012, 0., 0.), - Coor4D::raw(-455_673.814_189_040, 6_198_246.671_090_279, 0., 0.), - Coor4D::raw(-455_673.814_189_040,-6_198_246.671_090_279, 0., 0.) + Coor2D::raw( 691_875.632_139_661, 6_098_907.825_005_012), + Coor2D::raw( 691_875.632_139_661,-6_098_907.825_005_012), + Coor2D::raw(-455_673.814_189_040, 6_198_246.671_090_279), + Coor2D::raw(-455_673.814_189_040,-6_198_246.671_090_279) ]; let mut operands = geo.clone(); - ctx.apply(op, Fwd, &mut operands)?; + assert_eq!(ctx.apply(op, Fwd, &mut operands)?, 4); for i in 0..operands.len() { - assert!(operands[i].hypot2(&projected[i]) < 5e-3); + assert_float_eq!(operands[i].0, projected[i].0, abs_all <= 1e-8); } - ctx.apply(op, Inv, &mut operands)?; + assert_eq!(ctx.apply(op, Inv, &mut operands)?, 4); for i in 0..operands.len() { - assert!(operands[i].hypot2(&geo[i]) < 10e-8); + assert_float_eq!(operands[i].0, geo[i].0, abs_all <= 1e-12); } Ok(()) diff --git a/src/inner_op/webmerc.rs b/src/inner_op/webmerc.rs index cbe4186..0da3bde 100644 --- a/src/inner_op/webmerc.rs +++ b/src/inner_op/webmerc.rs @@ -6,7 +6,7 @@ use std::f64::consts::FRAC_PI_4; // ----- F O R W A R D ----------------------------------------------------------------- fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let a = ellps.semimajor_axis(); let mut successes = 0_usize; @@ -30,7 +30,7 @@ fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { // ----- I N V E R S E ----------------------------------------------------------------- fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize { - let ellps = op.params.ellps[0]; + let ellps = op.params.ellps(0); let a = ellps.semimajor_axis(); let mut successes = 0_usize; @@ -80,6 +80,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result #[cfg(test)] mod tests { use super::*; + use float_eq::assert_float_eq; #[test] fn webmerc() -> Result<(), Error> { @@ -102,14 +103,13 @@ mod tests { let mut operands = geo.clone(); ctx.apply(op, Fwd, &mut operands)?; for i in 0..operands.len() { - assert!(operands[i].hypot2(&projected[i]) < 1e-8); + assert_float_eq!(operands[i].0, projected[i].0, abs_all <= 1e-8); } // Roundtrip ctx.apply(op, Inv, &mut operands)?; for i in 0..operands.len() { - dbg!(operands[0].to_geo()); - assert!(operands[i].hypot2(&geo[i]) < 2e-9); + assert_float_eq!(operands[i].0, geo[i].0, abs_all <= 2e-9); } Ok(()) diff --git a/src/op/mod.rs b/src/op/mod.rs index 4689216..0a62e72 100644 --- a/src/op/mod.rs +++ b/src/op/mod.rs @@ -67,7 +67,19 @@ impl Op { _ctx: &dyn Context, ) -> Result { let def = parameters.definition.as_str(); - let params = ParsedParameters::new(parameters, gamut)?; + let mut params = ParsedParameters::new(parameters, gamut)?; + + // Convert lat_{0..4} and lon_{0..4} to radians + for i in ["lat_0", "lat_1", "lat_2", "lat_3"] { + let lat = *params.real.get(i).unwrap_or(&0.); + params.real.insert(i, lat); + } + + for i in ["lon_0", "lon_1", "lon_2", "lon_3"] { + let lon = *params.real.get(i).unwrap_or(&0.); + params.real.insert(i, lon); + } + let descriptor = OpDescriptor::new(def, fwd, inv); let steps = Vec::::new(); let id = OpHandle::new(); diff --git a/src/op/parameter.rs b/src/op/parameter.rs index 0344158..3ddc7dc 100644 --- a/src/op/parameter.rs +++ b/src/op/parameter.rs @@ -12,7 +12,7 @@ /// /// For a given operation, the union of the sets of its required and optional /// parameters is called the *gamut* of the operation. -#[derive(Debug)] +#[derive(Debug, PartialEq, Clone)] pub enum OpParameter { /// A flag is a boolean that is true if present, false if not Flag { key: &'static str }, diff --git a/src/op/parsed_parameters.rs b/src/op/parsed_parameters.rs index d7c57f8..b64daf4 100644 --- a/src/op/parsed_parameters.rs +++ b/src/op/parsed_parameters.rs @@ -11,18 +11,23 @@ use std::collections::BTreeSet; use super::*; +#[rustfmt::skip] +const ZERO_VALUED_IMPLICIT_GAMUT_ELEMENTS: [&str; 16] = [ + "x_0", "x_1", "x_2", "x_3", + "y_0", "y_1", "y_2", "y_3", + "lat_0", "lat_1", "lat_2", "lat_3", + "lon_0", "lon_1", "lon_2", "lon_3" +]; + +#[rustfmt::skip] +const UNIT_VALUED_IMPLICIT_GAMUT_ELEMENTS: [&str; 4] = [ + "k_0", "k_1", "k_2", "k_3" +]; + #[derive(Debug, Clone)] pub struct ParsedParameters { pub name: String, - // Commonly used options have hard-coded slots - pub ellps: [Ellipsoid; 2], - pub lat: [f64; 4], - pub lon: [f64; 4], - pub x: [f64; 4], - pub y: [f64; 4], - pub k: [f64; 4], - // Op-specific options are stored in B-Trees pub boolean: BTreeSet<&'static str>, pub natural: BTreeMap<&'static str, usize>, @@ -88,23 +93,34 @@ impl ParsedParameters { pub fn ignored(&self) -> Vec { self.ignored.clone() } - pub fn ellps(&self, index: usize) -> &Ellipsoid { - &self.ellps[index] + pub fn ellps(&self, index: usize) -> Ellipsoid { + // if 'ellps' was explicitly given, it will override 'ellps_0' + if index == 0 { + if let Some(e) = self.text.get("ellps") { + return Ellipsoid::named(e).unwrap(); + } + } + let key = format!("ellps_{index}"); + if let Some(e) = self.text.get(&key[..]) { + return Ellipsoid::named(e).unwrap(); + } + // If none of them existed, i.e. no defaults were given, we return the general default + Ellipsoid::default() + } + pub fn k(&self, index: usize) -> f64 { + *(self.real.get(&format!("k_{index}")[..]).unwrap_or(&1.)) } pub fn x(&self, index: usize) -> f64 { - self.x[index] + *(self.real.get(&format!("x_{index}")[..]).unwrap_or(&0.)) } pub fn y(&self, index: usize) -> f64 { - self.y[index] + *(self.real.get(&format!("y_{index}")[..]).unwrap_or(&0.)) } pub fn lat(&self, index: usize) -> f64 { - self.lat[index] + *self.real.get(&format!("lat_{index}")[..]).unwrap_or(&0.) } pub fn lon(&self, index: usize) -> f64 { - self.lon[index] - } - pub fn k(&self, index: usize) -> f64 { - self.k[index] + *self.real.get(&format!("lon_{index}")[..]).unwrap_or(&0.) } } @@ -127,18 +143,9 @@ impl ParsedParameters { let mut uuid = BTreeMap::<&'static str, uuid::Uuid>::new(); let fourier_coefficients = BTreeMap::<&'static str, FourierCoefficients>::new(); - // 'omit_fwd'/'omit_inv' are implicitly valid for all operators (including pipelines), - // so we append them to the gamut - let mut gamutt = Vec::new(); - for g in gamut { - gamutt.push(g); - } - gamutt.push(&OpParameter::Flag { key: "omit_fwd" }); - gamutt.push(&OpParameter::Flag { key: "omit_inv" }); - // Try to locate all accepted parameters, type check, and place them into // their proper bins - for p in gamutt { + for p in gamut { match *p { OpParameter::Flag { key } => { if let Some(value) = chase(globals, &locals, key)? { @@ -311,55 +318,30 @@ impl ParsedParameters { }; } - // Now handle the commonly used options with the hard-coded slots + // Default gamut elements - traditionally supported for all operators - let mut ellps = [Ellipsoid::default(), Ellipsoid::default()]; - let mut lat = [0.; 4]; - let mut lon = [0.; 4]; - let mut x = [0.; 4]; - let mut y = [0.; 4]; - let mut k = [1.; 4]; - - // ellps_{n} - for i in 0..2 { - let key = format!("ellps_{i}"); - if let Some(e) = text.get(&key[..]) { - ellps[i] = Ellipsoid::named(e)?; + // omit_fwd and omit_inv are implicitly valid for all ops + if let Some(value) = chase(globals, &locals, "omit_fwd")? { + if value.is_empty() || value.to_lowercase() == "true" { + boolean.insert("omit_fwd"); } } - // But `ellps` trumps `ellps_0` - if let Some(e) = text.get("ellps") { - ellps[0] = Ellipsoid::named(e)?; - } - - // lat_{n} - for i in 0..4 { - let key = format!("lat_{i}"); - lat[i] = (*real.get(&key[..]).unwrap_or(&0.)).to_radians(); - } - - // lon_{n} - for i in 0..4 { - let key = format!("lon_{i}"); - lon[i] = (*real.get(&key[..]).unwrap_or(&0.)).to_radians(); - } - - // x_{n} - for i in 0..4 { - let key = format!("x_{i}"); - x[i] = *real.get(&key[..]).unwrap_or(&0.); + if let Some(value) = chase(globals, &locals, "omit_inv")? { + if value.is_empty() || value.to_lowercase() == "true" { + boolean.insert("omit_inv"); + } } - // y_{n} - for i in 0..4 { - let key = format!("y_{i}"); - y[i] = *real.get(&key[..]).unwrap_or(&0.); + for k in ZERO_VALUED_IMPLICIT_GAMUT_ELEMENTS { + if !real.contains_key(k) { + real.insert(k, 0.); + } } - // k_{n} - for i in 0..4 { - let key = format!("k_{i}"); - k[i] = *real.get(&key[..]).unwrap_or(&0.); + for k in UNIT_VALUED_IMPLICIT_GAMUT_ELEMENTS { + if !real.contains_key(k) { + real.insert(k, 1.); + } } let name = locals @@ -374,12 +356,6 @@ impl ParsedParameters { let given = locals.clone(); let ignored: Vec = locals.into_keys().collect(); Ok(ParsedParameters { - ellps, - lat, - lon, - x, - y, - k, name, boolean, natural, @@ -522,7 +498,7 @@ mod tests { assert_eq!(*p.text.get("text").unwrap(), "text"); assert_eq!( - p.ellps[0].semimajor_axis(), + p.ellps(0).semimajor_axis(), Ellipsoid::new(123., 1. / 456.).semimajor_axis() );