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.ParameterDescriptor;
033: import org.opengis.parameter.ParameterDescriptorGroup;
034: import org.opengis.parameter.ParameterNotFoundException;
035: import org.opengis.parameter.ParameterValueGroup;
036: import org.opengis.referencing.operation.MathTransform;
037:
038: // Geotools dependencies
039: import org.geotools.resources.XMath;
040: import org.geotools.resources.i18n.Errors;
041: import org.geotools.resources.i18n.ErrorKeys;
042: import org.geotools.referencing.NamedIdentifier;
043: import org.geotools.metadata.iso.citation.Citations;
044:
045: /**
046: * Provides the transform equations for the Oblique Stereographic (EPSG code 9809).
047: * The formulas used below are not from the EPSG, but rather those of the
048: * "Oblique Stereographic Alternative" in the {@code libproj4} package
049: * written by Gerald Evenden. His work is acknowledged here and greatly appreciated.
050: * <p>
051: *
052: * The forward equations used in {@code libproj4} are the same as those given in the
053: * UNB reports for the Double Stereographic. The inverse equations are similar,
054: * but use different methods to iterate for the latitude.
055: * <p>
056: *
057: * <strong>References:</strong><ul>
058: * <li>{@code libproj4} is available at
059: * <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/">libproj4 Miscellanea</A><br>
060: * Relevent files are: {@code PJ_sterea.c}, {@code pj_gauss.c},
061: * {@code pj_fwd.c}, {@code pj_inv.c} and {@code lib_proj.h}</li>
062: * <li>Gerald Evenden. <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/sterea.pdf">
063: * "Supplementary PROJ.4 Notes - Oblique Stereographic Alternative"</A></li>
064: * <li>"Coordinate Conversions and Transformations including Formulas",
065: * EPSG Guidence Note Number 7, Version 19.</li>
066: * <li>Krakiwsky, E.J., D.B. Thomson, and R.R. Steeves. 1977. A Manual
067: * For Geodetic Coordinate Transformations in the Maritimes.
068: * Geodesy and Geomatics Engineering, UNB. Technical Report No. 48.</li>
069: * <li>Thomson, D.B., M.P. Mepham and R.R. Steeves. 1977.
070: * The Stereographic Double Projection.
071: * Surveying Engineering, University of New Brunswick. Technical Report No. 46.</li>
072: * </ul>
073: *
074: * @since 2.4
075: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/ObliqueStereographic.java $
076: * @version $Id: ObliqueStereographic.java 24576 2007-02-24 00:07:40Z desruisseaux $
077: * @author Rueben Schulz
078: */
079: public class ObliqueStereographic extends StereographicUSGS {
080: /*
081: * The tolerance used for the inverse iteration. This is smaller
082: * than the tolerance in the {@code StereographicOblique} superclass.
083: */
084: private static final double ITERATION_TOLERANCE = 1E-14;
085:
086: /**
087: * Maximum number of iterations for iterative computations.
088: */
089: private static final int MAXIMUM_ITERATIONS = 15;
090:
091: /**
092: * Maximum difference allowed when comparing real numbers.
093: */
094: private static final double EPSILON = 1E-6;
095:
096: /*
097: * Contstants used in the forward and inverse gauss methods.
098: */
099: private final double C, K, ratexp;
100:
101: /*
102: * Constants for the EPSG stereographic transform.
103: */
104: private final double phic0, cosc0, sinc0, R2;
105:
106: /**
107: * Constructs an oblique stereographic projection (EPSG equations).
108: *
109: * @param parameters The group of parameter values.
110: * @throws ParameterNotFoundException if a required parameter was not found.
111: */
112: protected ObliqueStereographic(final ParameterValueGroup parameters)
113: throws ParameterNotFoundException {
114: this (parameters, Provider.PARAMETERS);
115: }
116:
117: /**
118: * Constructs an oblique stereographic projection (EPSG equations).
119: *
120: * @param parameters The group of parameter values.
121: * @param descriptor The expected parameter descriptor.
122: * @throws ParameterNotFoundException if a required parameter was not found.
123: */
124: ObliqueStereographic(final ParameterValueGroup parameters,
125: final ParameterDescriptorGroup descriptor)
126: throws ParameterNotFoundException {
127: super (parameters, descriptor);
128:
129: // Compute constants
130: final double sphi = Math.sin(latitudeOfOrigin);
131: double cphi = Math.cos(latitudeOfOrigin);
132: cphi *= cphi;
133: R2 = 2.0 * Math.sqrt(1 - excentricitySquared)
134: / (1 - excentricitySquared * sphi * sphi);
135: C = Math.sqrt(1. + excentricitySquared * cphi * cphi
136: / (1. - excentricitySquared));
137: phic0 = Math.asin(sphi / C);
138: sinc0 = Math.sin(phic0);
139: cosc0 = Math.cos(phic0);
140: ratexp = 0.5 * C * excentricity;
141: K = Math.tan(0.5 * phic0 + Math.PI / 4)
142: / (Math.pow(Math.tan(0.5 * latitudeOfOrigin + Math.PI
143: / 4), C) * srat(excentricity * sphi, ratexp));
144: }
145:
146: /**
147: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
148: * (units in radians) and stores the result in {@code ptDst} (linear distance
149: * on a unit sphere).
150: */
151: protected Point2D transformNormalized(double x, double y,
152: Point2D ptDst) throws ProjectionException {
153: // Compute using USGS formulas, for comparaison later.
154: assert (ptDst = super .transformNormalized(x, y, ptDst)) != null;
155:
156: y = 2.0
157: * Math.atan(K
158: * Math.pow(Math.tan(0.5 * y + Math.PI / 4), C)
159: * srat(excentricity * Math.sin(y), ratexp))
160: - Math.PI / 2;
161: x *= C;
162: final double sinc = Math.sin(y);
163: final double cosc = Math.cos(y);
164: final double cosl = Math.cos(x);
165: final double k = R2
166: / (1.0 + sinc0 * sinc + cosc0 * cosc * cosl);
167: x = k * cosc * Math.sin(x);
168: y = k * (cosc0 * sinc - sinc0 * cosc * cosl);
169:
170: assert checkTransform(x, y, ptDst, 0.1);
171: if (ptDst != null) {
172: ptDst.setLocation(x, y);
173: return ptDst;
174: }
175: return new Point2D.Double(x, y);
176: }
177:
178: /**
179: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
180: * and stores the result in {@code ptDst}.
181: */
182: protected Point2D inverseTransformNormalized(double x, double y,
183: Point2D ptDst) throws ProjectionException {
184: // Compute using USGS formulas, for comparaison later.
185: assert (ptDst = super .inverseTransformNormalized(x, y, ptDst)) != null;
186:
187: final double rho = XMath.hypot(x, y);
188: if (Math.abs(rho) < EPSILON) {
189: x = 0.0;
190: y = phic0;
191: } else {
192: final double ce = 2.0 * Math.atan2(rho, R2);
193: final double sinc = Math.sin(ce);
194: final double cosc = Math.cos(ce);
195: x = Math.atan2(x * sinc, rho * cosc0 * cosc - y * sinc0
196: * sinc);
197: y = (cosc * sinc0) + (y * sinc * cosc0 / rho);
198:
199: if (Math.abs(y) >= 1.0) {
200: y = (y < 0.0) ? -Math.PI / 2.0 : Math.PI / 2.0;
201: } else {
202: y = Math.asin(y);
203: }
204: }
205:
206: // Begin pj_inv_gauss(...) method inlined
207: x /= C;
208: double num = Math.pow(Math.tan(0.5 * y + Math.PI / 4.0) / K,
209: 1.0 / C);
210: for (int i = MAXIMUM_ITERATIONS;;) {
211: double phi = 2.0
212: * Math.atan(num
213: * srat(excentricity * Math.sin(y), -0.5
214: * excentricity)) - Math.PI / 2.0;
215: if (Math.abs(phi - y) < ITERATION_TOLERANCE) {
216: break;
217: }
218: y = phi;
219: if (--i < 0) {
220: throw new ProjectionException(Errors
221: .format(ErrorKeys.NO_CONVERGENCE));
222: }
223: }
224: // End pj_inv_gauss(...) method inlined
225:
226: // TODO: the tolerance in the following assertion is quite large for
227: // an angle in radians. We should check if this is normal.
228: assert checkInverseTransform(x, y, ptDst, 0.01);
229: if (ptDst != null) {
230: ptDst.setLocation(x, y);
231: return ptDst;
232: }
233: return new Point2D.Double(x, y);
234: }
235:
236: /**
237: * A simple function used by the transforms.
238: */
239: private static double srat(double esinp, double exp) {
240: return Math.pow((1.0 - esinp) / (1.0 + esinp), exp);
241: }
242:
243: //////////////////////////////////////////////////////////////////////////////////////////
244: //////////////////////////////////////////////////////////////////////////////////////////
245: //////// ////////
246: //////// PROVIDERS ////////
247: //////// ////////
248: //////////////////////////////////////////////////////////////////////////////////////////
249: //////////////////////////////////////////////////////////////////////////////////////////
250:
251: /**
252: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
253: * provider} for a stereographic projection of any kind. The equations used are the one from
254: * EPSG.
255: *
256: * @since 2.4
257: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/ObliqueStereographic.java $
258: * @version $Id: ObliqueStereographic.java 24576 2007-02-24 00:07:40Z desruisseaux $
259: * @author Rueben Schulz
260: *
261: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
262: */
263: public static final class Provider extends Stereographic.Provider {
264: /**
265: * The parameters group.
266: */
267: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
268: new NamedIdentifier[] {
269: new NamedIdentifier(Citations.OGC,
270: "Oblique_Stereographic"),
271: new NamedIdentifier(Citations.EPSG,
272: "Oblique Stereographic"),
273: new NamedIdentifier(Citations.EPSG, "Roussilhe"),
274: new NamedIdentifier(Citations.EPSG, "9809"),
275: new NamedIdentifier(Citations.GEOTIFF,
276: "CT_ObliqueStereographic"),
277: new NamedIdentifier(Citations.ESRI,
278: "Double_Stereographic"),
279: new NamedIdentifier(Citations.GEOTOOLS, NAME) },
280: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
281: CENTRAL_MERIDIAN, LATITUDE_OF_ORIGIN,
282: SCALE_FACTOR, FALSE_EASTING, FALSE_NORTHING });
283:
284: /**
285: * Constructs a new provider.
286: */
287: public Provider() {
288: super (PARAMETERS);
289: }
290:
291: /**
292: * Creates the general case.
293: */
294: //@Override
295: MathTransform createMathTransform(
296: final ParameterValueGroup parameters,
297: final ParameterDescriptorGroup descriptor)
298: throws ParameterNotFoundException {
299: return new ObliqueStereographic(parameters, descriptor);
300: }
301: }
302: }
|