001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2003-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:
047: /**
048: * Albers Equal Area Projection (EPSG code 9822). This is a conic projection
049: * with parallels being unequally spaced arcs of concentric circles, more
050: * closely spaced at north and south edges of the map. Merideans
051: * are equally spaced radii of the same circles and intersect parallels at right
052: * angles. As the name implies, this projection minimizes distortion in areas.
053: * <p>
054: *
055: * The "standard_parallel_2" parameter is optional and will be given the
056: * same value as "standard_parallel_1" if not set (creating a 1 standard parallel
057: * projection).
058: * <p>
059: *
060: * NOTE: formulae used below are from a port, to java, of the
061: * 'proj4' package of the USGS survey. USGS work is acknowledged here.
062: * <p>
063: *
064: * <strong>References:</strong><ul>
065: * <li> Proj-4.4.7 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
066: * Relevent files are: PJ_aea.c, pj_fwd.c and pj_inv.c </li>
067: * <li> John P. Snyder (Map Projections - A Working Manual,
068: * U.S. Geological Survey Professional Paper 1395, 1987)</li>
069: * <li> "Coordinate Conversions and Transformations including Formulas",
070: * EPSG Guidence Note Number 7, Version 19.</li>
071: * </ul>
072: *
073: * @see <A HREF="http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html">Albers Equal-Area Conic Projection on MathWorld</A>
074: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/albers_equal_area_conic.html">"Albers_Conic_Equal_Area" on RemoteSensing.org</A>
075: * @see <A HREF="http://srmwww.gov.bc.ca/gis/bceprojection.html">British Columbia Albers Standard Projection</A>
076: *
077: * @since 2.1
078: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/AlbersEqualArea.java $
079: * @version $Id: AlbersEqualArea.java 24384 2007-02-14 00:23:05Z desruisseaux $
080: * @author Rueben Schulz
081: */
082: public class AlbersEqualArea extends MapProjection {
083: /**
084: * Maximum number of iterations for iterative computations.
085: */
086: private static final int MAXIMUM_ITERATIONS = 15;
087:
088: /**
089: * Difference allowed in iterative computations.
090: */
091: private static final double ITERATION_TOLERANCE = 1E-10;
092:
093: /**
094: * Maximum difference allowed when comparing real numbers.
095: */
096: private static final double EPSILON = 1E-6;
097:
098: /**
099: * Constants used by the spherical and elliptical Albers projection.
100: */
101: private final double n, c, rho0;
102:
103: /**
104: * An error condition indicating iteration will not converge for the
105: * inverse ellipse. See Snyder (14-20)
106: */
107: private final double ec;
108:
109: /**
110: * Standards parallel 1 in radians, for {@link #getParameterValues} implementation.
111: */
112: private final double phi1;
113:
114: /**
115: * Standards parallel 2 in radians, for {@link #getParameterValues} implementation.
116: */
117: private double phi2;
118:
119: /**
120: * Constructs a new map projection from the supplied parameters.
121: *
122: * @param parameters The parameter values in standard units.
123: * @throws ParameterNotFoundException if a mandatory parameter is missing.
124: */
125: protected AlbersEqualArea(final ParameterValueGroup parameters)
126: throws ParameterNotFoundException {
127: // Fetch parameters
128: super (parameters);
129: final Collection expected = getParameterDescriptors()
130: .descriptors();
131: phi1 = doubleValue(expected, Provider.STANDARD_PARALLEL_1,
132: parameters);
133: ensureLatitudeInRange(Provider.STANDARD_PARALLEL_1, phi1, true);
134: phi2 = doubleValue(expected, Provider.STANDARD_PARALLEL_2,
135: parameters);
136: if (Double.isNaN(phi2)) {
137: phi2 = phi1;
138: }
139: ensureLatitudeInRange(Provider.STANDARD_PARALLEL_2, phi2, true);
140:
141: // Compute Constants
142: if (Math.abs(phi1 + phi2) < EPSILON)
143: throw new IllegalArgumentException(Errors.format(
144: ErrorKeys.ANTIPODE_LATITUDES_$2, new Latitude(Math
145: .toDegrees(phi1)), new Latitude(Math
146: .toDegrees(phi2))));
147:
148: double sinphi = Math.sin(phi1);
149: double cosphi = Math.cos(phi1);
150: double n = sinphi;
151: boolean secant = (Math.abs(phi1 - phi2) >= EPSILON);
152: if (isSpherical) {
153: if (secant) {
154: n = 0.5 * (n + Math.sin(phi2));
155: }
156: c = cosphi * cosphi + n * 2 * sinphi;
157: rho0 = Math.sqrt(c - n * 2 * Math.sin(latitudeOfOrigin))
158: / n;
159: ec = Double.NaN;
160: } else {
161: double m1 = msfn(sinphi, cosphi);
162: double q1 = qsfn(sinphi);
163: if (secant) { /* secant cone */
164: sinphi = Math.sin(phi2);
165: cosphi = Math.cos(phi2);
166: double m2 = msfn(sinphi, cosphi);
167: double q2 = qsfn(sinphi);
168: n = (m1 * m1 - m2 * m2) / (q2 - q1);
169: }
170: c = m1 * m1 + n * q1;
171: rho0 = Math.sqrt(c - n * qsfn(Math.sin(latitudeOfOrigin)))
172: / n;
173: ec = 1.0
174: - .5
175: * (1.0 - excentricitySquared)
176: * Math.log((1.0 - excentricity)
177: / (1.0 + excentricity)) / excentricity;
178: }
179: this .n = n;
180: }
181:
182: /**
183: * {@inheritDoc}
184: */
185: public ParameterDescriptorGroup getParameterDescriptors() {
186: return Provider.PARAMETERS;
187: }
188:
189: /**
190: * {@inheritDoc}
191: */
192: public ParameterValueGroup getParameterValues() {
193: final ParameterValueGroup values = super .getParameterValues();
194: final Collection expected = getParameterDescriptors()
195: .descriptors();
196: set(expected, Provider.STANDARD_PARALLEL_1, values, phi1);
197: set(expected, Provider.STANDARD_PARALLEL_2, values, phi2);
198: return values;
199: }
200:
201: /**
202: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
203: * (units in radians) and stores the result in {@code ptDst} (linear distance
204: * on a unit sphere).
205: */
206: protected Point2D transformNormalized(double x, double y,
207: Point2D ptDst) throws ProjectionException {
208: x *= n;
209: double rho;
210: if (isSpherical) {
211: rho = c - n * 2 * Math.sin(y);
212: } else {
213: rho = c - n * qsfn(Math.sin(y));
214: }
215:
216: if (rho < 0.0) {
217: if (rho > -EPSILON) {
218: rho = 0.0;
219: } else {
220: throw new ProjectionException(Errors
221: .format(ErrorKeys.TOLERANCE_ERROR));
222: }
223: }
224: rho = Math.sqrt(rho) / n;
225: y = rho0 - rho * Math.cos(x);
226: x = rho * Math.sin(x);
227:
228: if (ptDst != null) {
229: ptDst.setLocation(x, y);
230: return ptDst;
231: }
232: return new Point2D.Double(x, y);
233: }
234:
235: /**
236: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
237: * and stores the result in {@code ptDst}.
238: */
239: protected Point2D inverseTransformNormalized(double x, double y,
240: Point2D ptDst) throws ProjectionException {
241: y = rho0 - y;
242: double rho = Math.sqrt(x * x + y * y);
243: if (rho > EPSILON) {
244: if (n < 0.0) {
245: rho = -rho;
246: x = -x;
247: y = -y;
248: }
249: x = Math.atan2(x, y) / n;
250: y = rho * n;
251: if (isSpherical) {
252: y = (c - y * y) / (n * 2);
253: if (Math.abs(y) <= 1.0) {
254: y = Math.asin(y);
255: } else {
256: y = (y < 0.0) ? -Math.PI / 2.0 : Math.PI / 2.0;
257: }
258: } else {
259: y = (c - y * y) / n;
260: if (Math.abs(ec - Math.abs(y)) > EPSILON) {
261: y = phi1(y);
262: } else {
263: y = (y < 0.0) ? -Math.PI / 2.0 : Math.PI / 2.0;
264: }
265: }
266: } else {
267: x = 0.0;
268: y = n > 0.0 ? Math.PI / 2.0 : -Math.PI / 2.0;
269: }
270:
271: if (ptDst != null) {
272: ptDst.setLocation(x, y);
273: return ptDst;
274: }
275: return new Point2D.Double(x, y);
276: }
277:
278: /**
279: * Iteratively solves equation (3-16) from Snyder.
280: *
281: * @param qs arcsin(q/2), used in the first step of iteration
282: * @return the latitude
283: */
284: private double phi1(final double qs) throws ProjectionException {
285: final double tone_es = 1 - excentricitySquared;
286: double phi = Math.asin(0.5 * qs);
287: if (excentricity < EPSILON) {
288: return phi;
289: }
290: for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
291: final double sinpi = Math.sin(phi);
292: final double cospi = Math.cos(phi);
293: final double con = excentricity * sinpi;
294: final double com = 1.0 - con * con;
295: final double dphi = 0.5
296: * com
297: * com
298: / cospi
299: * (qs / tone_es - sinpi / com + 0.5 / excentricity
300: * Math.log((1. - con) / (1. + con)));
301: phi += dphi;
302: if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
303: return phi;
304: }
305: }
306: throw new ProjectionException(Errors
307: .format(ErrorKeys.NO_CONVERGENCE));
308: }
309:
310: /**
311: * Calculates q, Snyder equation (3-12)
312: *
313: * @param sinphi sin of the latitude q is calculated for
314: * @return q from Snyder equation (3-12)
315: */
316: private double qsfn(final double sinphi) {
317: final double one_es = 1 - excentricitySquared;
318: if (excentricity >= EPSILON) {
319: final double con = excentricity * sinphi;
320: return (one_es * (sinphi / (1. - con * con) - (0.5 / excentricity)
321: * Math.log((1. - con) / (1. + con))));
322: } else {
323: return sinphi + sinphi;
324: }
325: }
326:
327: /**
328: * Returns a hash value for this projection.
329: */
330: public int hashCode() {
331: final long code = Double.doubleToLongBits(c);
332: return ((int) code ^ (int) (code >>> 32)) + 37
333: * super .hashCode();
334: }
335:
336: /**
337: * Compares the specified object with this map projection for equality.
338: */
339: public boolean equals(final Object object) {
340: if (object == this ) {
341: // Slight optimization
342: return true;
343: }
344: if (super .equals(object)) {
345: final AlbersEqualArea that = (AlbersEqualArea) object;
346: return equals(this .n, that.n) && equals(this .c, that.c)
347: && equals(this .rho0, that.rho0)
348: && equals(this .phi1, that.phi1)
349: && equals(this .phi2, that.phi2);
350: }
351: return false;
352: }
353:
354: //////////////////////////////////////////////////////////////////////////////////////////
355: //////////////////////////////////////////////////////////////////////////////////////////
356: //////// ////////
357: //////// PROVIDERS ////////
358: //////// ////////
359: //////////////////////////////////////////////////////////////////////////////////////////
360: //////////////////////////////////////////////////////////////////////////////////////////
361:
362: /**
363: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
364: * provider} for an {@linkplain AlbersEqualArea Albers Equal Area} projection (EPSG code 9822).
365: *
366: * @version $Id: AlbersEqualArea.java 24384 2007-02-14 00:23:05Z desruisseaux $
367: * @author Rueben Schulz
368: *
369: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
370: */
371: public static class Provider extends AbstractProvider {
372: /**
373: * The parameters group.
374: */
375: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
376: new NamedIdentifier[] {
377: new NamedIdentifier(Citations.OGC,
378: "Albers_Conic_Equal_Area"),
379: new NamedIdentifier(Citations.EPSG,
380: "Albers Equal Area"),
381: new NamedIdentifier(Citations.EPSG, "9822"),
382: new NamedIdentifier(Citations.GEOTIFF,
383: "CT_AlbersEqualArea"),
384: new NamedIdentifier(Citations.ESRI, "Albers"),
385: new NamedIdentifier(Citations.ESRI,
386: "Albers Equal Area Conic"),
387: new NamedIdentifier(
388: Citations.GEOTOOLS,
389: Vocabulary
390: .formatInternational(VocabularyKeys.ALBERS_EQUAL_AREA_PROJECTION)) },
391: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
392: CENTRAL_MERIDIAN, LATITUDE_OF_ORIGIN,
393: STANDARD_PARALLEL_1, STANDARD_PARALLEL_2,
394: FALSE_EASTING, FALSE_NORTHING });
395:
396: /**
397: * Constructs a new provider.
398: */
399: public Provider() {
400: super (PARAMETERS);
401: }
402:
403: /**
404: * Returns the operation type for this map projection.
405: */
406: public Class getOperationType() {
407: return ConicProjection.class;
408: }
409:
410: /**
411: * Creates a transform from the specified group of parameter values.
412: *
413: * @param parameters The group of parameter values.
414: * @return The created math transform.
415: * @throws ParameterNotFoundException if a required parameter was not found.
416: */
417: protected MathTransform createMathTransform(
418: final ParameterValueGroup parameters)
419: throws ParameterNotFoundException {
420: return new AlbersEqualArea(parameters);
421: }
422: }
423: }
|