001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: *
005: * (C) 2005-2006, Geotools Project Managment Committee (PMC)
006: * (C) 2003, Gerald I. Evenden
007: * (C) 2000, Frank Warmerdam
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; either
012: * version 2.1 of the License, or (at your option) any later version.
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.util.InternationalString;
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.parameter.InvalidParameterValueException;
037: import org.opengis.referencing.operation.CylindricalProjection;
038: import org.opengis.referencing.operation.MathTransform;
039:
040: // Geotools dependencies
041: import org.geotools.measure.Angle;
042: import org.geotools.measure.Latitude;
043: import org.geotools.metadata.iso.citation.Citations;
044: import org.geotools.referencing.NamedIdentifier;
045: import org.geotools.resources.i18n.Vocabulary;
046: import org.geotools.resources.i18n.VocabularyKeys;
047: import org.geotools.resources.i18n.Errors;
048: import org.geotools.resources.i18n.ErrorKeys;
049:
050: /**
051: * Oblique Mercator Projection. A conformal, oblique, cylindrical projection
052: * with the cylinder touching the ellipsoid (or sphere) along a great circle
053: * path (the central line). The {@linkplain Mercator} and {@linkplain TransverseMercator
054: * Transverse Mercator} projections can be thought of as special cases of the
055: * oblique mercator, where the central line is along the equator or a meridian,
056: * respectively. The Oblique Mercator projection has been used in Switzerland,
057: * Hungary, Madagascar, Malaysia, Borneo and the panhandle of Alaska.
058: * <p>
059: *
060: * The Oblique Mercator projection uses a (<var>U</var>,<var>V</var>) coordinate system,
061: * with the <var>U</var> axis along the central line. During the forward projection,
062: * coordinates from the ellipsoid are projected conformally to a sphere of constant total
063: * curvature, called the "aposphere", before being projected onto the plane. The projection
064: * coordinates are further convented to a (<var>X</var>,<var>Y</var>) coordinate system
065: * by rotating the calculated (<var>u</var>,<var>v</var>) coordinates to give output
066: * (<var>x</var>,<var>y</var>) coordinates. The rotation value is usually the same as
067: * the projection azimuth (the angle, east of north, of the central line), but some cases
068: * allow a separate rotation parameter.
069: * <p>
070: *
071: * There are two forms of the oblique mercator, differing in the origin of
072: * their grid coordinates. The {@linkplain HotineObliqueMercator Hotine Oblique Mercator}
073: * (EPSG code 9812) has grid coordinates start at the intersection of the central line
074: * and the equator of the aposphere. The {@linkplain ObliqueMercator Oblique Mercator}
075: * (EPSG code 9815) is the same, except the grid coordinates begin at the central point
076: * (where the latitude of center and central line intersect). ESRI separates these two
077: * case by appending {@code "Natural_Origin"} (for the {@code "Hotine_Oblique_Mercator"})
078: * and {@code "Center"} (for the {@code "Oblique_Mercator"}) to the projection names.
079: * <p>
080: *
081: * Two different methods are used to specify the central line for the oblique mercator:
082: * 1) a central point and an azimuth, east of north, describing the central line and
083: * 2) two points on the central line. The EPSG does not use the two point method,
084: * while ESRI separates the two cases by putting {@code "Azimuth"} and {@code "Two_Point"}
085: * in their projection names. Both cases use the point where the {@code "latitude_of_center"}
086: * parameter crosses the central line as the projection's central point.
087: * The {@linkplain #centralMeridian central meridian} is not a projection parameter,
088: * and is instead calculated as the intersection between the central line and the
089: * equator of the aposphere.
090: * <p>
091: *
092: * For the azimuth method, the central latitude cannot be ±90.0 degrees
093: * and the central line cannot be at a maximum or minimum latitude at the central point.
094: * In the two point method, the latitude of the first and second points cannot be
095: * equal. Also, the latitude of the first point and central point cannot be
096: * ±90.0 degrees. Furthermore, the latitude of the first point cannot be 0.0 and
097: * the latitude of the second point cannot be -90.0 degrees. A change of
098: * 10<sup>-7</sup> radians can allow calculation at these special cases. Snyder's restriction
099: * of the central latitude being 0.0 has been removed, since the equations appear
100: * to work correctly in this case.
101: * <p>
102: *
103: * Azimuth values of 0.0 and ±90.0 degrees are allowed (and used in Hungary
104: * and Switzerland), though these cases would usually use a Mercator or
105: * Transverse Mercator projection instead. Azimuth values > 90 degrees cause
106: * errors in the equations.
107: * <p>
108: *
109: * The oblique mercator is also called the "Rectified Skew Orthomorphic" (RSO). It appears
110: * is that the only difference from the oblique mercator is that the RSO allows the rotation
111: * from the (<var>U</var>,<var>V</var>) to (<var>X</var>,<var>Y</var>) coordinate system to
112: * be different from the azimuth. This separate parameter is called
113: * {@code "rectified_grid_angle"} (or {@code "XY_Plane_Rotation"} by ESRI) and is also
114: * included in the EPSG's parameters for the Oblique Mercator and Hotine Oblique Mercator.
115: * The rotation parameter is optional in all the non-two point projections and will be
116: * set to the azimuth if not specified.
117: * <p>
118: *
119: * Projection cases and aliases implemented by the {@link ObliqueMercator} are:
120: * <ul>
121: * <li>Oblique_Mercator (EPSG code 9815) - grid coordinates begin at the central point, has "rectified_grid_angle" parameter.</li>
122: * <li>Hotine_Oblique_Mercator_Azimuth_Center (ESRI) - grid coordinates begin at the central point.</li>
123: * <li>Rectified_Skew_Orthomorphic_Center (ESRI) - grid coordinates begin at the central point, has "rectified_grid_angle" parameter.</li>
124: *
125: * <li>Hotine_Oblique_Mercator (EPSG code 9812) - grid coordinates begin at the interseciton of the central line and aposphere equator, has "rectified_grid_angle" parameter.</li>
126: * <li>Hotine_Oblique_Mercator_Azimuth_Natural_Origin (ESRI) - grid coordinates begin at the interseciton of the central line and aposphere equator.</li>
127: * <li>Rectified_Skew_Orthomorphic_Natural_Origin (ESRI) - grid coordinates begin at the interseciton of the central line and aposphere equator, has "rectified_grid_angle" parameter.</li>
128: *
129: * <li>Hotine_Oblique_Mercator_Two_Point_Center (ESRI) - grid coordinates begin at the central point.</li>
130: * <li>Hotine_Oblique_Mercator_Two_Point_Natural_Origin (ESRI) - grid coordinates begin at the interseciton of the central line and aposphere equator.</li>
131: * </ul>
132: *
133: * <strong>References:</strong>
134: * <ul>
135: * <li>{@code libproj4} is available at
136: * <A HREF="http://members.bellatlantic.net/~vze2hc4d/proj4/">libproj4 Miscellanea</A><br>
137: * Relevent files are: {@code PJ_omerc.c}, {@code pj_tsfn.c},
138: * {@code pj_fwd.c}, {@code pj_inv.c} and {@code lib_proj.h}</li>
139: * <li> John P. Snyder (Map Projections - A Working Manual,
140: * U.S. Geological Survey Professional Paper 1395, 1987)</li>
141: * <li> "Coordinate Conversions and Transformations including Formulas",
142: * EPSG Guidence Note Number 7 part 2, Version 24.</li>
143: * <li>Gerald Evenden, 2004, <a href="http://members.verizon.net/~vze2hc4d/proj4/omerc.pdf">
144: * Documentation of revised Oblique Mercator</a></li>
145: * </ul>
146: *
147: * @see <A HREF="http://mathworld.wolfram.com/MercatorProjection.html">Oblique Mercator projection on MathWorld</A>
148: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/hotine_oblique_mercator.html">"hotine_oblique_mercator" on RemoteSensing.org</A>
149: * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/oblique_mercator.html">"oblique_mercator" on RemoteSensing.org</A>
150: *
151: * @since 2.1
152: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/ObliqueMercator.java $
153: * @version $Id: ObliqueMercator.java 24581 2007-02-26 01:36:46Z desruisseaux $
154: * @author Rueben Schulz
155: */
156: public class ObliqueMercator extends MapProjection {
157: /**
158: * Maximum difference allowed when comparing real numbers.
159: */
160: private static final double EPSILON = 1E-6;
161:
162: /**
163: * Maximum difference allowed when comparing latitudes.
164: */
165: private static final double EPSILON_LATITUDE = 1E-10;
166:
167: //////
168: ////// Map projection parameters. The following are NOT used by the transformation
169: ////// methods, but are stored in order to restitute them in WKT formatting. They
170: ////// are made visible ('protected' access) for documentation purpose and because
171: ////// they are user-supplied parameters, not derived coefficients.
172: //////
173:
174: /**
175: * Latitude of the projection centre. This is similar to the
176: * {@link #latitudeOfOrigin}, but the latitude of origin is the
177: * Earth equator on aposphere for the oblique mercator.
178: */
179: protected final double latitudeOfCentre;
180:
181: /**
182: * Longitude of the projection centre. This is <strong>NOT</strong> equal
183: * to the {@link #centralMeridian}, which is the meridian where the
184: * central line intersects the Earth equator on aposphere.
185: * <p>
186: * This parameter applies to the "azimuth" case only.
187: * It is set to {@link Double#NaN NaN} for the "two points" case.
188: */
189: protected final double longitudeOfCentre;
190:
191: /**
192: * The azimuth of the central line passing throught the centre of the projection, in radians.
193: */
194: protected final double azimuth;
195:
196: /**
197: * The rectified bearing of the central line, in radians. This is equals to the
198: * {@linkplain #azimuth} if the {@link Provider#RECTIFIED_GRID_ANGLE "rectified_grid_angle"}
199: * parameter value is not set.
200: */
201: protected final double rectifiedGridAngle;
202:
203: // The next parameters still private for now because I'm not
204: // sure if they should appear in some 'TwoPoint' subclass...
205: /**
206: * The latitude of the 1st point used to specify the central line, in radians.
207: * <p>
208: * This parameter applies to the "two points" case only.
209: * It is set to {@link Double#NaN NaN} for the "azimuth" case.
210: */
211: private final double latitudeOf1stPoint;
212:
213: /**
214: * The longitude of the 1st point used to specify the central line, in radians.
215: * <p>
216: * This parameter applies to the "two points" case only.
217: * It is set to {@link Double#NaN NaN} for the "azimuth" case.
218: */
219: private final double longitudeOf1stPoint;
220:
221: /**
222: * The latitude of the 2nd point used to specify the central line, in radians.
223: * <p>
224: * This parameter applies to the "two points" case only.
225: * It is set to {@link Double#NaN NaN} for the "azimuth" case.
226: */
227: private final double latitudeOf2ndPoint;
228:
229: /**
230: * The longitude of the 2nd point used to specify the central line, in radians.
231: * <p>
232: * This parameter applies to the "two points" case only.
233: * It is set to {@link Double#NaN NaN} for the "azimuth" case.
234: */
235: private final double longitudeOf2ndPoint;
236:
237: //////
238: ////// Map projection coefficients computed from the above parameters.
239: ////// They are the fields used for coordinate transformations.
240: //////
241:
242: /**
243: * Constants used in the transformation.
244: */
245: private final double B, A, E;
246:
247: /**
248: * Convenience values equal to {@link #A} / {@link #B},
249: * {@link #A}×{@link #B}, and {@link #B} / {@link #A}.
250: */
251: private final double ArB, AB, BrA;
252:
253: /**
254: * <var>v</var> values when the input latitude is a pole.
255: */
256: private final double v_pole_n, v_pole_s;
257:
258: /**
259: * Sine and Cosine values for gamma0 (the angle between the meridian
260: * and central line at the intersection between the central line and
261: * the Earth equator on aposphere).
262: */
263: private final double singamma0, cosgamma0;
264:
265: /**
266: * Sine and Cosine values for the rotation between (U,V) and
267: * (X,Y) coordinate systems
268: */
269: private final double sinrot, cosrot;
270:
271: /**
272: * <var>u</var> value (in (U,V) coordinate system) of the central point. Used in
273: * the oblique mercator case. The <var>v</var> value of the central point is 0.0.
274: */
275: private final double u_c;
276:
277: /**
278: * {@code true} if using two points on the central line to specify the azimuth.
279: */
280: final boolean twoPoint;
281:
282: /**
283: * Constructs a new map projection from the supplied parameters.
284: *
285: * @param parameters The parameter values in standard units.
286: * @throws ParameterNotFoundException if a mandatory parameter is missing.
287: *
288: * @since 2.4
289: *
290: * @todo Current implementation assumes a "azimuth" case. We may try to detect the
291: * "two points" case in a future version if needed.
292: */
293: protected ObliqueMercator(final ParameterValueGroup parameters)
294: throws ParameterNotFoundException {
295: this (parameters, Provider.PARAMETERS.descriptors(), false,
296: false);
297: }
298:
299: /**
300: * Constructs a new map projection from the supplied parameters.
301: *
302: * @param parameters The parameter values in standard units.
303: * @param expected The expected parameter descriptors.
304: * @param twoPoint {@code true} for the "two points" case, or {@code false} for the
305: * "azimuth" case. The former is used by ESRI but not EPSG.
306: * @param hotine {@code true} only if invoked by the {@link HotineObliqueMercator}
307: * constructor.
308: * @throws ParameterNotFoundException if a mandatory parameter is missing.
309: */
310: ObliqueMercator(final ParameterValueGroup parameters,
311: final Collection expected, final boolean twoPoint,
312: final boolean hotine) throws ParameterNotFoundException {
313: // Fetch parameters
314: super (parameters, expected);
315: this .twoPoint = twoPoint;
316: /*
317: * Sets the 'centralMeridian' and 'latitudeOfOrigin' fields to NaN for safety
318: * (they are not 'ObliqueMercator' parameters, so the super-class initialized
319: * them to 0.0 by default). The 'centralMeridian' value will be computed later.
320: */
321: centralMeridian = Double.NaN;
322: latitudeOfOrigin = Double.NaN;
323: latitudeOfCentre = doubleValue(expected,
324: Provider.LATITUDE_OF_CENTRE, parameters);
325: /*
326: * Checks that 'latitudeOfCentre' is not +- 90 degrees.
327: * Not checking if 'latitudeOfCentere' is 0, since equations behave correctly.
328: */
329: ensureLatitudeInRange(Provider.LATITUDE_OF_CENTRE,
330: latitudeOfCentre, false);
331: /*
332: * Computes common constants now. In a previous version of 'ObliqueMercator', those
333: * constants were computed a little bit later. We compute them now in order to have
334: * a single 'if (twoPoint) ... else' statement, which help us to keep every fields
335: * final and catch some potential errors at compile-time (e.g. unitialized fields).
336: */
337: final double com = Math.sqrt(1.0 - excentricitySquared);
338: final double sinphi0 = Math.sin(latitudeOfCentre);
339: final double cosphi0 = Math.cos(latitudeOfCentre);
340: double temp = cosphi0 * cosphi0;
341: B = Math.sqrt(1.0 + excentricitySquared * (temp * temp)
342: / (1.0 - excentricitySquared));
343: final double con = 1.0 - excentricitySquared * sinphi0
344: * sinphi0;
345: A = B * com / con;
346: final double D = B * com / (cosphi0 * Math.sqrt(con));
347: double F = D * D - 1.0;
348: if (F < 0.0) {
349: F = 0.0;
350: } else {
351: F = Math.sqrt(F);
352: if (latitudeOfCentre < 0.0) { // Taking sign of 'latitudeOfCentre'
353: F = -F;
354: }
355: }
356: F = F += D;
357: E = F * Math.pow(tsfn(latitudeOfCentre, sinphi0), B);
358: /*
359: * Computes the constants that depend on the "twoPoint" vs "azimuth" case. In the
360: * two points case, we compute them from (LAT_OF_1ST_POINT, LONG_OF_1ST_POINT) and
361: * (LAT_OF_2ND_POINT, LONG_OF_2ND_POINT). For the "azimuth" case, we compute them
362: * from LONGITUDE_OF_CENTRE and AZIMUTH.
363: */
364: final double gamma0;
365: if (twoPoint) {
366: longitudeOfCentre = Double.NaN; // This is for the "azimuth" case only.
367: latitudeOf1stPoint = doubleValue(expected,
368: Provider_TwoPoint.LAT_OF_1ST_POINT, parameters);
369: // Checks that latOf1stPoint is not +-90 degrees
370: ensureLatitudeInRange(Provider_TwoPoint.LAT_OF_1ST_POINT,
371: latitudeOf1stPoint, false);
372: longitudeOf1stPoint = doubleValue(expected,
373: Provider_TwoPoint.LONG_OF_1ST_POINT, parameters);
374: ensureLongitudeInRange(Provider_TwoPoint.LONG_OF_1ST_POINT,
375: longitudeOf1stPoint, true);
376: latitudeOf2ndPoint = doubleValue(expected,
377: Provider_TwoPoint.LAT_OF_2ND_POINT, parameters);
378: ensureLatitudeInRange(Provider_TwoPoint.LAT_OF_2ND_POINT,
379: latitudeOf2ndPoint, true);
380: double longitudeOf2ndPoint; // Will be assigned to the field later.
381: longitudeOf2ndPoint = doubleValue(expected,
382: Provider_TwoPoint.LONG_OF_2ND_POINT, parameters);
383: ensureLongitudeInRange(Provider_TwoPoint.LONG_OF_2ND_POINT,
384: longitudeOf2ndPoint, true);
385: /*
386: * Ensures that (phi1 != phi2), (phi1 != 0°) and (phi2 != -90°),
387: * as specified in class javadoc.
388: */
389: ParameterDescriptor desc = null;
390: Object value = null;
391: if (Math.abs(latitudeOf1stPoint - latitudeOf2ndPoint) < EPSILON_LATITUDE) {
392: desc = Provider_TwoPoint.LAT_OF_1ST_POINT;
393: value = Provider_TwoPoint.LAT_OF_2ND_POINT.getName()
394: .getCode();
395: // Exception will be thrown below.
396: }
397: if (Math.abs(latitudeOf1stPoint) < EPSILON_LATITUDE) {
398: desc = Provider_TwoPoint.LAT_OF_1ST_POINT;
399: value = new Latitude(latitudeOf1stPoint);
400: // Exception will be thrown below.
401: }
402: if (Math.abs(latitudeOf2ndPoint + Math.PI / 2.0) < EPSILON_LATITUDE) {
403: desc = Provider_TwoPoint.LAT_OF_2ND_POINT;
404: value = new Latitude(latitudeOf2ndPoint);
405: // Exception will be thrown below.
406: }
407: if (desc != null) {
408: final String name = desc.getName().getCode();
409: throw new InvalidParameterValueException(Errors.format(
410: ErrorKeys.ILLEGAL_ARGUMENT_$2, name, value),
411: name, value);
412: }
413: /*
414: * The coefficients for the "two points" case.
415: */
416: final double H = Math.pow(tsfn(latitudeOf1stPoint, Math
417: .sin(latitudeOf1stPoint)), B);
418: final double L = Math.pow(tsfn(latitudeOf2ndPoint, Math
419: .sin(latitudeOf2ndPoint)), B);
420: final double Fp = E / H;
421: final double P = (L - H) / (L + H);
422: double J = E * E;
423: J = (J - L * H) / (J + L * H);
424: double diff = longitudeOf1stPoint - longitudeOf2ndPoint;
425: if (diff < -Math.PI) {
426: longitudeOf2ndPoint -= 2.0 * Math.PI;
427: } else if (diff > Math.PI) {
428: longitudeOf2ndPoint += 2.0 * Math.PI;
429: }
430: this .longitudeOf2ndPoint = longitudeOf2ndPoint;
431: centralMeridian = rollLongitude(0.5
432: * (longitudeOf1stPoint + longitudeOf2ndPoint)
433: - Math
434: .atan(J
435: * Math
436: .tan(0.5
437: * B
438: * (longitudeOf1stPoint - longitudeOf2ndPoint))
439: / P) / B);
440: gamma0 = Math.atan(2.0
441: * Math.sin(B
442: * rollLongitude(longitudeOf1stPoint
443: - centralMeridian))
444: / (Fp - 1.0 / Fp));
445: azimuth = Math.asin(D * Math.sin(gamma0));
446: rectifiedGridAngle = azimuth;
447: } else {
448: /*
449: * Computes coefficients for the "azimuth" case. Set the 1st and 2nd points
450: * to (NaN,NaN) since they are specific to the "two points" case. They are
451: * involved in WKT formatting only, not in transformation calculation.
452: */
453: latitudeOf1stPoint = Double.NaN;
454: longitudeOf1stPoint = Double.NaN;
455: latitudeOf2ndPoint = Double.NaN;
456: longitudeOf2ndPoint = Double.NaN;
457:
458: longitudeOfCentre = doubleValue(expected,
459: Provider.LONGITUDE_OF_CENTRE, parameters);
460: ensureLongitudeInRange(Provider.LONGITUDE_OF_CENTRE,
461: longitudeOfCentre, true);
462:
463: azimuth = doubleValue(expected, Provider.AZIMUTH,
464: parameters);
465: // Already checked for +-360 deg. above.
466: if ((azimuth > -1.5 * Math.PI && azimuth < -0.5 * Math.PI)
467: || (azimuth > 0.5 * Math.PI && azimuth < 1.5 * Math.PI)) {
468: final String name = Provider.AZIMUTH.getName()
469: .getCode();
470: final Angle value = new Angle(Math.toDegrees(azimuth));
471: throw new InvalidParameterValueException(Errors.format(
472: ErrorKeys.ILLEGAL_ARGUMENT_$2, name, value),
473: name, value);
474: }
475:
476: temp = doubleValue(expected, Provider.RECTIFIED_GRID_ANGLE,
477: parameters);
478: if (Double.isNaN(temp)) {
479: temp = azimuth;
480: }
481: rectifiedGridAngle = temp;
482: gamma0 = Math.asin(Math.sin(azimuth) / D);
483: // Check for asin(+-1.00000001)
484: temp = 0.5 * (F - 1.0 / F) * Math.tan(gamma0);
485: if (Math.abs(temp) > 1.0) {
486: if (Math.abs(Math.abs(temp) - 1.0) > EPSILON) {
487: throw new IllegalArgumentException(Errors
488: .format(ErrorKeys.TOLERANCE_ERROR));
489: }
490: temp = (temp > 0) ? 1.0 : -1.0;
491: }
492: centralMeridian = longitudeOfCentre - Math.asin(temp) / B;
493: }
494: /*
495: * More coefficients common to all kind of oblique mercator.
496: */
497: singamma0 = Math.sin(gamma0);
498: cosgamma0 = Math.cos(gamma0);
499: sinrot = Math.sin(rectifiedGridAngle);
500: cosrot = Math.cos(rectifiedGridAngle);
501: ArB = A / B;
502: AB = A * B;
503: BrA = B / A;
504: v_pole_n = ArB
505: * Math.log(Math.tan(0.5 * (Math.PI / 2.0 - gamma0)));
506: v_pole_s = ArB
507: * Math.log(Math.tan(0.5 * (Math.PI / 2.0 + gamma0)));
508:
509: if (hotine) {
510: u_c = 0.0;
511: } else {
512: if (Math.abs(Math.abs(azimuth) - Math.PI / 2.0) < EPSILON_LATITUDE) {
513: // LongitudeOfCentre = NaN in twoPoint, but azimuth cannot be 90 here (lat1 != lat2)
514: u_c = A * (longitudeOfCentre - centralMeridian);
515: } else {
516: double u_c = Math.abs(ArB
517: * Math.atan2(Math.sqrt(D * D - 1.0), Math
518: .cos(azimuth)));
519: if (latitudeOfCentre < 0.0) {
520: u_c = -u_c;
521: }
522: this .u_c = u_c;
523: }
524: }
525: }
526:
527: /**
528: * {@inheritDoc}
529: */
530: public ParameterDescriptorGroup getParameterDescriptors() {
531: return (twoPoint) ? Provider_TwoPoint.PARAMETERS
532: : Provider.PARAMETERS;
533: }
534:
535: /**
536: * {@inheritDoc}
537: */
538: public ParameterValueGroup getParameterValues() {
539: final ParameterValueGroup values = super .getParameterValues();
540: final Collection expected = getParameterDescriptors()
541: .descriptors();
542: // Note: we don't need a "if (twoPoint) ... else" statement since
543: // the "set" method will actually set the value only if applicable.
544: set(expected, Provider.LATITUDE_OF_CENTRE, values,
545: latitudeOfCentre);
546: set(expected, Provider.LONGITUDE_OF_CENTRE, values,
547: longitudeOfCentre);
548: set(expected, Provider.AZIMUTH, values, azimuth);
549: set(expected, Provider.RECTIFIED_GRID_ANGLE, values,
550: rectifiedGridAngle);
551:
552: set(expected, Provider_TwoPoint.LAT_OF_1ST_POINT, values,
553: latitudeOf1stPoint);
554: set(expected, Provider_TwoPoint.LONG_OF_1ST_POINT, values,
555: longitudeOf1stPoint);
556: set(expected, Provider_TwoPoint.LAT_OF_2ND_POINT, values,
557: latitudeOf2ndPoint);
558: set(expected, Provider_TwoPoint.LONG_OF_2ND_POINT, values,
559: longitudeOf2ndPoint);
560: return values;
561: }
562:
563: /**
564: * {@inheritDoc}
565: */
566: protected Point2D transformNormalized(double x, double y,
567: Point2D ptDst) throws ProjectionException {
568: double u, v;
569: if (Math.abs(Math.abs(y) - Math.PI / 2.0) > EPSILON) {
570: double Q = E / Math.pow(tsfn(y, Math.sin(y)), B);
571: double temp = 1.0 / Q;
572: double S = 0.5 * (Q - temp);
573: double V = Math.sin(B * x);
574: double U = (S * singamma0 - V * cosgamma0)
575: / (0.5 * (Q + temp));
576: if (Math.abs(Math.abs(U) - 1.0) < EPSILON) {
577: throw new ProjectionException(Errors.format(
578: ErrorKeys.INFINITE_VALUE_$1, "v"));
579: }
580: v = 0.5 * ArB * Math.log((1.0 - U) / (1.0 + U));
581: temp = Math.cos(B * x);
582: if (Math.abs(temp) < EPSILON_LATITUDE) {
583: u = AB * x;
584: } else {
585: u = ArB
586: * Math.atan2((S * cosgamma0 + V * singamma0),
587: temp);
588: }
589: } else {
590: v = y > 0 ? v_pole_n : v_pole_s;
591: u = ArB * y;
592: }
593:
594: u -= u_c;
595: x = v * cosrot + u * sinrot;
596: y = u * cosrot - v * sinrot;
597:
598: if (ptDst != null) {
599: ptDst.setLocation(x, y);
600: return ptDst;
601: }
602: return new Point2D.Double(x, y);
603: }
604:
605: /**
606: * {@inheritDoc}
607: */
608: protected Point2D inverseTransformNormalized(double x, double y,
609: Point2D ptDst) throws ProjectionException {
610: double v = x * cosrot - y * sinrot;
611: double u = y * cosrot + x * sinrot + u_c;
612:
613: double Qp = Math.exp(-BrA * v);
614: double temp = 1.0 / Qp;
615: double Sp = 0.5 * (Qp - temp);
616: double Vp = Math.sin(BrA * u);
617: double Up = (Vp * cosgamma0 + Sp * singamma0)
618: / (0.5 * (Qp + temp));
619: if (Math.abs(Math.abs(Up) - 1.0) < EPSILON) {
620: x = 0.0;
621: y = Up < 0.0 ? -Math.PI / 2.0 : Math.PI / 2.0;
622: } else {
623: y = Math.pow(E / Math.sqrt((1. + Up) / (1. - Up)), 1.0 / B); //calculate t
624: y = cphi2(y);
625: x = -Math.atan2((Sp * cosgamma0 - Vp * singamma0), Math
626: .cos(BrA * u))
627: / B;
628: }
629:
630: if (ptDst != null) {
631: ptDst.setLocation(x, y);
632: return ptDst;
633: }
634: return new Point2D.Double(x, y);
635: }
636:
637: /**
638: * Maximal error (in metres) tolerated for assertion, if enabled.
639: *
640: * @param longitude The longitude in decimal degrees.
641: * @param latitude The latitude in decimal degrees.
642: * @return The tolerance level for assertions, in meters.
643: */
644: protected double getToleranceForAssertions(final double longitude,
645: final double latitude) {
646: if (Math.abs(longitude - centralMeridian) / 2
647: + Math.abs(latitude - latitudeOfCentre) > 10) {
648: // When far from the valid area, use a larger tolerance.
649: return 1;
650: }
651: return super .getToleranceForAssertions(longitude, latitude);
652: }
653:
654: /**
655: * Returns a hash value for this projection.
656: */
657: public int hashCode() {
658: long code = Double.doubleToLongBits(latitudeOfCentre);
659: code = code * 37 + Double.doubleToLongBits(longitudeOfCentre);
660: code = code * 37 + Double.doubleToLongBits(azimuth);
661: code = code * 37 + Double.doubleToLongBits(rectifiedGridAngle);
662: code = code * 37 + Double.doubleToLongBits(latitudeOf1stPoint);
663: code = code * 37 + Double.doubleToLongBits(latitudeOf2ndPoint);
664: return ((int) code ^ (int) (code >>> 32)) + 37
665: * super .hashCode();
666: }
667:
668: /**
669: * Compares the specified object with this map projection for equality.
670: */
671: public boolean equals(final Object object) {
672: if (object == this ) {
673: // Slight optimization
674: return true;
675: }
676: if (super .equals(object)) {
677: final ObliqueMercator that = (ObliqueMercator) object;
678: return this .twoPoint == that.twoPoint
679: && equals(this .latitudeOfCentre,
680: that.latitudeOfCentre)
681: && equals(this .longitudeOfCentre,
682: that.longitudeOfCentre)
683: && equals(this .azimuth, that.azimuth)
684: && equals(this .rectifiedGridAngle,
685: that.rectifiedGridAngle)
686: && equals(this .latitudeOf1stPoint,
687: that.latitudeOf1stPoint)
688: && equals(this .longitudeOf1stPoint,
689: that.longitudeOf1stPoint)
690: && equals(this .latitudeOf2ndPoint,
691: that.latitudeOf2ndPoint)
692: && equals(this .longitudeOf2ndPoint,
693: that.longitudeOf2ndPoint)
694: && equals(this .u_c, that.u_c);
695: // Note: "u_c" is a derived parameter, so in theory we don't need to compare it.
696: // However we still compare it as a safety, because it takes a different
697: // value in the "hotine" case.
698: }
699: return false;
700: }
701:
702: //////////////////////////////////////////////////////////////////////////////////////////
703: //////////////////////////////////////////////////////////////////////////////////////////
704: //////// ////////
705: //////// PROVIDERS ////////
706: //////// ////////
707: //////////////////////////////////////////////////////////////////////////////////////////
708: //////////////////////////////////////////////////////////////////////////////////////////
709:
710: /**
711: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
712: * provider} for an {@linkplain ObliqueMercator Oblique Mercator} projection (EPSG code 9815).
713: *
714: * @since 2.1
715: * @version $Id: ObliqueMercator.java 24581 2007-02-26 01:36:46Z desruisseaux $
716: * @author Rueben Schulz
717: *
718: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
719: */
720: public static class Provider extends AbstractProvider {
721: /**
722: * The operation parameter descriptor for the {@link #latitudeOfCentre latitudeOfCentre}
723: * parameter value. Valid values range is from -90 to 90°. Default value is 0.
724: */
725: public static final ParameterDescriptor LATITUDE_OF_CENTRE = createDescriptor(
726: new NamedIdentifier[] {
727: new NamedIdentifier(Citations.OGC,
728: "latitude_of_center"),
729: new NamedIdentifier(Citations.EPSG,
730: "Latitude of projection centre"),
731: new NamedIdentifier(Citations.ESRI,
732: "Latitude_Of_Center"),
733: new NamedIdentifier(Citations.GEOTIFF,
734: "CenterLat") }, 0, -90, 90,
735: NonSI.DEGREE_ANGLE);
736:
737: /**
738: * @deprecated Renamed {@link #LATITUDE_OF_CENTRE}.
739: */
740: public static final ParameterDescriptor LAT_OF_CENTRE = LATITUDE_OF_CENTRE;
741:
742: /**
743: * The operation parameter descriptor for the {@link #longitudeOfCentre longitudeOfCentre}
744: * parameter value. Valid values range is from -180 to 180°. Default value is 0.
745: */
746: public static final ParameterDescriptor LONGITUDE_OF_CENTRE = createDescriptor(
747: new NamedIdentifier[] {
748: new NamedIdentifier(Citations.OGC,
749: "longitude_of_center"),
750: new NamedIdentifier(Citations.EPSG,
751: "Longitude of projection centre"),
752: new NamedIdentifier(Citations.ESRI,
753: "Longitude_Of_Center"),
754: new NamedIdentifier(Citations.GEOTIFF,
755: "CenterLong") }, 0, -180, 180,
756: NonSI.DEGREE_ANGLE);
757:
758: /**
759: * @deprecated Renamed {@link #LONGITUDE_OF_CENTRE}.
760: */
761: public static final ParameterDescriptor LONG_OF_CENTRE = LONGITUDE_OF_CENTRE;
762:
763: /**
764: * The operation parameter descriptor for the {@link #azimuth azimuth}
765: * parameter value. Valid values range is from -360 to -270, -90 to 90,
766: * and 270 to 360 degrees. Default value is 0.
767: */
768: public static final ParameterDescriptor AZIMUTH = createDescriptor(
769: new NamedIdentifier[] {
770: new NamedIdentifier(Citations.OGC, "azimuth"),
771: new NamedIdentifier(Citations.ESRI, "Azimuth"),
772: new NamedIdentifier(Citations.EPSG,
773: "Azimuth of initial line"),
774: new NamedIdentifier(Citations.GEOTIFF,
775: "AzimuthAngle") }, 0, -360, 360,
776: NonSI.DEGREE_ANGLE);
777:
778: /**
779: * The operation parameter descriptor for the {@link #rectifiedGridAngle
780: * rectifiedGridAngle} parameter value. It is an optional parameter with
781: * valid values ranging from -360 to 360°. Default value is {@link #azimuth azimuth}.
782: */
783: public static final ParameterDescriptor RECTIFIED_GRID_ANGLE = createOptionalDescriptor(
784: new NamedIdentifier[] {
785: new NamedIdentifier(Citations.OGC,
786: "rectified_grid_angle"),
787: new NamedIdentifier(Citations.EPSG,
788: "Angle from Rectified to Skew Grid"),
789: new NamedIdentifier(Citations.ESRI,
790: "XY_Plane_Rotation"),
791: new NamedIdentifier(Citations.GEOTIFF,
792: "RectifiedGridAngle") }, -360, 360,
793: NonSI.DEGREE_ANGLE);
794:
795: /**
796: * The localized name for "Oblique Mercator".
797: */
798: static final InternationalString NAME = Vocabulary
799: .formatInternational(VocabularyKeys.OBLIQUE_MERCATOR_PROJECTION);
800:
801: /**
802: * The parameters group.
803: */
804: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
805: new NamedIdentifier[] {
806: new NamedIdentifier(Citations.OGC,
807: "Oblique_Mercator"),
808: new NamedIdentifier(Citations.EPSG,
809: "Oblique Mercator"),
810: new NamedIdentifier(Citations.EPSG, "9815"),
811: new NamedIdentifier(Citations.GEOTIFF,
812: "CT_ObliqueMercator"),
813: new NamedIdentifier(Citations.ESRI,
814: "Hotine_Oblique_Mercator_Azimuth_Center"),
815: new NamedIdentifier(Citations.ESRI,
816: "Rectified_Skew_Orthomorphic_Center"),
817: new NamedIdentifier(Citations.GEOTOOLS, NAME) },
818: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
819: LONGITUDE_OF_CENTRE, LATITUDE_OF_CENTRE,
820: AZIMUTH, RECTIFIED_GRID_ANGLE, SCALE_FACTOR,
821: FALSE_EASTING, FALSE_NORTHING });
822:
823: /**
824: * Constructs a new provider.
825: */
826: public Provider() {
827: super (PARAMETERS);
828: }
829:
830: /**
831: * Constructs a new provider.
832: */
833: protected Provider(final ParameterDescriptorGroup params) {
834: super (params);
835: }
836:
837: /**
838: * Returns the operation type for this map projection.
839: */
840: public Class getOperationType() {
841: return CylindricalProjection.class;
842: }
843:
844: /**
845: * Creates a transform from the specified group of parameter values.
846: *
847: * @param parameters The group of parameter values.
848: * @return The created math transform.
849: * @throws ParameterNotFoundException if a required parameter was not found.
850: */
851: protected MathTransform createMathTransform(
852: final ParameterValueGroup parameters)
853: throws ParameterNotFoundException {
854: final Collection descriptors = PARAMETERS.descriptors();
855: return new ObliqueMercator(parameters, descriptors, false,
856: false);
857: }
858: }
859:
860: /**
861: * The {@linkplain org.geotools.referencing.operation.MathTransformProvider math transform
862: * provider} for a {@linkplain ObliqueMercator Oblique Mercator} projection, specified
863: * with two points on the central line (instead of a central point and azimuth).
864: *
865: * @since 2.1
866: * @version $Id: ObliqueMercator.java 24581 2007-02-26 01:36:46Z desruisseaux $
867: * @author Rueben Schulz
868: *
869: * @see org.geotools.referencing.operation.DefaultMathTransformFactory
870: */
871: public static class Provider_TwoPoint extends Provider {
872: /**
873: * The operation parameter descriptor for the {@code latitudeOf1stPoint}
874: * parameter value. Valid values range is from -90 to 90°. Default value is 0.
875: */
876: public static final ParameterDescriptor LAT_OF_1ST_POINT = createDescriptor(
877: new NamedIdentifier[] { new NamedIdentifier(
878: Citations.ESRI, "Latitude_Of_1st_Point") }, 0,
879: -90, 90, NonSI.DEGREE_ANGLE);
880:
881: /**
882: * The operation parameter descriptor for the {@code longitudeOf1stPoint}
883: * parameter value. Valid values range is from -180 to 180°. Default value is 0.
884: */
885: public static final ParameterDescriptor LONG_OF_1ST_POINT = createDescriptor(
886: new NamedIdentifier[] { new NamedIdentifier(
887: Citations.ESRI, "Longitude_Of_1st_Point") }, 0,
888: -180, 180, NonSI.DEGREE_ANGLE);
889:
890: /**
891: * The operation parameter descriptor for the {@code latitudeOf2ndPoint}
892: * parameter value. Valid values range is from -90 to 90°. Default value is 0.
893: */
894: public static final ParameterDescriptor LAT_OF_2ND_POINT = createDescriptor(
895: new NamedIdentifier[] { new NamedIdentifier(
896: Citations.ESRI, "Latitude_Of_2nd_Point") }, 0,
897: -90, 90, NonSI.DEGREE_ANGLE);
898:
899: /**
900: * The operation parameter descriptor for the {@code longitudeOf2ndPoint}
901: * parameter value. Valid values range is from -180 to 180°. Default value is 0.
902: */
903: public static final ParameterDescriptor LONG_OF_2ND_POINT = createDescriptor(
904: new NamedIdentifier[] { new NamedIdentifier(
905: Citations.ESRI, "Longitude_Of_2nd_Point") }, 0,
906: -180, 180, NonSI.DEGREE_ANGLE);
907:
908: /**
909: * The parameters group.
910: */
911: static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(
912: new NamedIdentifier[] {
913: new NamedIdentifier(Citations.ESRI,
914: "Hotine_Oblique_Mercator_Two_Point_Center"),
915: new NamedIdentifier(Citations.GEOTOOLS, NAME) },
916: new ParameterDescriptor[] { SEMI_MAJOR, SEMI_MINOR,
917: LAT_OF_1ST_POINT, LONG_OF_1ST_POINT,
918: LAT_OF_2ND_POINT, LONG_OF_2ND_POINT,
919: LATITUDE_OF_CENTRE, SCALE_FACTOR,
920: FALSE_EASTING, FALSE_NORTHING });
921:
922: /**
923: * Constructs a new provider.
924: */
925: public Provider_TwoPoint() {
926: super (PARAMETERS);
927: }
928:
929: /**
930: * Constructs a new provider.
931: */
932: protected Provider_TwoPoint(
933: final ParameterDescriptorGroup params) {
934: super (params);
935: }
936:
937: /**
938: * Creates a transform from the specified group of parameter values.
939: *
940: * @param parameters The group of parameter values.
941: * @return The created math transform.
942: * @throws ParameterNotFoundException if a required parameter was not found.
943: */
944: protected MathTransform createMathTransform(
945: final ParameterValueGroup parameters)
946: throws ParameterNotFoundException {
947: final Collection descriptors = PARAMETERS.descriptors();
948: return new ObliqueMercator(parameters, descriptors, true,
949: false);
950: }
951: }
952: }
|