Skip to content

Commit

Permalink
Simplify ParsedParameters (#58)
Browse files Browse the repository at this point in the history 
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
  • Loading branch information
@busstoptaktik
busstoptaktik authored Aug 11, 2023
1 parent 3ffc6cf commit dc5bcd6
Show file tree
Hide file tree
Showing 17 changed files with 248 additions and 259 deletions.
20 changes: 10 additions & 10 deletions src/context/minimal.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,24 +136,24 @@ 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.);

ctx.apply(op, Inv, &mut data)?;
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.);
Expand Down Expand Up @@ -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")?;
Expand Down Expand Up @@ -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,
Expand Down
38 changes: 20 additions & 18 deletions src/inner_op/btmerc.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -134,22 +134,24 @@ pub fn utm(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
}

// 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)));
Expand Down
14 changes: 8 additions & 6 deletions src/inner_op/cart.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand All @@ -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;
Expand Down
7 changes: 6 additions & 1 deletion src/inner_op/curvature.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -116,6 +116,11 @@ pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result<Op, Error> {
));
}

// 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)
}

Expand Down
4 changes: 2 additions & 2 deletions src/inner_op/geodesic.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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();
Expand Down
29 changes: 18 additions & 11 deletions src/inner_op/laea.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -191,7 +191,8 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
let def = &parameters.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()));
Expand All @@ -216,7 +217,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>

// --- 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();
Expand Down Expand Up @@ -266,6 +267,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
#[cfg(test)]
mod tests {
use super::*;
use float_eq::assert_float_eq;

#[test]
fn laea_oblique() -> Result<(), Error> {
Expand All @@ -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)?;
Expand Down
6 changes: 3 additions & 3 deletions src/inner_op/latitude.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -137,7 +137,7 @@ pub const GAMUT: [OpParameter; 8] = [

pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result<Op, Error> {
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") {
Expand Down
54 changes: 32 additions & 22 deletions src/inner_op/lcc.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 };
Expand Down Expand Up @@ -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 };
Expand Down Expand Up @@ -119,18 +121,25 @@ pub const GAMUT: [OpParameter; 9] = [
pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
let def = &parameters.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 {
Expand All @@ -140,8 +149,9 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>

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(
Expand Down Expand Up @@ -189,7 +199,7 @@ pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error>
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::<Op>::new();
Expand Down
Loading

0 comments on commit dc5bcd6

Please sign in to comment.