gistools/proj/project/
gnom.rs

1use crate::proj::{
2    CoordinateStep, EPS10, GeodGeodesic, GeodGeodesicline, GeodMask, Proj, ProjMethod, ProjMode,
3    ProjectCoordinates, TransformCoordinates, geod_geninverse, geod_genposition, geod_init,
4    geod_lineinit,
5};
6use alloc::rc::Rc;
7use core::{cell::RefCell, f64::consts::FRAC_PI_2};
8use libm::{asin, atan, atan2, cos, fabs, hypot, sin, sqrt};
9
10/// Gnomonic Variables
11#[derive(Debug, Default, Clone, PartialEq)]
12pub struct GnomData {
13    sinph0: f64,
14    cosph0: f64,
15    mode: ProjMode,
16    g: GeodGeodesic,
17}
18
19/// # Gnomonic (gnom)
20///
21/// For a sphere, the gnomonic projection is a projection from the center of
22/// the sphere onto a plane tangent to the center point of the projection.
23/// This projects great circles to straight lines.  For an ellipsoid, it is
24/// the limit of a doubly azimuthal projection, a projection where the
25/// azimuths from 2 points are preserved, as the two points merge into the
26/// center point.  In this case, geodesics project to approximately straight
27/// lines (these are exactly straight if the geodesic includes the center
28/// point).  For details, see Section 8 of :cite:`Karney2013`.
29///
30/// **Classification**: Azimuthal
31///
32/// **Available forms**: Forward and inverse, spherical and ellipsoidal
33///
34/// **Defined area**: Within a quarter circumference of the center point
35///
36/// **Alias**: gnom
37///
38/// **Domain**: 2D
39///
40/// **Input type**: Geodetic coordinates
41///
42/// **Output type**: Projected coordinates
43///
44/// ## Projection String
45/// ```ini
46/// +proj=gnom +lat_0=90 +lon_0=-50 +R=6.4e6
47/// ```
48///
49/// ## Required Parameters
50/// - None, all parameters are optional for this projection.
51///
52/// ## Optional Parameters
53/// - `+lon_0`: Longitude of origin (central meridian).
54/// - `+lat_0`: Latitude of origin.
55/// - `+x_0`: False easting.
56/// - `+y_0`: False northing.
57/// - `+ellps`: Ellipsoid.
58/// - `+R`: Earth radius.
59///
60/// Reference:
61/// Wolfram Mathworld "Gnomonic Projection"
62/// <http://mathworld.wolfram.com/GnomonicProjection.html>
63/// Accessed: 12th November 2009
64///
65/// ![Gnomonic](https://github.com/Open-S2/gis-tools/blob/master/assets/proj4/projections/images/gnom.png?raw=true)
66#[derive(Debug, Clone, PartialEq)]
67pub struct GnomonicProjection {
68    proj: Rc<RefCell<Proj>>,
69    store: RefCell<GnomData>,
70    method: ProjMethod,
71}
72impl ProjectCoordinates for GnomonicProjection {
73    fn code(&self) -> i64 {
74        -1
75    }
76    fn name(&self) -> &'static str {
77        "Gnomonic"
78    }
79    fn names() -> &'static [&'static str] {
80        &["Gnomonic", "gnom"]
81    }
82}
83impl CoordinateStep for GnomonicProjection {
84    fn new(proj: Rc<RefCell<Proj>>) -> Self {
85        let mut store = GnomData::default();
86        let method: ProjMethod;
87        {
88            let proj = &mut proj.borrow_mut();
89
90            method = if proj.es == 0. {
91                if fabs(fabs(proj.phi0) - FRAC_PI_2) < EPS10 {
92                    store.mode = if proj.phi0 < 0. { ProjMode::SPole } else { ProjMode::NPole };
93                } else if fabs(proj.phi0) < EPS10 {
94                    store.mode = ProjMode::Equit;
95                } else {
96                    store.mode = ProjMode::Obliq;
97                    store.sinph0 = sin(proj.phi0);
98                    store.cosph0 = cos(proj.phi0);
99                }
100                ProjMethod::Spheroidal
101            } else {
102                geod_init(&mut store.g, 1., proj.f);
103                ProjMethod::Ellipsoidal
104            };
105            proj.es = 0.;
106        }
107        GnomonicProjection { proj, store: store.into(), method }
108    }
109    fn forward<P: TransformCoordinates>(&self, p: &mut P) {
110        if self.method == ProjMethod::Ellipsoidal {
111            gnom_e_forward(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
112        } else {
113            gnom_s_forward(&mut self.store.borrow_mut(), p);
114        }
115    }
116    fn inverse<P: TransformCoordinates>(&self, p: &mut P) {
117        if self.method == ProjMethod::Ellipsoidal {
118            gnom_e_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
119        } else {
120            gnom_s_inverse(&mut self.store.borrow_mut(), &self.proj.borrow(), p);
121        }
122    }
123}
124
125/// Gnomonic Spheroidal Forward
126pub fn gnom_s_forward<P: TransformCoordinates>(gnom: &mut GnomData, p: &mut P) {
127    let mut y;
128
129    let sinphi = sin(p.phi());
130    let cosphi = cos(p.phi());
131    let mut coslam = cos(p.lam());
132
133    match gnom.mode {
134        ProjMode::Equit => {
135            y = cosphi * coslam;
136        }
137        ProjMode::Obliq => {
138            y = gnom.sinph0 * sinphi + gnom.cosph0 * cosphi * coslam;
139        }
140        ProjMode::SPole => {
141            y = -sinphi;
142        }
143        ProjMode::NPole => {
144            y = sinphi;
145        }
146    }
147
148    if y <= EPS10 {
149        panic!("Coordinate outside projection domain");
150    }
151
152    y = 1. / y;
153    let x = (y) * cosphi * sin(p.lam());
154    match gnom.mode {
155        ProjMode::Equit => {
156            y *= sinphi;
157        }
158        ProjMode::Obliq => {
159            y *= gnom.cosph0 * sinphi - gnom.sinph0 * cosphi * coslam;
160        }
161        ProjMode::NPole => {
162            coslam = -coslam;
163            y *= cosphi * coslam;
164        }
165        ProjMode::SPole => {
166            y *= cosphi * coslam;
167        }
168    }
169
170    p.set_x(x);
171    p.set_y(y);
172}
173
174/// Gnomonic Spheroidal inverse
175pub fn gnom_s_inverse<P: TransformCoordinates>(gnom: &mut GnomData, proj: &Proj, p: &mut P) {
176    let mut x = p.x();
177    let mut y = p.y();
178
179    let rh = hypot(x, y);
180    let mut phi = atan(rh);
181    let lam;
182    let sinz = sin(phi);
183    let cosz = sqrt(1. - sinz * sinz);
184
185    if fabs(rh) <= EPS10 {
186        phi = proj.phi0;
187        lam = 0.;
188    } else {
189        match gnom.mode {
190            ProjMode::Obliq => {
191                phi = cosz * gnom.sinph0 + y * sinz * gnom.cosph0 / rh;
192                if fabs(phi) >= 1. {
193                    phi = if phi > 0. { FRAC_PI_2 } else { -FRAC_PI_2 };
194                } else {
195                    phi = asin(phi);
196                }
197                y = (cosz - gnom.sinph0 * sin(phi)) * rh;
198                x *= sinz * gnom.cosph0;
199            }
200            ProjMode::Equit => {
201                phi = y * sinz / rh;
202                if fabs(phi) >= 1. {
203                    phi = if phi > 0. { FRAC_PI_2 } else { -FRAC_PI_2 };
204                } else {
205                    phi = asin(phi);
206                }
207                y = cosz * rh;
208                x *= sinz;
209            }
210            ProjMode::SPole => {
211                phi -= FRAC_PI_2;
212            }
213            ProjMode::NPole => {
214                phi = FRAC_PI_2 - phi;
215                y = -y;
216            }
217        }
218        lam = atan2(x, y);
219    }
220
221    p.set_phi(phi);
222    p.set_lam(lam);
223}
224
225/// Gnomonic Ellipsoidal Forward
226pub fn gnom_e_forward<P: TransformCoordinates>(gnom: &mut GnomData, proj: &Proj, p: &mut P) {
227    let lat0 = proj.phi0.to_degrees();
228    let lon0 = 0.;
229    let lat1 = p.phi().to_degrees();
230    let lon1 = p.lam().to_degrees();
231    let mut azi0 = 0.;
232    let mut m = 0.;
233    let mut _m = 0.;
234
235    geod_geninverse(
236        &mut gnom.g,
237        lat0,
238        lon0,
239        lat1,
240        lon1,
241        &mut 0.,
242        &mut azi0,
243        &mut 0.,
244        &mut m,
245        &mut _m,
246        &mut 0.,
247        &mut 0.,
248    );
249    if _m <= 0. {
250        panic!("Coordinate outside projection domain {_m}");
251    } else {
252        let rho = m / _m;
253        azi0 = azi0.to_radians();
254        p.set_x(rho * sin(azi0));
255        p.set_y(rho * cos(azi0));
256    }
257}
258
259/// Gnomonic Ellipsoidal inverse
260pub fn gnom_e_inverse<P: TransformCoordinates>(gnom: &mut GnomData, proj: &Proj, p: &mut P) {
261    let numit_ = 10;
262    let eps_ = 0.01 * sqrt(f64::EPSILON);
263    let lat0 = proj.phi0.to_degrees();
264    let lon0 = 0.;
265    let azi0 = atan2(p.x(), p.y()).to_degrees();
266    let mut rho = hypot(p.x(), p.y());
267    let mut s = atan(rho);
268    let little = rho <= 1.;
269    if !little {
270        rho = 1. / rho;
271    }
272    let mut l = GeodGeodesicline::default();
273    geod_lineinit(
274        &mut l,
275        &gnom.g,
276        lat0,
277        lon0,
278        azi0,
279        GeodMask::GeodLatitude as u32
280            | GeodMask::GeodLongitude as u32
281            | GeodMask::GeodDistanceIn as u32
282            | GeodMask::GeodReducedlength as u32
283            | GeodMask::GeodGeodesicScale as u32,
284    );
285
286    let mut lat1 = 0.;
287    let mut lon1 = 0.;
288    let mut count = numit_;
289    let mut trip = 0;
290    while count != 0 {
291        count -= 1;
292        // double m, M;
293        let mut m = 0.;
294        let mut _m = 0.;
295        geod_genposition(
296            &l, 0, s, &mut lat1, &mut lon1, &mut 0., &mut s, &mut m, &mut _m, &mut 0., &mut 0.,
297        );
298        if trip != 0 {
299            break;
300        }
301        // If little, solve rho(s) = rho with drho(s)/ds = 1/M^2
302        // else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2
303        let ds = if little { (m - rho * _m) * _m } else { (rho * m - _m) * m };
304        s -= ds;
305        // Reversed test to allow escape with NaNs
306        if fabs(ds) < eps_ {
307            trip += 1;
308        }
309    }
310    if trip != 0 {
311        p.set_phi(lat1.to_radians());
312        p.set_lam(lon1.to_radians());
313    } else {
314        panic!("Coordinate outside projection domain");
315    }
316}
gistools/proj/project/gnom.rs

gistools/proj/project/
gnom.rs
