From 05a2da231bcc92f763ef66d64296c672c9e7aef7 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Fri, 8 Nov 2024 12:08:33 +0100 Subject: [PATCH] Support conversion between permanent tide systems --- README.md | 2 +- ruminations/002-rumination.md | 61 +++++++++++++ src/bibliography/mod.rs | 20 +++++ src/inner_op/mod.rs | 4 +- src/inner_op/permtide.rs | 159 ++++++++++++++++++++++++++++++++++ 5 files changed, 244 insertions(+), 2 deletions(-) create mode 100644 src/inner_op/permtide.rs diff --git a/README.md b/README.md index 2a24e0e..571e6c1 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ Rust Geodesy provides a number of **features** to support a number of **objectiv The most important **features** are -- a set of more than 20 geodetic **transformation primitives** +- a set of more than 30 geodetic **transformation primitives** - a set of more than 40 primitives for **operations on the ellipsoid** - a means for **composing** these primitives into more complex operations. diff --git a/ruminations/002-rumination.md b/ruminations/002-rumination.md index 022dcab..a47bc80 100644 --- a/ruminations/002-rumination.md +++ b/ruminations/002-rumination.md @@ -42,6 +42,8 @@ $ echo 553036. -124509 | kp "dms:in | geo:out" - [`molodensky`](#operator-molodensky): The full and abridged Molodensky transformations - [`noop`](#operator-noop): The no-operation - [`omerc`](#operator-omerc): The oblique Mercator projection +- [`permtide`](#operator-permtide): + Convert geoid undulations between different permanent tide systems - [`pop`](#operator-pop): Pop a dimension from the stack into the operands - [`push`](#operator-push): Push a dimension from the operands onto the stack - [`stack`](#operator-stack): Push/pop/swap dimensions from the operands onto the stack @@ -861,6 +863,65 @@ and RG does not support PROJ's "indirectly given azimuth" case. --- +### Operator `permtide` + +**Purpose:** Convert geoid undulations between different permanent tide systems + +**Description:** + +Since the orbits of the sun and the moon (as observed from the earth) +are concentrated at lower latitudes, the mean tidal effect of these +celestial bodies do not vanish, but results in a non-zero mean potential. + +Hence, if we compute the long term mean of a series of repeated +levellings between two fixed points at different latitudes, then +we will eliminate the time-varying parts of the lunar and solar tidal +potentials, but the non-vanishing long term mean will still blend into +our attempt to measure the geopotential difference between the two +points. This is known as *the mean tide* case. + +If correcting for the mean as well, we formally obtain a more pure +*geo*-potential. This is known as the *zero-tide* case, and is +the equivalent to formally moving all external gravitating masses to +infinity. + +But since the permanent tide not only influences the potential, but +also the shape of the earth's crust, there is a secondary effect from +the external gravitating bodies due to the deformation. When we +formally remove that as well, we are left with what is known as the +*non-tidal* or *tide free* case. + +In height systems, we must discern between *mean tide*, +*zero tide*, and *tide free* conventions, and adapt the corresponding +geoid model to fit with the convention. Hence, this operator uses +geoid-centric terminology and sign conventions. + + + + +| Argument | Description | +|----------|-------------| +| `inv` | swap forward and inverse operations | +| `ellps=name` | Use ellipsoid `name` for the conversion | +| `k` | zero frequency Love number. Defaults to $0.3$ | +| `from=system` | Convert from either `mean`, `zero` or `free` tide system | +| `to=system` | Convert to either `mean`, `zero` or `free` tide system | + +**Example**: Convert a geoid model using the zero-tide convention to a +corresponding model using the mean-tide convention + +```geodesy +permtide from=zero to=mean ellps=GRS80 +``` + +**See also:** +[Martin Losch and Verena Seufer, 2003:](https://mitgcm.org/~mlosch/geoidcookbook.pdf) +*How to Compute Geoid Undulations (Geoid Height Relative to a Given +Reference Ellipsoid) from Spherical Harmonic Coefficients for +Satellite Altimetry Applications* + +--- + ### Operator `pop` **DEPRECATED!** Use [`stack`](#operator-stack) diff --git a/src/bibliography/mod.rs b/src/bibliography/mod.rs index 261f058..b25babd 100644 --- a/src/bibliography/mod.rs +++ b/src/bibliography/mod.rs @@ -103,6 +103,20 @@ pub enum Bibliography { /// [pdf](https://gfzpublic.gfz-potsdam.de/rest/items/item_8827_5/component/file_130038/content). Kru12, + /// Martin Losch and Verena Seufer, 2003: + /// *How to Compute Geoid Undulations (Geoid Height Relative to a Given Reference Ellipsoid) + /// from Spherical Harmonic Coefficients for Satellite Altimetry Applications* + /// Unpublished manuscript, 10 pp. + /// [html](https://mitgcm.org/~mlosch/geoidcookbook/geoidcookbook.html) + /// [pdf](https://mitgcm.org/~mlosch/geoidcookbook.pdf). + Los03, + + /// Gérard Petit and Brian Luzum (eds), 2010: + /// *IERS conventions (2010)*. + /// IERS technical note 36. + /// Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main, 179 pp. + /// [pdf](https://iers-conventions.obspm.fr/content/tn36.pdf) + /// Knud Poder and Karsten Engsager, 1998. /// *Some Conformal Mappings and Transformations for Geodesy and Topographic Cartography*. /// Copenhagen, Denmark: Geodetic Division, KMS, @@ -116,6 +130,12 @@ pub enum Bibliography { /// [DOI](https://doi.org/10.1080/00396265.2016.1191748). Ruf16, + /// Dru A. Smith, 1998: + /// *There is no such thing as "The" EGM96 geoid: Subtle points on the use of a global geopotential model*, + /// IGeS Bulletin No. 8, International Geoid Service, Milan, Italy, pp. 17-28, + /// [URL](https://www.ngs.noaa.gov/PUBS_LIB/EGM96_GEOID_PAPER/egm96_geoid_paper.html). + Dru98, + /// T. Vincenty (1975) *Direct and Inverse Solutions of Geodesics on the Ellipsoid /// with application of nested equations*. /// Survey Review, 23(176): 88-93. diff --git a/src/inner_op/mod.rs b/src/inner_op/mod.rs index 1f768be..450cea6 100644 --- a/src/inner_op/mod.rs +++ b/src/inner_op/mod.rs @@ -25,6 +25,7 @@ mod merc; mod molodensky; mod noop; mod omerc; +mod permtide; pub(crate) mod pipeline; // Needed by Op for instantiation mod pushpop; mod somerc; @@ -35,7 +36,7 @@ mod units; mod webmerc; #[rustfmt::skip] -const BUILTIN_OPERATORS: [(&str, OpConstructor); 35] = [ +const BUILTIN_OPERATORS: [(&str, OpConstructor); 36] = [ ("adapt", OpConstructor(adapt::new)), ("addone", OpConstructor(addone::new)), ("axisswap", OpConstructor(axisswap::new)), @@ -58,6 +59,7 @@ const BUILTIN_OPERATORS: [(&str, OpConstructor); 35] = [ ("webmerc", OpConstructor(webmerc::new)), ("molodensky", OpConstructor(molodensky::new)), ("omerc", OpConstructor(omerc::new)), + ("permtide", OpConstructor(permtide::new)), ("somerc", OpConstructor(somerc::new)), ("tmerc", OpConstructor(tmerc::new)), ("unitconvert", OpConstructor(unitconvert::new)), diff --git a/src/inner_op/permtide.rs b/src/inner_op/permtide.rs new file mode 100644 index 0000000..a2d2fcc --- /dev/null +++ b/src/inner_op/permtide.rs @@ -0,0 +1,159 @@ +/// Permanent tide systems +use crate::authoring::*; + +// ----- F O R W A R D ----------------------------------------------------------------- + +fn fwd(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 Ok(coefficient) = op.params.real("coefficient") else { + return successes; + }; + + for i in 0..n { + let mut coord = operands.get_coord(i); + let phibar = ellps.latitude_geographic_to_geocentric(coord[1]); + let s = phibar.sin(); + coord[2] += coefficient * (-0.198) * (1.5 * s * s - 0.5); + 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 { + let mut successes = 0_usize; + let n = operands.len(); + let ellps = op.params.ellps(0); + let Ok(coefficient) = op.params.real("coefficient") else { + return successes; + }; + + for i in 0..n { + let mut coord = operands.get_coord(i); + let phibar = ellps.latitude_geographic_to_geocentric(coord[1]); + let s = phibar.sin(); + coord[2] -= coefficient * (-0.198) * (1.5 * s * s - 0.5); + operands.set_coord(i, &coord); + successes += 1; + } + + successes +} + +// ----- C O N S T R U C T O R --------------------------------------------------------- + +// Example... +#[rustfmt::skip] +pub const GAMUT: [OpParameter; 5] = [ + OpParameter::Flag { key: "inv" }, + OpParameter::Real { key: "k", default: Some(0.3) }, + OpParameter::Text { key: "ellps", default: Some("GRS80") }, + OpParameter::Text { key: "from", default: None }, + OpParameter::Text { key: "to", default: None } +]; + +pub fn new(parameters: &RawParameters, ctx: &dyn Context) -> Result { + let mut op = Op::plain(parameters, InnerOp(fwd), Some(InnerOp(inv)), &GAMUT, ctx)?; + let k = op.params.real("k")?; + + let Ok(to) = op.params.text("to") else { + return Err(Error::MissingParam( + "permtide: must specify 'to=' as exactly one of {'mean', 'zero', 'free'}".to_string(), + )); + }; + let Ok(from) = op.params.text("from") else { + return Err(Error::MissingParam( + "permtide: must specify 'from=' as exactly one of {'mean', 'zero', 'free'}".to_string(), + )); + }; + + let coefficient = match (to.as_str(), from.as_str()) { + ("mean", "mean") => 0.0, + ("mean", "zero") => 1.0, + ("mean", "free") => 1.0 + k, + ("zero", "zero") => 0.0, + ("zero", "mean") => -1.0, + ("zero", "free") => k, + ("free", "free") => 0.0, + ("free", "mean") => -(1.0 + k), + ("free", "zero") => -k, + _ => f64::NAN, + }; + + if coefficient.is_nan() { + return Err(Error::BadParam( + "'to=' or 'from='".to_string(), + "must be one of {'mean', 'zero', 'free'}".to_string(), + )); + } + + op.params.real.insert("coefficient", coefficient); + Ok(op) +} + +// ----- T E S T S --------------------------------------------------------------------- + +#[cfg(test)] +mod tests { + use super::*; + use float_eq::assert_float_eq; + + #[test] + fn permtide() -> Result<(), Error> { + let mut ctx = Minimal::default(); + + // A test point near Copenhagen + let pnt = Coor4D::geo(55., 12., 0., 0.); + + // Mean -> zero + let op = ctx.op("permtide from=mean to=zero ellps=GRS80")?; + let mut operands = [pnt]; + ctx.apply(op, Fwd, &mut operands)?; + assert_float_eq!(operands[0][2], 0.099407199, abs_all <= 1e-8); + ctx.apply(op, Inv, &mut operands)?; + assert_float_eq!(operands[0][2], pnt[2], abs_all <= 1e-12); + + // Mean -> free + let op = ctx.op("permtide from=mean to=free ellps=GRS80")?; + let mut operands = [pnt]; + ctx.apply(op, Fwd, &mut operands)?; + assert_float_eq!(operands[0][2], 0.1292293587824579, abs_all <= 1e-8); + ctx.apply(op, Inv, &mut operands)?; + assert_float_eq!(operands[0][2], pnt[2], abs_all <= 1e-12); + + // Inversion + let fwd_op = ctx.op("permtide from=mean to=zero ellps=GRS80")?; + let inv_op = ctx.op("permtide from=zero to=mean ellps=GRS80 inv")?; + let mut operands = [pnt]; + ctx.apply(fwd_op, Fwd, &mut operands)?; + let fwd_h = operands[0][2]; + + let mut operands = [pnt]; + ctx.apply(inv_op, Fwd, &mut operands)?; + let inv_h = operands[0][2]; + assert_float_eq!(fwd_h, inv_h, abs_all <= 1e-20); + + // Bad tide system names + assert!(matches!( + ctx.op("permtide from=cheese to=zero ellps=GRS80"), + Err(Error::BadParam(_, _)) + )); + + // Missing tide system names + assert!(matches!(ctx.op("permtide"), Err(Error::MissingParam(_)))); + assert!(matches!( + ctx.op("permtide to=zero ellps=GRS80"), + Err(Error::MissingParam(_)) + )); + assert!(matches!( + ctx.op("permtide from=mean ellps=GRS80"), + Err(Error::MissingParam(_)) + )); + Ok(()) + } +}