Skip to content

Commit

Permalink
Support conversion between permanent tide systems
Browse files Browse the repository at this point in the history
  • Loading branch information
busstoptaktik committed Nov 8, 2024
1 parent 867d70f commit 05a2da2
Show file tree
Hide file tree
Showing 5 changed files with 244 additions and 2 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
61 changes: 61 additions & 0 deletions ruminations/002-rumination.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
20 changes: 20 additions & 0 deletions src/bibliography/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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.
Expand Down
4 changes: 3 additions & 1 deletion src/inner_op/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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)),
Expand All @@ -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)),
Expand Down
159 changes: 159 additions & 0 deletions src/inner_op/permtide.rs
Original file line number Diff line number Diff line change
@@ -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<Op, Error> {
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(())
}
}

0 comments on commit 05a2da2

Please sign in to comment.