Support axisswap

With the recent merge of @Rennzie's implementation of the PROJ
unitconvert operator, we just need the (somewhat simpler) axisswap
to support PROJ's way of doing what `adapt` otherwise does for RG.

ruminations/002-rumination.md

@@ -22,6 +22,7 @@ $ echo 553036. -124509 | kp "dms:in | geo:out"
 - [Prologue](#prologue)
 - [A brief `kp` HOWTO](#a-brief-kp-howto)
 - [`adapt`](#operator-adapt): The order-and-unit adaptor
+- [`axisswap`](#operator-axisswap): The axis order adaptor
 - [`cart`](#operator-cart): The geographical-to-cartesian converter
 - [`curvature`](#operator-curvature): Radii of curvature
 - [`deformation`](#operator-deformation): Kinematic datum shift using a
@@ -138,6 +139,46 @@ where `geo:out` could be defined as `geo:in inv`.
 
 ---
 
+### Operator `axisswap`
+
+**Purpose:** Swap the order of coordinate elements in a coordinate tuple
+
+**Description:** In the `axisswap` model, the coordinate axes are numbered 1,2,3,4 and the axis swapping process is specified through the `order` argument, by providing a comma separated list of the reorganized order e.g.:
+
+```txt
+order=2,1,3,4
+```
+
+for swapping the first two axes.
+
+Axis indices may be prefixed by a minus sign, `-` to indicate a 180 degree swapping of the axis in question:
+
+```txt
+order=2,-1,3,4
+```
+
+which will make the second axis of the output equal to the negative of the first axis of the input.
+
+Postfix nonconsequential axis indices may be left out so:
+
+```txt
+order=2,-1
+```
+
+will give the same result as the previous example.
+
+**Usage:** Typically, `axisswap` (like `adapt` and `unitconvert`) is used in one or both ends of a pipeline, to match data between the RG internal representation and the requirements of the external coordinate representation:
+
+```txt
+axisswap order=2,1 | utm zone=32 | axisswap order=2,1
+```
+
+**Note:** This is an attempt to replicate Kristian Evers' PROJ operator of the [same name](https://proj.org/en/9.3/operations/conversions/axisswap.html), and any discrepancies should, as a general rule, be interpreted as errors in this implementation. Exceptions to this rule are all functionality related to PROJ's continued (but deprecated and undocumented) support of the classsical PROJ.4 syntax `axis=enu`, etc.
+
+**See also:** The documentation for the corresponding [PROJ operator](https://proj.org/en/9.3/operations/conversions/axisswap.html)
+
+---
+
 ### Operator `cart`
 
 **Purpose:** Convert from geographic coordinates + ellipsoidal height to geocentric cartesian coordinates
@@ -878,7 +919,7 @@ unitconvert xy_in=deg xy_out=rad
 ```
 
 **See also:** [PROJ documentation](https://proj.org/en/9.2/operations/conversions/unitconvert.html): *Unit Conversion*.
-A noticeable difference from PROJ is that time unit conversions are not yet supported. 
+A noticeable difference from PROJ is that time unit conversions are not yet supported.
 
 ---
 

src/inner_op/axisswap.rs

@@ -0,0 +1,219 @@
+use super::*;
+
+// ----- F O R W A R D -----------------------------------------------------------------
+
+fn fwd(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
+    let n = operands.len();
+
+    // We default to order=1,2,3,4, so if order is not given, we are done already
+    let Ok(order) = op.params.series("order") else {
+        return n;
+    };
+
+    let dimensionality = order.len();
+
+    let mut pos = [0_usize, 1, 2, 3];
+    let mut sgn = [1., 1., 1., 1.];
+    for (index, value) in order.iter().enumerate() {
+        pos[index] = (value.abs() - 1.) as usize;
+        sgn[index] = 1_f64.copysign(*value);
+    }
+
+    let mut successes = 0_usize;
+    for i in 0..n {
+        let inp = operands.get_coord(i);
+        let mut out = inp;
+        for index in 0..dimensionality {
+            out[index] = inp[pos[index]] * sgn[index];
+        }
+        operands.set_coord(i, &out);
+        successes += 1;
+    }
+
+    successes
+}
+
+// ----- I N V E R S E -----------------------------------------------------------------
+
+fn inv(op: &Op, _ctx: &dyn Context, operands: &mut dyn CoordinateSet) -> usize {
+    let n = operands.len();
+
+    // We default to order=1,2,3,4, so if order is not given, we are done already
+    let Ok(order) = op.params.series("order") else {
+        return n;
+    };
+
+    let dimensionality = order.len();
+
+    let mut pos = [0_usize, 1, 2, 3];
+    let mut sgn = [1., 1., 1., 1.];
+    for (index, value) in order.iter().enumerate() {
+        pos[index] = (value.abs() - 1.) as usize;
+        sgn[index] = 1_f64.copysign(*value);
+    }
+
+    let mut successes = 0_usize;
+    for i in 0..n {
+        let inp = operands.get_coord(i);
+        let mut out = inp;
+        for index in 0..dimensionality {
+            out[pos[index]] = inp[index] * sgn[index];
+        }
+        operands.set_coord(i, &out);
+        successes += 1;
+    }
+
+    successes
+}
+
+// ----- C O N S T R U C T O R ---------------------------------------------------------
+
+// Example...
+#[rustfmt::skip]
+pub const GAMUT: [OpParameter; 2] = [
+    OpParameter::Flag { key: "inv" },
+    OpParameter::Series { key: "order", default: Some("1,2,3,4") },
+];
+
+pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result<Op, Error> {
+    let op = Op::plain(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT, ctx)?;
+
+    // We default to order=1,2,3,4, so if order is not given, all is OK
+    let Ok(order) = op.params.series("order") else {
+        return Ok(op);
+    };
+
+    if order.len() > 4 {
+        return Err(Error::BadParam(
+            "order".to_string(),
+            "More than 4 indices given".to_string(),
+        ));
+    }
+
+    // While the Series type returns a Vec<f64>, the elements must be convertible to i64
+    // and further to (x.abs() as usize) for use as array indices
+    for &o in order {
+        let i = o as i64;
+        if (i as f64) != o || i == 0 || (i.unsigned_abs() as usize) > order.len() {
+            return Err(Error::BadParam("order".to_string(), o.to_string()));
+        }
+    }
+
+    // PROJ does not allow duplicate axes, presumably for a well considered reason,
+    // so neither do we
+    for o in 1_u64..5 {
+        if order.iter().filter(|x| (x.abs() as u64) == o).count() > 1 {
+            return Err(Error::BadParam(
+                "order".to_string(),
+                "duplicate axis specified".to_string(),
+            ));
+        }
+    }
+
+    Ok(op)
+}
+
+// ----- T E S T S ---------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn four_dim() -> Result<(), Error> {
+        let mut ctx = Minimal::default();
+        let op = ctx.op("axisswap order=2,1,-3,-4")?;
+
+        let mut operands = [Coor4D([1., 2., 3., 4.])];
+
+        // Forward
+        ctx.apply(op, Fwd, &mut operands)?;
+        assert_eq!(operands[0][0], 2.);
+        assert_eq!(operands[0][1], 1.);
+        assert_eq!(operands[0][2], -3.);
+        assert_eq!(operands[0][3], -4.);
+
+        // Inverse + roundtrip
+        ctx.apply(op, Inv, &mut operands)?;
+        assert_eq!(operands[0][0], 1.);
+        assert_eq!(operands[0][1], 2.);
+        assert_eq!(operands[0][2], 3.);
+        assert_eq!(operands[0][3], 4.);
+
+        Ok(())
+    }
+
+    #[test]
+    fn default_order() -> Result<(), Error> {
+        let mut ctx = Minimal::default();
+        let op = ctx.op("axisswap")?;
+
+        let mut operands = [Coor4D([1., 2., 3., 4.])];
+
+        // Forward
+        ctx.apply(op, Fwd, &mut operands)?;
+        assert_eq!(operands[0][0], 1.);
+        assert_eq!(operands[0][1], 2.);
+        assert_eq!(operands[0][2], 3.);
+        assert_eq!(operands[0][3], 4.);
+
+        // Inverse + roundtrip
+        ctx.apply(op, Inv, &mut operands)?;
+        assert_eq!(operands[0][0], 1.);
+        assert_eq!(operands[0][1], 2.);
+        assert_eq!(operands[0][2], 3.);
+        assert_eq!(operands[0][3], 4.);
+
+        Ok(())
+    }
+
+    #[test]
+    fn two_dim() -> Result<(), Error> {
+        let mut ctx = Minimal::default();
+
+        let op = ctx.op("axisswap order=2,-1")?;
+        let mut operands = [Coor4D([1., 2., 3., 4.])];
+
+        // Forward
+        ctx.apply(op, Fwd, &mut operands)?;
+        assert_eq!(operands[0][0], 2.);
+        assert_eq!(operands[0][1], -1.);
+        assert_eq!(operands[0][2], 3.);
+        assert_eq!(operands[0][3], 4.);
+
+        // Inverse + roundtrip
+        ctx.apply(op, Inv, &mut operands)?;
+        assert_eq!(operands[0][0], 1.);
+        assert_eq!(operands[0][1], 2.);
+        assert_eq!(operands[0][2], 3.);
+        assert_eq!(operands[0][3], 4.);
+        Ok(())
+    }
+
+    #[test]
+    fn bad_parameters() -> Result<(), Error> {
+        let mut ctx = Minimal::default();
+
+        // Too many indices
+        let op = ctx.op("axisswap order=4,4,4,2,-1");
+        assert!(matches!(op, Err(Error::BadParam(_, _))));
+
+        // Repeated indices
+        let op = ctx.op("axisswap order=4,-4,2,-1");
+        assert!(matches!(op, Err(Error::BadParam(_, _))));
+
+        // Index exceeding dimensionality
+        let op = ctx.op("axisswap order=2,3");
+        assert!(matches!(op, Err(Error::BadParam(_, _))));
+
+        // Missing indices ('order' becomes a flag)
+        let op = ctx.op("axisswap order");
+        assert!(matches!(op, Err(Error::BadParam(_, _))));
+
+        // Missing all args: axisswap succeeds and becomes a no-op
+        let op = ctx.op("axisswap");
+        assert!(op.is_ok());
+
+        Ok(())
+    }
+}

src/inner_op/mod.rs

@@ -7,6 +7,7 @@ use crate::authoring::*;
 
 mod adapt;
 mod addone;
+mod axisswap;
 mod btmerc;
 mod cart;
 mod curvature;
@@ -30,9 +31,10 @@ mod units;
 mod webmerc;
 
 #[rustfmt::skip]
-const BUILTIN_OPERATORS: [(&str, OpConstructor); 31] = [
+const BUILTIN_OPERATORS: [(&str, OpConstructor); 32] = [
     ("adapt",        OpConstructor(adapt::new)),
     ("addone",       OpConstructor(addone::new)),
+    ("axisswap",     OpConstructor(axisswap::new)),
     ("btmerc",       OpConstructor(btmerc::new)),
     ("butm",         OpConstructor(btmerc::utm)),
     ("cart",         OpConstructor(cart::new)),
@@ -49,17 +51,19 @@ const BUILTIN_OPERATORS: [(&str, OpConstructor); 31] = [
     ("merc",         OpConstructor(merc::new)),
     ("webmerc",      OpConstructor(webmerc::new)),
     ("molodensky",   OpConstructor(molodensky::new)),
-    ("noop",         OpConstructor(noop::new)),
     ("omerc",        OpConstructor(omerc::new)),
     ("somerc",       OpConstructor(somerc::new)),
     ("tmerc",        OpConstructor(tmerc::new)),
+    ("unitconvert",  OpConstructor(unitconvert::new)),
     ("utm",          OpConstructor(tmerc::utm)),
-    ("unitconvert", OpConstructor(unitconvert::new)),
+
+    // Pipeline handlers
     ("pipeline",     OpConstructor(pipeline::new)),
     ("pop",          OpConstructor(pipeline::pop)),
     ("push",         OpConstructor(pipeline::push)),
 
     // Some commonly used noop-aliases
+    ("noop",         OpConstructor(noop::new)),
     ("longlat",      OpConstructor(noop::new)),
     ("latlon",       OpConstructor(noop::new)),
     ("latlong",      OpConstructor(noop::new)),