001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
006: * (C) 2003, Gerald I. Evenden
007: * (C) 2001, Institut de Recherche pour le Développement
008: * (C) 2000, Frank Warmerdam
009: * (C) 1999, Fisheries and Oceans Canada
010: *
011: * This library is free software; you can redistribute it and/or
012: * modify it under the terms of the GNU Lesser General Public
013: * License as published by the Free Software Foundation; either
014: * version 2.1 of the License, or (at your option) any later version.
015: *
016: * This library is distributed in the hope that it will be useful,
017: * but WITHOUT ANY WARRANTY; without even the implied warranty of
018: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
019: * Lesser General Public License for more details.
020: *
021: * This package contains formulas from the PROJ package of USGS.
022: * USGS's work is fully acknowledged here. This derived work has
023: * been relicensed under LGPL with Frank Warmerdam's permission.
024: */
025: package org.geotools.referencing.operation.projection;
026:
027: // J2SE dependencies and extensions
028: import java.awt.geom.Point2D;
029:
030: // OpenGIS dependencies
031: import org.opengis.parameter.ParameterValueGroup;
032: import org.opengis.parameter.ParameterDescriptorGroup;
033: import org.opengis.parameter.ParameterNotFoundException;
034:
035: // Geotools dependencies
036: import org.geotools.resources.XMath;
037: import org.geotools.resources.i18n.Errors;
038: import org.geotools.resources.i18n.ErrorKeys;
039:
040: /**
041: * The USGS oblique/equatorial case of the Stereographic projection. This is similar but
042: * <strong>NOT</strong> equal to EPSG code 9809 ({@code "Oblique_Stereographic"} EPSG name).
043: * The later is rather implemented by {@link ObliqueStereographic}.
044: * <p>
045: * This class is not public in order to keep names that closely match the ones in common usage
046: * (i.e. this projection is called just "Stereographic" in ESRI). Furthermore, the "USGS" name
047: * is not really accurate for a class to be extended by {@link ObliqueStereographic}.
048: *
049: * @since 2.4
050: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/StereographicUSGS.java $
051: * @version $Id: StereographicUSGS.java 24576 2007-02-24 00:07:40Z desruisseaux $
052: * @author André Gosselin
053: * @author Martin Desruisseaux
054: * @author Rueben Schulz
055: */
056: class StereographicUSGS extends Stereographic {
057: /**
058: * Maximum number of iterations for iterative computations.
059: */
060: private static final int MAXIMUM_ITERATIONS = 15;
061:
062: /**
063: * Difference allowed in iterative computations.
064: */
065: private static final double ITERATION_TOLERANCE = 1E-10;
066:
067: /**
068: * Maximum difference allowed when comparing real numbers.
069: */
070: private static final double EPSILON = 1E-6;
071:
072: /**
073: * Constants used for the oblique projections. All those constants are completly determined by
074: * {@link #latitudeOfOrigin}. Concequently, there is no need to test them in {@link #hashCode}
075: * or {@link #equals} methods.
076: */
077: final double k0, sinphi0, cosphi0, chi1, sinChi1, cosChi1;
078:
079: /**
080: * Constructs an oblique stereographic projection (USGS equations).
081: *
082: * @param parameters The group of parameter values.
083: * @throws ParameterNotFoundException if a required parameter was not found.
084: */
085: protected StereographicUSGS(final ParameterValueGroup parameters)
086: throws ParameterNotFoundException {
087: this (parameters, Provider.PARAMETERS);
088: }
089:
090: /**
091: * Constructs an oblique stereographic projection (USGS equations).
092: *
093: * @param parameters The group of parameter values.
094: * @param descriptor The expected parameter descriptor.
095: * @throws ParameterNotFoundException if a required parameter was not found.
096: */
097: StereographicUSGS(final ParameterValueGroup parameters,
098: final ParameterDescriptorGroup descriptor)
099: throws ParameterNotFoundException {
100: super (parameters, descriptor);
101: if (Math.abs(latitudeOfOrigin) < EPSILON) { // Equatorial
102: latitudeOfOrigin = 0;
103: cosphi0 = 1.0;
104: sinphi0 = 0.0;
105: chi1 = 0.0;
106: cosChi1 = 1.0;
107: sinChi1 = 0.0;
108: } else { // Oblique
109: cosphi0 = Math.cos(latitudeOfOrigin);
110: sinphi0 = Math.sin(latitudeOfOrigin);
111: chi1 = 2.0 * Math.atan(ssfn(latitudeOfOrigin, sinphi0))
112: - (Math.PI / 2);
113: cosChi1 = Math.cos(chi1);
114: sinChi1 = Math.sin(chi1);
115: }
116: // part of (14 - 15)
117: k0 = 2.0 * msfn(sinphi0, cosphi0);
118: }
119:
120: /**
121: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
122: * (units in radians) and stores the result in {@code ptDst} (linear distance
123: * on a unit sphere).
124: */
125: protected Point2D transformNormalized(double x, double y,
126: Point2D ptDst) throws ProjectionException {
127: final double chi = 2.0 * Math.atan(ssfn(y, Math.sin(y)))
128: - (Math.PI / 2);
129: final double sinChi = Math.sin(chi);
130: final double cosChi = Math.cos(chi);
131: final double cosChi_cosLon = cosChi * Math.cos(x);
132: final double A = k0 / cosChi1
133: / (1 + sinChi1 * sinChi + cosChi1 * cosChi_cosLon);
134: x = A * cosChi * Math.sin(x);
135: y = A * (cosChi1 * sinChi - sinChi1 * cosChi_cosLon);
136:
137: if (ptDst != null) {
138: ptDst.setLocation(x, y);
139: return ptDst;
140: }
141: return new Point2D.Double(x, y);
142: }
143:
144: /**
145: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
146: * and stores the result in {@code ptDst}.
147: */
148: protected Point2D inverseTransformNormalized(double x, double y,
149: Point2D ptDst) throws ProjectionException {
150: final double rho = Math.sqrt(x * x + y * y);
151: final double ce = 2.0 * Math.atan2(rho * cosChi1, k0);
152: final double cosce = Math.cos(ce);
153: final double since = Math.sin(ce);
154: final boolean rhoIs0 = Math.abs(rho) < EPSILON;
155: final double chi = rhoIs0 ? chi1 : Math.asin(cosce * sinChi1
156: + (y * since * cosChi1 / rho));
157: final double tp = Math.tan(Math.PI / 4.0 + chi / 2.0);
158:
159: // parts of (21-36) used to calculate longitude
160: final double t = x * since;
161: final double ct = rho * cosChi1 * cosce - y * sinChi1 * since;
162:
163: // Compute latitude using iterative technique (3-4)
164: final double halfe = excentricity / 2.0;
165: double phi0 = chi;
166: for (int i = MAXIMUM_ITERATIONS;;) {
167: final double esinphi = excentricity * Math.sin(phi0);
168: final double phi = 2
169: * Math.atan(tp
170: * Math.pow((1 + esinphi) / (1 - esinphi),
171: halfe)) - (Math.PI / 2);
172: if (Math.abs(phi - phi0) < ITERATION_TOLERANCE) {
173: // TODO: checking rho may be redundant
174: x = rhoIs0
175: || (Math.abs(t) < EPSILON && Math.abs(ct) < EPSILON) ? 0.0
176: : Math.atan2(t, ct);
177: y = phi;
178: break;
179: }
180: phi0 = phi;
181: if (--i < 0) {
182: throw new ProjectionException(Errors
183: .format(ErrorKeys.NO_CONVERGENCE));
184: }
185: }
186:
187: if (ptDst != null) {
188: ptDst.setLocation(x, y);
189: return ptDst;
190: }
191: return new Point2D.Double(x, y);
192: }
193:
194: /**
195: * Maximal error (in metres) tolerated for assertions, if enabled.
196: */
197: //@Override
198: protected double getToleranceForAssertions(final double longitude,
199: final double latitude) {
200: final double delta = Math.abs(longitude - centralMeridian) / 2
201: + Math.abs(latitude - latitudeOfOrigin);
202: if (delta > 40) {
203: return 0.5;
204: }
205: if (delta > 15) {
206: return 0.1;
207: }
208: return super .getToleranceForAssertions(longitude, latitude);
209: }
210:
211: /**
212: * Computes part of function (3-1) from Snyder.
213: */
214: final double ssfn(double phi, double sinphi) {
215: sinphi *= excentricity;
216: return Math.tan((Math.PI / 4.0) + phi / 2.0)
217: * Math.pow((1 - sinphi) / (1 + sinphi),
218: excentricity / 2.0);
219: }
220:
221: /**
222: * Provides the transform equations for the spherical case of the
223: * Stereographic projection.
224: *
225: * @version $Id: StereographicUSGS.java 24576 2007-02-24 00:07:40Z desruisseaux $
226: * @author Martin Desruisseaux
227: * @author Rueben Schulz
228: */
229: static final class Spherical extends StereographicUSGS {
230: /**
231: * A constant used in the transformations. This constant hides the {@code k0}
232: * constant from the ellipsoidal case. The spherical and ellipsoidal {@code k0}
233: * are not computed in the same way, and we preserve the ellipsoidal {@code k0}
234: * in {@link Stereographic} in order to allow assertions to work.
235: */
236: private static final double k0 = 2;
237:
238: /**
239: * Constructs a spherical oblique stereographic projection.
240: *
241: * @param parameters The group of parameter values.
242: * @param descriptor The expected parameter descriptor.
243: * @throws ParameterNotFoundException if a required parameter was not found.
244: */
245: Spherical(final ParameterValueGroup parameters,
246: final ParameterDescriptorGroup descriptor)
247: throws ParameterNotFoundException {
248: super (parameters, descriptor);
249: ensureSpherical();
250: }
251:
252: /**
253: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
254: * (units in radians) and stores the result in {@code ptDst} (linear distance
255: * on a unit sphere).
256: */
257: protected Point2D transformNormalized(double x, double y,
258: Point2D ptDst) throws ProjectionException {
259: // Compute using ellipsoidal formulas, for comparaison later.
260: assert (ptDst = super .transformNormalized(x, y, ptDst)) != null;
261:
262: final double coslat = Math.cos(y);
263: final double sinlat = Math.sin(y);
264: final double coslon = Math.cos(x);
265: double f = 1.0 + sinphi0 * sinlat + cosphi0 * coslat
266: * coslon; // (21-4)
267: if (f < EPSILON) {
268: throw new ProjectionException(Errors
269: .format(ErrorKeys.VALUE_TEND_TOWARD_INFINITY));
270: }
271: f = k0 / f;
272: x = f * coslat * Math.sin(x); // (21-2)
273: y = f * (cosphi0 * sinlat - sinphi0 * coslat * coslon); // (21-3)
274:
275: assert checkTransform(x, y, ptDst);
276: if (ptDst != null) {
277: ptDst.setLocation(x, y);
278: return ptDst;
279: }
280: return new Point2D.Double(x, y);
281: }
282:
283: /**
284: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
285: * and stores the result in {@code ptDst}.
286: */
287: protected Point2D inverseTransformNormalized(double x,
288: double y, Point2D ptDst) throws ProjectionException {
289: // Compute using ellipsoidal formulas, for comparaison later.
290: assert (ptDst = super .inverseTransformNormalized(x, y,
291: ptDst)) != null;
292:
293: final double rho = Math.sqrt(x * x + y * y);
294: if (Math.abs(rho) < EPSILON) {
295: y = latitudeOfOrigin;
296: x = 0.0;
297: } else {
298: final double c = 2.0 * Math.atan(rho / k0);
299: final double cosc = Math.cos(c);
300: final double sinc = Math.sin(c);
301: final double ct = rho * cosphi0 * cosc - y * sinphi0
302: * sinc; // (20-15)
303: final double t = x * sinc; // (20-15)
304: y = Math
305: .asin(cosc * sinphi0 + y * sinc * cosphi0 / rho); // (20-14)
306: x = (Math.abs(ct) < EPSILON && Math.abs(t) < EPSILON) ? 0.0
307: : Math.atan2(t, ct);
308: }
309:
310: assert checkInverseTransform(x, y, ptDst);
311: if (ptDst != null) {
312: ptDst.setLocation(x, y);
313: return ptDst;
314: }
315: return new Point2D.Double(x, y);
316: }
317: }
318: }
|