source: osm/applications/editors/josm/plugins/opendata/includes/org/geotools/referencing/operation/projection/LambertAzimuthalEqualArea.java@ 28055

Last change on this file since 28055 was 28055, checked in by donvip, 12 years ago

opendata: download improvements
File size: 22.5 KB
Line 
1/*
2 * GeoTools - The Open Source Java GIS Toolkit
3 * http://geotools.org
4 *
5 * (C) 2006-2008, Open Source Geospatial Foundation (OSGeo)
6 *
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation;
10 * version 2.1 of the License.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * This package contains formulas from the PROJ package of USGS.
18 * USGS's work is fully acknowledged here. This derived work has
19 * been relicensed under LGPL with Frank Warmerdam's permission.
20 */
21package org.geotools.referencing.operation.projection;
22
23import java.awt.geom.Point2D;
24import java.util.Collection;
25import javax.measure.unit.NonSI;
26import org.opengis.parameter.GeneralParameterDescriptor;
27import org.opengis.parameter.ParameterDescriptor;
28import org.opengis.parameter.ParameterDescriptorGroup;
29import org.opengis.parameter.ParameterNotFoundException;
30import org.opengis.parameter.ParameterValueGroup;
31import org.opengis.referencing.operation.MathTransform;
32import org.geotools.metadata.iso.citation.Citations;
33import org.geotools.referencing.NamedIdentifier;
34import org.geotools.resources.i18n.ErrorKeys;
35
36import static java.lang.Math.*;
37
38
39/**
40 * Lambert Azimuthal Equal Area (EPSG code 9820).
41 * <p>
42 * <b>References:</b>
43 * <ul>
44 * <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li>
45 * <li> John P. Snyder (Map Projections - A Working Manual,
46 * U.S. Geological Survey Professional Paper 1395)</li>
47 * </ul>
48 *
49 * @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection on MathWorld</A>
50 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area" on RemoteSensing.org</A>
51 *
52 * @since 2.4
53 * @version $Id: LambertAzimuthalEqualArea.java 37299 2011-05-25 05:21:24Z mbedward $
54 *
55 * @source $URL: http://svn.osgeo.org/geotools/branches/2.7.x/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/LambertAzimuthalEqualArea.java $
56 * @author Gerald Evenden (for original code in Proj4)
57 * @author Beate Stollberg
58 * @author Martin Desruisseaux
59 */
60public class LambertAzimuthalEqualArea extends MapProjection {
61 /** For cross-version compatibility. */
62 private static final long serialVersionUID = 1639914708790574760L;
63
64 /** Maximum difference allowed when comparing real numbers. */
65 private static final double EPSILON = 1E-7;
66
67 /** Epsilon for the comparaison of small quantities. */
68 private static final double FINE_EPSILON = 1E-10;
69
70 /** Epsilon for the comparaison of latitudes. */
71 private static final double EPSILON_LATITUDE = 1E-10;
72
73 /** Constants for authalic latitude. */
74 private static final double P00 = 0.33333333333333333333,
75 P01 = 0.17222222222222222222,
76 P02 = 0.10257936507936507936,
77 P10 = 0.06388888888888888888,
78 P11 = 0.06640211640211640211,
79 P20 = 0.01641501294219154443;
80
81 /** The projection mode. */
82 static final int OBLIQUE=0, EQUATORIAL=1, NORTH_POLE=2, SOUTH_POLE=3;
83
84 /** The projection mode for this particular instance. */
85 final int mode;
86
87 /** Constant parameters. */
88 final double sinb1, cosb1, xmf, ymf, mmf, qp, dd, rq;
89
90 /** Coefficients for authalic latitude. */
91 private final double APA0, APA1, APA2;
92
93 /**
94 * Constructs a new map projection from the supplied parameters.
95 *
96 * @param parameters The parameter values in standard units.
97 * @throws ParameterNotFoundException if a mandatory parameter is missing.
98 */
99 protected LambertAzimuthalEqualArea(final ParameterValueGroup parameters)
100 throws ParameterNotFoundException
101 {
102 // Fetch parameters
103 super(parameters);
104 final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
105 latitudeOfOrigin = doubleValue(expected, Provider.LATITUDE_OF_CENTRE, parameters);
106 centralMeridian = doubleValue(expected, Provider.LONGITUDE_OF_CENTRE, parameters);
107 ensureLatitudeInRange (Provider.LATITUDE_OF_CENTRE, latitudeOfOrigin, true);
108 ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTRE, centralMeridian, true);
109 /*
110 * Detects the mode (oblique, etc.).
111 */
112 final double t = abs(latitudeOfOrigin);
113 if (abs(t - PI/2) < EPSILON_LATITUDE) {
114 mode = latitudeOfOrigin < 0.0 ? SOUTH_POLE : NORTH_POLE;
115 } else if (abs(t) < EPSILON_LATITUDE) {
116 mode = EQUATORIAL;
117 } else {
118 mode = OBLIQUE;
119 }
120 /*
121 * Computes the constants for authalic latitude.
122 */
123 final double es2 = excentricitySquared * excentricitySquared;
124 final double es3 = excentricitySquared * es2;
125 APA0 = P02 * es3 + P01 * es2 + P00 * excentricitySquared;
126 APA1 = P11 * es3 + P10 * es2;
127 APA2 = P20 * es3;
128
129 final double sinphi;
130 qp = qsfn(1);
131 rq = sqrt(0.5 * qp);
132 mmf = 0.5 / (1 - excentricitySquared);
133 sinphi = sin(latitudeOfOrigin);
134 if (isSpherical) {
135 sinb1 = sin(latitudeOfOrigin);
136 cosb1 = cos(latitudeOfOrigin);
137 } else {
138 sinb1 = qsfn(sinphi) / qp;
139 cosb1 = sqrt(1.0 - sinb1 * sinb1);
140 }
141 switch (mode) {
142 case NORTH_POLE: // Fall through
143 case SOUTH_POLE: {
144 dd = 1.0;
145 xmf = ymf = rq;
146 break;
147 }
148 case EQUATORIAL: {
149 dd = 1.0 / rq;
150 xmf = 1.0;
151 ymf = 0.5 * qp;
152 break;
153 }
154 case OBLIQUE: {
155 dd = cos(latitudeOfOrigin) /
156 (sqrt(1.0 - excentricitySquared * sinphi * sinphi) * rq * cosb1);
157 xmf = rq * dd;
158 ymf = rq / dd;
159 break;
160 }
161 default: {
162 throw new AssertionError(mode);
163 }
164 }
165 }
166
167 /**
168 * {@inheritDoc}
169 */
170 public ParameterDescriptorGroup getParameterDescriptors() {
171 return Provider.PARAMETERS;
172 }
173
174 /**
175 * {@inheritDoc}
176 */
177 @Override
178 public ParameterValueGroup getParameterValues() {
179 final ParameterValueGroup values = super.getParameterValues();
180 final Collection<GeneralParameterDescriptor> expected = getParameterDescriptors().descriptors();
181 set(expected, Provider.LATITUDE_OF_CENTRE, values, latitudeOfOrigin);
182 set(expected, Provider.LONGITUDE_OF_CENTRE, values, centralMeridian);
183 return values;
184 }
185
186 /**
187 * Transforms the specified (<var>&lambda;</var>,<var>&phi;</var>) coordinates
188 * (units in radians) and stores the result in {@code ptDst} (linear distance
189 * on a unit sphere).
190 */
191 protected Point2D transformNormalized(final double lambda, final double phi, Point2D ptDst)
192 throws ProjectionException
193 {
194 final double coslam = cos(lambda);
195 final double sinlam = sin(lambda);
196 final double sinphi = sin(phi);
197 double q = qsfn(sinphi);
198 final double sinb, cosb, b, c, x, y;
199 switch (mode) {
200 case OBLIQUE: {
201 sinb = q / qp;
202 cosb = sqrt(1.0 - sinb * sinb);
203 c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam;
204 b = sqrt(2.0 / c);
205 y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam);
206 x = xmf * b * cosb * sinlam;
207 break;
208 }
209 case EQUATORIAL: {
210 sinb = q / qp;
211 cosb = sqrt(1.0 - sinb * sinb);
212 c = 1.0 + cosb * coslam;
213 b = sqrt(2.0 / c);
214 y = ymf * b * sinb;
215 x = xmf * b * cosb * sinlam;
216 break;
217 }
218 case NORTH_POLE: {
219 c = (PI / 2) + phi;
220 q = qp - q;
221 if (q >= 0.0) {
222 b = sqrt(q);
223 x = b * sinlam;
224 y = coslam * -b;
225 } else {
226 x = y = 0.;
227 }
228 break;
229 }
230 case SOUTH_POLE: {
231 c = phi - (PI / 2);
232 q = qp + q;
233 if (q >= 0.0) {
234 b = sqrt(q);
235 x = b * sinlam;
236 y = coslam * +b;
237 } else {
238 x = y = 0.;
239 }
240 break;
241 }
242 default: {
243 throw new AssertionError(mode);
244 }
245 }
246 if (abs(c) < EPSILON_LATITUDE) {
247 throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
248 }
249 if (ptDst != null) {
250 ptDst.setLocation(x,y);
251 return ptDst;
252 }
253 return new Point2D.Double(x,y);
254 }
255
256 /**
257 * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
258 * and stores the result in {@code ptDst}.
259 */
260 @Override
261 @SuppressWarnings("fallthrough")
262 protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
263 throws ProjectionException
264 {
265 final double lambda, phi;
266 switch (mode) {
267 case EQUATORIAL: // Fall through
268 case OBLIQUE: {
269 x /= dd;
270 y *= dd;
271 final double rho = hypot(x, y);
272 if (rho < FINE_EPSILON) {
273 lambda = 0.0;
274 phi = latitudeOfOrigin;
275 } else {
276 double sCe, cCe, /*q,*/ ab;
277 sCe = 2.0 * asin(0.5 * rho / rq);
278 cCe = cos(sCe);
279 sCe = sin(sCe);
280 x *= sCe;
281 if (mode == OBLIQUE) {
282 ab = cCe * sinb1 + y * sCe * cosb1 / rho;
283 //q = qp * ab;
284 y = rho * cosb1 * cCe - y * sinb1 * sCe;
285 } else {
286 ab = y * sCe / rho;
287 //q = qp * ab;
288 y = rho * cCe;
289 }
290 lambda = atan2(x, y);
291 phi = authlat(asin(ab));
292 }
293 break;
294 }
295 case NORTH_POLE: {
296 y = -y;
297 // Fall through
298 }
299 case SOUTH_POLE: {
300 final double q = x*x + y*y;
301 if (q == 0) {
302 lambda = 0.;
303 phi = latitudeOfOrigin;
304 } else {
305 double ab = 1.0 - q / qp;
306 if (mode == SOUTH_POLE) {
307 ab = -ab;
308 }
309 lambda = atan2(x, y);
310 phi = authlat(asin(ab));
311 }
312 break;
313 }
314 default: {
315 throw new AssertionError(mode);
316 }
317 }
318 if (ptDst != null) {
319 ptDst.setLocation(lambda, phi);
320 return ptDst;
321 }
322 return new Point2D.Double(lambda, phi);
323 }
324
325
326 /**
327 * Provides the transform equations for the spherical case.
328 *
329 * @version $Id: LambertAzimuthalEqualArea.java 37299 2011-05-25 05:21:24Z mbedward $
330 * @author Martin Desruisseaux
331 */
332 private static final class Spherical extends LambertAzimuthalEqualArea {
333 /**
334 * For cross-version compatibility.
335 */
336 private static final long serialVersionUID = 2091431369806844342L;
337
338 /**
339 * Constructs a new map projection from the suplied parameters.
340 *
341 * @param parameters The parameter values in standard units.
342 * @throws ParameterNotFoundException if a mandatory parameter is missing.
343 */
344 protected Spherical(final ParameterValueGroup parameters)
345 throws ParameterNotFoundException
346 {
347 super(parameters);
348 ensureSpherical();
349 }
350
351 /**
352 * Transforms the specified (<var>&lambda;</var>,<var>&phi;</var>) coordinates
353 * (units in radians) and stores the result in {@code ptDst} (linear distance
354 * on a unit sphere).
355 */
356 @Override
357 protected Point2D transformNormalized(final double lambda, final double phi, Point2D ptDst)
358 throws ProjectionException
359 {
360 // Compute using ellipsoidal formulas, for comparaison later.
361 assert (ptDst = super.transformNormalized(lambda, phi, ptDst)) != null;
362
363 final double sinphi = sin(phi);
364 final double cosphi = cos(phi);
365 final double coslam = cos(lambda);
366 double x,y;
367 switch (mode) {
368 case EQUATORIAL: {
369 y = 1.0 + cosphi * coslam;
370 if (y <= FINE_EPSILON) {
371 throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
372 }
373 y = sqrt(2.0 / y);
374 x = y * cosphi * sin(lambda);
375 y *= sinphi;
376 break;
377 }
378 case OBLIQUE: {
379 y = 1.0 + sinb1 * sinphi + cosb1 * cosphi * coslam;
380 if (y <= FINE_EPSILON) {
381 throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
382 }
383 y = sqrt(2.0 / y);
384 x = y * cosphi * sin(lambda);
385 y *= cosb1 * sinphi - sinb1 * cosphi * coslam;
386 break;
387 }
388 case NORTH_POLE: {
389 if (abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
390 throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
391 }
392 y = (PI/4) - phi * 0.5;
393 y = 2.0 * sin(y);
394 x = y * sin(lambda);
395 y *= -coslam;
396 break;
397 }
398 case SOUTH_POLE: {
399 if (abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
400 throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
401 }
402 y = (PI/4) - phi * 0.5;
403 y = 2.0 * cos(y);
404 x = y * sin(lambda);
405 y *= +coslam;
406 break;
407 }
408 default: {
409 throw new AssertionError(mode);
410 }
411 }
412 assert checkTransform(x, y, ptDst);
413 if (ptDst != null) {
414 ptDst.setLocation(x,y);
415 return ptDst;
416 }
417 return new Point2D.Double(x,y);
418 }
419
420 /**
421 * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
422 * and stores the result in {@code ptDst} using equations for a sphere.
423 */
424 @Override
425 protected Point2D inverseTransformNormalized(double x, double y, Point2D ptDst)
426 throws ProjectionException
427 {
428 // Compute using ellipsoidal formulas, for comparaison later.
429 assert (ptDst = super.inverseTransformNormalized(x, y, ptDst)) != null;
430
431 double lambda, phi;
432 final double rh = hypot(x, y);
433 phi = rh * 0.5;
434 if (phi > 1.0) {
435 throw new ProjectionException(ErrorKeys.TOLERANCE_ERROR);
436 }
437 phi = 2.0 * asin(phi);
438 switch (mode) {
439 case EQUATORIAL: {
440 final double sinz = sin(phi);
441 final double cosz = cos(phi);
442 phi = abs(rh) <= FINE_EPSILON ? 0.0 : asin(y * sinz / rh);
443 x *= sinz;
444 y = cosz * rh;
445 lambda = (y == 0) ? 0.0 : atan2(x, y);
446 break;
447 }
448 case OBLIQUE: {
449 final double sinz = sin(phi);
450 final double cosz = cos(phi);
451 phi = abs(rh) <= FINE_EPSILON ? latitudeOfOrigin :
452 asin(cosz * sinb1 + y * sinz * cosb1 / rh);
453 x *= sinz * cosb1;
454 y = (cosz - sin(phi) * sinb1) * rh;
455 lambda = (y == 0) ? 0.0 : atan2(x, y);
456 break;
457 }
458 case NORTH_POLE: {
459 phi = (PI / 2) - phi;
460 lambda = atan2(x, -y);
461 break;
462 }
463 case SOUTH_POLE: {
464 phi -= (PI / 2);
465 lambda = atan2(x, y);
466 break;
467 }
468 default: {
469 throw new AssertionError(mode);
470 }
471 }
472 assert checkInverseTransform(lambda, phi, ptDst);
473 if (ptDst != null) {
474 ptDst.setLocation(lambda, phi);
475 return ptDst;
476 }
477 return new Point2D.Double(lambda, phi);
478 }
479 }
480
481 /**
482 * Calculates <var>q</var>, Snyder equation (3-12)
483 *
484 * @param sinphi sin of the latitude <var>q</var> is calculated for.
485 * @return <var>q</var> from Snyder equation (3-12).
486 */
487 private double qsfn(final double sinphi) {
488 if (excentricity >= EPSILON) {
489 final double con = excentricity * sinphi;
490 return ((1.0 - excentricitySquared) * (sinphi / (1.0 - con*con) -
491 (0.5 / excentricity) * log((1.0 - con) / (1.0 + con))));
492 } else {
493 return sinphi + sinphi;
494 }
495 }
496
497 /**
498 * Determines latitude from authalic latitude.
499 */
500 private double authlat(final double beta) {
501 final double t = beta + beta;
502 return beta + APA0 * sin(t) + APA1 * sin(t+t) + APA2 * sin(t+t+t);
503 }
504
505
506
507
508 //////////////////////////////////////////////////////////////////////////////////////////
509 //////////////////////////////////////////////////////////////////////////////////////////
510 //////// ////////
511 //////// PROVIDERS ////////
512 //////// ////////
513 //////////////////////////////////////////////////////////////////////////////////////////
514 //////////////////////////////////////////////////////////////////////////////////////////
515
516 /**
517 * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
518 * provider} for an {@linkplain LambertAzimuthalEqualArea Lambert Equal Area} projection
519 * (EPSG code 9820).
520 *
521 * @since 2.4
522 * @version $Id: LambertAzimuthalEqualArea.java 37299 2011-05-25 05:21:24Z mbedward $
523 * @author Beate Stollberg
524 *
525 * @see org.geotools.referencing.operation.DefaultMathTransformFactory
526 */
527 public static class Provider extends AbstractProvider {
528 /**
529 * For cross-version compatibility.
530 */
531 private static final long serialVersionUID = 3877793025552244132L;
532
533 /**
534 * The operation parameter descriptor for the {@link #latitudeOfOrigin}
535 * parameter value. Valid values range is from -90 to 90°. Default value is 0.
536 */
537 public static final ParameterDescriptor LATITUDE_OF_CENTRE = createDescriptor(
538 new NamedIdentifier[] {
539 new NamedIdentifier(Citations.OGC, "latitude_of_center"),
540 new NamedIdentifier(Citations.EPSG, "Latitude of natural origin"),
541 new NamedIdentifier(Citations.EPSG, "Spherical latitude of origin"),
542 new NamedIdentifier(Citations.ESRI, "Latitude_Of_Origin"),
543 new NamedIdentifier(Citations.GEOTIFF, "ProjCenterLat")
544 },
545 0, -90, 90, NonSI.DEGREE_ANGLE);
546
547 /**
548 * The operation parameter descriptor for the {@link #centralMeridian}
549 * parameter value. Valid values range is from -180 to 180°. Default value is 0.
550 */
551 public static final ParameterDescriptor LONGITUDE_OF_CENTRE = createDescriptor(
552 new NamedIdentifier[] {
553 new NamedIdentifier(Citations.OGC, "longitude_of_center"),
554 new NamedIdentifier(Citations.EPSG, "Longitude of natural origin"),
555 new NamedIdentifier(Citations.EPSG, "Spherical longitude of origin"),
556 new NamedIdentifier(Citations.ESRI, "Central_Meridian"),
557 new NamedIdentifier(Citations.GEOTIFF, "ProjCenterLong")
558 },
559 0, -180, 180, NonSI.DEGREE_ANGLE);
560
561 /**
562 * The parameters group.
563 */
564 static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
565 new NamedIdentifier(Citations.OGC, "Lambert_Azimuthal_Equal_Area"),
566 new NamedIdentifier(Citations.EPSG, "Lambert Azimuthal Equal Area"),
567 new NamedIdentifier(Citations.EPSG, "Lambert Azimuthal Equal Area (Spherical)"),
568 new NamedIdentifier(Citations.GEOTIFF, "CT_LambertAzimEqualArea"),
569 new NamedIdentifier(Citations.EPSG, "9820"),
570 }, new ParameterDescriptor[] {
571 SEMI_MAJOR, SEMI_MINOR,
572 LATITUDE_OF_CENTRE, LONGITUDE_OF_CENTRE,
573 FALSE_EASTING, FALSE_NORTHING
574 });
575
576 /**
577 * Constructs a new provider.
578 */
579 public Provider() {
580 super(PARAMETERS);
581 }
582
583 /**
584 * Creates a transform from the specified group of parameter values.
585 *
586 * @param parameters The group of parameter values.
587 * @return The created math transform.
588 * @throws ParameterNotFoundException if a required parameter was not found.
589 */
590 public MathTransform createMathTransform(final ParameterValueGroup parameters)
591 throws ParameterNotFoundException
592 {
593 return isSpherical(parameters) ? new Spherical(parameters) :
594 new LambertAzimuthalEqualArea(parameters);
595 }
596 }
597}
Note: See TracBrowser for help on using the repository browser.