mappers/projections/azimuthal_equidistant.rs
1//! The azimuthal equidistant projection is an azimuthal map projection.
2//! It has the useful properties that all points on the map are at proportionally
3//! correct distances from the center point, and that all points on the map are at the
4//! correct azimuth (direction) from the center point. A useful application for this
5//! type of projection is a polar projection which shows all meridians (lines of longitude) as straight,
6//! with distances from the pole represented correctly [(Wikipedia, 2022)](https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection).
7//!
8//! This projection uses Geodesic computation (defined by [C. F. F. Karney (2013)](https://doi.org/10.1007/s00190-012-0578-z))
9//! to compute distances and azimuths between projected point and origin. So it might be slower than some other projections.
10//!
11//! Summary by [Snyder (1987)](https://pubs.er.usgs.gov/publication/pp1395):
12//!
13//! - Azimuthal.
14//! - Distances measured from the center are true.
15//! - Distances not measured along radii from the center are not correct.
16//! - The center of projection is the only point without distortion.
17//! - Directions from the center are true (except on some oblique and equatorial ellipsoidal forms).
18//! - Neither equal-area nor conformal.
19//! - All meridians on the polar aspect, the central meridian on other aspects, and the Equator on the equatorial aspect are straight lines.
20//! - Parallels on the polar projection are circles spaced at true intervals (equidistant for the sphere).
21//! - The outer meridian of a hemisphere on the equatorial aspect (for the sphere) is a circle.
22//! - All other meridians and parallels are complex curves.
23//! - Not a perspective projection.
24//! - Point opposite the center is shown as a circle (for the sphere) surrounding the map.
25//! - Used in the polar aspect for world maps and maps of polar hemispheres.
26//! - Used in the oblique aspect for atlas maps of continents and world maps for aviation and radio use.
27//! - Known for many centuries in the polar aspect.
28
29use crate::{
30 Ellipsoid, Projection, ProjectionError,
31 errors::{ensure_finite, ensure_within_range, unpack_required_parameter},
32};
33use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic};
34
35#[cfg(feature = "tracing")]
36use tracing::instrument;
37
38/// Main projection struct that is constructed from [`AzimuthalEquidistantBuilder`] and used for computations.
39#[derive(Copy, Clone, PartialEq, PartialOrd, Debug)]
40pub struct AzimuthalEquidistant {
41 lon_0: f64,
42 lat_0: f64,
43 geod: Geodesic,
44}
45
46impl AzimuthalEquidistant {
47 /// Initializes builder with default values.
48 /// Projection parameters can be set with builder methods,
49 /// refer to the documentation of those methods to check which parmeters are required
50 /// and default values for optional arguments.
51 #[must_use]
52 pub fn builder() -> AzimuthalEquidistantBuilder {
53 AzimuthalEquidistantBuilder::default()
54 }
55}
56
57/// Builder struct which allows to construct [`AzimuthalEquidistant`] projection.
58/// Refer to the documentation of this struct's methods to check which parmeters are required
59/// and default values for optional arguments.
60#[derive(Copy, Clone, PartialEq, PartialOrd, Debug)]
61pub struct AzimuthalEquidistantBuilder {
62 ref_lon: Option<f64>,
63 ref_lat: Option<f64>,
64 ellipsoid: Ellipsoid,
65}
66
67impl Default for AzimuthalEquidistantBuilder {
68 fn default() -> Self {
69 Self {
70 ref_lon: None,
71 ref_lat: None,
72 ellipsoid: Ellipsoid::WGS84,
73 }
74 }
75}
76
77impl AzimuthalEquidistantBuilder {
78 /// *(required)* Sets reference longitude and latitude. Point (0, 0) on the map will be at this coordinates.
79 pub const fn ref_lonlat(&mut self, lon: f64, lat: f64) -> &mut Self {
80 self.ref_lon = Some(lon);
81 self.ref_lat = Some(lat);
82 self
83 }
84
85 /// *(optional)* Sets reference [`Ellipsoid`], defaults to [`WGS84`](Ellipsoid::WGS84).
86 pub const fn ellipsoid(&mut self, ellps: Ellipsoid) -> &mut Self {
87 self.ellipsoid = ellps;
88 self
89 }
90
91 /// AEQD projection constructor.
92 ///
93 /// To reduce computational overhead of projection functions this
94 /// constructor is non-trivial and tries to do as much projection computations as possible.
95 /// Thus creating a new structure can involve a significant computational overhead.
96 /// When projecting multiple coordinates only one instance of the structure should be created
97 /// and copied/borrowed as needed.
98 ///
99 ///
100 /// # Errors
101 ///
102 /// Returns [`ProjectionError`] with additional information when:
103 ///
104 /// - one or more longitudes are not within -180..180 range.
105 /// - one or more latitudes are not within -90..90 range.
106 /// - one or more arguments are not finite.
107 pub fn initialize_projection(&self) -> Result<AzimuthalEquidistant, ProjectionError> {
108 let ref_lon = unpack_required_parameter!(self, ref_lon);
109 let ref_lat = unpack_required_parameter!(self, ref_lat);
110 let ellps = self.ellipsoid;
111 ensure_finite!(ref_lon, ref_lat);
112
113 ensure_within_range!(ref_lon, -180.0..180.0);
114 ensure_within_range!(ref_lat, -90.0..90.0);
115
116 Ok(AzimuthalEquidistant {
117 lon_0: ref_lon,
118 lat_0: ref_lat,
119 geod: ellps.into(),
120 })
121 }
122}
123
124impl Projection for AzimuthalEquidistant {
125 #[cfg_attr(feature = "tracing", instrument(level = "trace"))]
126 fn project_unchecked(&self, lon: f64, lat: f64) -> (f64, f64) {
127 let (s12, azi1, _, _) = self.geod.inverse(self.lat_0, self.lon_0, lat, lon);
128
129 let x = s12 * azi1.to_radians().sin();
130 let y = s12 * azi1.to_radians().cos();
131
132 (x, y)
133 }
134
135 #[cfg_attr(feature = "tracing", instrument(level = "trace"))]
136 fn inverse_project_unchecked(&self, x: f64, y: f64) -> (f64, f64) {
137 let azi1 = x.atan2(y).to_degrees();
138 let s12 = x.hypot(y);
139
140 let (lat, lon) = self.geod.direct(self.lat_0, self.lon_0, azi1, s12);
141
142 (lon, lat)
143 }
144}