gistools/proj/project/
ortho.rs

1use crate::proj::{
2    AZIMUTH_PROJECTION_CENTRE, CoordinateStep, Coords, EPS10, ORTHOGRAPHIC, Proj, ProjMethod,
3    ProjMode, ProjValue, ProjectCoordinates, TransformCoordinates, adjlon,
4};
5use alloc::rc::Rc;
6use core::{
7    cell::RefCell,
8    f64::consts::{FRAC_PI_2, PI},
9};
10use libm::{acos, asin, atan2, cos, fabs, hypot, sin, sqrt};
11
12/// Orthographic variables
13#[derive(Debug, Default, Clone, PartialEq)]
14pub struct OrthoData {
15    sinph0: f64,
16    cosph0: f64,
17    nu0: f64,
18    y_shift: f64,
19    y_scale: f64,
20    mode: ProjMode,
21    sinalpha: f64,
22    cosalpha: f64,
23}
24
25/// # Orthographic
26///
27/// The orthographic projection is a perspective azimuthal projection centered
28/// around a given latitude and longitude.
29///
30/// **Classification**: Azimuthal
31///
32/// **Available forms**: Forward and inverse, spherical and ellipsoidal
33///
34/// **Defined area**: Global, although only one hemisphere can be seen at a time
35///
36/// **Alias**: ortho
37///
38/// **Domain**: 2D
39///
40/// **Input type**: Geodetic coordinates
41///
42/// **Output type**: Projected coordinates
43///
44/// ## Projection String
45/// ```ini
46/// +proj=ortho
47/// ```
48///
49/// ## Required Parameters
50/// - None
51///
52/// ## Optional Parameters
53/// - `+alpha=<value>`: Azimuth clockwise from north at the center of projection. *Defaults to 0.0.* (added in PROJ 9.5.0)
54/// - `+k_0=<value>`: Scale factor. Determines scale factor used in the projection. *Defaults to 1.0.* (added in PROJ 9.5.0)
55/// - `+lon_0=<value>`
56/// - `+lat_0=<value>`
57/// - `+ellps=<value>`
58/// - `+R=<value>`
59/// - `+x_0=<value>`
60/// - `+y_0=<value>`
61///
62/// ## Notes
63/// - Before PROJ 7.2, only the spherical formulation was implemented. To replicate PROJ < 7.2 results with newer versions, force the ellipsoid to a sphere using `+f=0`.
64/// - This projection method corresponds to `EPSG:9840` (or `EPSG:1130` with `k_0` or `alpha`).
65///
66/// ![Orthographic](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/ortho.png?raw=true)
67#[derive(Debug, Clone, PartialEq)]
68pub struct OrthographicProjection {
69    proj: Rc<RefCell<Proj>>,
70    store: RefCell<OrthoData>,
71    method: ProjMethod,
72}
73impl ProjectCoordinates for OrthographicProjection {
74    fn code(&self) -> i64 {
75        ORTHOGRAPHIC
76    }
77    fn name(&self) -> &'static str {
78        "Orthographic"
79    }
80    fn names() -> &'static [&'static str] {
81        &["Orthographic", "ortho"]
82    }
83}
84impl CoordinateStep for OrthographicProjection {
85    fn new(proj: Rc<RefCell<Proj>>) -> Self {
86        let mut store = OrthoData::default();
87        let method: ProjMethod;
88        {
89            let proj = &mut proj.borrow_mut();
90            store.sinph0 = sin(proj.phi0);
91            store.cosph0 = cos(proj.phi0);
92            if fabs(fabs(proj.phi0) - FRAC_PI_2) <= EPS10 {
93                store.mode = if proj.phi0 < 0. { ProjMode::SPole } else { ProjMode::NPole };
94            } else if fabs(proj.phi0) > EPS10 {
95                store.mode = ProjMode::Obliq;
96            } else {
97                store.mode = ProjMode::Equit;
98            }
99            method = if proj.es == 0. {
100                ProjMethod::Spheroidal
101            } else {
102                store.nu0 = 1.0 / sqrt(1.0 - proj.es * store.sinph0 * store.sinph0);
103                store.y_shift = proj.es * store.nu0 * store.sinph0 * store.cosph0;
104                store.y_scale = 1.0 / sqrt(1.0 - proj.es * store.cosph0 * store.cosph0);
105                ProjMethod::Ellipsoidal
106            };
107
108            let alpha =
109                proj.params.get(&AZIMUTH_PROJECTION_CENTRE).unwrap_or(&ProjValue::default()).f64();
110            store.sinalpha = sin(alpha);
111            store.cosalpha = cos(alpha);
112        }
113        OrthographicProjection { proj, store: store.into(), method }
114    }
115    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
116        if self.method == ProjMethod::Spheroidal {
117            ortho_s_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
118        } else {
119            ortho_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
120        }
121    }
122    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
123        if self.method == ProjMethod::Spheroidal {
124            ortho_s_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
125        } else {
126            ortho_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
127        }
128    }
129}
130
131fn throw_error() {
132    panic!("Coordinate outside projection domain");
133}
134
135/// Equal Earth Spheroidal forward project
136pub fn ortho_s_forward<P: TransformCoordinates>(ortho: &mut OrthoData, proj: &Proj, p: &mut P) {
137    let mut y;
138
139    let cosphi = cos(p.phi());
140    let mut coslam = cos(p.lam());
141    match ortho.mode {
142        ProjMode::Equit => {
143            if cosphi * coslam < -EPS10 {
144                throw_error();
145            }
146            y = sin(p.phi());
147        }
148        ProjMode::Obliq => {
149            let sinphi = sin(p.phi());
150            // Is the point visible from the projection plane ?
151            // From
152            // https://lists.osgeo.org/pipermail/proj/2020-September/009831.html
153            // this is the dot product of the normal of the ellipsoid at the center
154            // of the projection and at the point considered for projection.
155            // [cos(phi)*cos(lambda), cos(phi)*sin(lambda), sin(phi)]
156            // Also from Snyder's Map Projection - A working manual, equation (5-3),
157            // page 149
158            if ortho.sinph0 * sinphi + ortho.cosph0 * cosphi * coslam < -EPS10 {
159                throw_error();
160            }
161            y = ortho.cosph0 * sinphi - ortho.sinph0 * cosphi * coslam;
162        }
163        ProjMode::NPole => {
164            coslam = -coslam;
165            if fabs(p.phi() - proj.phi0) - EPS10 > FRAC_PI_2 {
166                throw_error();
167            }
168            y = cosphi * coslam;
169        }
170        ProjMode::SPole => {
171            if fabs(p.phi() - proj.phi0) - EPS10 > FRAC_PI_2 {
172                throw_error();
173            }
174            y = cosphi * coslam;
175        }
176    }
177    let mut x = cosphi * sin(p.lam());
178
179    let xp = x;
180    let yp = y;
181    x = (xp * ortho.cosalpha - yp * ortho.sinalpha) * proj.k0;
182    y = (xp * ortho.sinalpha + yp * ortho.cosalpha) * proj.k0;
183
184    p.set_x(x);
185    p.set_y(y);
186}
187
188/// Equal Earth Spheroidal inverse project
189pub fn ortho_s_inverse<P: TransformCoordinates>(ortho: &mut OrthoData, proj: &Proj, p: &mut P) {
190    let lam;
191    let mut phi;
192
193    let xf = p.x();
194    let yf = p.y();
195    let mut x = (ortho.cosalpha * xf + ortho.sinalpha * yf) / proj.k0;
196    let mut y = (-ortho.sinalpha * xf + ortho.cosalpha * yf) / proj.k0;
197
198    let rh = hypot(x, y);
199    let mut sinc = rh;
200    if sinc > 1. {
201        if (sinc - 1.) > EPS10 {
202            throw_error();
203        }
204        sinc = 1.;
205    }
206    let cosc = sqrt(1. - sinc * sinc); /* in this range OK */
207    if fabs(rh) <= EPS10 {
208        phi = proj.phi0;
209        lam = 0.0;
210    } else {
211        match ortho.mode {
212            ProjMode::NPole => {
213                y = -y;
214                phi = acos(sinc);
215            }
216            ProjMode::SPole => {
217                phi = -acos(sinc);
218            }
219            ProjMode::Equit => {
220                phi = y * sinc / rh;
221                x *= sinc;
222                y = cosc * rh;
223                // goto sinchk;
224                if fabs(phi) >= 1. {
225                    phi = if phi < 0. { -FRAC_PI_2 } else { FRAC_PI_2 };
226                } else {
227                    phi = asin(phi);
228                }
229            }
230            ProjMode::Obliq => {
231                phi = cosc * ortho.sinph0 + y * sinc * ortho.cosph0 / rh;
232                y = (cosc - ortho.sinph0 * phi) * rh;
233                x *= sinc * ortho.cosph0;
234                // goto sinchk;
235                if fabs(phi) >= 1. {
236                    phi = if phi < 0. { -FRAC_PI_2 } else { FRAC_PI_2 };
237                } else {
238                    phi = asin(phi);
239                }
240            }
241        }
242        lam = if y == 0. && (ortho.mode == ProjMode::Obliq || ortho.mode == ProjMode::Equit) {
243            if x == 0. {
244                0.
245            } else if x < 0. {
246                -FRAC_PI_2
247            } else {
248                FRAC_PI_2
249            }
250        } else {
251            atan2(x, y)
252        };
253    }
254
255    p.set_lam(lam);
256    p.set_phi(phi);
257}
258
259/// Equal Earth Ellipsoidal forward project
260pub fn ortho_e_forward<P: TransformCoordinates>(ortho: &mut OrthoData, proj: &Proj, p: &mut P) {
261    // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic
262    let cosphi = cos(p.phi());
263    let sinphi = sin(p.phi());
264    let coslam = cos(p.lam());
265    let sinlam = sin(p.lam());
266
267    // Is the point visible from the projection plane ?
268    // Same condition as in spherical case
269    if ortho.sinph0 * sinphi + ortho.cosph0 * cosphi * coslam < -EPS10 {
270        throw_error();
271    }
272
273    let nu = 1.0 / sqrt(1.0 - proj.es * sinphi * sinphi);
274    let xp = nu * cosphi * sinlam;
275    let yp = nu * (sinphi * ortho.cosph0 - cosphi * ortho.sinph0 * coslam)
276        + proj.es * (ortho.nu0 * ortho.sinph0 - nu * sinphi) * ortho.cosph0;
277    p.set_x((ortho.cosalpha * xp - ortho.sinalpha * yp) * proj.k0);
278    p.set_y((ortho.sinalpha * xp + ortho.cosalpha * yp) * proj.k0);
279}
280
281/// Equal Earth Ellipsoidal inverse project
282pub fn ortho_e_inverse<P: TransformCoordinates>(ortho: &mut OrthoData, proj: &Proj, p: &mut P) {
283    let sq = |x: f64| -> f64 { x * x };
284
285    let xf = p.x();
286    let yf = p.y();
287    let x = (ortho.cosalpha * xf + ortho.sinalpha * yf) / proj.k0;
288    let y = (-ortho.sinalpha * xf + ortho.cosalpha * yf) / proj.k0;
289
290    if ortho.mode == ProjMode::NPole || ortho.mode == ProjMode::SPole {
291        // Polar case. Forward case equations can be simplified as:
292        // x = nu * cosphi * sinlam
293        // y = nu * -cosphi * coslam * sign(phi0)
294        // ==> lam = atan2(x, -y * sign(phi0))
295        // ==> x^2 + y^2 = nu^2 * cosphi^2
296        //                rh^2 = cosphi^2 / (1 - es * sinphi^2)
297        // ==>  cosphi^2 = rh^2 * (1 - es) / (1 - es * rh^2)
298
299        let rh2 = sq(x) + sq(y);
300        if rh2 >= 1. - 1e-15 {
301            if (rh2 - 1.) > EPS10 {
302                throw_error();
303            }
304            p.set_phi(0.);
305        } else {
306            p.set_phi(
307                acos(sqrt(rh2 * proj.one_es / (1. - proj.es * rh2)))
308                    * (if ortho.mode == ProjMode::NPole { 1. } else { -1. }),
309            );
310        }
311        p.set_lam(atan2(x, y * (if ortho.mode == ProjMode::NPole { -1. } else { 1. })));
312        return;
313    }
314
315    if ortho.mode == ProjMode::Equit {
316        // Equatorial case. Forward case equations can be simplified as:
317        // x = nu * cosphi * sinlam
318        // y  = nu * sinphi * (1 - proj.es)
319        // x^2 * (1 - es * sinphi^2) = (1 - sinphi^2) * sinlam^2
320        // y^2 / ((1 - es)^2 + y^2 * es) = sinphi^2
321
322        // Equation of the ellipse
323        if sq(x) + sq(y * (proj.a / proj.b)) > 1. + 1.0e-11 {
324            throw_error();
325        }
326
327        let sinphi2 = if y == 0. { 0. } else { 1.0 / (sq((1. - proj.es) / y) + proj.es) };
328        if sinphi2 > 1. - 1e-11 {
329            p.set_phi(FRAC_PI_2 * (if y > 0. { 1. } else { -1. }));
330            p.set_lam(0.);
331            return;
332        }
333        p.set_phi(asin(sqrt(sinphi2)) * (if y > 0. { 1. } else { -1. }));
334        let sinlam = x * sqrt((1. - proj.es * sinphi2) / (1. - sinphi2));
335        if fabs(sinlam) - 1. > -1e-15 {
336            p.set_lam(FRAC_PI_2 * (if x > 0. { 1. } else { -1. }));
337        } else {
338            p.set_lam(asin(sinlam));
339        }
340        return;
341    }
342
343    // Using ortho.sinph0 * sinphi + ortho.cosph0 * cosphi * coslam == 0 (visibity
344    // condition of the forward case) in the forward equations, and a lot of
345    // substitution games...
346    // PJ_XY xy_recentered;
347    let mut xy_recentered = Coords::default();
348    xy_recentered.set_x(x);
349    xy_recentered.set_y((y - ortho.y_shift) / ortho.y_scale);
350    if sq(x) + sq(xy_recentered.y()) > 1. + 1e-11 {
351        throw_error();
352    }
353
354    // From EPSG guidance note 7.2, March 2020, §3.3.5 Orthographic
355
356    // It suggests as initial guess:
357    // lp.lam = 0;
358    // lp.phi = proj.phi0;
359    // But for poles, this will not converge well. Better use:
360    ortho_s_inverse(ortho, proj, &mut xy_recentered);
361    p.set_x(xy_recentered.x());
362    p.set_y(xy_recentered.y() * ortho.y_scale + ortho.y_shift);
363
364    for _ in 0..20 {
365        let cosphi = cos(p.phi());
366        let sinphi = sin(p.phi());
367        let coslam = cos(p.lam());
368        let sinlam = sin(p.lam());
369        let one_minus_es_sinphi2 = 1.0 - proj.es * sinphi * sinphi;
370        let nu = 1.0 / sqrt(one_minus_es_sinphi2);
371        let mut xy_new = Coords::default();
372        xy_new.set_x(nu * cosphi * sinlam);
373        xy_new.set_y(
374            nu * (sinphi * ortho.cosph0 - cosphi * ortho.sinph0 * coslam)
375                + proj.es * (ortho.nu0 * ortho.sinph0 - nu * sinphi) * ortho.cosph0,
376        );
377        let rho = (1.0 - proj.es) * nu / one_minus_es_sinphi2;
378        let j11 = -rho * sinphi * sinlam;
379        let j12 = nu * cosphi * coslam;
380        let j21 = rho * (cosphi * ortho.cosph0 + sinphi * ortho.sinph0 * coslam);
381        let j22 = nu * ortho.sinph0 * cosphi * sinlam;
382        let d = j11 * j22 - j12 * j21;
383        let dx = x - xy_new.x();
384        let dy = y - xy_new.y();
385        let dphi = (j22 * dx - j12 * dy) / d;
386        let dlam = (-j21 * dx + j11 * dy) / d;
387        p.set_phi(p.phi() + dphi);
388        if p.phi() > FRAC_PI_2 {
389            p.set_phi(FRAC_PI_2 - (p.phi() - FRAC_PI_2));
390            p.set_lam(adjlon(p.lam() + PI));
391        } else if p.phi() < -FRAC_PI_2 {
392            p.set_phi(-FRAC_PI_2 + (-FRAC_PI_2 - p.phi()));
393            p.set_lam(adjlon(p.lam() + PI));
394        }
395        p.set_lam(p.lam() + dlam);
396        if fabs(dphi) < 1e-12 && fabs(dlam) < 1e-12 {
397            return;
398        }
399    }
400    throw_error();
401}
gistools/proj/project/ortho.rs

gistools/proj/project/
ortho.rs
