source: josm/trunk/src/org/openstreetmap/josm/data/projection/proj/LambertAzimuthalEqualArea.java@ 13632

Last change on this file since 13632 was 12013, checked in by bastiK, 7 years ago

see #11889 - backport improved version of Math.toDegrees and Math.toRadians from Java 9
  • Property svn:eol-style set to native
File size: 9.0 KB
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.data.projection.proj;
3
4import static org.openstreetmap.josm.tools.I18n.tr;
5
6import org.openstreetmap.josm.data.Bounds;
7import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
8import org.openstreetmap.josm.tools.Utils;
9
10/**
11 * Lambert Azimuthal Equal Area (EPSG code 9820).
12 * <p>
13 * This class has been derived from the implementation of the Geotools project;
14 * git 8cbf52d, org.geotools.referencing.operation.projection.LambertAzimuthalEqualArea
15 * at the time of migration.
16 * <p>
17 * <b>References:</b>
18 * <ul>
19 * <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li>
20 * <li> John P. Snyder (Map Projections - A Working Manual,
21 * U.S. Geological Survey Professional Paper 1395)</li>
22 * </ul>
23 *
24 * @author Gerald Evenden (for original code in Proj4)
25 * @author Beate Stollberg
26 * @author Martin Desruisseaux
27 *
28 * @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection</A>
29 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area"</A>
30 */
31public class LambertAzimuthalEqualArea extends AbstractProj {
32
33 /** Maximum difference allowed when comparing real numbers. */
34 private static final double EPSILON = 1E-7;
35
36 /** Epsilon for the comparison of small quantities. */
37 private static final double FINE_EPSILON = 1E-10;
38
39 /** Epsilon for the comparison of latitudes. */
40 private static final double EPSILON_LATITUDE = 1E-10;
41
42 /** Constants for authalic latitude. */
43 private static final double P00 = 0.33333333333333333333;
44 private static final double P01 = 0.17222222222222222222;
45 private static final double P02 = 0.10257936507936507936;
46 private static final double P10 = 0.06388888888888888888;
47 private static final double P11 = 0.06640211640211640211;
48 private static final double P20 = 0.01641501294219154443;
49
50 /** The projection mode. */
51 private enum Mode { OBLIQUE, EQUATORIAL, NORTH_POLE, SOUTH_POLE }
52
53 /** The projection mode for this particular instance. */
54 private Mode mode;
55
56 /** Constant parameters. */
57 private double sinb1, cosb1, xmf, ymf, qp, dd, rq;
58
59 /** Coefficients for authalic latitude. */
60 private double aPA0, aPA1, aPA2;
61
62 private double latitudeOfOrigin;
63
64 @Override
65 public String getName() {
66 return tr("Lambert Azimuthal Equal Area");
67 }
68
69 @Override
70 public String getProj4Id() {
71 return "laea";
72 }
73
74 @Override
75 public void initialize(ProjParameters params) throws ProjectionConfigurationException {
76 super.initialize(params);
77
78 if (params.lat0 == null)
79 throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", "lat_0"));
80
81 latitudeOfOrigin = Utils.toRadians(params.lat0);
82 /*
83 * Detects the mode (oblique, etc.).
84 */
85 final double t = Math.abs(latitudeOfOrigin);
86 if (Math.abs(t - Math.PI/2) < EPSILON_LATITUDE) {
87 mode = latitudeOfOrigin < 0.0 ? Mode.SOUTH_POLE : Mode.NORTH_POLE;
88 } else if (Math.abs(t) < EPSILON_LATITUDE) {
89 mode = Mode.EQUATORIAL;
90 } else {
91 mode = Mode.OBLIQUE;
92 }
93 /*
94 * Computes the constants for authalic latitude.
95 */
96 final double es2 = e2 * e2;
97 final double es3 = e2 * es2;
98 aPA0 = P02 * es3 + P01 * es2 + P00 * e2;
99 aPA1 = P11 * es3 + P10 * es2;
100 aPA2 = P20 * es3;
101
102 final double sinphi;
103 qp = qsfn(1);
104 rq = Math.sqrt(0.5 * qp);
105 sinphi = Math.sin(latitudeOfOrigin);
106 sinb1 = qsfn(sinphi) / qp;
107 cosb1 = Math.sqrt(1.0 - sinb1 * sinb1);
108 switch (mode) {
109 case NORTH_POLE: // Fall through
110 case SOUTH_POLE:
111 dd = 1.0;
112 xmf = ymf = rq;
113 break;
114 case EQUATORIAL:
115 dd = 1.0 / rq;
116 xmf = 1.0;
117 ymf = 0.5 * qp;
118 break;
119 case OBLIQUE:
120 dd = Math.cos(latitudeOfOrigin) / (Math.sqrt(1.0 - e2 * sinphi * sinphi) * rq * cosb1);
121 xmf = rq * dd;
122 ymf = rq / dd;
123 break;
124 default:
125 throw new AssertionError(mode);
126 }
127 }
128
129 @Override
130 public double[] project(final double phi, final double lambda) {
131 final double coslam = Math.cos(lambda);
132 final double sinlam = Math.sin(lambda);
133 final double sinphi = Math.sin(phi);
134 double q = qsfn(sinphi);
135 final double sinb, cosb, b, c, x, y;
136 switch (mode) {
137 case OBLIQUE:
138 sinb = q / qp;
139 cosb = Math.sqrt(1.0 - sinb * sinb);
140 c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam;
141 b = Math.sqrt(2.0 / c);
142 y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam);
143 x = xmf * b * cosb * sinlam;
144 break;
145 case EQUATORIAL:
146 sinb = q / qp;
147 cosb = Math.sqrt(1.0 - sinb * sinb);
148 c = 1.0 + cosb * coslam;
149 b = Math.sqrt(2.0 / c);
150 y = ymf * b * sinb;
151 x = xmf * b * cosb * sinlam;
152 break;
153 case NORTH_POLE:
154 c = (Math.PI / 2) + phi;
155 q = qp - q;
156 if (q >= 0.0) {
157 b = Math.sqrt(q);
158 x = b * sinlam;
159 y = coslam * -b;
160 } else {
161 x = y = 0.;
162 }
163 break;
164 case SOUTH_POLE:
165 c = phi - (Math.PI / 2);
166 q = qp + q;
167 if (q >= 0.0) {
168 b = Math.sqrt(q);
169 x = b * sinlam;
170 y = coslam * +b;
171 } else {
172 x = y = 0.;
173 }
174 break;
175 default:
176 throw new AssertionError(mode);
177 }
178 if (Math.abs(c) < EPSILON_LATITUDE) {
179 return new double[] {0, 0}; // this is an error, we should handle it somehow
180 }
181 return new double[] {x, y};
182 }
183
184 @Override
185 public double[] invproject(double x, double y) {
186 switch (mode) {
187 case EQUATORIAL: // Fall through
188 case OBLIQUE:
189 return invprojectEO(x, y);
190 case NORTH_POLE:
191 return invprojectNS(x, -y);
192 case SOUTH_POLE:
193 return invprojectNS(x, y);
194 default:
195 throw new AssertionError(mode);
196 }
197 }
198
199 private double[] invprojectEO(double x, double y) {
200 final double lambda;
201 final double phi;
202 x /= dd;
203 y *= dd;
204 final double rho = Math.hypot(x, y);
205 if (rho < FINE_EPSILON) {
206 lambda = 0.0;
207 phi = latitudeOfOrigin;
208 } else {
209 final double ab;
210 double sCe = 2.0 * Math.asin(0.5 * rho / rq);
211 double cCe = Math.cos(sCe);
212 sCe = Math.sin(sCe);
213 x *= sCe;
214 if (mode == Mode.OBLIQUE) {
215 ab = cCe * sinb1 + y * sCe * cosb1 / rho;
216 y = rho * cosb1 * cCe - y * sinb1 * sCe;
217 } else {
218 ab = y * sCe / rho;
219 y = rho * cCe;
220 }
221 lambda = Math.atan2(x, y);
222 phi = authlat(Math.asin(ab));
223 }
224 return new double[] {phi, lambda};
225 }
226
227 private double[] invprojectNS(double x, double y) {
228 final double lambda;
229 final double phi;
230 final double q = x*x + y*y;
231 if (q == 0) {
232 lambda = 0.;
233 phi = latitudeOfOrigin;
234 } else {
235 double ab = 1.0 - q / qp;
236 if (mode == Mode.SOUTH_POLE) {
237 ab = -ab;
238 }
239 lambda = Math.atan2(x, y);
240 phi = authlat(Math.asin(ab));
241 }
242 return new double[] {phi, lambda};
243 }
244
245 /**
246 * Calculates <var>q</var>, Snyder equation (3-12)
247 *
248 * @param sinphi sin of the latitude <var>q</var> is calculated for.
249 * @return <var>q</var> from Snyder equation (3-12).
250 */
251 private double qsfn(final double sinphi) {
252 if (e >= EPSILON) {
253 final double con = e * sinphi;
254 return (1.0 - e2) * (sinphi / (1.0 - con*con) -
255 (0.5 / e) * Math.log((1.0 - con) / (1.0 + con)));
256 } else {
257 return sinphi + sinphi;
258 }
259 }
260
261 /**
262 * Determines latitude from authalic latitude.
263 * @param beta authalic latitude
264 * @return corresponding latitude
265 */
266 private double authlat(final double beta) {
267 final double t = beta + beta;
268 return beta + aPA0 * Math.sin(t) + aPA1 * Math.sin(t+t) + aPA2 * Math.sin(t+t+t);
269 }
270
271 @Override
272 public Bounds getAlgorithmBounds() {
273 return new Bounds(-89, -174, 89, 174, false);
274 }
275}
Note: See TracBrowser for help on using the repository browser.