1 | // License: GPL. For details, see LICENSE file.
|
---|
2 | package org.openstreetmap.josm.data.projection.proj;
|
---|
3 |
|
---|
4 | import static org.openstreetmap.josm.tools.I18n.tr;
|
---|
5 |
|
---|
6 | import org.openstreetmap.josm.data.Bounds;
|
---|
7 | import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
|
---|
8 | import org.openstreetmap.josm.tools.JosmRuntimeException;
|
---|
9 |
|
---|
10 | /**
|
---|
11 | * Azimuthal Equidistant projection.
|
---|
12 | * <p>
|
---|
13 | * This implementation does not include the Guam or Micronesia variants.
|
---|
14 | *
|
---|
15 | * @author Gerald Evenden (original PROJ.4 implementation in C)
|
---|
16 | * @author Ben Caradoc-Davies (Transient Software Limited)
|
---|
17 | * @see <a href="https://pubs.er.usgs.gov/publication/pp1395"><em>Map Projections: A Working Manual</em>, Snyder (1987), pages 191-202</a>
|
---|
18 | * @see <a href="http://geotiff.maptools.org/proj_list/azimuthal_equidistant.html">PROJ.4 notes on parameters</a>
|
---|
19 | * @see <a href="https://github.com/OSGeo/proj.4/blob/master/src/PJ_aeqd.c">PROJ.4 implemention in C</a>
|
---|
20 | * @see <a href="https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection">Wikipedia</a>
|
---|
21 | * @see <a href="http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html">Wolfram Alpha</a>
|
---|
22 | * @since 13598
|
---|
23 | */
|
---|
24 | public class AzimuthalEquidistant extends AbstractProj {
|
---|
25 |
|
---|
26 | /**
|
---|
27 | * Less strict tolerance.
|
---|
28 | */
|
---|
29 | public static final double EPS10 = 1.e-10;
|
---|
30 |
|
---|
31 | /**
|
---|
32 | * Stricter tolerance.
|
---|
33 | */
|
---|
34 | public static final double TOL = 1.e-14;
|
---|
35 |
|
---|
36 | /**
|
---|
37 | * Half of π.
|
---|
38 | */
|
---|
39 | public static final double HALF_PI = Math.PI / 2;
|
---|
40 |
|
---|
41 | /**
|
---|
42 | * The four possible modes or aspects of the projection.
|
---|
43 | */
|
---|
44 | public enum Mode {
|
---|
45 | /** North pole */
|
---|
46 | NORTH_POLAR,
|
---|
47 | /** South pole */
|
---|
48 | SOUTH_POLAR,
|
---|
49 | /** Equator */
|
---|
50 | EQUATORIAL,
|
---|
51 | /** Oblique */
|
---|
52 | OBLIQUE;
|
---|
53 | }
|
---|
54 |
|
---|
55 | /**
|
---|
56 | * Length of semi-major axis, in metres. This is named '<var>a</var>' or '<var>R</var>'
|
---|
57 | * (Radius in spherical cases) in Snyder.
|
---|
58 | */
|
---|
59 | protected double semiMajor;
|
---|
60 |
|
---|
61 | /**
|
---|
62 | * Length of semi-minor axis, in metres. This is named '<var>b</var>' in Snyder.
|
---|
63 | */
|
---|
64 | protected double semiMinor;
|
---|
65 |
|
---|
66 | /**
|
---|
67 | * Central longitude in <u>radians</u>. Default value is 0, the Greenwich meridian.
|
---|
68 | * This is called '<var>lambda0</var>' in Snyder.
|
---|
69 | */
|
---|
70 | protected double centralMeridian;
|
---|
71 |
|
---|
72 | /**
|
---|
73 | * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
|
---|
74 | * This is called '<var>phi0</var>' in Snyder.
|
---|
75 | */
|
---|
76 | protected double latitudeOfOrigin;
|
---|
77 |
|
---|
78 | /**
|
---|
79 | * The mode or aspect of the projection.
|
---|
80 | */
|
---|
81 | protected Mode mode;
|
---|
82 |
|
---|
83 | /**
|
---|
84 | * Geodesic calculator used for this projection. Not used and set to null for polar projections.
|
---|
85 | */
|
---|
86 | //protected Geodesic geodesic; // See https://josm.openstreetmap.de/ticket/16129#comment:21
|
---|
87 |
|
---|
88 | /**
|
---|
89 | * The sine of the central latitude of the projection.
|
---|
90 | */
|
---|
91 | protected double sinph0;
|
---|
92 |
|
---|
93 | /**
|
---|
94 | * The cosine of the central latitude of the projection.
|
---|
95 | */
|
---|
96 | protected double cosph0;
|
---|
97 |
|
---|
98 | /**
|
---|
99 | * Meridian distance from the equator to the pole. Not used and set to NaN for non-polar projections.
|
---|
100 | */
|
---|
101 | protected double mp;
|
---|
102 |
|
---|
103 | @Override
|
---|
104 | public String getName() {
|
---|
105 | return tr("Azimuthal Equidistant");
|
---|
106 | }
|
---|
107 |
|
---|
108 | @Override
|
---|
109 | public String getProj4Id() {
|
---|
110 | return "aeqd";
|
---|
111 | }
|
---|
112 |
|
---|
113 | @Override
|
---|
114 | public void initialize(ProjParameters params) throws ProjectionConfigurationException {
|
---|
115 | super.initialize(params);
|
---|
116 | if (params.lon0 == null)
|
---|
117 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lon_0"));
|
---|
118 | if (params.lat0 == null)
|
---|
119 | throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
|
---|
120 | centralMeridian = Math.toRadians(params.lon0);
|
---|
121 | latitudeOfOrigin = Math.toRadians(params.lat0);
|
---|
122 | semiMajor = params.ellps.a;
|
---|
123 | semiMinor = params.ellps.b;
|
---|
124 | if (Math.abs(latitudeOfOrigin - HALF_PI) < EPS10) {
|
---|
125 | mode = Mode.NORTH_POLAR;
|
---|
126 | mp = mlfn(HALF_PI, 1, 0);
|
---|
127 | sinph0 = 1;
|
---|
128 | cosph0 = 0;
|
---|
129 | } else if (Math.abs(latitudeOfOrigin + HALF_PI) < EPS10) {
|
---|
130 | mode = Mode.SOUTH_POLAR;
|
---|
131 | mp = mlfn(-HALF_PI, -1, 0);
|
---|
132 | sinph0 = -1;
|
---|
133 | cosph0 = 0;
|
---|
134 | } else if (Math.abs(latitudeOfOrigin) < EPS10) {
|
---|
135 | mode = Mode.EQUATORIAL;
|
---|
136 | mp = Double.NaN;
|
---|
137 | sinph0 = 0;
|
---|
138 | cosph0 = 1;
|
---|
139 | //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor);
|
---|
140 | throw new ProjectionConfigurationException("Equatorial AzimuthalEquidistant not yet supported");
|
---|
141 | } else {
|
---|
142 | mode = Mode.OBLIQUE;
|
---|
143 | mp = Double.NaN;
|
---|
144 | sinph0 = Math.sin(latitudeOfOrigin);
|
---|
145 | cosph0 = Math.cos(latitudeOfOrigin);
|
---|
146 | //geodesic = new Geodesic(semiMajor, (semiMajor - semiMinor) / semiMajor);
|
---|
147 | throw new ProjectionConfigurationException("Oblique AzimuthalEquidistant not yet supported");
|
---|
148 | }
|
---|
149 | }
|
---|
150 |
|
---|
151 | @Override
|
---|
152 | public Bounds getAlgorithmBounds() {
|
---|
153 | return new Bounds(-89, -180, 89, 180, false);
|
---|
154 | }
|
---|
155 |
|
---|
156 | @Override
|
---|
157 | public double[] project(double latRad, double lonRad) {
|
---|
158 | return spherical ? projectSpherical(latRad, lonRad) : projectEllipsoidal(latRad, lonRad);
|
---|
159 | }
|
---|
160 |
|
---|
161 | @Override
|
---|
162 | public double[] invproject(double east, double north) {
|
---|
163 | return spherical ? invprojectSpherical(east, north) : invprojectEllipsoidal(east, north);
|
---|
164 | }
|
---|
165 |
|
---|
166 | double[] projectSpherical(double latRad, double lonRad) {
|
---|
167 | double x = 0;
|
---|
168 | double y = 0;
|
---|
169 | double sinphi = Math.sin(latRad);
|
---|
170 | double cosphi = Math.cos(latRad);
|
---|
171 | double coslam = Math.cos(lonRad);
|
---|
172 | switch (mode) {
|
---|
173 | case EQUATORIAL:
|
---|
174 | case OBLIQUE:
|
---|
175 | if (mode == Mode.EQUATORIAL) {
|
---|
176 | y = cosphi * coslam;
|
---|
177 | } else { // Oblique
|
---|
178 | y = sinph0 * sinphi + cosph0 * cosphi * coslam;
|
---|
179 | }
|
---|
180 | if (Math.abs(Math.abs(y) - 1) < TOL) {
|
---|
181 | if (y < 0) {
|
---|
182 | throw new JosmRuntimeException("TOLERANCE_ERROR");
|
---|
183 | } else {
|
---|
184 | x = 0;
|
---|
185 | y = 0;
|
---|
186 | }
|
---|
187 | } else {
|
---|
188 | y = Math.acos(y);
|
---|
189 | y /= Math.sin(y);
|
---|
190 | x = y * cosphi * Math.sin(lonRad);
|
---|
191 | y *= (mode == Mode.EQUATORIAL) ? sinphi
|
---|
192 | : (cosph0 * sinphi - sinph0 * cosphi * coslam);
|
---|
193 | }
|
---|
194 | break;
|
---|
195 | case NORTH_POLAR:
|
---|
196 | latRad = -latRad;
|
---|
197 | coslam = -coslam;
|
---|
198 | // fall through
|
---|
199 | case SOUTH_POLAR:
|
---|
200 | if (Math.abs(latRad - HALF_PI) < EPS10) {
|
---|
201 | throw new JosmRuntimeException("TOLERANCE_ERROR");
|
---|
202 | }
|
---|
203 | y = HALF_PI + latRad;
|
---|
204 | x = y * Math.sin(lonRad);
|
---|
205 | y *= coslam;
|
---|
206 | break;
|
---|
207 | }
|
---|
208 | return new double[] {x, y};
|
---|
209 | }
|
---|
210 |
|
---|
211 | double[] invprojectSpherical(double east, double north) {
|
---|
212 | double x = east;
|
---|
213 | double y = north;
|
---|
214 | double lambda = 0;
|
---|
215 | double phi = 0;
|
---|
216 | double c_rh = Math.hypot(x, y);
|
---|
217 | if (c_rh > Math.PI) {
|
---|
218 | if (c_rh - EPS10 > Math.PI) {
|
---|
219 | throw new JosmRuntimeException("TOLERANCE_ERROR");
|
---|
220 | }
|
---|
221 | c_rh = Math.PI;
|
---|
222 | } else if (c_rh < EPS10) {
|
---|
223 | phi = latitudeOfOrigin;
|
---|
224 | lambda = 0.;
|
---|
225 | } else {
|
---|
226 | if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) {
|
---|
227 | double sinc = Math.sin(c_rh);
|
---|
228 | double cosc = Math.cos(c_rh);
|
---|
229 | if (mode == Mode.EQUATORIAL) {
|
---|
230 | phi = aasin(y * sinc / c_rh);
|
---|
231 | x *= sinc;
|
---|
232 | y = cosc * c_rh;
|
---|
233 | } else { // Oblique
|
---|
234 | phi = aasin(cosc * sinph0 + y * sinc * cosph0 / c_rh);
|
---|
235 | y = (cosc - sinph0 * Math.sin(phi)) * c_rh;
|
---|
236 | x *= sinc * cosph0;
|
---|
237 | }
|
---|
238 | lambda = (y == 0) ? 0 : Math.atan2(x, y);
|
---|
239 | } else if (mode == Mode.NORTH_POLAR) {
|
---|
240 | phi = HALF_PI - c_rh;
|
---|
241 | lambda = Math.atan2(x, -y);
|
---|
242 | } else { // South Polar
|
---|
243 | phi = c_rh - HALF_PI;
|
---|
244 | lambda = Math.atan2(x, y);
|
---|
245 | }
|
---|
246 | }
|
---|
247 | return new double[] {phi, lambda};
|
---|
248 | }
|
---|
249 |
|
---|
250 | double[] projectEllipsoidal(double latRad, double lonRad) {
|
---|
251 | double x = 0;
|
---|
252 | double y = 0;
|
---|
253 | double coslam = Math.cos(lonRad);
|
---|
254 | double cosphi = Math.cos(latRad);
|
---|
255 | double sinphi = Math.sin(latRad);
|
---|
256 | switch (mode) {
|
---|
257 | case NORTH_POLAR:
|
---|
258 | coslam = -coslam;
|
---|
259 | // fall through
|
---|
260 | case SOUTH_POLAR:
|
---|
261 | double rho = Math.abs(mp - mlfn(latRad, sinphi, cosphi));
|
---|
262 | x = rho * Math.sin(lonRad);
|
---|
263 | y = rho * coslam;
|
---|
264 | break;
|
---|
265 | case EQUATORIAL:
|
---|
266 | case OBLIQUE:
|
---|
267 | if (Math.abs(lonRad) < EPS10 && Math.abs(latRad - latitudeOfOrigin) < EPS10) {
|
---|
268 | x = 0;
|
---|
269 | y = 0;
|
---|
270 | break;
|
---|
271 | }
|
---|
272 | /*GeodesicData g = geodesic.Inverse(Math.toDegrees(latitudeOfOrigin),
|
---|
273 | Math.toDegrees(centralMeridian), Math.toDegrees(latRad),
|
---|
274 | Math.toDegrees(lonRad + centralMeridian));
|
---|
275 | double azi1 = Math.toRadians(g.azi1);
|
---|
276 | x = g.s12 * Math.sin(azi1) / semiMajor;
|
---|
277 | y = g.s12 * Math.cos(azi1) / semiMajor;*/
|
---|
278 | break;
|
---|
279 | }
|
---|
280 | return new double[] {x, y};
|
---|
281 | }
|
---|
282 |
|
---|
283 | double[] invprojectEllipsoidal(double east, double north) {
|
---|
284 | double x = east;
|
---|
285 | double y = north;
|
---|
286 | double lambda = 0;
|
---|
287 | double phi = 0;
|
---|
288 | double c = Math.hypot(x, y);
|
---|
289 | if (c < EPS10) {
|
---|
290 | phi = latitudeOfOrigin;
|
---|
291 | lambda = 0;
|
---|
292 | } else {
|
---|
293 | if (mode == Mode.OBLIQUE || mode == Mode.EQUATORIAL) {
|
---|
294 | /*double x2 = x * semiMajor;
|
---|
295 | double y2 = y * semiMajor;
|
---|
296 | double azi1 = Math.atan2(x2, y2);
|
---|
297 | double s12 = Math.sqrt(x2 * x2 + y2 * y2);
|
---|
298 | GeodesicData g = geodesic.Direct(Math.toDegrees(latitudeOfOrigin),
|
---|
299 | Math.toDegrees(centralMeridian), Math.toDegrees(azi1), s12);
|
---|
300 | phi = Math.toRadians(g.lat2);
|
---|
301 | lambda = Math.toRadians(g.lon2);*/
|
---|
302 | lambda -= centralMeridian;
|
---|
303 | } else { // Polar
|
---|
304 | phi = invMlfn((mode == Mode.NORTH_POLAR) ? (mp - c) : (mp + c));
|
---|
305 | lambda = Math.atan2(x, (mode == Mode.NORTH_POLAR) ? -y : y);
|
---|
306 | }
|
---|
307 | }
|
---|
308 | return new double[] {phi, lambda};
|
---|
309 | }
|
---|
310 | }
|
---|