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

Last change on this file since 11601 was 11601, checked in by Don-vip, 8 years ago

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