miniproj_ops/ops/
popvis_pseudo_mercator.rs

1//This file is licensed under EUPL v1.2 as part of the Digital Earth Viewer
2
3use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};
4
5use crate::{
6    ellipsoid::Ellipsoid, traits::GetterContstruct, DbContstruct, Projection, PseudoSerialize,
7};
8
9#[derive(Copy, Clone, Debug)]
10pub struct PopVisPseudoMercatorParams {
11    /// longitude of natural origin
12    lon_orig: f64,
13    /// latitude of natural origin
14    lat_orig: f64,
15    /// false easting
16    false_e: f64,
17    /// false northing
18    false_n: f64,
19}
20
21impl PopVisPseudoMercatorParams {
22    pub const fn new(lon_orig: f64, lat_orig: f64, false_e: f64, false_n: f64) -> Self {
23        Self {
24            lat_orig,
25            lon_orig,
26            false_e,
27            false_n,
28        }
29    }
30
31    /// Get longitude of natural origin, radians.
32    pub fn lon_orig(&self) -> f64 {
33        self.lon_orig
34    }
35
36    /// Get latitude of natural origin, radians.
37    pub fn lat_orig(&self) -> f64 {
38        self.lat_orig
39    }
40
41    /// Get false easting.
42    pub fn false_e(&self) -> f64 {
43        self.false_e
44    }
45
46    /// Get false northing.
47    pub fn false_n(&self) -> f64 {
48        self.false_n
49    }
50}
51
52/// Transverse Mercator coordinate operation (EPSG:9807).
53#[allow(non_snake_case)]
54#[derive(Copy, Clone, Debug)]
55pub struct PopVisPseudoMercatorProjection {
56    pub false_e: f64,
57    pub false_n: f64,
58    pub ellipsoid_a: f64,
59    pub lon_orig: f64,
60}
61
62impl PopVisPseudoMercatorProjection {
63    #[allow(non_snake_case)]
64    pub fn new(ell: &Ellipsoid, params: &PopVisPseudoMercatorParams) -> Self {
65        Self {
66            ellipsoid_a: ell.a(),
67            lon_orig: params.lon_orig(),
68            false_e: params.false_e(),
69            false_n: params.false_n(),
70        }
71    }
72}
73
74impl Projection for PopVisPseudoMercatorProjection {
75    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
76    /// longitude & latitude in radians
77    #[allow(non_snake_case)]
78    fn rad_to_projected(&self, longitude: f64, latitude: f64) -> (f64, f64) {
79        (
80            self.false_e + self.ellipsoid_a * (longitude - self.lon_orig),
81            self.false_n + self.ellipsoid_a * (FRAC_PI_4 + latitude / 2f64).tan().ln(),
82        )
83    }
84
85    /// as per IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – March 2020
86    /// longitude & latitude in radians
87    #[allow(non_snake_case)]
88    fn projected_to_rad(&self, easting: f64, northing: f64) -> (f64, f64) {
89        let D = (self.false_n - northing) / self.ellipsoid_a;
90        (
91            ((easting - self.false_e) / self.ellipsoid_a) + self.lon_orig,
92            FRAC_PI_2 - 2.0 * D.exp().atan(),
93        )
94    }
95}
96
97impl PseudoSerialize for PopVisPseudoMercatorProjection {
98    fn to_constructed(&self) -> String {
99        format!(
100            r"PopVisPseudoMercatorProjection{{
101    ellipsoid_a: {}f64,
102    lon_orig: {}f64,
103    false_e: {}f64,
104    false_n: {}f64,
105}}",
106            self.ellipsoid_a, self.lon_orig, self.false_e, self.false_n,
107        )
108    }
109}
110
111impl DbContstruct for PopVisPseudoMercatorProjection {
112    fn from_database_params(params: &[(u32, f64)], ellipsoid: &Ellipsoid) -> Self {
113        let params = PopVisPseudoMercatorParams::new(
114            params
115                .iter()
116                .find_map(|(c, v)| if *c == 8802 { Some(*v) } else { None })
117                .unwrap(),
118            params
119                .iter()
120                .find_map(|(c, v)| if *c == 8801 { Some(*v) } else { None })
121                .unwrap(),
122            params
123                .iter()
124                .find_map(|(c, v)| if *c == 8806 { Some(*v) } else { None })
125                .unwrap(),
126            params
127                .iter()
128                .find_map(|(c, v)| if *c == 8807 { Some(*v) } else { None })
129                .unwrap(),
130        );
131        Self::new(ellipsoid, &params)
132    }
133}
134
135impl GetterContstruct for PopVisPseudoMercatorProjection {
136    fn with_db_getter<G>(mut getter: G, ellipsoid: &Ellipsoid) -> Option<Self>
137    where
138        G: FnMut(u32) -> Option<f64>,
139    {
140        let params = PopVisPseudoMercatorParams::new(
141            getter(8802)?,
142            getter(8801)?,
143            getter(8806)?,
144            getter(8807)?,
145        );
146        Some(Self::new(ellipsoid, &params))
147    }
148}
149
150pub fn direct_projection(params: &[(u32, f64)], ell: Ellipsoid) -> String {
151    PopVisPseudoMercatorProjection::from_database_params(params, &ell).to_constructed()
152}
153#[cfg(test)]
154mod tests {
155
156    use crate::ellipsoid::Ellipsoid;
157    use crate::traits::*;
158
159    use super::PopVisPseudoMercatorParams;
160    use super::PopVisPseudoMercatorProjection;
161
162    #[test]
163    fn popvis_mercator_consistency() {
164        let ell = Ellipsoid::from_a_f_inv(6378137.0, 298.257223563);
165        let params =
166            PopVisPseudoMercatorParams::new(0.0f64.to_radians(), 0.0f64.to_radians(), 0.0, 0.0);
167
168        let projection = PopVisPseudoMercatorProjection::new(&ell, &params);
169        let easting_goal = -11169055.58;
170        let northing_goal = 2800000.00;
171        let (lon, lat) = projection.projected_to_deg(easting_goal, northing_goal);
172        let (easting, northing) = projection.deg_to_projected(lon, lat);
173
174        eprintln!("easting: {easting_goal} - {easting}");
175        eprintln!("northing: {northing_goal} - {northing}");
176
177        assert!((easting - easting_goal).abs() < 0.01);
178
179        assert!((northing - northing_goal).abs() < 0.01);
180    }
181}
miniproj_ops/ops/popvis_pseudo_mercator.rs

miniproj_ops/ops/
popvis_pseudo_mercator.rs
