gistools/proj/project/
gstmerc.rs

1use crate::proj::{CoordinateStep, Proj, ProjectCoordinates, TransformCoordinates, phi2, tsfn};
2use alloc::rc::Rc;
3use core::cell::RefCell;
4use libm::{asin, atan, cos, cosh, exp, log, pow, sin, sinh, sqrt};
5
6/// # Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)
7///
8/// **Classification**: Conformal
9///
10/// **Available forms**: Forward and inverse, spherical projection
11///
12/// **Defined area**: Global
13///
14/// **Alias**: gstmerc
15///
16/// **Domain**: 2D
17///
18/// **Input type**: Geodetic coordinates
19///
20/// **Output type**: Projected coordinates
21///
22/// ## Projection String
23/// ```ini
24/// +proj=gstmerc
25/// ```
26///
27/// ## Optional Parameters
28/// - `+k_0=<value>`: Scale factor at the central meridian.
29/// - `+lon_0=<value>`: Longitude of the central meridian.
30/// - `+lat_0=<value>`: Latitude of origin.
31/// - `+ellps=<value>`: Ellipsoid name (e.g., GRS80, WGS84).
32/// - `+R=<value>`: Radius of the sphere (used in spherical projections).
33/// - `+x_0=<value>`: False easting.
34/// - `+y_0=<value>`: False northing.
35///
36/// ## Usage Example
37/// ```bash
38/// echo 12 55 | proj +proj=gstmerc +ellps=WGS84
39/// echo 12 55 | proj +proj=gstmerc +k_0=1 +lon_0=0 +x_0=500000 +y_0=0
40/// ```
41///
42/// ![Gauss-Schreiber Transverse Mercator](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/gstmerc.png?raw=true)
43#[derive(Debug, Default, Clone, PartialEq)]
44pub struct GstmercData {
45    lamc: f64,
46    phic: f64,
47    c: f64,
48    n1: f64,
49    n2: f64,
50    xs: f64,
51    ys: f64,
52}
53
54/// Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion) Projection
55#[derive(Debug, Clone, PartialEq)]
56pub struct GaussSchreiberTransverseMercatorProjection {
57    proj: Rc<RefCell<Proj>>,
58    store: RefCell<GstmercData>,
59}
60impl ProjectCoordinates for GaussSchreiberTransverseMercatorProjection {
61    fn code(&self) -> i64 {
62        -1
63    }
64    fn name(&self) -> &'static str {
65        "Gauss-Schreiber Transverse Mercator"
66    }
67    fn names() -> &'static [&'static str] {
68        &[
69            "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)",
70            "Gauss-Schreiber Transverse Mercator",
71            "gstmerc",
72        ]
73    }
74}
75impl CoordinateStep for GaussSchreiberTransverseMercatorProjection {
76    fn new(proj: Rc<RefCell<Proj>>) -> Self {
77        let mut store = GstmercData::default();
78        {
79            let proj = proj.borrow();
80            store.lamc = proj.lam0;
81            store.n1 = sqrt(1. + proj.es * pow(cos(proj.phi0), 4.0) / (1. - proj.es));
82            store.phic = asin(sin(proj.phi0) / store.n1);
83            store.c = log(tsfn(-store.phic, -sin(proj.phi0) / store.n1, 0.0))
84                - store.n1 * log(tsfn(-proj.phi0, -sin(proj.phi0), proj.e));
85            store.n2 = proj.k0 * proj.a * sqrt(1. - proj.es)
86                / (1. - proj.es * sin(proj.phi0) * sin(proj.phi0));
87            store.xs = 0.;
88            store.ys = -store.n2 * store.phic;
89        }
90        GaussSchreiberTransverseMercatorProjection { proj, store: store.into() }
91    }
92    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
93        gstmerc_s_forward(&self.store.borrow(), &self.proj.borrow(), p);
94    }
95    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
96        gstmerc_s_inverse(&self.store.borrow(), &self.proj.borrow(), p);
97    }
98}
99
100/// Gauss-Schreiber Transverse Mercator Spheroidal forward project
101pub fn gstmerc_s_forward<P: TransformCoordinates>(gstmerc: &GstmercData, proj: &Proj, p: &mut P) {
102    let l = gstmerc.n1 * p.lam();
103    let ls = gstmerc.c + gstmerc.n1 * log(tsfn(-p.phi(), -sin(p.phi()), proj.e));
104    let sin_ls1 = sin(l) / cosh(ls);
105    let ls1 = log(tsfn(-asin(sin_ls1), -sin_ls1, 0.0));
106    p.set_x((gstmerc.xs + gstmerc.n2 * ls1) * proj.ra);
107    p.set_y((gstmerc.ys + gstmerc.n2 * atan(sinh(ls) / cos(l))) * proj.ra);
108}
109
110/// Gauss-Schreiber Transverse Mercator Spheroidal inverse project
111pub fn gstmerc_s_inverse<P: TransformCoordinates>(gstmerc: &GstmercData, proj: &Proj, p: &mut P) {
112    let l = atan(
113        sinh((p.x() * proj.a - gstmerc.xs) / gstmerc.n2)
114            / cos((p.y() * proj.a - gstmerc.ys) / gstmerc.n2),
115    );
116    let sin_c = sin((p.y() * proj.a - gstmerc.ys) / gstmerc.n2)
117        / cosh((p.x() * proj.a - gstmerc.xs) / gstmerc.n2);
118    let lc = log(tsfn(-asin(sin_c), -sin_c, 0.0));
119    p.set_lam(l / gstmerc.n1);
120    p.set_phi(-phi2(exp((lc - gstmerc.c) / gstmerc.n1), proj.e));
121}
gistools/proj/project/gstmerc.rs

gistools/proj/project/
gstmerc.rs
