mappers/projections/
modified_azimuthal_equidistant.rs

1//! This is a modified version of the [`Azimuthal Equidistant`](crate::projections::azimuthal_equidistant) projection,
2//! defined for islands of Micronesia and described by [Snyder (1987)](https://pubs.er.usgs.gov/publication/pp1395).
3//! It does not use Geodesic calculations so it is faster than the original projection, but significantly
4//! diverges from the AEQD projection at bigger scales.
5//!
6//! The azimuthal equidistant projection is an azimuthal map projection.
7//! It has the useful properties that all points on the map are at proportionally
8//! correct distances from the center point, and that all points on the map are at the
9//! correct azimuth (direction) from the center point. A useful application for this
10//! type of projection is a polar projection which shows all meridians (lines of longitude) as straight,
11//! with distances from the pole represented correctly [(Wikipedia, 2022)](https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection).
12//!
13//! Summary by [Snyder (1987)](https://pubs.er.usgs.gov/publication/pp1395):
14//!
15//! - Azimuthal.
16//! - Distances measured from the center are true.
17//! - Distances not measured along radii from the center are not correct.
18//! - The center of projection is the only point without distortion.
19//! - Directions from the center are true (except on some oblique and equatorial ellipsoidal forms).
20//! - Neither equal-area nor conformal.
21//! - All meridians on the polar aspect, the central meridian on other aspects, and the Equator on the equatorial aspect are straight lines.
22//! - Parallels on the polar projection are circles spaced at true intervals (equidistant for the sphere).
23//! - The outer meridian of a hemisphere on the equatorial aspect (for the sphere) is a circle.
24//! - All other meridians and parallels are complex curves.
25//! - Not a perspective projection.
26//! - Point opposite the center is shown as a circle (for the sphere) surrounding the map.
27//! - Used in the polar aspect for world maps and maps of polar hemispheres.
28//! - Used in the oblique aspect for atlas maps of continents and world maps for aviation and radio use.
29//! - Known for many centuries in the polar aspect.
30
31use float_cmp::approx_eq;
32
33use crate::{
34    Projection, ProjectionError,
35    ellipsoids::Ellipsoid,
36    errors::{ensure_finite, ensure_within_range, unpack_required_parameter},
37};
38
39#[cfg(feature = "tracing")]
40use tracing::instrument;
41
42/// Main projection struct that is constructed from [`ModifiedAzimuthalEquidistantBuilder`] and used for computations.
43#[derive(Copy, Clone, PartialEq, PartialOrd, Debug)]
44pub struct ModifiedAzimuthalEquidistant {
45    lon_0: f64,
46    lat_0: f64,
47    n_1: f64,
48    g: f64,
49    ellps: Ellipsoid,
50}
51
52impl ModifiedAzimuthalEquidistant {
53    /// Initializes builder with default values.
54    /// Projection parameters can be set with builder methods,
55    /// refer to the documentation of those methods to check which parmeters are required
56    /// and default values for optional arguments.
57    #[must_use]
58    pub fn builder() -> ModifiedAzimuthalEquidistantBuilder {
59        ModifiedAzimuthalEquidistantBuilder::default()
60    }
61}
62
63/// Builder struct which allows to construct [`ModifiedAzimuthalEquidistant`] projection.
64/// Refer to the documentation of this struct's methods to check which parmeters are required
65/// and default values for optional arguments.
66#[derive(Copy, Clone, PartialEq, PartialOrd, Debug)]
67pub struct ModifiedAzimuthalEquidistantBuilder {
68    ref_lon: Option<f64>,
69    ref_lat: Option<f64>,
70    ellipsoid: Ellipsoid,
71}
72
73impl Default for ModifiedAzimuthalEquidistantBuilder {
74    fn default() -> Self {
75        Self {
76            ref_lon: None,
77            ref_lat: None,
78            ellipsoid: Ellipsoid::WGS84,
79        }
80    }
81}
82
83impl ModifiedAzimuthalEquidistantBuilder {
84    /// *(required)* Sets reference longitude and latitude. Point (0, 0) on the map will be at this coordinates.
85    pub const fn ref_lonlat(&mut self, lon: f64, lat: f64) -> &mut Self {
86        self.ref_lon = Some(lon);
87        self.ref_lat = Some(lat);
88        self
89    }
90
91    /// *(optional)* Sets reference [`Ellipsoid`], defaults to [`WGS84`](Ellipsoid::WGS84).
92    pub const fn ellipsoid(&mut self, ellps: Ellipsoid) -> &mut Self {
93        self.ellipsoid = ellps;
94        self
95    }
96
97    /// `ModifiedAzimuthalEquidistant` projection constructor.
98    ///
99    /// To reduce computational overhead of projection functions this
100    /// constructor is non-trivial and tries to do as much projection computations as possible.
101    /// Thus creating a new structure can involve a significant computational overhead.
102    /// When projecting multiple coordinates only one instance of the structure should be created
103    /// and copied/borrowed as needed.
104    ///
105    /// # Errors
106    ///
107    /// Returns [`ProjectionError`] with additional information when:
108    ///
109    /// - one or more longitudes are not within -180..180 range.
110    /// - one or more latitudes are not within -90..90 range.
111    /// - one or more arguments are not finite.
112    pub fn initialize_projection(&self) -> Result<ModifiedAzimuthalEquidistant, ProjectionError> {
113        let ref_lon = unpack_required_parameter!(self, ref_lon);
114        let ref_lat = unpack_required_parameter!(self, ref_lat);
115        let ellps = self.ellipsoid;
116        ensure_finite!(ref_lon, ref_lat);
117
118        ensure_within_range!(ref_lon, -180.0..180.0);
119        ensure_within_range!(ref_lat, -90.0..90.0);
120
121        let lon_0 = ref_lon.to_radians();
122        let lat_0 = ref_lat.to_radians();
123
124        let n_1 = ellps.A / ellps.E.powi(2).mul_add(-(lat_0.sin()).powi(2), 1.0).sqrt();
125        let g = ellps.E * lat_0.sin() / ellps.E.mul_add(-ellps.E, 1.0).sqrt();
126
127        Ok(ModifiedAzimuthalEquidistant {
128            lon_0,
129            lat_0,
130            n_1,
131            g,
132            ellps,
133        })
134    }
135}
136
137impl Projection for ModifiedAzimuthalEquidistant {
138    #[inline]
139    #[allow(clippy::many_single_char_names)]
140    #[cfg_attr(feature = "tracing", instrument(level = "trace"))]
141    fn project_unchecked(&self, lon: f64, lat: f64) -> (f64, f64) {
142        let lon = lon.to_radians();
143        let lat = lat.to_radians();
144
145        let n = self.ellps.A / self.ellps.E.powi(2).mul_add(-(lat.sin()).powi(2), 1.0).sqrt();
146
147        let psi = self.ellps.E.mul_add(-self.ellps.E, 1.0).mul_add(lat.tan(), (self.ellps.E.powi(2) * self.n_1 * self.lat_0.sin()) / (n * lat.cos()))
148        .atan();
149
150        let az = ((lon - self.lon_0).sin())
151            .atan2(self.lat_0.cos().mul_add(psi.tan(), -(self.lat_0.sin() * (lon - self.lon_0).cos())));
152
153        let s = if approx_eq!(f64, az.sin(), 0.0) {
154            self.lat_0.cos().mul_add(psi.sin(), -(self.lat_0.sin() * psi.cos()))
155                .asin()
156                .abs()
157                * az.cos().signum()
158        } else {
159            (((lon - self.lon_0).sin() * psi.cos()) / (az.sin())).asin()
160        };
161
162        let h = self.ellps.E * self.lat_0.cos() * az.cos() / self.ellps.E.mul_add(-self.ellps.E, 1.0).sqrt();
163
164        let c = (self.n_1 * s)
165            * ((s.powi(5) / 48.0) * self.g).mul_add(-h, (s.powi(4) / 120.0).mul_add(h.powi(2).mul_add(7.0f64.mul_add(-h.powi(2), 4.0), -(3.0 * self.g.powi(2) * 7.0f64.mul_add(-h.powi(2), 1.0))), ((s.powi(3) / 8.0) * self.g * h).mul_add(2.0f64.mul_add(-h.powi(2), 1.0), 1.0 - (s.powi(2) * h.powi(2) * h.mul_add(-h, 1.0) / 6.0))));
166
167        let x = c * az.sin();
168        let y = c * az.cos();
169
170        (x, y)
171    }
172
173    #[inline]
174    #[cfg_attr(feature = "tracing", instrument(level = "trace"))]
175    fn inverse_project_unchecked(&self, x: f64, y: f64) -> (f64, f64) {
176        let c = x.hypot(y);
177        let az = x.atan2(y);
178
179        let big_a =
180            -self.ellps.E * self.ellps.E * ((self.lat_0.cos()).powi(2)) * ((az.cos()).powi(2))
181                / self.ellps.E.mul_add(-self.ellps.E, 1.0);
182        let big_b = 3.0
183            * self.ellps.E
184            * self.ellps.E
185            * (1.0 - big_a)
186            * self.lat_0.sin()
187            * self.lat_0.cos()
188            * az.cos()
189            / self.ellps.E.mul_add(-self.ellps.E, 1.0);
190        let big_d = c / self.n_1;
191        let big_e = big_d
192            - (big_a * (1.0 + big_a) * big_d.powi(3) / 6.0)
193            - (big_b * 3.0f64.mul_add(big_a, 1.0) * big_d.powi(4) / 24.0);
194        let big_f = 1.0 - (big_a * big_e * big_e / 2.0) - (big_b * big_e.powi(3) / 6.0);
195
196        let psi =
197            self.lat_0.sin().mul_add(big_e.cos(), self.lat_0.cos() * big_e.sin() * az.cos()).asin();
198
199        let lon = self.lon_0 + (az.sin() * big_e.sin() / psi.cos()).asin();
200        let lat = ((1.0 - (self.ellps.E * self.ellps.E * big_f * self.lat_0.sin() / psi.sin()))
201            * psi.tan()
202            / self.ellps.E.mul_add(-self.ellps.E, 1.0))
203            .atan();
204
205        let lon = lon.to_degrees();
206        let lat = lat.to_degrees();
207
208        (lon, lat)
209    }
210}
mappers/projections/modified_azimuthal_equidistant.rs

mappers/projections/
modified_azimuthal_equidistant.rs
