001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2003, 2004 Geotools Project Managment Committee (PMC)
006: * (C) 2001, Institut de Recherche pour le Développement
007: * (C) 2000, Frank Warmerdam
008: * (C) 1999, Fisheries and Oceans Canada
009: *
010: * This library is free software; you can redistribute it and/or
011: * modify it under the terms of the GNU Lesser General Public
012: * License as published by the Free Software Foundation; either
013: * version 2.1 of the License, or (at your option) any later version.
014: *
015: * This library is distributed in the hope that it will be useful,
016: * but WITHOUT ANY WARRANTY; without even the implied warranty of
017: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
018: * Lesser General Public License for more details.
019: *
020: * This package contains formulas from the PROJ package of USGS.
021: * USGS's work is fully acknowledged here. This derived work has
022: * been relicensed under LGPL with Frank Warmerdam's permission.
023: */
024: package org.geotools.referencing.operation.projection;
025:
026: // J2SE dependencies
027: import java.awt.geom.Point2D;
028: import java.util.Collection;
029:
030: // OpenGIS dependencies
031: import org.opengis.parameter.ParameterValueGroup;
032: import org.opengis.parameter.ParameterNotFoundException;
033:
034: // Geotools dependencies
035: import org.geotools.measure.Latitude;
036: import org.geotools.resources.i18n.Errors;
037: import org.geotools.resources.i18n.ErrorKeys;
038:
039: /**
040: * Lambert Conical Conformal Projection. Areas and shapes are deformed
041: * as one moves away from standard parallels. The angles are true in
042: * a limited area. This projection is used for the charts of North America.
043: * <p>
044: *
045: * This implementation provides transforms for three cases of the lambert conic
046: * conformal projection:
047: * <ul>
048: * <li>{@code Lambert_Conformal_Conic_1SP} (EPSG code 9801)</li>
049: * <li>{@code Lambert_Conformal_Conic_2SP} (EPSG code 9802)</li>
050: * <li>{@code Lambert_Conic_Conformal_2SP_Belgium} (EPSG code 9803)</li>
051: * <li>{@code Lambert_Conformal_Conic} - An alias for the ESRI 2SP case
052: * that includes a scale_factor parameter</li>
053: * </ul>
054: *
055: * For the 1SP case the latitude of origin is used as the standard parallel (SP).
056: * To use 1SP with a latitude of origin different from the SP, use the 2SP
057: * and set the SP1 to the single SP. The "standard_parallel_2"
058: * parameter is optional and will be given the same value as "standard_parallel_1"
059: * if not set (creating a 1 standard parallel projection).
060: * <p>
061: *
062: * <strong>References:</strong><ul>
063: * <li>John P. Snyder (Map Projections - A Working Manual,<br>
064: * U.S. Geological Survey Professional Paper 1395, 1987)</li>
065: * <li>"Coordinate Conversions and Transformations including Formulas",<br>
066: * EPSG Guidence Note Number 7, Version 19.</li>
067: * </ul>
068: *
069: * @see <A HREF="http://mathworld.wolfram.com/LambertConformalConicProjection.html">Lambert conformal conic projection on MathWorld</A>
070: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html">lambert_conic_conformal_1sp</A>
071: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html">lambert_conic_conformal_2sp</A>
072: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp_belgium.html">lambert_conic_conformal_2sp_belgium</A>
073: *
074: * @since 2.1
075: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/LambertConformal.java $
076: * @version $Id: LambertConformal.java 25485 2007-05-11 19:12:35Z desruisseaux $
077: * @author André Gosselin
078: * @author Martin Desruisseaux
079: * @author Rueben Schulz
080: */
081: public abstract class LambertConformal extends MapProjection {
082: /**
083: * Maximum difference allowed when comparing real numbers.
084: */
085: private static final double EPSILON = 1E-6;
086:
087: /**
088: * Constant for the belgium 2SP case. This is 29.2985 seconds, given here in radians.
089: */
090: private static final double BELGE_A = 0.00014204313635987700;
091:
092: /**
093: * Standards parallel 1 in radians, for {@link #getParameterValues} implementation.
094: */
095: private final double phi1;
096:
097: /**
098: * Standards parallel 2 in radians, for {@link #getParameterValues} implementation.
099: */
100: final private double phi2;
101:
102: /**
103: * Internal variables for computation.
104: */
105: private final double n, F, rho0;
106:
107: /**
108: * {@code true} for Belgium 2SP.
109: */
110: private final boolean belgium;
111:
112: /**
113: * Constructs a new map projection from the supplied parameters.
114: *
115: * @param parameters The parameter values in standard units.
116: * @throws ParameterNotFoundException if a mandatory parameter is missing.
117: */
118: protected LambertConformal(final ParameterValueGroup parameters)
119: throws ParameterNotFoundException {
120: this (parameters, false);
121: }
122:
123: /**
124: * Constructs a new map projection from the supplied parameters.
125: *
126: * @param parameters The parameter values in standard units.
127: * @param belgium {@code true} for the Belgium 2SP case.
128: * @throws ParameterNotFoundException if a mandatory parameter is missing.
129: */
130: LambertConformal(final ParameterValueGroup parameters,
131: final boolean belgium) throws ParameterNotFoundException {
132: //Fetch parameters
133: super (parameters);
134: final Collection expected = getParameterDescriptors()
135: .descriptors();
136: final boolean sp2 = expected
137: .contains(AbstractProvider.STANDARD_PARALLEL_2);
138: this .belgium = belgium;
139: if (sp2) {
140: double phi2;
141: phi1 = doubleValue(expected,
142: AbstractProvider.STANDARD_PARALLEL_1, parameters);
143: ensureLatitudeInRange(AbstractProvider.STANDARD_PARALLEL_1,
144: phi1, true);
145: phi2 = doubleValue(expected,
146: AbstractProvider.STANDARD_PARALLEL_2, parameters);
147: if (Double.isNaN(phi2)) {
148: phi2 = phi1;
149: }
150: this .phi2 = phi2;
151: ensureLatitudeInRange(AbstractProvider.STANDARD_PARALLEL_2,
152: phi2, true);
153: } else {
154: if (belgium) {
155: throw new IllegalArgumentException();
156: }
157: // EPSG says the 1SP case uses the latitude of origin as the SP
158: phi1 = phi2 = latitudeOfOrigin;
159: }
160: // Compute constants
161: if (Math.abs(phi1 + phi2) < EPSILON) {
162: throw new IllegalArgumentException(Errors.format(
163: ErrorKeys.ANTIPODE_LATITUDES_$2, new Latitude(Math
164: .toDegrees(phi1)), new Latitude(Math
165: .toDegrees(phi2))));
166: }
167: final double cosphi1 = Math.cos(phi1);
168: final double sinphi1 = Math.sin(phi1);
169: final boolean secant = Math.abs(phi1 - phi2) > EPSILON; // Should be 'true' for 2SP case.
170: if (isSpherical) {
171: if (secant) {
172: n = Math.log(cosphi1 / Math.cos(phi2))
173: / Math.log(Math.tan((Math.PI / 4) + 0.5 * phi2)
174: / Math.tan((Math.PI / 4) + 0.5 * phi1));
175: } else {
176: n = sinphi1;
177: }
178: F = cosphi1
179: * Math.pow(Math.tan((Math.PI / 4) + 0.5 * phi1), n)
180: / n;
181: if (Math.abs(Math.abs(latitudeOfOrigin) - (Math.PI / 2)) >= EPSILON) {
182: rho0 = F
183: * Math.pow(Math.tan((Math.PI / 4) + 0.5
184: * latitudeOfOrigin), -n);
185: } else {
186: rho0 = 0.0;
187: }
188: } else {
189: final double m1 = msfn(sinphi1, cosphi1);
190: final double t1 = tsfn(phi1, sinphi1);
191: if (secant) {
192: final double sinphi2 = Math.sin(phi2);
193: final double m2 = msfn(sinphi2, Math.cos(phi2));
194: final double t2 = tsfn(phi2, sinphi2);
195: n = Math.log(m1 / m2) / Math.log(t1 / t2);
196: } else {
197: n = sinphi1;
198: }
199: F = m1 * Math.pow(t1, -n) / n;
200: if (Math.abs(Math.abs(latitudeOfOrigin) - (Math.PI / 2)) >= EPSILON) {
201: rho0 = F
202: * Math.pow(tsfn(latitudeOfOrigin, Math
203: .sin(latitudeOfOrigin)), n);
204: } else {
205: rho0 = 0.0;
206: }
207: }
208: }
209:
210: /**
211: * {@inheritDoc}
212: */
213: public ParameterValueGroup getParameterValues() {
214: final ParameterValueGroup values = super .getParameterValues();
215: final Collection expected = getParameterDescriptors()
216: .descriptors();
217: set(expected, AbstractProvider.STANDARD_PARALLEL_1, values,
218: phi1);
219: set(expected, AbstractProvider.STANDARD_PARALLEL_2, values,
220: phi2);
221: return values;
222: }
223:
224: /**
225: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
226: * (units in radians) and stores the result in {@code ptDst} (linear distance
227: * on a unit sphere).
228: */
229: protected Point2D transformNormalized(double x, double y,
230: Point2D ptDst) throws ProjectionException {
231: double rho;
232: //Snyder p. 108
233: if (Math.abs(Math.abs(y) - (Math.PI / 2)) < EPSILON) {
234: if (y * n <= 0) {
235: throw new ProjectionException(Errors.format(
236: ErrorKeys.POLE_PROJECTION_$1, new Latitude(Math
237: .toDegrees(y))));
238: } else {
239: rho = 0;
240: }
241: } else if (isSpherical) {
242: rho = F * Math.pow(Math.tan((Math.PI / 4) + 0.5 * y), -n);
243: } else {
244: rho = F * Math.pow(tsfn(y, Math.sin(y)), n);
245: }
246:
247: x *= n;
248: if (belgium) {
249: x -= BELGE_A;
250: }
251: y = rho0 - rho * Math.cos(x);
252: x = rho * Math.sin(x);
253:
254: if (ptDst != null) {
255: ptDst.setLocation(x, y);
256: return ptDst;
257: }
258: return new Point2D.Double(x, y);
259: }
260:
261: /**
262: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
263: * and stores the result in {@code ptDst}.
264: */
265: protected Point2D inverseTransformNormalized(double x, double y,
266: Point2D ptDst) throws ProjectionException {
267: double theta;
268: y = rho0 - y;
269: double rho = Math.sqrt(x * x + y * y); // Zero when the latitude is 90 degrees.
270: if (rho > EPSILON) {
271: if (n < 0) {
272: rho = -rho;
273: x = -x;
274: y = -y;
275: }
276: theta = Math.atan2(x, y);
277: if (belgium) {
278: theta += BELGE_A;
279: }
280: x = theta / n;
281: if (isSpherical) {
282: y = 2.0 * Math.atan(Math.pow(F / rho, 1.0 / n))
283: - (Math.PI / 2);
284: } else {
285: y = cphi2(Math.pow(rho / F, 1.0 / n));
286: }
287: } else {
288: x = 0.0;
289: y = n < 0 ? -(Math.PI / 2) : (Math.PI / 2);
290: }
291: if (ptDst != null) {
292: ptDst.setLocation(x, y);
293: return ptDst;
294: }
295: return new Point2D.Double(x, y);
296: }
297:
298: /**
299: * Returns a hash value for this projection.
300: */
301: public int hashCode() {
302: /*
303: * This code should be computed fast. Consequently, we do not use all fields
304: * in this object. Two {@code LambertConformal} objects with different
305: * {@link #phi1} and {@link #phi2} should compute a F value different enough.
306: */
307: final long code = Double.doubleToLongBits(F);
308: return ((int) code ^ (int) (code >>> 32)) + 37
309: * super .hashCode();
310: }
311:
312: /**
313: * Compares the specified object with this map projection for equality.
314: */
315: public boolean equals(final Object object) {
316: if (object == this ) {
317: // Slight optimization
318: return true;
319: }
320: if (super .equals(object)) {
321: final LambertConformal that = (LambertConformal) object;
322: return (this .belgium == that.belgium)
323: && equals(this .n, that.n) && equals(this .F, that.F)
324: && equals(this .rho0, that.rho0)
325: && equals(this .phi1, that.phi1)
326: && equals(this .phi2, that.phi2);
327: }
328: return false;
329: }
330: }
|