001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2006, Geotools Project Managment Committee (PMC)
006: * (C) 2000, Frank Warmerdam
007: * (C) 1995, Gerald Evenden
008: *
009: * This library is free software; you can redistribute it and/or
010: * modify it under the terms of the GNU Lesser General Public
011: * License as published by the Free Software Foundation;
012: * version 2.1 of the License.
013: *
014: * This library is distributed in the hope that it will be useful,
015: * but WITHOUT ANY WARRANTY; without even the implied warranty of
016: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
017: * Lesser General Public License for more details.
018: *
019: * This package contains formulas from the PROJ package of USGS.
020: * USGS's work is fully acknowledged here. This derived work has
021: * been relicensed under LGPL with Frank Warmerdam's permission.
022: */
023: package org.geotools.referencing.operation.projection;
024:
025: // J2SE dependencies and extensions
026: import java.awt.geom.Point2D;
027: import java.util.Collection;
028: import javax.units.NonSI;
029:
030: // OpenGIS dependencies
031: import org.opengis.parameter.ParameterDescriptor;
032: import org.opengis.parameter.ParameterDescriptorGroup;
033: import org.opengis.parameter.ParameterNotFoundException;
034: import org.opengis.parameter.ParameterValueGroup;
035: import org.opengis.referencing.operation.ConicProjection;
036: import org.opengis.referencing.operation.MathTransform;
037:
038: // Geotools dependencies
039: import org.geotools.measure.Latitude;
040: import org.geotools.metadata.iso.citation.Citations;
041: import org.geotools.referencing.NamedIdentifier;
042: import org.geotools.resources.i18n.VocabularyKeys;
043: import org.geotools.resources.i18n.Vocabulary;
044: import org.geotools.resources.i18n.ErrorKeys;
045: import org.geotools.resources.i18n.Errors;
046: import org.geotools.resources.XMath;
047:
048: /**
049: * Lambert Azimuthal Equal Area (EPSG code 9820).
050: *
051: * <strong>References:</strong><ul>
052: * <li> A. Annoni, C. Luzet, E.Gubler and J. Ihde - Map Projections for Europe</li>
053: * <li> John P. Snyder (Map Projections - A Working Manual,
054: * U.S. Geological Survey Professional Paper 1395)</li>
055: * </ul>
056: *
057: * @see <A HREF="http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html">Lambert Azimuthal Equal-Area Projection on MathWorld</A>
058: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_azimuthal_equal_area.html">"Lambert_Azimuthal_Equal_Area" on RemoteSensing.org</A>
059: *
060: * @since 2.4
061: * @version $Id: LambertAzimuthalEqualArea.java 24333 2007-02-10 00:16:22Z desruisseaux $
062: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/LambertAzimuthalEqualArea.java $
063: * @author Gerald Evenden
064: * @author Beate Stollberg
065: * @author Martin Desruisseaux
066: */
067: public class LambertAzimuthalEqualArea extends MapProjection {
068: /** Maximum difference allowed when comparing real numbers. */
069: private static final double EPSILON = 1E-7;
070:
071: /** Epsilon for the comparaison of small quantities. */
072: private static final double FINE_EPSILON = 1E-10;
073:
074: /** Epsilon for the comparaison of latitudes. */
075: private static final double EPSILON_LATITUDE = 1E-10;
076:
077: /** Constants for authalic latitude. */
078: private static final double P00 = 0.33333333333333333333,
079: P01 = 0.17222222222222222222, P02 = 0.10257936507936507936,
080: P10 = 0.06388888888888888888, P11 = 0.06640211640211640211,
081: P20 = 0.01641501294219154443;
082:
083: /** The projection mode. */
084: static final int OBLIQUE = 0, EQUATORIAL = 1, NORTH_POLE = 2,
085: SOUTH_POLE = 3;
086:
087: /** The projection mode for this particular instance. */
088: final int mode;
089:
090: /** Constant parameters. */
091: final double sinb1, cosb1, xmf, ymf, mmf, qp, dd, rq;
092:
093: /** Coefficients for authalic latitude. */
094: private final double APA0, APA1, APA2;
095:
096: /**
097: * Constructs a new map projection from the supplied parameters.
098: *
099: * @param parameters The parameter values in standard units.
100: * @throws ParameterNotFoundException if a mandatory parameter is missing.
101: */
102: protected LambertAzimuthalEqualArea(
103: final ParameterValueGroup parameters)
104: throws ParameterNotFoundException {
105: // Fetch parameters
106: super (parameters);
107: final Collection expected = getParameterDescriptors()
108: .descriptors();
109: latitudeOfOrigin = doubleValue(expected,
110: Provider.LATITUDE_OF_CENTRE, parameters);
111: centralMeridian = doubleValue(expected,
112: Provider.LONGITUDE_OF_CENTRE, parameters);
113: ensureLatitudeInRange(Provider.LATITUDE_OF_CENTRE,
114: latitudeOfOrigin, true);
115: ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTRE,
116: centralMeridian, true);
117: /*
118: * Detects the mode (oblique, etc.).
119: */
120: final double t = Math.abs(latitudeOfOrigin);
121: if (Math.abs(t - Math.PI / 2) < EPSILON_LATITUDE) {
122: mode = latitudeOfOrigin < 0.0 ? SOUTH_POLE : NORTH_POLE;
123: } else if (Math.abs(t) < EPSILON_LATITUDE) {
124: mode = EQUATORIAL;
125: } else {
126: mode = OBLIQUE;
127: }
128: /*
129: * Computes the constants for authalic latitude.
130: */
131: final double es2 = excentricitySquared * excentricitySquared;
132: final double es3 = excentricitySquared * es2;
133: APA0 = P02 * es3 + P01 * es2 + P00 * excentricitySquared;
134: APA1 = P11 * es3 + P10 * es2;
135: APA2 = P20 * es3;
136:
137: final double sinphi;
138: qp = qsfn(1);
139: rq = Math.sqrt(0.5 * qp);
140: mmf = 0.5 / (1 - excentricitySquared);
141: sinphi = Math.sin(latitudeOfOrigin);
142: if (isSpherical) {
143: sinb1 = Math.sin(latitudeOfOrigin);
144: cosb1 = Math.cos(latitudeOfOrigin);
145: } else {
146: sinb1 = qsfn(sinphi) / qp;
147: cosb1 = Math.sqrt(1.0 - sinb1 * sinb1);
148: }
149: switch (mode) {
150: case NORTH_POLE: // Fall through
151: case SOUTH_POLE: {
152: dd = 1.0;
153: xmf = ymf = rq;
154: break;
155: }
156: case EQUATORIAL: {
157: dd = 1.0 / rq;
158: xmf = 1.0;
159: ymf = 0.5 * qp;
160: break;
161: }
162: case OBLIQUE: {
163: dd = Math.cos(latitudeOfOrigin)
164: / (Math.sqrt(1.0 - excentricitySquared * sinphi
165: * sinphi)
166: * rq * cosb1);
167: xmf = rq * dd;
168: ymf = rq / dd;
169: break;
170: }
171: default: {
172: throw new AssertionError(mode);
173: }
174: }
175: }
176:
177: /**
178: * {@inheritDoc}
179: */
180: public ParameterDescriptorGroup getParameterDescriptors() {
181: return Provider.PARAMETERS;
182: }
183:
184: /**
185: * {@inheritDoc}
186: */
187: public ParameterValueGroup getParameterValues() {
188: final ParameterValueGroup values = super .getParameterValues();
189: final Collection expected = getParameterDescriptors()
190: .descriptors();
191: set(expected, Provider.LATITUDE_OF_CENTRE, values,
192: latitudeOfOrigin);
193: set(expected, Provider.LONGITUDE_OF_CENTRE, values,
194: centralMeridian);
195: return values;
196: }
197:
198: /**
199: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
200: * (units in radians) and stores the result in {@code ptDst} (linear distance
201: * on a unit sphere).
202: */
203: protected Point2D transformNormalized(final double lambda,
204: final double phi, Point2D ptDst) throws ProjectionException {
205: final double coslam = Math.cos(lambda);
206: final double sinlam = Math.sin(lambda);
207: final double sinphi = Math.sin(phi);
208: double q = qsfn(sinphi);
209: final double sinb, cosb, b, c, x, y;
210: switch (mode) {
211: case OBLIQUE: {
212: sinb = q / qp;
213: cosb = Math.sqrt(1.0 - sinb * sinb);
214: c = 1.0 + sinb1 * sinb + cosb1 * cosb * coslam;
215: b = Math.sqrt(2.0 / c);
216: y = ymf * b * (cosb1 * sinb - sinb1 * cosb * coslam);
217: x = xmf * b * cosb * sinlam;
218: break;
219: }
220: case EQUATORIAL: {
221: sinb = q / qp;
222: cosb = Math.sqrt(1.0 - sinb * sinb);
223: c = 1.0 + cosb * coslam;
224: b = Math.sqrt(2.0 / c);
225: y = ymf * b * sinb;
226: x = xmf * b * cosb * sinlam;
227: break;
228: }
229: case NORTH_POLE: {
230: c = (Math.PI / 2) + phi;
231: q = qp - q;
232: if (q >= 0.0) {
233: b = Math.sqrt(q);
234: x = b * sinlam;
235: y = coslam * -b;
236: } else {
237: x = y = 0.;
238: }
239: break;
240: }
241: case SOUTH_POLE: {
242: c = phi - (Math.PI / 2);
243: q = qp + q;
244: if (q >= 0.0) {
245: b = Math.sqrt(q);
246: x = b * sinlam;
247: y = coslam * +b;
248: } else {
249: x = y = 0.;
250: }
251: break;
252: }
253: default: {
254: throw new AssertionError(mode);
255: }
256: }
257: if (Math.abs(c) < EPSILON_LATITUDE) {
258: throw toleranceError();
259: }
260: if (ptDst != null) {
261: ptDst.setLocation(x, y);
262: return ptDst;
263: }
264: return new Point2D.Double(x, y);
265: }
266:
267: /**
268: * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
269: * and stores the result in {@code ptDst}.
270: */
271: protected Point2D inverseTransformNormalized(double x, double y,
272: Point2D ptDst) throws ProjectionException {
273: final double lambda, phi;
274: switch (mode) {
275: case EQUATORIAL: // Fall through
276: case OBLIQUE: {
277: x /= dd;
278: y *= dd;
279: final double rho = XMath.hypot(x, y);
280: if (rho < FINE_EPSILON) {
281: lambda = 0.0;
282: phi = latitudeOfOrigin;
283: } else {
284: double sCe, cCe, q, ab;
285: sCe = 2.0 * Math.asin(0.5 * rho / rq);
286: cCe = Math.cos(sCe);
287: sCe = Math.sin(sCe);
288: x *= sCe;
289: if (mode == OBLIQUE) {
290: ab = cCe * sinb1 + y * sCe * cosb1 / rho;
291: q = qp * ab;
292: y = rho * cosb1 * cCe - y * sinb1 * sCe;
293: } else {
294: ab = y * sCe / rho;
295: q = qp * ab;
296: y = rho * cCe;
297: }
298: lambda = Math.atan2(x, y);
299: phi = authlat(Math.asin(ab));
300: }
301: break;
302: }
303: case NORTH_POLE: {
304: y = -y;
305: // Fall through
306: }
307: case SOUTH_POLE: {
308: final double q = x * x + y * y;
309: if (q == 0) {
310: lambda = 0.;
311: phi = latitudeOfOrigin;
312: } else {
313: double ab = 1.0 - q / qp;
314: if (mode == SOUTH_POLE) {
315: ab = -ab;
316: }
317: lambda = Math.atan2(x, y);
318: phi = authlat(Math.asin(ab));
319: }
320: break;
321: }
322: default: {
323: throw new AssertionError(mode);
324: }
325: }
326: if (ptDst != null) {
327: ptDst.setLocation(lambda, phi);
328: return ptDst;
329: }
330: return new Point2D.Double(lambda, phi);
331: }
332:
333: /**
334: * Provides the transform equations for the spherical case.
335: *
336: * @version $Id: LambertAzimuthalEqualArea.java 24333 2007-02-10 00:16:22Z desruisseaux $
337: * @author Martin Desruisseaux
338: */
339: private static final class Spherical extends
340: LambertAzimuthalEqualArea {
341: /**
342: * Constructs a new map projection from the suplied parameters.
343: *
344: * @param parameters The parameter values in standard units.
345: * @throws ParameterNotFoundException if a mandatory parameter is missing.
346: */
347: protected Spherical(final ParameterValueGroup parameters)
348: throws ParameterNotFoundException {
349: super (parameters);
350: ensureSpherical();
351: }
352:
353: /**
354: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
355: * (units in radians) and stores the result in {@code ptDst} (linear distance
356: * on a unit sphere).
357: */
358: protected Point2D transformNormalized(final double lambda,
359: final double phi, Point2D ptDst)
360: throws ProjectionException {
361: // Compute using ellipsoidal formulas, for comparaison later.
362: assert (ptDst = super .transformNormalized(lambda, phi,
363: ptDst)) != null;
364:
365: final double sinphi = Math.sin(phi);
366: final double cosphi = Math.cos(phi);
367: final double coslam = Math.cos(lambda);
368: double x, y;
369: switch (mode) {
370: case EQUATORIAL: {
371: y = 1.0 + cosphi * coslam;
372: if (y <= FINE_EPSILON) {
373: throw toleranceError();
374: }
375: y = Math.sqrt(2.0 / y);
376: x = y * cosphi * Math.sin(lambda);
377: y *= sinphi;
378: break;
379: }
380: case OBLIQUE: {
381: y = 1.0 + sinb1 * sinphi + cosb1 * cosphi * coslam;
382: if (y <= FINE_EPSILON) {
383: throw toleranceError();
384: }
385: y = Math.sqrt(2.0 / y);
386: x = y * cosphi * Math.sin(lambda);
387: y *= cosb1 * sinphi - sinb1 * cosphi * coslam;
388: break;
389: }
390: case NORTH_POLE: {
391: if (Math.abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
392: throw toleranceError();
393: }
394: y = (Math.PI / 4) - phi * 0.5;
395: y = 2.0 * Math.sin(y);
396: x = y * Math.sin(lambda);
397: y *= -coslam;
398: break;
399: }
400: case SOUTH_POLE: {
401: if (Math.abs(phi + latitudeOfOrigin) < EPSILON_LATITUDE) {
402: throw toleranceError();
403: }
404: y = (Math.PI / 4) - phi * 0.5;
405: y = 2.0 * Math.cos(y);
406: x = y * Math.sin(lambda);
407: y *= +coslam;
408: break;
409: }
410: default: {
411: throw new AssertionError(mode);
412: }
413: }
414: assert checkTransform(x, y, ptDst);
415: if (ptDst != null) {
416: ptDst.setLocation(x, y);
417: return ptDst;
418: }
419: return new Point2D.Double(x, y);
420: }
421:
422: /**
423: * Transforms the specified (<var>x</var>,<var>y</var>) coordinate
424: * and stores the result in {@code ptDst} using equations for a sphere.
425: */
426: protected Point2D inverseTransformNormalized(double x,
427: double y, Point2D ptDst) throws ProjectionException {
428: // Compute using ellipsoidal formulas, for comparaison later.
429: assert (ptDst = super .inverseTransformNormalized(x, y,
430: ptDst)) != null;
431:
432: double lambda, phi;
433: final double rh = XMath.hypot(x, y);
434: phi = rh * 0.5;
435: if (phi > 1.0) {
436: throw toleranceError();
437: }
438: phi = 2.0 * Math.asin(phi);
439: switch (mode) {
440: case EQUATORIAL: {
441: final double sinz = Math.sin(phi);
442: final double cosz = Math.cos(phi);
443: phi = Math.abs(rh) <= FINE_EPSILON ? 0.0 : Math.asin(y
444: * sinz / rh);
445: x *= sinz;
446: y = cosz * rh;
447: lambda = (y == 0) ? 0.0 : Math.atan2(x, y);
448: break;
449: }
450: case OBLIQUE: {
451: final double sinz = Math.sin(phi);
452: final double cosz = Math.cos(phi);
453: phi = Math.abs(rh) <= FINE_EPSILON ? latitudeOfOrigin
454: : Math.asin(cosz * sinb1 + y * sinz * cosb1
455: / rh);
456: x *= sinz * cosb1;
457: y = (cosz - Math.sin(phi) * sinb1) * rh;
458: lambda = (y == 0) ? 0.0 : Math.atan2(x, y);
459: break;
460: }
461: case NORTH_POLE: {
462: phi = (Math.PI / 2) - phi;
463: lambda = Math.atan2(x, -y);
464: break;
465: }
466: case SOUTH_POLE: {
467: phi -= (Math.PI / 2);
468: lambda = Math.atan2(x, y);
469: break;
470: }
471: default: {
472: throw new AssertionError(mode);
473: }
474: }
475: assert checkInverseTransform(lambda, phi, ptDst);
476: if (ptDst != null) {
477: ptDst.setLocation(lambda, phi);
478: return ptDst;
479: }
480: return new Point2D.Double(lambda, phi);
481: }
482: }
483:
484: /**
485: * Calculates <var>q</var>, Snyder equation (3-12)
486: *
487: * @param sinphi sin of the latitude <var>q</var> is calculated for.
488: * @return <var>q</var> from Snyder equation (3-12).
489: */
490: private double qsfn(final double sinphi) {
491: if (excentricity >= EPSILON) {
492: final double con = excentricity * sinphi;
493: return ((1.0 - excentricitySquared) * (sinphi
494: / (1.0 - con * con) - (0.5 / excentricity)
495: * Math.log((1.0 - con) / (1.0 + con))));
496: } else {
497: return sinphi + sinphi;
498: }
499: }
500:
501: /**
502: * Determines latitude from authalic latitude.
503: */
504: private double authlat(final double beta) {
505: final double t = beta + beta;
506: return beta + APA0 * Math.sin(t) + APA1 * Math.sin(t + t)
507: + APA2 * Math.sin(t + t + t);
508: }
509:
510: /**
511: * Returns an exception for a tolerance error (error code -20 in Proj4).
512: */
513: private static ProjectionException toleranceError() {
514: return new ProjectionException(Errors
515: .format(ErrorKeys.TOLERANCE_ERROR));
516: }
517:
518: //////////////////////////////////////////////////////////////////////////////////////////
519: //////////////////////////////////////////////////////////////////////////////////////////
520: //////// ////////
521: //////// PROVIDERS ////////
522: //////// ////////
523: //////////////////////////////////////////////////////////////////////////////////////////
524: //////////////////////////////////////////////////////////////////////////////////////////
525:
526: /**
527: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
528: * provider} for an {@linkplain LambertAzimuthalEqualArea Lambert Equal Area} projection
529: * (EPSG code 9820).
530: *
531: * @since 2.4
532: * @version $Id: LambertAzimuthalEqualArea.java 24333 2007-02-10 00:16:22Z desruisseaux $
533: * @author Beate Stollberg
534: *
535: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
536: */
537: public static class Provider extends AbstractProvider {
538: /**
539: * The operation parameter descriptor for the {@link #latitudeOfOrigin}
540: * parameter value. Valid values range is from -90 to 90°. Default value is 0.
541: */
542: public static final ParameterDescriptor LATITUDE_OF_CENTRE = createDescriptor(
543: new NamedIdentifier[] {
544: new NamedIdentifier(Citations.OGC,
545: "latitude_of_center"),
546: new NamedIdentifier(Citations.EPSG,
547: "Latitude of natural origin"),
548: new NamedIdentifier(Citations.ESRI,
549: "latitude_of_origin"),
550: new NamedIdentifier(Citations.GEOTIFF,
551: "ProjCenterLat") }, 0, -90, 90,
552: NonSI.DEGREE_ANGLE);
553:
554: /**
555: * The operation parameter descriptor for the {@link #centralMeridian}
556: * parameter value. Valid values range is from -180 to 180°. Default value is 0.
557: */
558: public static final ParameterDescriptor LONGITUDE_OF_CENTRE = createDescriptor(
559: new NamedIdentifier[] {
560: new NamedIdentifier(Citations.OGC,
561: "longitude_of_center"),
562: new NamedIdentifier(Citations.EPSG,
563: "Longitude of natural origin"),
564: new NamedIdentifier(Citations.ESRI,
565: "central_meridian"),
566: new NamedIdentifier(Citations.GEOTIFF,
567: "ProjCenterLong") }, 0, -180, 180,
568: NonSI.DEGREE_ANGLE);
569:
570: /**
571: * The parameters group.
572: */
573: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
574: new NamedIdentifier[] {
575: new NamedIdentifier(Citations.OGC,
576: "Lambert_Azimuthal_Equal_Area"),
577: new NamedIdentifier(Citations.EPSG,
578: "Lambert Azimuthal Equal Area"),
579: new NamedIdentifier(Citations.GEOTIFF,
580: "CT_LambertAzimEqualArea"),
581: new NamedIdentifier(Citations.EPSG, "9820"), },
582: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
583: LATITUDE_OF_CENTRE, LONGITUDE_OF_CENTRE,
584: FALSE_EASTING, FALSE_NORTHING });
585:
586: /**
587: * Constructs a new provider.
588: */
589: public Provider() {
590: super (PARAMETERS);
591: }
592:
593: /**
594: * Creates a transform from the specified group of parameter values.
595: *
596: * @param parameters The group of parameter values.
597: * @return The created math transform.
598: * @throws ParameterNotFoundException if a required parameter was not found.
599: */
600: public MathTransform createMathTransform(
601: final ParameterValueGroup parameters)
602: throws ParameterNotFoundException {
603: return isSpherical(parameters) ? new Spherical(parameters)
604: : new LambertAzimuthalEqualArea(parameters);
605: }
606: }
607: }
|