001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
006: * (C) 2002, Centre for Computational Geography
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.AffineTransform;
029: import java.awt.geom.Point2D;
030: import java.util.Collection;
031:
032: // OpenGIS dependencies
033: import org.opengis.parameter.ParameterDescriptor;
034: import org.opengis.parameter.ParameterDescriptorGroup;
035: import org.opengis.parameter.ParameterNotFoundException;
036: import org.opengis.parameter.ParameterValueGroup;
037: import org.opengis.referencing.operation.CylindricalProjection;
038: import org.opengis.referencing.operation.MathTransform;
039: import org.opengis.referencing.ReferenceIdentifier;
040:
041: // Geotools dependencies
042: import org.geotools.metadata.iso.citation.Citations;
043: import org.geotools.referencing.NamedIdentifier;
044: import org.geotools.referencing.operation.transform.ConcatenatedTransform;
045: import org.geotools.referencing.operation.transform.ProjectiveTransform;
046: import org.geotools.resources.i18n.VocabularyKeys;
047: import org.geotools.resources.i18n.Vocabulary;
048: import org.geotools.resources.i18n.ErrorKeys;
049: import org.geotools.resources.i18n.Errors;
050:
051: /**
052: * Transverse Mercator Projection (EPSG code 9807). This
053: * is a cylindrical projection, in which the cylinder has been rotated 90°.
054: * Instead of being tangent to the equator (or to an other standard latitude),
055: * it is tangent to a central meridian. Deformation are more important as we
056: * are going futher from the central meridian. The Transverse Mercator
057: * projection is appropriate for region wich have a greater extent north-south
058: * than east-west.
059: * <p>
060: *
061: * The elliptical equations used here are series approximations, and their accuracy
062: * decreases as points move farther from the central meridian of the projection.
063: * The forward equations here are accurate to a less than a mm ±10 degrees from
064: * the central meridian, a few mm ±15 degrees from the
065: * central meridian and a few cm ±20 degrees from the central meridian.
066: * The spherical equations are not approximations and should always give the
067: * correct values.
068: * <p>
069: *
070: * There are a number of versions of the transverse mercator projection
071: * including the Universal (UTM) and Modified (MTM) Transverses Mercator
072: * projections. In these cases the earth is divided into zones. For the UTM
073: * the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from
074: * 180 degrees longitude, and between lats 84 degrees North and 80
075: * degrees South. The central meridian is taken as the center of the zone
076: * and the latitude of origin is the equator. A scale factor of 0.9996 and
077: * false easting of 500000m is used for all zones and a false northing of 10000000m
078: * is used for zones in the southern hemisphere.
079: * <p>
080: *
081: * NOTE: formulas used below are not those of Snyder, but rather those
082: * from the {@code proj4} package of the USGS survey, which
083: * have been reproduced verbatim. USGS work is acknowledged here.
084: * <p>
085: *
086: * <strong>References:</strong><ul>
087: * <li> Proj-4.4.6 available at <A HREF="http://www.remotesensing.org/proj">www.remotesensing.org/proj</A><br>
088: * Relevent files are: {@code PJ_tmerc.c}, {@code pj_mlfn.c}, {@code pj_fwd.c} and {@code pj_inv.c}.</li>
089: * <li> John P. Snyder (Map Projections - A Working Manual,
090: * U.S. Geological Survey Professional Paper 1395, 1987).</li>
091: * <li> "Coordinate Conversions and Transformations including Formulas",
092: * EPSG Guidence Note Number 7, Version 19.</li>
093: * </ul>
094: *
095: * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Transverse Mercator projection on MathWorld</A>
096: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/transverse_mercator.html">"Transverse_Mercator" on RemoteSensing.org</A>
097: *
098: * @since 2.1
099: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/TransverseMercator.java $
100: * @version $Id: TransverseMercator.java 25262 2007-04-23 21:11:16Z desruisseaux $
101: * @author André Gosselin
102: * @author Martin Desruisseaux
103: * @author Rueben Schulz
104: */
105: public class TransverseMercator extends MapProjection {
106: /**
107: * Maximum number of iterations for iterative computations.
108: */
109: private static final int MAXIMUM_ITERATIONS = 15;
110:
111: /**
112: * Relative iteration precision used in the {@code mlfn} method.
113: * This overrides the value in the {@link MapProjection} class.
114: */
115: private static final double ITERATION_TOLERANCE = 1E-11;
116:
117: /**
118: * Maximum difference allowed when comparing real numbers.
119: */
120: private static final double EPSILON = 1E-6;
121:
122: /**
123: * Maximum difference allowed when comparing latitudes.
124: */
125: private static final double EPSILON_LATITUDE = 1E-10;
126:
127: /**
128: * A derived quantity of excentricity, computed by <code>e'˛ = (a˛-b˛)/b˛ = es/(1-es)</code>
129: * where <var>a</var> is the semi-major axis length and <var>b</bar> is the semi-minor axis
130: * length.
131: */
132: private final double esp;
133:
134: /**
135: * Meridian distance at the {@code latitudeOfOrigin}.
136: * Used for calculations for the ellipsoid.
137: */
138: private final double ml0;
139:
140: /**
141: * Constant needed for the <code>mlfn<code> method.
142: * Setup at construction time.
143: */
144: private final double en0, en1, en2, en3, en4;
145:
146: /**
147: * Constants used to calculate {@link #en0}, {@link #en1},
148: * {@link #en2}, {@link #en3}, {@link #en4}.
149: */
150: private static final double C00 = 1.0, C02 = 0.25, C04 = 0.046875,
151: C06 = 0.01953125, C08 = 0.01068115234375, C22 = 0.75,
152: C44 = 0.46875, C46 = 0.01302083333333333333,
153: C48 = 0.00712076822916666666, C66 = 0.36458333333333333333,
154: C68 = 0.00569661458333333333, C88 = 0.3076171875;
155:
156: /**
157: * Contants used for the forward and inverse transform for the eliptical
158: * case of the Transverse Mercator.
159: */
160: private static final double FC1 = 1.00000000000000000000000, // 1/1
161: FC2 = 0.50000000000000000000000, // 1/2
162: FC3 = 0.16666666666666666666666, // 1/6
163: FC4 = 0.08333333333333333333333, // 1/12
164: FC5 = 0.05000000000000000000000, // 1/20
165: FC6 = 0.03333333333333333333333, // 1/30
166: FC7 = 0.02380952380952380952380, // 1/42
167: FC8 = 0.01785714285714285714285; // 1/56
168:
169: /**
170: * Constructs a new map projection from the supplied parameters.
171: *
172: * @param parameters The parameter values in standard units.
173: * @throws ParameterNotFoundException if a mandatory parameter is missing.
174: */
175: protected TransverseMercator(final ParameterValueGroup parameters)
176: throws ParameterNotFoundException {
177: //Fetch parameters
178: super (parameters);
179:
180: // Compute constants
181: esp = excentricitySquared / (1.0 - excentricitySquared);
182:
183: double t;
184: en0 = C00
185: - excentricitySquared
186: * (C02 + excentricitySquared
187: * (C04 + excentricitySquared
188: * (C06 + excentricitySquared * C08)));
189: en1 = excentricitySquared
190: * (C22 - excentricitySquared
191: * (C04 + excentricitySquared
192: * (C06 + excentricitySquared * C08)));
193: en2 = (t = excentricitySquared * excentricitySquared)
194: * (C44 - excentricitySquared
195: * (C46 + excentricitySquared * C48));
196: en3 = (t *= excentricitySquared)
197: * (C66 - excentricitySquared * C68);
198: en4 = t * excentricitySquared * C88;
199: ml0 = mlfn(latitudeOfOrigin, Math.sin(latitudeOfOrigin), Math
200: .cos(latitudeOfOrigin));
201: }
202:
203: /**
204: * {@inheritDoc}
205: */
206: public ParameterDescriptorGroup getParameterDescriptors() {
207: return Provider.PARAMETERS;
208: }
209:
210: /**
211: * Transforms the specified (<var>λ</var>,<var>φ</var>) coordinates
212: * (units in radians) and stores the result in {@code ptDst} (linear distance
213: * on a unit sphere).
214: */
215: protected Point2D transformNormalized(double x, double y,
216: Point2D ptDst) throws ProjectionException {
217: double sinphi = Math.sin(y);
218: double cosphi = Math.cos(y);
219:
220: double t = (Math.abs(cosphi) > EPSILON) ? sinphi / cosphi : 0;
221: t *= t;
222: double al = cosphi * x;
223: double als = al * al;
224: al /= Math.sqrt(1.0 - excentricitySquared * sinphi * sinphi);
225: double n = esp * cosphi * cosphi;
226:
227: /* NOTE: meridinal distance at latitudeOfOrigin is always 0 */
228: y = (mlfn(y, sinphi, cosphi) - ml0 + sinphi
229: * al
230: * x
231: * FC2
232: * (1.0 + FC4
233: * als
234: * (5.0 - t + n * (9.0 + 4.0 * n) + FC6
235: * als
236: * (61.0 + t * (t - 58.0) + n
237: * (270.0 - 330.0 * t) + FC8
238: * als
239: * (1385.0 + t
240: * (t * (543.0 - t) - 3111.0))))));
241:
242: x = al
243: * (FC1 + FC3
244: * als
245: * (1.0 - t + n + FC5
246: * als
247: * (5.0 + t * (t - 18.0) + n
248: * (14.0 - 58.0 * t) + FC7
249: * als
250: * (61.0 + t
251: * (t * (179.0 - t) - 479.0)))));
252:
253: if (ptDst != null) {
254: ptDst.setLocation(x, y);
255: return ptDst;
256: }
257: return new Point2D.Double(x, y);
258: }
259:
260: /**
261: * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
262: * and stores the result in {@code ptDst}.
263: */
264: protected Point2D inverseTransformNormalized(double x, double y,
265: Point2D ptDst) throws ProjectionException {
266: double phi = inv_mlfn(ml0 + y);
267:
268: if (Math.abs(phi) >= (Math.PI / 2)) {
269: y = y < 0.0 ? -(Math.PI / 2) : (Math.PI / 2);
270: x = 0.0;
271: } else {
272: double sinphi = Math.sin(phi);
273: double cosphi = Math.cos(phi);
274: double t = (Math.abs(cosphi) > EPSILON) ? sinphi / cosphi
275: : 0.0;
276: double n = esp * cosphi * cosphi;
277: double con = 1.0 - excentricitySquared * sinphi * sinphi;
278: double d = x * Math.sqrt(con);
279: con *= t;
280: t *= t;
281: double ds = d * d;
282:
283: y = phi
284: - (con * ds / (1.0 - excentricitySquared))
285: * FC2
286: * (1.0 - ds
287: * FC4
288: * (5.0 + t * (3.0 - 9.0 * n) + n
289: * (1.0 - 4 * n) - ds
290: * FC6
291: * (61.0
292: + t
293: * (90.0 - 252.0 * n + 45.0 * t)
294: + 46.0 * n - ds
295: * FC8
296: * (1385.0 + t
297: * (3633.0 + t
298: * (4095.0 + 1574.0 * t))))));
299:
300: x = d
301: * (FC1 - ds
302: * FC3
303: * (1.0 + 2.0 * t + n - ds
304: * FC5
305: * (5.0 + t
306: * (28.0 + 24 * t + 8.0 * n)
307: + 6.0 * n - ds
308: * FC7
309: * (61.0 + t
310: * (662.0 + t
311: * (1320.0 + 720.0 * t))))))
312: / cosphi;
313: }
314:
315: if (ptDst != null) {
316: ptDst.setLocation(x, y);
317: return ptDst;
318: }
319: return new Point2D.Double(x, y);
320: }
321:
322: /**
323: * {@inheritDoc}
324: */
325: protected double getToleranceForAssertions(final double longitude,
326: final double latitude) {
327: if (Math.abs(longitude - centralMeridian) > 0.26) { // 15 degrees
328: // When far from the valid area, use a larger tolerance.
329: return 2.5;
330: } else if (Math.abs(longitude - centralMeridian) > 0.22) { // 12.5 degrees
331: return 1.0;
332: } else if (Math.abs(longitude - centralMeridian) > 0.17) { // 10 degrees
333: return 0.5;
334: }
335: // a normal tolerance
336: return 1E-6;
337: }
338:
339: /**
340: * Provides the transform equations for the spherical case of the
341: * TransverseMercator projection.
342: *
343: * @version $Id: TransverseMercator.java 25262 2007-04-23 21:11:16Z desruisseaux $
344: * @author André Gosselin
345: * @author Martin Desruisseaux
346: * @author Rueben Schulz
347: */
348: private static final class Spherical extends TransverseMercator {
349: /**
350: * Constructs a new map projection from the suplied parameters.
351: *
352: * @param parameters The parameter values in standard units.
353: * @throws ParameterNotFoundException if a mandatory parameter is missing.
354: */
355: protected Spherical(final ParameterValueGroup parameters)
356: throws ParameterNotFoundException {
357: super (parameters);
358: ensureSpherical();
359: }
360:
361: /**
362: * {@inheritDoc}
363: */
364: protected Point2D transformNormalized(double x, double y,
365: Point2D ptDst) throws ProjectionException {
366: // Compute using ellipsoidal formulas, for comparaison later.
367: final double normalizedLongitude = x;
368: assert (ptDst = super .transformNormalized(x, y, ptDst)) != null;
369:
370: double cosphi = Math.cos(y);
371: double b = cosphi * Math.sin(x);
372: if (Math.abs(Math.abs(b) - 1.0) <= EPSILON) {
373: throw new ProjectionException(Errors
374: .format(ErrorKeys.VALUE_TEND_TOWARD_INFINITY));
375: }
376:
377: //Using Snyder's equation for calculating y, instead of the one used in Proj4
378: //poential problems when y and x = 90 degrees, but behaves ok in tests
379: y = Math.atan2(Math.tan(y), Math.cos(x)) - latitudeOfOrigin; /* Snyder 8-3 */
380: x = 0.5 * Math.log((1.0 + b) / (1.0 - b)); /* Snyder 8-1 */
381:
382: assert checkTransform(
383: x,
384: y,
385: ptDst,
386: getToleranceForSphereAssertions(normalizedLongitude));
387: if (ptDst != null) {
388: ptDst.setLocation(x, y);
389: return ptDst;
390: }
391: return new Point2D.Double(x, y);
392: }
393:
394: /**
395: * {@inheritDoc}
396: */
397: protected Point2D inverseTransformNormalized(double x,
398: double y, Point2D ptDst) throws ProjectionException {
399: // Compute using ellipsoidal formulas, for comparaison later.
400: assert (ptDst = super .inverseTransformNormalized(x, y,
401: ptDst)) != null;
402:
403: double t = Math.exp(x);
404: double sinhX = 0.5 * (t - 1.0 / t); // sinh(x)
405: double cosD = Math.cos(latitudeOfOrigin + y);
406: double phi = Math.asin(Math.sqrt((1.0 - cosD * cosD)
407: / (1.0 + sinhX * sinhX)));
408: // correct for the fact that we made everything positive using sqrt(x*x)
409: y = ((y + latitudeOfOrigin) < 0.0) ? -phi : phi;
410: x = (Math.abs(sinhX) <= EPSILON && Math.abs(cosD) <= EPSILON) ? 0.0
411: : Math.atan2(sinhX, cosD);
412:
413: assert checkInverseTransform(x, y, ptDst,
414: getToleranceForSphereAssertions(x));
415: if (ptDst != null) {
416: ptDst.setLocation(x, y);
417: return ptDst;
418: }
419: return new Point2D.Double(x, y);
420: }
421:
422: /**
423: * Maximal error tolerated for assertions in the spherical case. When assertions
424: * are enabled, every projection using spherical formulas is followed by a projection
425: * using the ellipsical formulas, and the results are compared. If a distance greater
426: * than the tolerance level is found, then an {@link AssertionError} will be thrown.
427: * Subclasses may override this method if they need to relax the tolerance level.
428: * <p>
429: * The default implementation allows a larger tolerance when comparing spherical to
430: * elliptical calculations when we are far from the central meridian (elliptical
431: * calculations are not valid here).
432: *
433: * @param longitude The longitude standardized (in radians) with
434: * {@linkplain #centralMeridian} already removed from it.
435: * @return The tolerance level for assertions, in meters.
436: */
437: protected double getToleranceForSphereAssertions(
438: final double longitude) {
439: if (Math.abs(Math.abs(longitude) - Math.PI / 2.0) < EPSILON_LATITUDE) { // 90 degrees
440: // elliptical equations are at their worst here
441: return 1E+18;
442: }
443: if (Math.abs(longitude) > 0.26) { // 15 degrees
444: // When far from the valid area, use a very larger tolerance.
445: return 1E+6;
446: }
447: // a normal tolerance
448: return 1E-6;
449: }
450: }
451:
452: /**
453: * Calculates the meridian distance. This is the distance along the central
454: * meridian from the equator to {@code phi}. Accurate to < 1e-5 meters
455: * when used in conjuction with typical major axis values.
456: *
457: * @param phi latitude to calculate meridian distance for.
458: * @param sphi sin(phi).
459: * @param cphi cos(phi).
460: * @return meridian distance for the given latitude.
461: */
462: private final double mlfn(final double phi, double sphi, double cphi) {
463: cphi *= sphi;
464: sphi *= sphi;
465: return en0 * phi - cphi
466: * (en1 + sphi * (en2 + sphi * (en3 + sphi * (en4))));
467: }
468:
469: /**
470: * Calculates the latitude ({@code phi}) from a meridian distance.
471: * Determines phi to {@link #ITERATION_TOLERANCE} radians, about
472: * 1e-6 seconds.
473: *
474: * @param arg meridian distance to calulate latitude for.
475: * @return the latitude of the meridian distance.
476: * @throws ProjectionException if the iteration does not converge.
477: */
478: private final double inv_mlfn(double arg)
479: throws ProjectionException {
480: double s, t, phi, k = 1.0 / (1.0 - excentricitySquared);
481: int i;
482: phi = arg;
483: for (i = MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
484: if (--i < 0) {
485: throw new ProjectionException(Errors
486: .format(ErrorKeys.NO_CONVERGENCE));
487: }
488: s = Math.sin(phi);
489: t = 1.0 - excentricitySquared * s * s;
490: t = (mlfn(phi, s, Math.cos(phi)) - arg)
491: * (t * Math.sqrt(t)) * k;
492: phi -= t;
493: if (Math.abs(t) < ITERATION_TOLERANCE) {
494: return phi;
495: }
496: }
497: }
498:
499: /**
500: * Convenience method computing the zone code from the central meridian.
501: * Information about zones convention must be specified in argument. Two
502: * widely set of arguments are of Universal Transverse Mercator (UTM) and
503: * Modified Transverse Mercator (MTM) projections:<br>
504: * <p>
505: *
506: * UTM projection (zones numbered from 1 to 60):<br>
507: * <p>
508: * {@code getZone(-177, 6);}<br>
509: * <p>
510: * MTM projection (zones numbered from 1 to 120):<br>
511: * <p>
512: * {@code getZone(-52.5, -3);}<br>
513: *
514: * @param centralLongitudeZone1 Longitude in the middle of zone 1, in decimal degrees
515: * relative to Greenwich. Positive longitudes are toward east, and negative
516: * longitudes toward west.
517: * @param zoneWidth Number of degrees of longitudes in one zone. A positive value
518: * means that zones are numbered from west to east (i.e. in the direction of
519: * positive longitudes). A negative value means that zones are numbered from
520: * east to west.
521: * @return The zone number. First zone is numbered 1.
522: */
523: private int getZone(final double centralLongitudeZone1,
524: final double zoneWidth) {
525: final double zoneCount = Math.abs(360 / zoneWidth);
526: double t;
527: t = centralLongitudeZone1 - 0.5 * zoneWidth; // Longitude at the beginning of the first zone.
528: t = Math.toDegrees(centralMeridian) - t; // Degrees of longitude between the central longitude and longitude 1.
529: t = Math.floor(t / zoneWidth + EPSILON); // Number of zones between the central longitude and longitude 1.
530: t -= zoneCount * Math.floor(t / zoneCount); // If negative, bring back to the interval 0 to (zoneCount-1).
531: return ((int) t) + 1;
532: }
533:
534: /**
535: * Convenience method returning the meridian in the middle of current zone. This meridian is
536: * typically the central meridian. This method may be invoked to make sure that the central
537: * meridian is correctly set.
538: *
539: * @param centralLongitudeZone1 Longitude in the middle of zone 1, in decimal degrees
540: * relative to Greenwich. Positive longitudes are toward east, and negative
541: * longitudes toward west.
542: * @param zoneWidth Number of degrees of longitudes in one zone. A positive value
543: * means that zones are numbered from west to east (i.e. in the direction of
544: * positive longitudes). A negative value means that zones are numbered from
545: * east to west.
546: * @return The central meridian.
547: */
548: private double getCentralMedirian(
549: final double centralLongitudeZone1, final double zoneWidth) {
550: double t;
551: t = centralLongitudeZone1
552: + (getZone(centralLongitudeZone1, zoneWidth) - 1)
553: * zoneWidth;
554: t -= 360 * Math.floor((t + 180) / 360); // Bring back into [-180..+180] range.
555: return t;
556: }
557:
558: /**
559: * Convenience method computing the zone code from the central meridian.
560: *
561: * @return The zone number, using the scalefactor and false easting
562: * to decide if this is a UTM or MTM case. Returns 0 if the
563: * case of the projection cannot be determined.
564: */
565: public int getZone() {
566: // UTM
567: if (scaleFactor == 0.9996 && falseEasting == 500000.0) {
568: return getZone(-177, 6);
569: }
570: // MTM
571: if (scaleFactor == 0.9999 && falseEasting == 304800.0) {
572: return getZone(-52.5, -3);
573: }
574: // unknown
575: throw new IllegalStateException(Errors
576: .format(ErrorKeys.UNKNOW_PROJECTION_TYPE));
577: }
578:
579: /**
580: * Convenience method returning the meridian in the middle of current zone. This meridian is
581: * typically the central meridian. This method may be invoked to make sure that the central
582: * meridian is correctly set.
583: *
584: * @return The central meridian, using the scalefactor and false easting
585: * to decide if this is a UTM or MTM case. Returns {@link Double#NaN}
586: * if the case of the projection cannot be determined.
587: */
588: public double getCentralMeridian() {
589: // UTM
590: if (scaleFactor == 0.9996 && falseEasting == 500000.0) {
591: return getCentralMedirian(-177, 6);
592: }
593: // MTM
594: if (scaleFactor == 0.9999 && falseEasting == 304800.0) {
595: return getCentralMedirian(-52.5, -3);
596: }
597: // unknown
598: throw new IllegalStateException(Errors
599: .format(ErrorKeys.UNKNOW_PROJECTION_TYPE));
600: }
601:
602: /**
603: * Returns a hash value for this projection.
604: */
605: public int hashCode() {
606: final long code = Double.doubleToLongBits(ml0);
607: return ((int) code ^ (int) (code >>> 32)) + 37
608: * super .hashCode();
609: }
610:
611: /**
612: * Compares the specified object with this map projection for equality.
613: */
614: public boolean equals(final Object object) {
615: if (object == this ) {
616: // Slight optimization
617: return true;
618: }
619: // Relevant parameters are already compared in MapProjection
620: return super .equals(object);
621: }
622:
623: //////////////////////////////////////////////////////////////////////////////////////////
624: //////////////////////////////////////////////////////////////////////////////////////////
625: //////// ////////
626: //////// PROVIDERS ////////
627: //////// ////////
628: //////////////////////////////////////////////////////////////////////////////////////////
629: //////////////////////////////////////////////////////////////////////////////////////////
630:
631: /**
632: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
633: * provider} for a {@linkplain TransverseMercator Transverse Mercator} projection (EPSG code
634: * 9807).
635: *
636: * @since 2.1
637: * @version $Id: TransverseMercator.java 25262 2007-04-23 21:11:16Z desruisseaux $
638: * @author Martin Desruisseaux
639: * @author Rueben Schulz
640: *
641: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
642: */
643: public static class Provider extends AbstractProvider {
644: /**
645: * Returns a descriptor group for the specified parameters.
646: */
647: static ParameterDescriptorGroup createDescriptorGroup(
648: final ReferenceIdentifier[] identifiers) {
649: return createDescriptorGroup(
650: identifiers,
651: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
652: CENTRAL_MERIDIAN, LATITUDE_OF_ORIGIN,
653: SCALE_FACTOR, FALSE_EASTING, FALSE_NORTHING });
654: }
655:
656: /**
657: * The parameters group.
658: */
659: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new NamedIdentifier[] {
660: new NamedIdentifier(Citations.OGC,
661: "Transverse_Mercator"),
662: new NamedIdentifier(Citations.EPSG,
663: "Transverse Mercator"),
664: new NamedIdentifier(Citations.EPSG, "Gauss-Kruger"),
665: new NamedIdentifier(Citations.EPSG, "9807"),
666: new NamedIdentifier(Citations.GEOTIFF,
667: "CT_TransverseMercator"),
668: new NamedIdentifier(Citations.ESRI,
669: "Transverse_Mercator"),
670: new NamedIdentifier(Citations.ESRI, "Gauss_Kruger"),
671: new NamedIdentifier(
672: Citations.GEOTOOLS,
673: Vocabulary
674: .formatInternational(VocabularyKeys.TRANSVERSE_MERCATOR_PROJECTION)) });
675:
676: /**
677: * Constructs a new provider.
678: */
679: public Provider() {
680: super (PARAMETERS);
681: }
682:
683: /**
684: * Constructs a new provider with the specified parameters.
685: */
686: Provider(final ParameterDescriptorGroup descriptor) {
687: super (descriptor);
688: }
689:
690: /**
691: * Returns the operation type for this map projection.
692: */
693: public Class getOperationType() {
694: return CylindricalProjection.class;
695: }
696:
697: /**
698: * Creates a transform from the specified group of parameter values.
699: *
700: * @param parameters The group of parameter values.
701: * @return The created math transform.
702: * @throws ParameterNotFoundException if a required parameter was not found.
703: */
704: protected MathTransform createMathTransform(
705: final ParameterValueGroup parameters)
706: throws ParameterNotFoundException {
707: if (isSpherical(parameters)) {
708: return new Spherical(parameters);
709: } else {
710: return new TransverseMercator(parameters);
711: }
712: }
713: }
714:
715: /**
716: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
717: * provider} for a South Orientated {@linkplain TransverseMercator Transverse Mercator}
718: * projection (EPSG code 9808). Note that at the contrary of what this class name suggest,
719: * the projected coordinates are still increasing toward North. This is because all
720: * {@linkplain MapProjection map projections} in Geotools must complies with standard axis
721: * orientations. The real axis flip is performed by the CRS framework outside this package.
722: * See "<cite>Axis units and orientation</cite>" in
723: * {@linkplain org.geotools.referencing.operation.projection package description} for details.
724: * <p>
725: * The usual Transverse Mercator formulas are:
726: * <p>
727: * <ul>
728: * <li><var>easting</var> = <var>{@linkplain #falseEasting false easting}</var> + <var>px</var></li>
729: * <li><var>northing</var> = <var>{@linkplain #falseNorthing false northing}</var> + <var>py</var></li>
730: * </ul>
731: * <p>
732: * The Transverse Mercator South Orientated Projection formulas are:
733: * <p>
734: * <ul>
735: * <li><var>westing</var> = <var>{@linkplain #falseEasting false easting}</var> - <var>px</var></li>
736: * <li><var>southing</var> = <var>{@linkplain #falseNorthing false northing}</var> - <var>py</var></li>
737: * </ul>
738: * <p>
739: * Where the <var>px</var> and <var>py</var> terms are the same in both cases.
740: * Transforms created by this provider actually computes
741: * (<var>easting</var>,<var>northing</var>) = (-<var>westing</var>,-<var>southing</var>).
742: * This is equivalent to a {@link TransverseMercator} projection with
743: * {@link #falseEasting falseEasting} and {@link #falseNorthing falseNorthing} sign reverted.
744: * This operation is implemented as a concatenation of a North-orientated transverse mercator
745: * projection with an affine transform for (<var>false easting</var>,<var>false northing</var>)
746: * correction.
747: *
748: * @since 2.2
749: * @version $Id: TransverseMercator.java 25262 2007-04-23 21:11:16Z desruisseaux $
750: * @author Martin Desruisseaux
751: *
752: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
753: */
754: public static class Provider_SouthOrientated extends Provider {
755: /**
756: * Constructs a new provider.
757: */
758: public Provider_SouthOrientated() {
759: super (createDescriptorGroup(new NamedIdentifier[] {
760: new NamedIdentifier(Citations.EPSG,
761: "Transverse Mercator (South Orientated)"),
762: new NamedIdentifier(Citations.EPSG, "9808") }));
763: }
764:
765: /**
766: * Creates a transform from the specified group of parameter values.
767: *
768: * @param parameters The group of parameter values.
769: * @return The created math transform.
770: * @throws ParameterNotFoundException if a required parameter was not found.
771: */
772: public MathTransform createMathTransform(
773: final ParameterValueGroup parameters)
774: throws ParameterNotFoundException {
775: final MapProjection projection = (MapProjection) super
776: .createMathTransform(parameters);
777: if (projection.falseEasting == 0
778: && projection.falseNorthing == 0) {
779: return projection;
780: }
781: final AffineTransform step = AffineTransform
782: .getTranslateInstance(-2 * projection.falseEasting,
783: -2 * projection.falseNorthing);
784: return ConcatenatedTransform.create(ProjectiveTransform
785: .create(step), projection);
786: }
787: }
788: }
|