Peano
EllipsoidProjection.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "../Constants.h"
4 
5 #include <cmath>
6 #include <utility>
7 
8 namespace projection {
9  namespace albers_equal_area {
10  namespace ellipsoid {
11  namespace constants = ::applications::exahype2::swe;
12  /*
13  * @brief
14  * Projects geographic coordinates (longitude, latitude) onto a planar surface
15  * using the Albers Equal-Area Conic projection for the WGS_1984 ellipsoid.
16  *
17  * This implementation is based on the equations provided in USGS Professional Paper 1395.
18  * It supports standard parallels and is suitable for regional-scale mapping where
19  * area preservation is important.
20  *
21  * @see https://doi.org/10.3133/pp1395
22  */
23  inline auto forward(const auto lon, const auto lat) {
24  constexpr auto a = constants::aeaEllipsoidA; // Halbachse des WGS_1984-Ellipsoids (Meter)
25  constexpr auto e = constants::aeaEllipsoidE;
26  constexpr auto stdParallel1 = constants::aeaStdParallel1; // Standard Parallel 1 (Grad)
27  constexpr auto stdParallel2 = constants::aeaStdParallel2; // Standard Parallel 2 (Grad)
28  constexpr auto lon0 = constants::aeaLon0; // Central Meridian (Grad)
29  constexpr auto lat0 = constants::aeaLat0; // Latitude of Origin (Grad)
30 
31  constexpr auto toRad = M_PI / 180.0;
32 
33  constexpr auto stdParallel1Rad = stdParallel1 * toRad;
34  constexpr auto stdParallel2Rad = stdParallel2 * toRad;
35  constexpr auto lat0Rad = lat0 * toRad;
36  constexpr auto lon0Rad = lon0 * toRad;
37 
38  const auto lonRad = lon * toRad;
39  const auto latRad = lat * toRad;
40 
41  constexpr auto e2 = e * e;
42  constexpr auto oneMinusE2 = 1.0 - (e * e);
43 
44  auto m = [](const auto phi) { return std::cos(phi) / std::sqrt(1.0 - (e2 * std::sin(phi) * std::sin(phi))); };
45 
46  static const auto m1 = m(stdParallel1Rad);
47  static const auto m2 = m(stdParallel2Rad);
48 
49  auto q = [](const auto phi) {
50  const auto sinPhi = std::sin(phi);
51  const auto sinPhi2 = sinPhi * sinPhi;
52  return oneMinusE2
53  * ((std::sin(phi) / (1.0 - (e2 * sinPhi2))) - ((1.0 / (2.0 * e)) * std::log((1.0 - (e * std::sin(phi))) / (1.0 + (e * std::sin(phi))))));
54  };
55 
56  static const auto q0 = q(lat0Rad);
57  static const auto q1 = q(stdParallel1Rad);
58  static const auto q2 = q(stdParallel2Rad);
59 
60 
61  const auto qLat = q(latRad);
62 
63  static const auto n = ((m1 * m1) - (m2 * m2)) / (q2 - q1);
64  static const auto C = (m1 * m1) + (n * q1);
65  static const auto p0 = a * std::sqrt(C - (n * q0)) / n;
66 
67  const auto p = a * std::sqrt(C - (n * qLat)) / n;
68  const auto theta = n * (lonRad - lon0Rad);
69 
70  const auto x = p * std::sin(theta);
71  const auto y = p0 - (p * std::cos(theta));
72 
73  return std::pair{x, y};
74  };
75 
86  inline auto inverse(const auto x, const auto y) {
87  constexpr size_t maxIter = 10;
88  constexpr auto tol = 1e-12;
89 
90  constexpr auto a = constants::aeaEllipsoidA; // Halbachse des WGS_1984-Ellipsoids (Meter)
91  constexpr auto e = constants::aeaEllipsoidE;
92  constexpr auto stdParallel1 = constants::aeaStdParallel1; // Standard Parallel 1 (Grad)
93  constexpr auto stdParallel2 = constants::aeaStdParallel2; // Standard Parallel 2 (Grad)
94  constexpr auto lon0 = constants::aeaLon0; // Central Meridian (Grad)
95  constexpr auto lat0 = constants::aeaLat0; // Latitude of Origin (Grad)
96 
97  constexpr auto toRad = M_PI / 180.0;
98 
99  constexpr auto stdParallel1Rad = stdParallel1 * toRad;
100  constexpr auto stdParallel2Rad = stdParallel2 * toRad;
101  constexpr auto lat0Rad = lat0 * toRad;
102  constexpr auto lon0Rad = lon0 * toRad;
103 
104  constexpr auto e2 = e * e;
105  constexpr auto oneMinusE2 = 1.0 - (e * e);
106 
107  auto m = [](const auto phi) { return std::cos(phi) / std::sqrt(1.0 - (e2 * std::sin(phi) * std::sin(phi))); };
108 
109  static const auto m1 = m(stdParallel1Rad);
110  static const auto m2 = m(stdParallel2Rad);
111 
112  auto q = [](const auto phi) {
113  const auto sinPhi = std::sin(phi);
114  const auto sinPhi2 = sinPhi * sinPhi;
115  return oneMinusE2
116  * ((std::sin(phi) / (1.0 - (e2 * sinPhi2))) - ((1.0 / (2.0 * e)) * std::log((1.0 - (e * std::sin(phi))) / (1.0 + (e * std::sin(phi))))));
117  };
118 
119  static const auto q0 = q(lat0Rad);
120  static const auto q1 = q(stdParallel1Rad);
121  static const auto q2 = q(stdParallel2Rad);
122 
123  static const auto n = ((m1 * m1) - (m2 * m2)) / (q2 - q1);
124  static const auto C = (m1 * m1) + (n * q1);
125  static const auto p0 = a * std::sqrt(C - (n * q0)) / n;
126 
127  const auto p = std::sqrt((x * x) + ((p0 - y) * (p0 - y)));
128  const auto qnew = (C - (p * p * n * n / (a * a))) / n;
129 
130  auto phi = std::asin(qnew / 2);
131 
132 #pragma unroll
133  for (int i = 0; i < maxIter; ++i) {
134  const auto sinPhi = std::sin(phi);
135  const auto cosPhi = std::cos(phi);
136  const auto sinPhi2 = sinPhi * sinPhi;
137  const auto denom = 1.0 - (e2 * sinPhi2);
138 
139  const auto num1 = qnew / oneMinusE2;
140  const auto num2 = sinPhi / denom;
141  const auto num3 = (1.0 / (2.0 * e)) * std::log((1.0 - e * sinPhi) / (1.0 + e * sinPhi));
142 
143  const auto deltaPhi = (1.0 - (e2 * sinPhi2)) * (1.0 - (e2 * sinPhi2)) / (2.0 * cosPhi) * (num1 - num2 + num3);
144  const auto phiNew = phi + deltaPhi;
145 
146  if (std::abs(phiNew - phi) < tol) {
147  break;
148  }
149  phi = phiNew;
150  }
151 
152 
153  const auto theta = std::atan2(x, p0 - y);
154  const auto lon = lon0Rad + (theta / n);
155 
156  return std::pair{
157  lon / toRad,
158  phi / toRad,
159  };
160  };
161  } // namespace ellipsoid
162  } // namespace albers_equal_area
163 } // namespace projection
csv_to_plot.y
list y
Definition: csv_to_plot.py:98
projection::albers_equal_area::ellipsoid::inverse
auto inverse(const auto x, const auto y)
Inverts the Albers Equal-Area Conic projection for WGS_1984 and returns the geographic coordinates (l...
Definition: EllipsoidProjection.h:86
projection::albers_equal_area::ellipsoid::forward
auto forward(const auto lon, const auto lat)
Definition: EllipsoidProjection.h:23