gistools/proj/project/
omerc.rs

1use crate::proj::{
2    ANGLE_RECTIFIED_TO_SKEW_GRID, AZIMUTH_PROJECTION_CENTRE, CoordinateStep,
3    HOTINE_OBLIQUE_MERCATOR_VARIANT_A, HOTINE_OBLIQUE_MERCATOR_VARIANT_B, LATITUDE_OF_FIRST_POINT,
4    LATITUDE_OF_SECOND_POINT, LONGITUDE_OF_FIRST_POINT, LONGITUDE_OF_PROJECTION_CENTRE,
5    LONGITUDE_OF_SECOND_POINT, NO_OFF, NO_ROTATION, NO_UOFF, Proj, ProjValue, ProjectCoordinates,
6    TransformCoordinates, aasin, adjlon, phi2, tsfn,
7};
8use alloc::rc::Rc;
9use core::{
10    cell::RefCell,
11    f64::consts::{FRAC_PI_2, FRAC_PI_4, PI, TAU},
12};
13use libm::{atan, atan2, cos, exp, fabs, log, pow, sin, sqrt, tan};
14
15/// Oblique Mercator Variables
16/// INTERNA
17#[derive(Debug, Default, Clone, PartialEq)]
18pub struct OmercData {
19    a: f64,
20    b: f64,
21    e: f64,
22    ab: f64,
23    ar_b: f64,
24    br_a: f64,
25    r_b: f64,
26    singam: f64,
27    cosgam: f64,
28    sinrot: f64,
29    cosrot: f64,
30    v_pole_n: f64,
31    v_pole_s: f64,
32    u_0: f64,
33    no_rot: bool,
34}
35
36const TOL: f64 = 1e-7;
37const EPS: f64 = 1e-10;
38
39/// Hotine Oblique Mercator (variant A) Projection
40///
41/// EPSG Codes Used by Hotine Oblique Mercator (variant A): 8811, 8812, 8813, 8814, 8815, 8806, 8807
42///
43/// See [`ObliqueMercatorProjection`] for full documentation.
44pub type HotineObliqueMercatorVariantAProjection =
45    ObliqueMercatorProjection<HOTINE_OBLIQUE_MERCATOR_VARIANT_A>;
46/// Hotine Oblique Mercator (variant B) Projection
47///
48/// EPSG Codes Used by Hotine Oblique Mercator (variant B): 8811, 8812, 8813, 8814, 8815, 8816, 8817
49///
50/// See [`ObliqueMercatorProjection`] for full documentation.
51pub type HotineObliqueMercatorVariantBProjection =
52    ObliqueMercatorProjection<HOTINE_OBLIQUE_MERCATOR_VARIANT_B>;
53
54/// # Oblique Mercator
55///
56/// **Classification**: Conformal cylindrical
57///
58/// **Available forms**: Forward and inverse, spherical and ellipsoidal
59///
60/// **Defined area**: Global, but reasonably accurate only within 15 degrees of the oblique central line
61///
62/// **Alias**: omerc
63///
64/// **Domain**: 2D
65///
66/// **Input type**: Geodetic coordinates
67///
68/// **Output type**: Projected coordinates
69///
70/// ## Projection String
71/// ```ini
72/// +proj=omerc +lat_1=45 +lat_2=55
73/// ```
74///
75/// ## Required Parameters
76/// - `+lat_1=<value>`: Latitude of the first point on the central line.
77/// - `+lat_2=<value>`: Latitude of the second point on the central line.
78///
79/// ## Optional Parameters
80/// - `+alpha=<value>`: Azimuth of the centerline clockwise from north at the center point of the line.
81/// - `+gamma=<value>`: Azimuth of the centerline clockwise from north of the rectified bearing of the centerline.
82/// - `+lonc=<value>`: Longitude of the projection center (overrides `+lon_0`).
83/// - `+lat_0=<value>`: Latitude of the projection center.
84/// - `+no_rot`: Disables rectification (historical reason).
85/// - `+no_off`: Disables origin offset to the center of projection.
86/// - `+k_0=<value>`: Scale factor at the central line.
87/// - `+x_0=<value>`: False easting.
88/// - `+y_0=<value>`: False northing.
89///
90/// ## Usage Example
91/// ```bash
92/// echo 12 55 | proj +proj=omerc +alpha=90 +ellps=GRS80
93/// echo 12 55 | proj +proj=omerc +alpha=0 +R=6400000
94/// echo 12 55 | proj +proj=omerc +lon_1=0 +lat_1=-1 +lon_2=0 +lat_2=0 +R=6400000
95/// echo 12 55 | proj +proj=tmerc +R=6400000
96/// echo 12 55 | proj +proj=omerc +lon_1=-1 +lat_1=1 +lon_2=0 +lat_2=0 +ellps=GRS80
97/// echo 10.536498003 56.229892362 | cs2cs +proj=longlat +ellps=GRS80 +to +proj=omerc +axis=wnu +lonc=9.46 +lat_0=56.13333333 +x_0=-266906.229 +y_0=189617.957 +k=0.9999537 +alpha=-0.76324 +gamma=0 +ellps=GRS80
98/// ```
99///
100/// ## Caveats
101/// The two-point method with no rectification is probably only marginally useful.
102///
103/// ![Oblique Mercator](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/omerc.png?raw=true)
104#[derive(Debug, Clone, PartialEq)]
105pub struct ObliqueMercatorProjection<const C: i64> {
106    proj: Rc<RefCell<Proj>>,
107    store: RefCell<OmercData>,
108}
109impl<const C: i64> ProjectCoordinates for ObliqueMercatorProjection<C> {
110    fn code(&self) -> i64 {
111        C
112    }
113    fn name(&self) -> &'static str {
114        "Oblique Mercator"
115    }
116    fn names() -> &'static [&'static str] {
117        &[
118            "Hotine_Oblique_Mercator",
119            "Hotine Oblique Mercator",
120            "Hotine_Oblique_Mercator_Azimuth_Natural_Origin",
121            "Hotine Oblique Mercator Azimuth Natural Origin",
122            "Hotine_Oblique_Mercator_Two_Point_Natural_Origin",
123            "Hotine Oblique Mercator Two Point Natural Origin",
124            "Hotine_Oblique_Mercator_Azimuth_Center",
125            "Hotine Oblique Mercator Azimuth Center",
126            "Hotine Oblique Mercator (variant A)",
127            "Hotine Oblique Mercator (variant B)",
128            "Oblique_Mercator",
129            "Oblique Mercator",
130            "omerc",
131        ]
132    }
133}
134impl<const C: i64> CoordinateStep for ObliqueMercatorProjection<C> {
135    fn new(proj: Rc<RefCell<Proj>>) -> Self {
136        let mut store = OmercData::default();
137        {
138            let proj = &mut proj.borrow_mut();
139
140            let mut no_off = false;
141            let mut lam1 = 0.;
142            let mut lam2 = 0.;
143            let mut phi1 = 0.;
144            let mut phi2 = 0.;
145            let gamma0;
146            let mut con;
147            let mut _f;
148            let _d;
149            let mut lamc = 0.;
150
151            store.no_rot = proj.params.get(&NO_ROTATION).unwrap_or(&ProjValue::default()).bool();
152            let mut alpha_c =
153                proj.params.get(&AZIMUTH_PROJECTION_CENTRE).unwrap_or(&ProjValue::default()).f64();
154            let alp = alpha_c != 0.;
155            let mut gamma = proj
156                .params
157                .get(&ANGLE_RECTIFIED_TO_SKEW_GRID)
158                .unwrap_or(&ProjValue::default())
159                .f64();
160            let gam = gamma != 0.;
161            if alp || gam {
162                lamc = proj
163                    .params
164                    .get(&LONGITUDE_OF_PROJECTION_CENTRE)
165                    .unwrap_or(&ProjValue::default())
166                    .f64();
167                // For libproj4 compatibility
168                no_off = proj.params.get(&NO_OFF).unwrap_or(&ProjValue::default()).bool()
169                    // for backward compatibility
170                    || proj.params.get(&NO_UOFF).unwrap_or(&ProjValue::default()).bool();
171            } else {
172                lam1 = proj
173                    .params
174                    .get(&LONGITUDE_OF_FIRST_POINT)
175                    .unwrap_or(&ProjValue::default())
176                    .f64();
177                phi1 = proj
178                    .params
179                    .get(&LATITUDE_OF_FIRST_POINT)
180                    .unwrap_or(&ProjValue::default())
181                    .f64();
182                lam2 = proj
183                    .params
184                    .get(&LONGITUDE_OF_SECOND_POINT)
185                    .unwrap_or(&ProjValue::default())
186                    .f64();
187                phi2 = proj
188                    .params
189                    .get(&LATITUDE_OF_SECOND_POINT)
190                    .unwrap_or(&ProjValue::default())
191                    .f64();
192                con = fabs(phi1);
193
194                if fabs(phi1) > FRAC_PI_2 - TOL {
195                    panic!("Invalid value for lat_1: |lat_1| should be < 90°");
196                }
197                if fabs(phi2) > FRAC_PI_2 - TOL {
198                    panic!("Invalid value for lat_2: |lat_2| should be < 90°");
199                }
200                if fabs(phi1 - phi2) <= TOL {
201                    panic!("Invalid value for lat_1/lat_2: lat_1 should be different from lat_2");
202                }
203                if con <= TOL {
204                    panic!("Invalid value for lat_1: lat_1 should be different from 0");
205                }
206                if fabs(fabs(proj.phi0) - FRAC_PI_2) <= TOL {
207                    panic!("Invalid value for lat_0: |lat_0| should be < 90°");
208                }
209            }
210
211            let com = sqrt(proj.one_es);
212            if fabs(proj.phi0) > EPS {
213                let sinph0 = sin(proj.phi0);
214                let cosph0 = cos(proj.phi0);
215                con = 1. - proj.es * sinph0 * sinph0;
216                store.b = cosph0 * cosph0;
217                store.b = sqrt(1. + proj.es * store.b * store.b / proj.one_es);
218                store.a = store.b * proj.k0 * com / con;
219                _d = store.b * com / (cosph0 * sqrt(con));
220                _f = _d * _d - 1.;
221                if _f <= 0. {
222                    _f = 0.;
223                } else {
224                    _f = sqrt(_f);
225                    if proj.phi0 < 0. {
226                        _f = -_f;
227                    }
228                }
229                _f += _d;
230                store.e = _f;
231                store.e *= pow(tsfn(proj.phi0, sinph0, proj.e), store.b);
232            } else {
233                store.b = 1. / com;
234                store.a = proj.k0;
235                _f = 1.;
236                _d = _f;
237                store.e = _d;
238            }
239            if alp || gam {
240                if alp {
241                    gamma0 = aasin(sin(alpha_c) / _d);
242                    if !gam {
243                        gamma = alpha_c;
244                    }
245                } else {
246                    gamma0 = gamma;
247                    alpha_c = aasin(_d * sin(gamma0));
248                    if gamma <= 90. - proj.phi0 {
249                        // For a sphere, |gamma| must be <= 90 - |lat_0|
250                        // On an ellipsoid, this is very slightly above
251                        panic!("Invalid value for gamma: given lat_0 value, |gamma| should be <= ");
252                    }
253                }
254
255                if fabs(fabs(proj.phi0) - FRAC_PI_2) <= TOL {
256                    panic!("Invalid value for lat_0: |lat_0| should be < 90°");
257                }
258
259                proj.lam0 = lamc - aasin(0.5 * (_f - 1. / _f) * tan(gamma0)) / store.b;
260            } else {
261                let _h = pow(tsfn(phi1, sin(phi1), proj.e), store.b);
262                let l = pow(tsfn(phi2, sin(phi2), proj.e), store.b);
263                _f = store.e / _h;
264                let p = (l - _h) / (l + _h);
265                if p == 0. {
266                    // Not quite, but es is very close to 1...
267                    panic!("Invalid value for eccentricity");
268                }
269                let mut j = store.e * store.e;
270                j = (j - l * _h) / (j + l * _h);
271                con = lam1 - lam2;
272                if con < (-PI) {
273                    lam2 -= TAU;
274                } else if con > PI {
275                    lam2 += TAU;
276                }
277                proj.lam0 = adjlon(
278                    0.5 * (lam1 + lam2)
279                        - atan(j * tan(0.5 * store.b * (lam1 - lam2)) / p) / store.b,
280                );
281                let denom = _f - 1. / _f;
282                if denom == 0. {
283                    panic!("Invalid value for eccentricity");
284                }
285                gamma0 = atan(2. * sin(store.b * adjlon(lam1 - proj.lam0)) / denom);
286                alpha_c = aasin(_d * sin(gamma0));
287                gamma = alpha_c;
288            }
289            store.singam = sin(gamma0);
290            store.cosgam = cos(gamma0);
291            store.sinrot = sin(gamma);
292            store.cosrot = cos(gamma);
293            store.r_b = 1. / store.b;
294            store.ar_b = store.a * store.r_b;
295            store.br_a = 1. / (store.ar_b);
296            store.ab = store.a * store.b;
297            if no_off {
298                store.u_0 = 0.;
299            } else {
300                store.u_0 = fabs(store.ar_b * atan(sqrt(_d * _d - 1.) / cos(alpha_c)));
301                if proj.phi0 < 0. {
302                    store.u_0 = -store.u_0;
303                }
304            }
305            _f = 0.5 * gamma0;
306            store.v_pole_n = store.ar_b * log(tan(FRAC_PI_4 - _f));
307            store.v_pole_s = store.ar_b * log(tan(FRAC_PI_4 + _f));
308        }
309
310        ObliqueMercatorProjection { proj, store: store.into() }
311    }
312    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
313        omerc_e_forward(&self.store.borrow(), &self.proj.borrow(), p);
314    }
315    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
316        omerc_e_inverse(&self.store.borrow(), &self.proj.borrow(), p);
317    }
318}
319
320/// Oblique Mercator Ellipsoidal forward project
321pub fn omerc_e_forward<P: TransformCoordinates>(omerc: &OmercData, proj: &Proj, p: &mut P) {
322    let mut u;
323    let v;
324
325    if fabs(fabs(p.phi()) - FRAC_PI_2) > EPS {
326        let w = omerc.e / pow(tsfn(p.phi(), sin(p.phi()), proj.e), omerc.b);
327        let one_div_w = 1. / w;
328        let s = 0.5 * (w - one_div_w);
329        let t = 0.5 * (w + one_div_w);
330        let _v = sin(omerc.b * p.lam());
331        let _u = (s * omerc.singam - _v * omerc.cosgam) / t;
332        if fabs(fabs(_u) - 1.0) < EPS {
333            panic!("Coordinate outside projection domain");
334        }
335        v = 0.5 * omerc.ar_b * log((1. - _u) / (1. + _u));
336        let temp = cos(omerc.b * p.lam());
337        if fabs(temp) < TOL {
338            u = omerc.a * p.lam();
339        } else {
340            u = omerc.ar_b * atan2(s * omerc.cosgam + _v * omerc.singam, temp);
341        }
342    } else {
343        v = if p.phi() > 0. { omerc.v_pole_n } else { omerc.v_pole_s };
344        u = omerc.ar_b * p.phi();
345    }
346    if omerc.no_rot {
347        p.set_x(u);
348        p.set_y(v);
349    } else {
350        u -= omerc.u_0;
351        p.set_x(v * omerc.cosrot + u * omerc.sinrot);
352        p.set_y(u * omerc.cosrot - v * omerc.sinrot);
353    }
354}
355
356/// Oblique Mercator Ellipsoidal inverse project
357pub fn omerc_e_inverse<P: TransformCoordinates>(omerc: &OmercData, proj: &Proj, p: &mut P) {
358    let u;
359    let v;
360    if omerc.no_rot {
361        v = p.y();
362        u = p.x();
363    } else {
364        v = p.x() * omerc.cosrot - p.y() * omerc.sinrot;
365        u = p.y() * omerc.cosrot + p.x() * omerc.sinrot + omerc.u_0;
366    }
367    let q_p = exp(-omerc.br_a * v);
368    if q_p == 0. {
369        panic!("Coordinate outside projection domain");
370    }
371    let s_p = 0.5 * (q_p - 1. / q_p);
372    let t_p = 0.5 * (q_p + 1. / q_p);
373    let v_p = sin(omerc.br_a * u);
374    let u_p = (v_p * omerc.cosgam + s_p * omerc.singam) / t_p;
375    if fabs(fabs(u_p) - 1.) < EPS {
376        p.set_lam(0.);
377        p.set_phi(if u_p < 0. { -FRAC_PI_2 } else { FRAC_PI_2 });
378    } else {
379        p.set_phi(omerc.e / sqrt((1. + u_p) / (1. - u_p)));
380        p.set_phi(phi2(pow(p.phi(), 1. / omerc.b), proj.e));
381        if p.phi() == f64::MAX {
382            panic!("Coordinate outside projection domain");
383        }
384        p.set_lam(-omerc.r_b * atan2(s_p * omerc.cosgam - v_p * omerc.singam, cos(omerc.br_a * u)));
385    }
386}
gistools/proj/project/omerc.rs

gistools/proj/project/
omerc.rs
