38 add senmerc operator (#44)

* Add senmerc operator

* update change log

CHANGELOG.md

@@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Support for 2D, 3D and 4D coordinates as tuples ([number, number]) or objects ({x, y, z, t}) to the wrapper
 - Unstable support for loading a grid from the network
 - Support for the `axisswap` operator via [Geodesy:#84](https://github.com/busstoptaktik/geodesy/pull/84)
+- `senmerc` operator [#38](https://github.com/Rennzie/geodesy-wasm/issues/38)
 
 ### Changed
 

src/geodesy/operators/senmerc.rs

@@ -1,55 +1,123 @@
-/// An as yet undefined `Sensat Mercator` operator.
-/// Webmercator with a scaling applied to the z coordinate.
+//! Sensat Mercator
+//! A bastardisation for 3D where the Z value is scaled by 1 / cos(lat) in the forward direction
 use geodesy_rs::authoring::*;
-
-// ----- C O M M O N -------------------------------------------------------------------
+use std::f64::consts::FRAC_PI_2;
+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 {
-    operands.len()
+fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
+    let ellps = op.params.ellps(0);
+    let a = ellps.semimajor_axis();
+
+    let mut successes = 0_usize;
+    let length = operands.len();
+    for i in 0..length {
+        let mut coord = operands.get_coord(i);
+        // Easting
+        coord[0] *= a;
+
+        // Northing
+        let lat = coord[1];
+        coord[1] = a * (FRAC_PI_4 + lat / 2.0).tan().ln();
+
+        // Altitude
+        coord[2] *= 1.0 / lat.cos();
+
+        operands.set_coord(i, &coord);
+        successes += 1;
+    }
+
+    successes
 }
 
 // ----- I N V E R S E -----------------------------------------------------------------
 
-fn inv(_op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
-    operands.len()
+fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
+    let ellps = op.params.ellps(0);
+    let a = ellps.semimajor_axis();
+
+    let mut successes = 0_usize;
+    let length = operands.len();
+    for i in 0..length {
+        let mut coord = operands.get_coord(i);
+
+        // Easting -> Longitude
+        coord[0] /= a;
+
+        // Northing -> Latitude
+        coord[1] = FRAC_PI_2 - 2.0 * (-coord[1] / a).exp().atan();
+
+        // Altitude
+        coord[2] *= coord[1].cos();
+
+        operands.set_coord(i, &coord);
+        successes += 1;
+    }
+
+    successes
 }
 
 // ----- C O N S T R U C T O R ---------------------------------------------------------
 
 #[rustfmt::skip]
-pub const GAMUT: [OpParameter; 0] = [
+pub const GAMUT: [OpParameter; 2] = [
+    OpParameter::Flag { key: "inv" },
+    OpParameter::Text { key: "ellps",  default: Some("WGS84") },
 ];
 
 pub fn new(parameters: &RawParameters, _ctx: &dyn Context) -> Result<Op, Error> {
-    Op::plain(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT, _ctx)
+    let def = &parameters.definition;
+    let params = ParsedParameters::new(parameters, &GAMUT)?;
+
+    let descriptor = OpDescriptor::new(def, InnerOp(fwd), Some(InnerOp(inv)));
+    let steps = Vec::<Op>::new();
+    let id = OpHandle::new();
+
+    Ok(Op {
+        descriptor,
+        params,
+        steps,
+        id,
+    })
 }
 
 // ----- T E S T S ---------------------------------------------------------------------
 
 #[cfg(test)]
 mod tests {
     use super::*;
-    const GDA94: Coor4D = Coor4D([-4052051.7643, 4212836.2017, -2545106.0245, 0.0]);
+    use float_eq::assert_float_eq;
 
     #[test]
-    // Borrowed from Rust Geodesy `noop` inner_op test
-    fn no_change() -> Result<(), Error> {
+    fn senmerc() -> Result<(), Error> {
         let mut ctx = Minimal::default();
-        ctx.register_op("noop", OpConstructor(new));
-        let op = ctx.op("noop xy_in=us-ft z_in=us-ft")?;
+        ctx.register_op("senmerc", OpConstructor(new));
+        let definition = "senmerc";
+        let op = ctx.op(definition)?;
 
-        // EPSG:1134 - 3 parameter, ED50/WGS84, s = sqrt(27) m
-        let mut operands = [GDA94];
+        let ldn = [Coor4D::geo(51.505, -0.09, 30., 0.)];
+
+        let projected = [Coor4D::raw(
+            -10_018.754_171_394_621,
+            6_711_113.243_704_713,
+            48.196_925_789_267_21,
+            0.,
+        )];
 
         // Forward
+        let mut operands = ldn;
         ctx.apply(op, Fwd, &mut operands)?;
-        assert_eq!(operands[0], GDA94);
+        for i in 0..operands.len() {
+            assert_float_eq!(operands[i].0, projected[i].0, abs_all <= 1e-8);
+        }
 
-        // Inverse + roundtrip
+        // Roundtrip
         ctx.apply(op, Inv, &mut operands)?;
-        assert_eq!(operands[0], GDA94);
+        for i in 0..operands.len() {
+            assert_float_eq!(operands[i].0, ldn[i].0, abs_all <= 2e-9);
+        }
+
         Ok(())
     }
 }