0001: /*
0002: * GeoTools - OpenSource mapping toolkit
0003: * http://geotools.org
0004: *
0005: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
0006: * (C) 2001, Institut de Recherche pour le Développement
0007: * (C) 2000, Frank Warmerdam
0008: * (C) 1999, Fisheries and Oceans Canada
0009: *
0010: * This library is free software; you can redistribute it and/or
0011: * modify it under the terms of the GNU Lesser General Public
0012: * License as published by the Free Software Foundation; either
0013: * version 2.1 of the License, or (at your option) any later version.
0014: *
0015: * This library is distributed in the hope that it will be useful,
0016: * but WITHOUT ANY WARRANTY; without even the implied warranty of
0017: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
0018: * Lesser General Public License for more details.
0019: *
0020: * This package contains formulas from the PROJ package of USGS.
0021: * USGS's work is fully acknowledged here. This derived work has
0022: * been relicensed under LGPL with Frank Warmerdam's permission.
0023: */
0024: package org.geotools.referencing.operation.projection;
0025:
0026: // J2SE dependencies and extensions
0027: import java.awt.geom.Point2D;
0028: import java.io.Serializable;
0029: import java.util.Collection;
0030: import javax.units.NonSI;
0031: import javax.units.SI;
0032: import javax.units.Unit;
0033:
0034: // OpenGIS dependencies
0035: import org.opengis.parameter.InvalidParameterValueException;
0036: import org.opengis.parameter.ParameterDescriptor;
0037: import org.opengis.parameter.ParameterDescriptorGroup;
0038: import org.opengis.parameter.ParameterNotFoundException;
0039: import org.opengis.parameter.ParameterValue;
0040: import org.opengis.parameter.ParameterValueGroup;
0041: import org.opengis.referencing.operation.MathTransform;
0042: import org.opengis.referencing.operation.MathTransform2D;
0043: import org.opengis.referencing.operation.Projection;
0044: import org.opengis.referencing.operation.TransformException;
0045:
0046: // Geotools dependencies
0047: import org.geotools.measure.Latitude;
0048: import org.geotools.measure.Longitude;
0049: import org.geotools.metadata.iso.citation.Citations;
0050: import org.geotools.referencing.NamedIdentifier;
0051: import org.geotools.referencing.operation.MathTransformProvider;
0052: import org.geotools.referencing.operation.transform.AbstractMathTransform;
0053: import org.geotools.resources.XMath;
0054: import org.geotools.resources.i18n.Errors;
0055: import org.geotools.resources.i18n.ErrorKeys;
0056:
0057: /**
0058: * Base class for transformation services between ellipsoidal and cartographic projections.
0059: * This base class provides the basic feature needed for all methods (no need to overrides
0060: * methods). Subclasses must "only" implements the following methods:
0061: * <ul>
0062: * <li>{@link #getParameterValues}</li>
0063: * <li>{@link #transformNormalized}</li>
0064: * <li>{@link #inverseTransformNormalized}</li>
0065: * </ul>
0066: * <p>
0067: * <strong>NOTE:</strong>Serialization of this class is appropriate for short-term storage
0068: * or RMI use, but will probably not be compatible with future version. For long term storage,
0069: * WKT (Well Know Text) or XML (not yet implemented) are more appropriate.
0070: *
0071: * @since 2.0
0072: * @version $Id: MapProjection.java 26518 2007-08-10 17:05:01Z desruisseaux $
0073: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/MapProjection.java $
0074: * @author André Gosselin
0075: * @author Martin Desruisseaux
0076: * @author Rueben Schulz
0077: *
0078: * @see <A HREF="http://mathworld.wolfram.com/MapProjection.html">Map projections on MathWorld</A>
0079: * @see <A HREF="http://atlas.gc.ca/site/english/learningresources/carto_corner/map_projections.html">Map projections on the atlas of Canada</A>
0080: * @tutorial http://www.geotools.org/display/GEOTOOLS/How+to+add+new+projections
0081: */
0082: public abstract class MapProjection extends AbstractMathTransform
0083: implements MathTransform2D, Serializable {
0084: /**
0085: * Maximum difference allowed when comparing real numbers.
0086: */
0087: private static final double EPSILON = 1E-6;
0088:
0089: /**
0090: * Maximum difference allowed when comparing longitudes or latitudes in degrees.
0091: * This tolerance do not apply to angle in radians. A tolerance of 1E-4 is about
0092: * 10 kilometers.
0093: */
0094: private static final double ANGLE_TOLERANCE = 1E-4;
0095:
0096: /**
0097: * Difference allowed in iterative computations.
0098: */
0099: private static final double ITERATION_TOLERANCE = 1E-10;
0100:
0101: /**
0102: * Maximum number of iterations for iterative computations.
0103: */
0104: private static final int MAXIMUM_ITERATIONS = 15;
0105:
0106: /**
0107: * Ellipsoid excentricity, equals to <code>sqrt({@link #excentricitySquared})</code>.
0108: * Value 0 means that the ellipsoid is spherical.
0109: *
0110: * @see #excentricitySquared
0111: * @see #isSpherical
0112: */
0113: protected final double excentricity;
0114:
0115: /**
0116: * The square of excentricity: e² = (a²-b²)/a² where
0117: * <var>e</var> is the {@linkplain #excentricity excentricity},
0118: * <var>a</var> is the {@linkplain #semiMajor semi major} axis length and
0119: * <var>b</var> is the {@linkplain #semiMinor semi minor} axis length.
0120: *
0121: * @see #excentricity
0122: * @see #semiMajor
0123: * @see #semiMinor
0124: * @see #isSpherical
0125: */
0126: protected final double excentricitySquared;
0127:
0128: /**
0129: * {@code true} if this projection is spherical. Spherical model has identical
0130: * {@linkplain #semiMajor semi major} and {@linkplain #semiMinor semi minor} axis
0131: * length, and an {@linkplain #excentricity excentricity} zero.
0132: *
0133: * @see #excentricity
0134: * @see #semiMajor
0135: * @see #semiMinor
0136: */
0137: protected final boolean isSpherical;
0138:
0139: /**
0140: * Length of semi-major axis, in metres. This is named '<var>a</var>' or '<var>R</var>'
0141: * (Radius in spherical cases) in Snyder.
0142: *
0143: * @see #excentricity
0144: * @see #semiMinor
0145: */
0146: protected final double semiMajor;
0147:
0148: /**
0149: * Length of semi-minor axis, in metres. This is named '<var>b</var>' in Snyder.
0150: *
0151: * @see #excentricity
0152: * @see #semiMajor
0153: */
0154: protected final double semiMinor;
0155:
0156: /**
0157: * Central longitude in <u>radians</u>. Default value is 0, the Greenwich meridian.
0158: * This is called '<var>lambda0</var>' in Snyder.
0159: *
0160: * <strong>Consider this field as final</strong>. It is not final only
0161: * because some classes need to modify it at construction time.
0162: */
0163: protected double centralMeridian;
0164:
0165: /**
0166: * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
0167: * This is called '<var>phi0</var>' in Snyder.
0168: *
0169: * <strong>Consider this field as final</strong>. It is not final only
0170: * because some classes need to modify it at construction time.
0171: */
0172: protected double latitudeOfOrigin;
0173:
0174: /**
0175: * The scale factor. Default value is 1. Named '<var>k</var>' in Snyder.
0176: *
0177: * <strong>Consider this field as final</strong>. It is not final only
0178: * because some classes need to modify it at construction time.
0179: */
0180: protected double scaleFactor;
0181:
0182: /**
0183: * False easting, in metres. Default value is 0.
0184: */
0185: protected final double falseEasting;
0186:
0187: /**
0188: * False northing, in metres. Default value is 0.
0189: */
0190: protected final double falseNorthing;
0191:
0192: /**
0193: * Global scale factor. Default value {@code globalScale} is equal
0194: * to {@link #semiMajor}×{@link #scaleFactor}.
0195: *
0196: * <strong>Consider this field as final</strong>. It is not final only
0197: * because some classes need to modify it at construction time.
0198: */
0199: protected double globalScale;
0200:
0201: /**
0202: * The inverse of this map projection. Will be created only when needed.
0203: */
0204: private transient MathTransform inverse;
0205:
0206: /**
0207: * Constructs a new map projection from the suplied parameters.
0208: *
0209: * @param values The parameter values in standard units.
0210: * The following parameter are recognized:
0211: * <ul>
0212: * <li>"semi_major" (mandatory: no default)</li>
0213: * <li>"semi_minor" (mandatory: no default)</li>
0214: * <li>"central_meridian" (default to 0°)</li>
0215: * <li>"latitude_of_origin" (default to 0°)</li>
0216: * <li>"scale_factor" (default to 1 )</li>
0217: * <li>"false_easting" (default to 0 )</li>
0218: * <li>"false_northing" (default to 0 )</li>
0219: * </ul>
0220: * @throws ParameterNotFoundException if a mandatory parameter is missing.
0221: */
0222: protected MapProjection(final ParameterValueGroup values)
0223: throws ParameterNotFoundException {
0224: this (values, null);
0225: }
0226:
0227: /**
0228: * Constructor invoked by sub-classes when we can't rely on
0229: * {@link #getParameterDescriptors} before the construction
0230: * is completed. This is the case when the later method depends on
0231: * the value of some class's attribute, which has not yet been set.
0232: * An example is {@link Mercator#getParameterDescriptors}.
0233: *
0234: * This method is not public because it is not a very elegant hack, and
0235: * a work around exists. For example Mercator_1SP and Mercator_2SP could
0236: * be implemented by two separated classes, in which case {@link #getParameterDescriptors}
0237: * returns a constant and can be safely invoked in a constructor. We do
0238: * not always use this cleaner way in the projection package because it
0239: * is going to contains a lot of.. well... projections, and we will try
0240: * to reduce the amount of class declarations.
0241: */
0242: MapProjection(final ParameterValueGroup values, Collection expected)
0243: throws ParameterNotFoundException {
0244: if (expected == null) {
0245: expected = getParameterDescriptors().descriptors();
0246: }
0247: semiMajor = doubleValue(expected, AbstractProvider.SEMI_MAJOR,
0248: values);
0249: semiMinor = doubleValue(expected, AbstractProvider.SEMI_MINOR,
0250: values);
0251: centralMeridian = doubleValue(expected,
0252: AbstractProvider.CENTRAL_MERIDIAN, values);
0253: latitudeOfOrigin = doubleValue(expected,
0254: AbstractProvider.LATITUDE_OF_ORIGIN, values);
0255: scaleFactor = doubleValue(expected,
0256: AbstractProvider.SCALE_FACTOR, values);
0257: falseEasting = doubleValue(expected,
0258: AbstractProvider.FALSE_EASTING, values);
0259: falseNorthing = doubleValue(expected,
0260: AbstractProvider.FALSE_NORTHING, values);
0261: isSpherical = (semiMajor == semiMinor);
0262: excentricitySquared = 1.0 - (semiMinor * semiMinor)
0263: / (semiMajor * semiMajor);
0264: excentricity = Math.sqrt(excentricitySquared);
0265: globalScale = scaleFactor * semiMajor;
0266: ensureLongitudeInRange(AbstractProvider.CENTRAL_MERIDIAN,
0267: centralMeridian, true);
0268: ensureLatitudeInRange(AbstractProvider.LATITUDE_OF_ORIGIN,
0269: latitudeOfOrigin, true);
0270: }
0271:
0272: /**
0273: * Returns {@code true} if the specified parameter can apply to this map projection.
0274: * The set of expected parameters must be supplied. The default implementation just
0275: * invokes {@code expected.contains(param)}. Some subclasses will override this method
0276: * in order to handle {@link ModifiedParameterDescriptor} in a special way.
0277: *
0278: * @see #doubleValue
0279: * @see #set
0280: */
0281: boolean isExpectedParameter(final Collection expected,
0282: final ParameterDescriptor param) {
0283: return expected.contains(param);
0284: }
0285:
0286: /**
0287: * Returns the parameter value for the specified operation parameter. Values are
0288: * automatically converted into the standard units specified by the supplied
0289: * {@code param} argument, except {@link NonSI#DEGREE_ANGLE degrees} which
0290: * are converted to {@link SI#RADIAN radians}.
0291: *
0292: * @param expected The value returned by {@code getParameterDescriptors().descriptors()}.
0293: * @param param The parameter to look for.
0294: * @param group The parameter value group to search into.
0295: * @return The requested parameter value, or {@code NaN} if {@code param} is
0296: * {@linkplain MathTransformProvider#createOptionalDescriptor optional}
0297: * and the user didn't provided any value.
0298: * @throws ParameterNotFoundException if the parameter is not found.
0299: *
0300: * @see MathTransformProvider#doubleValue
0301: */
0302: final double doubleValue(final Collection expected,
0303: final ParameterDescriptor param,
0304: final ParameterValueGroup group)
0305: throws ParameterNotFoundException {
0306: if (isExpectedParameter(expected, param)) {
0307: /*
0308: * Gets the value supplied by the user. The conversion from decimal
0309: * degrees to radians (if needed) is performed by AbstractProvider.
0310: */
0311: return AbstractProvider.doubleValue(param, group);
0312: }
0313: /*
0314: * The constructor asked for a parameter value that do not apply to the type of the
0315: * projection to be created. Returns a default value common to all projection types,
0316: * but this value should not be used in projection computations.
0317: */
0318: double v;
0319: final Object value = param.getDefaultValue();
0320: if (value instanceof Number) {
0321: v = ((Number) value).doubleValue();
0322: if (NonSI.DEGREE_ANGLE.equals(param.getUnit())) {
0323: v = Math.toRadians(v);
0324: }
0325: } else {
0326: v = Double.NaN;
0327: }
0328: return v;
0329: }
0330:
0331: /**
0332: * Ensures that this projection has equals semi-major and semi-minor axis. This method is
0333: * invoked by constructors of classes implementing only spherical formulas.
0334: */
0335: final void ensureSpherical() throws IllegalArgumentException {
0336: if (!isSpherical) {
0337: throw new IllegalArgumentException(Errors
0338: .format(ErrorKeys.ELLIPTICAL_NOT_SUPPORTED));
0339: }
0340: }
0341:
0342: /**
0343: * Ensures that the absolute value of a latitude is equals to the specified value, up to
0344: * the tolerance value. The expected value is usually either {@code 0} or {@code PI/2}
0345: * (the equator or a pole).
0346: *
0347: * @param y Latitude to check, in radians.
0348: * @param expected The expected value, in radians.
0349: * @throws IllegalArgumentException if the latitude is not the expected one.
0350: */
0351: static void ensureLatitudeEquals(final ParameterDescriptor name,
0352: double y, double expected) throws IllegalArgumentException {
0353: if (!(Math.abs(Math.abs(y) - expected) < EPSILON)) {
0354: y = Math.toDegrees(y);
0355: final String n = name.getName().getCode();
0356: throw new InvalidParameterValueException(Errors.format(
0357: ErrorKeys.ILLEGAL_ARGUMENT_$2, n, new Latitude(y)),
0358: n, y);
0359: }
0360: }
0361:
0362: /**
0363: * Ensures that the latitude is within allowed limits (±π/2).
0364: * This method is useful to check the validity of projection parameters,
0365: * like {@link #latitudeOfOrigin}.
0366: *
0367: * @param y Latitude to check, in radians.
0368: * @param edge {@code true} to accept latitudes of ±π/2.
0369: * @throws IllegalArgumentException if the latitude is out of range.
0370: */
0371: static void ensureLatitudeInRange(final ParameterDescriptor name,
0372: double y, final boolean edge)
0373: throws IllegalArgumentException {
0374: if (edge ? (y >= Latitude.MIN_VALUE * Math.PI / 180 && y <= Latitude.MAX_VALUE
0375: * Math.PI / 180)
0376: : (y > Latitude.MIN_VALUE * Math.PI / 180 && y < Latitude.MAX_VALUE
0377: * Math.PI / 180)) {
0378: return;
0379: }
0380: y = Math.toDegrees(y);
0381: throw new InvalidParameterValueException(Errors.format(
0382: ErrorKeys.LATITUDE_OUT_OF_RANGE_$1, new Latitude(y)),
0383: name.getName().getCode(), y);
0384: }
0385:
0386: /**
0387: * Ensures that the longitue is within allowed limits (±π).
0388: * This method is used to check the validity of projection parameters,
0389: * like {@link #centralMeridian}.
0390: *
0391: * @param x Longitude to verify, in radians.
0392: * @param edge {@code true} for accepting longitudes of ±π.
0393: * @throws IllegalArgumentException if the longitude is out of range.
0394: */
0395: static void ensureLongitudeInRange(final ParameterDescriptor name,
0396: double x, final boolean edge)
0397: throws IllegalArgumentException {
0398: if (edge ? (x >= Longitude.MIN_VALUE * Math.PI / 180 && x <= Longitude.MAX_VALUE
0399: * Math.PI / 180)
0400: : (x > Longitude.MIN_VALUE * Math.PI / 180 && x < Longitude.MAX_VALUE
0401: * Math.PI / 180)) {
0402: return;
0403: }
0404: x = Math.toDegrees(x);
0405: throw new InvalidParameterValueException(Errors.format(
0406: ErrorKeys.LONGITUDE_OUT_OF_RANGE_$1, new Longitude(x)),
0407: name.getName().getCode(), x);
0408: }
0409:
0410: /**
0411: * Set the value in a parameter group. This convenience method is used
0412: * by subclasses for {@link #getParameterValues} implementation. Values
0413: * are automatically converted from radians to decimal degrees if needed.
0414: *
0415: * @param expected The value returned by {@code getParameterDescriptors().descriptors()}.
0416: * @param param One of the {@link AbstractProvider} constants.
0417: * @param group The group in which to set the value.
0418: * @param value The value to set.
0419: */
0420: final void set(final Collection expected,
0421: final ParameterDescriptor param,
0422: final ParameterValueGroup group, double value) {
0423: if (isExpectedParameter(expected, param)) {
0424: if (NonSI.DEGREE_ANGLE.equals(param.getUnit())) {
0425: /*
0426: * Converts radians to degrees and try to fix rounding error
0427: * (e.g. -61.500000000000014 --> -61.5). This is necessary
0428: * in order to avoid a bias when formatting a transform and
0429: * parsing it again.
0430: */
0431: value = Math.toDegrees(value);
0432: double old = value;
0433: value = XMath.fixRoundingError(value, 12);
0434: if (value == old) {
0435: /*
0436: * The attempt to fix rounding error failed. Try again with the
0437: * assumption that the true value is a multiple of 1/3 of angle
0438: * (e.g. 51.166666666666664 --> 51.166666666666666), which is
0439: * common in the EPSG database.
0440: */
0441: old *= 3;
0442: final double test = XMath.fixRoundingError(old, 12);
0443: if (test != old) {
0444: value = test / 3;
0445: }
0446: }
0447: }
0448: group.parameter(param.getName().getCode()).setValue(value);
0449: }
0450: }
0451:
0452: /**
0453: * Returns the parameter descriptors for this map projection.
0454: * This is used for a providing a default implementation of
0455: * {@link #getParameterValues}, as well as arguments checking.
0456: */
0457: public abstract ParameterDescriptorGroup getParameterDescriptors();
0458:
0459: /**
0460: * Returns the parameter values for this map projection.
0461: *
0462: * @return A copy of the parameter values for this map projection.
0463: */
0464: public ParameterValueGroup getParameterValues() {
0465: final ParameterDescriptorGroup descriptor = getParameterDescriptors();
0466: final Collection expected = descriptor.descriptors();
0467: // TODO: remove the cast below once we will be allowed to use J2SE 1.5.
0468: final ParameterValueGroup values = (ParameterValueGroup) descriptor
0469: .createValue();
0470: set(expected, AbstractProvider.SEMI_MAJOR, values, semiMajor);
0471: set(expected, AbstractProvider.SEMI_MINOR, values, semiMinor);
0472: set(expected, AbstractProvider.CENTRAL_MERIDIAN, values,
0473: centralMeridian);
0474: set(expected, AbstractProvider.LATITUDE_OF_ORIGIN, values,
0475: latitudeOfOrigin);
0476: set(expected, AbstractProvider.SCALE_FACTOR, values,
0477: scaleFactor);
0478: set(expected, AbstractProvider.FALSE_EASTING, values,
0479: falseEasting);
0480: set(expected, AbstractProvider.FALSE_NORTHING, values,
0481: falseNorthing);
0482: return values;
0483: }
0484:
0485: /**
0486: * Returns the dimension of input points.
0487: */
0488: public final int getSourceDimensions() {
0489: return 2;
0490: }
0491:
0492: /**
0493: * Returns the dimension of output points.
0494: */
0495: public final int getTargetDimensions() {
0496: return 2;
0497: }
0498:
0499: //////////////////////////////////////////////////////////////////////////////////////////
0500: //////// ////////
0501: //////// TRANSFORMATION METHODS ////////
0502: //////// Includes an inner class for inverse projections. ////////
0503: //////// ////////
0504: //////////////////////////////////////////////////////////////////////////////////////////
0505:
0506: /**
0507: * Returns the orthodromic distance between the two specified points using a spherical
0508: * approximation. This is used for assertions only.
0509: */
0510: private double orthodromicDistance(final Point2D source,
0511: final Point2D target) {
0512: final double y1 = Math.toRadians(source.getY());
0513: final double y2 = Math.toRadians(target.getY());
0514: final double dx = Math.toRadians(Math.abs(target.getX()
0515: - source.getX()) % 360);
0516: double rho = Math.sin(y1) * Math.sin(y2) + Math.cos(y1)
0517: * Math.cos(y2) * Math.cos(dx);
0518: if (rho > +1) {
0519: assert rho <= +(1 + EPSILON) : rho;
0520: rho = +1;
0521: }
0522: if (rho < -1) {
0523: assert rho >= -(1 + EPSILON) : rho;
0524: rho = -1;
0525: }
0526: return Math.acos(rho) * semiMajor;
0527: }
0528:
0529: /**
0530: * Check point for private use by {@link #checkReciprocal}. This class is necessary in order
0531: * to avoid never-ending loop in {@code assert} statements (when an {@code assert}
0532: * calls {@code transform(...)}, which calls {@code inverse.transform(...)}, which
0533: * calls {@code transform(...)}, etc.).
0534: */
0535: private static final class CheckPoint extends Point2D.Double {
0536: public CheckPoint(final Point2D point) {
0537: super (point.getX(), point.getY());
0538: }
0539: }
0540:
0541: /**
0542: * Check if the transform of {@code point} is close enough to {@code target}.
0543: * "Close enough" means that the two points are separated by a distance shorter than
0544: * {@link #getToleranceForAssertions}. This method is used for assertions with J2SE 1.4.
0545: *
0546: * @param point Point to transform, in decimal degrees if {@code inverse} is {@code false}.
0547: * @param target Point to compare to, in metres if {@code inverse} is {@code false}.
0548: * @param inverse {@code true} for an inverse transform instead of a direct one.
0549: * @return {@code true} if the two points are close enough.
0550: */
0551: private boolean checkReciprocal(Point2D point,
0552: final Point2D target, final boolean inverse)
0553: throws ProjectionException {
0554: if (!(point instanceof CheckPoint))
0555: try {
0556: point = new CheckPoint(point);
0557: final double longitude;
0558: final double latitude;
0559: final double distance;
0560: if (inverse) {
0561: // Computes orthodromic distance (spherical model) in metres.
0562: point = ((MathTransform2D) inverse()).transform(
0563: point, point);
0564: distance = orthodromicDistance(point, target);
0565: longitude = point.getX();
0566: latitude = point.getY();
0567: } else {
0568: // Computes cartesian distance in metres.
0569: longitude = point.getX();
0570: latitude = point.getY();
0571: point = transform(point, point);
0572: distance = point.distance(target);
0573: }
0574: if (distance > getToleranceForAssertions(longitude,
0575: latitude)) {
0576: /*
0577: * Do not fail for NaN values. For other cases we must throw a ProjectionException,
0578: * not an AssertionError, because some code like CRS.transform(CoordinateOperation,
0579: * ...) will project points that are know to be suspicious by surrounding them in
0580: * "try ... catch" statements. Failure are normal in their case and we want to let
0581: * them handle the exception the way they are used to.
0582: */
0583: throw new ProjectionException(
0584: Errors
0585: .format(
0586: ErrorKeys.PROJECTION_CHECK_FAILED_$4,
0587: new Double(distance),
0588: new Longitude(
0589: longitude
0590: - Math
0591: .toDegrees(centralMeridian)),
0592: new Latitude(
0593: latitude
0594: - Math
0595: .toDegrees(latitudeOfOrigin)),
0596: getParameterDescriptors()
0597: .getName()
0598: .getCode()));
0599: }
0600: } catch (TransformException exception) {
0601: final ProjectionException error = new ProjectionException(
0602: exception.getLocalizedMessage());
0603: error.initCause(exception);
0604: throw error;
0605: }
0606: return true;
0607: }
0608:
0609: /**
0610: * Checks if transform using spherical formulas produces the same result
0611: * than ellipsoidal formulas. This method is invoked during assertions only.
0612: *
0613: * @param x The easting computed by spherical formulas, in metres.
0614: * @param y The northing computed by spherical formulas, in metres.
0615: * @param expected The (easting,northing) computed by ellipsoidal formulas.
0616: * @param tolerance The tolerance (optional).
0617: */
0618: static boolean checkTransform(final double x, final double y,
0619: final Point2D expected, final double tolerance) {
0620: compare("x", expected.getX(), x, tolerance);
0621: compare("y", expected.getY(), y, tolerance);
0622: return tolerance < Double.POSITIVE_INFINITY;
0623: }
0624:
0625: /**
0626: * Default version of {@link #checkTransform(double,double,Point2D,double)}.
0627: */
0628: static boolean checkTransform(final double x, final double y,
0629: final Point2D expected) {
0630: return checkTransform(x, y, expected, EPSILON);
0631: }
0632:
0633: /**
0634: * Checks if inverse transform using spherical formulas produces the same result
0635: * than ellipsoidal formulas. This method is invoked during assertions only.
0636: * <p>
0637: * <strong>Note:</strong> this method ignores the longitude if the latitude is
0638: * at a pole, because in such case the longitude is meanless.
0639: *
0640: * @param longitude The longitude computed by spherical formulas, in radians.
0641: * @param latitude The latitude computed by spherical formulas, in radians.
0642: * @param expected The (longitude,latitude) computed by ellipsoidal formulas.
0643: * @param tolerance The tolerance (optional).
0644: */
0645: static boolean checkInverseTransform(final double longitude,
0646: final double latitude, final Point2D expected,
0647: final double tolerance) {
0648: compare("latitude", expected.getY(), latitude, tolerance);
0649: if (Math.abs(Math.PI / 2 - Math.abs(latitude)) > EPSILON) {
0650: compare("longitude", expected.getX(), longitude, tolerance);
0651: }
0652: return tolerance < Double.POSITIVE_INFINITY;
0653: }
0654:
0655: /**
0656: * Default version of {@link #checkInverseTransform(double,double,Point2D,double)}.
0657: */
0658: static boolean checkInverseTransform(double longitude,
0659: double latitude, Point2D expected) {
0660: return checkInverseTransform(longitude, latitude, expected,
0661: EPSILON);
0662: }
0663:
0664: /**
0665: * Compares two value for equality up to some tolerance threshold. This is used during
0666: * assertions only. The comparaison do not fails if at least one value to compare is
0667: * {@link Double#NaN}.
0668: * <p>
0669: * <strong>Hack:</strong> if the {@code variable} name starts by lower-case {@code L}
0670: * (as in "longitude" and "latitude"), then the value is assumed to be an angle in
0671: * radians. This is used for formatting an error message, if needed.
0672: */
0673: private static void compare(String variable, double expected,
0674: double actual, double tolerance) {
0675: if (Math.abs(expected - actual) > tolerance) {
0676: if (variable.charAt(0) == 'l') {
0677: actual = Math.toDegrees(actual);
0678: expected = Math.toDegrees(expected);
0679: }
0680: throw new AssertionError(Errors.format(
0681: ErrorKeys.TEST_FAILURE_$3, variable, new Double(
0682: expected), new Double(actual)));
0683: }
0684: }
0685:
0686: /**
0687: * Transforms the specified coordinate and stores the result in {@code ptDst}.
0688: * This method returns longitude as <var>x</var> values in the range {@code [-PI..PI]}
0689: * and latitude as <var>y</var> values in the range {@code [-PI/2..PI/2]}. It will be
0690: * checked by the caller, so this method doesn't need to performs this check.
0691: * <p>
0692: *
0693: * Input coordinates are also guarenteed to have the {@link #falseEasting}
0694: * and {@link #falseNorthing} removed and be divided by {@link #globalScale}
0695: * before this method is invoked. After this method is invoked, the
0696: * {@link #centralMeridian} is added to the {@code x} results
0697: * in {@code ptDst}. This means that projections that implement this method
0698: * are performed on an ellipse (or sphere) with a semiMajor axis of 1.0.
0699: * <p>
0700: *
0701: * In <A HREF="http://www.remotesensing.org/proj/">PROJ.4</A>, the same
0702: * standardization, described above, is handled by {@code pj_inv.c}.
0703: * Therefore when porting projections from PROJ.4, the inverse transform
0704: * equations can be used directly here with minimal change.
0705: * In the equations of Snyder, {@link #falseEasting}, {@link #falseNorthing}
0706: * and {@link #scaleFactor} are usually not given.
0707: * When implementing these equations here, you will not
0708: * need to add the {@link #centralMeridian} to the output longitude or remove the
0709: * {@link #semiMajor} ('<var>a</var>' or '<var>R</var>').
0710: *
0711: * @param x The easting of the coordinate, linear distance on a unit sphere or ellipse.
0712: * @param y The northing of the coordinate, linear distance on a unit sphere or ellipse.
0713: * @param ptDst the specified coordinate point that stores the result of transforming
0714: * {@code ptSrc}, or {@code null}. Ordinates will be in
0715: * <strong>radians</strong>.
0716: * @return the coordinate point after transforming {@code x}, {@code y}
0717: * and storing the result in {@code ptDst}.
0718: * @throws ProjectionException if the point can't be transformed.
0719: *
0720: * @todo The {@code ptDst} argument will be removed and the return type changed if
0721: * RFE <A HREF="http://developer.java.sun.com/developer/bugParade/bugs/4222792.html">4222792</A>
0722: * is implemented efficiently in a future J2SE release (maybe J2SE 1.6?).
0723: */
0724: protected abstract Point2D inverseTransformNormalized(double x,
0725: double y, final Point2D ptDst) throws ProjectionException;
0726:
0727: /**
0728: * Transforms the specified coordinate and stores the result in {@code ptDst}.
0729: * This method is guaranteed to be invoked with values of <var>x</var> in the range
0730: * {@code [-PI..PI]} and values of <var>y</var> in the range {@code [-PI/2..PI/2]}.
0731: * <p>
0732: *
0733: * Coordinates are also guaranteed to have the {@link #centralMeridian}
0734: * removed from <var>x</var> before this method is invoked. After this method
0735: * is invoked, the results in {@code ptDst} are multiplied by {@link #globalScale},
0736: * and the {@link #falseEasting} and {@link #falseNorthing} are added.
0737: * This means that projections that implement this method are performed on an
0738: * ellipse (or sphere) with a semiMajor axis of 1.0.
0739: * <p>
0740: *
0741: * In <A HREF="http://www.remotesensing.org/proj/">PROJ.4</A>, the same
0742: * standardization, described above, is handled by {@code pj_fwd.c}.
0743: * Therefore when porting projections from PROJ.4, the forward transform equations can
0744: * be used directly here with minimal change. In the equations of Snyder,
0745: * {@link #falseEasting}, {@link #falseNorthing} and {@link #scaleFactor}
0746: * are usually not given. When implementing these equations here, you will not
0747: * need to remove the {@link #centralMeridian} from <var>x</var> or apply the
0748: * {@link #semiMajor} ('<var>a</var>' or '<var>R</var>').
0749: *
0750: * @param x The longitude of the coordinate, in <strong>radians</strong>.
0751: * @param y The latitude of the coordinate, in <strong>radians</strong>.
0752: * @param ptDst the specified coordinate point that stores the result of transforming
0753: * {@code ptSrc}, or {@code null}. Ordinates will be in a
0754: * dimensionless unit, as a linear distance on a unit sphere or ellipse.
0755: * @return the coordinate point after transforming {@code x}, {@code y}
0756: * and storing the result in {@code ptDst}.
0757: * @throws ProjectionException if the point can't be transformed.
0758: *
0759: * @todo The {@code ptDst} argument will be removed and the return type changed if
0760: * RFE <A HREF="http://developer.java.sun.com/developer/bugParade/bugs/4222792.html">4222792</A>
0761: * is implemented efficiently in a future J2SE release (maybe J2SE 1.6?).
0762: */
0763: protected abstract Point2D transformNormalized(double x, double y,
0764: final Point2D ptDst) throws ProjectionException;
0765:
0766: /**
0767: * Transforms the specified {@code ptSrc} and stores the result in {@code ptDst}.
0768: * <p>
0769: *
0770: * This method standardizes the source {@code x} coordinate
0771: * by removing the {@link #centralMeridian}, before invoking
0772: * <code>{@link #transformNormalized transformNormalized}(x, y, ptDst)</code>.
0773: * It also multiplies by {@link #globalScale} and adds the {@link #falseEasting} and
0774: * {@link #falseNorthing} to the point returned by the {@code transformNormalized(...)}
0775: * call.
0776: *
0777: * @param ptSrc the specified coordinate point to be transformed.
0778: * Ordinates must be in decimal degrees.
0779: * @param ptDst the specified coordinate point that stores the result of transforming
0780: * {@code ptSrc}, or {@code null}. Ordinates will be in metres.
0781: * @return the coordinate point after transforming {@code ptSrc} and storing
0782: * the result in {@code ptDst}.
0783: * @throws ProjectionException if the point can't be transformed.
0784: */
0785: public final Point2D transform(final Point2D ptSrc, Point2D ptDst)
0786: throws ProjectionException {
0787: final double x = ptSrc.getX();
0788: final double y = ptSrc.getY();
0789: // Note: the following tests should not fails for NaN values.
0790: if (x < (Longitude.MIN_VALUE - ANGLE_TOLERANCE)
0791: || x > (Longitude.MAX_VALUE + ANGLE_TOLERANCE)) {
0792: throw new PointOutsideEnvelopeException(Errors.format(
0793: ErrorKeys.LONGITUDE_OUT_OF_RANGE_$1, new Longitude(
0794: x)));
0795: }
0796: if (y < (Latitude.MIN_VALUE - ANGLE_TOLERANCE)
0797: || y > (Latitude.MAX_VALUE + ANGLE_TOLERANCE)) {
0798: throw new PointOutsideEnvelopeException(Errors
0799: .format(ErrorKeys.LATITUDE_OUT_OF_RANGE_$1,
0800: new Latitude(y)));
0801: }
0802: /*
0803: * Makes sure that the longitude before conversion stay within +/- PI radians. As a
0804: * special case, we do not check the range if no rotation were applied on the longitude.
0805: * This is because the user may have a big area ranging from -180° to +180°. With the
0806: * slight rounding errors related to map projections, the 180° longitude may be slightly
0807: * over the limit. Rolling the longitude would changes its sign. For example a bounding
0808: * box from 30° to +180° would become 30° to -180°, which is probably not what the user
0809: * wanted.
0810: */
0811: ptDst = transformNormalized(
0812: centralMeridian != 0 ? rollLongitude(Math.toRadians(x)
0813: - centralMeridian) : Math.toRadians(x), Math
0814: .toRadians(y), ptDst);
0815: ptDst.setLocation(globalScale * ptDst.getX() + falseEasting,
0816: globalScale * ptDst.getY() + falseNorthing);
0817:
0818: assert checkReciprocal(ptDst, (ptSrc != ptDst) ? ptSrc
0819: : new Point2D.Double(x, y), true);
0820: return ptDst;
0821: }
0822:
0823: /**
0824: * Transforms a list of coordinate point ordinal values. Ordinates must be
0825: * (<var>longitude</var>,<var>latitude</var>) pairs in decimal degrees.
0826: *
0827: * @throws ProjectionException if a point can't be transformed. This method try
0828: * to transform every points even if some of them can't be transformed.
0829: * Non-transformable points will have value {@link Double#NaN}. If more
0830: * than one point can't be transformed, then this exception may be about
0831: * an arbitrary point.
0832: */
0833: public final void transform(final double[] src, int srcOffset,
0834: final double[] dest, int dstOffset, int numPts)
0835: throws ProjectionException {
0836: /*
0837: * Vérifie s'il faudra parcourir le tableau en sens inverse.
0838: * Ce sera le cas si les tableaux source et destination se
0839: * chevauchent et que la destination est après la source.
0840: */
0841: final boolean reverse = (src == dest && srcOffset < dstOffset && srcOffset
0842: + (2 * numPts) > dstOffset);
0843: if (reverse) {
0844: srcOffset += 2 * numPts;
0845: dstOffset += 2 * numPts;
0846: }
0847: final Point2D.Double point = new Point2D.Double();
0848: ProjectionException firstException = null;
0849: while (--numPts >= 0) {
0850: try {
0851: point.x = src[srcOffset++];
0852: point.y = src[srcOffset++];
0853: transform(point, point);
0854: dest[dstOffset++] = point.x;
0855: dest[dstOffset++] = point.y;
0856: } catch (ProjectionException exception) {
0857: dest[dstOffset++] = Double.NaN;
0858: dest[dstOffset++] = Double.NaN;
0859: if (firstException == null) {
0860: firstException = exception;
0861: }
0862: }
0863: if (reverse) {
0864: srcOffset -= 4;
0865: dstOffset -= 4;
0866: }
0867: }
0868: if (firstException != null) {
0869: throw firstException;
0870: }
0871: }
0872:
0873: /**
0874: * Transforms a list of coordinate point ordinal values. Ordinates must be
0875: * (<var>longitude</var>,<var>latitude</var>) pairs in decimal degrees.
0876: *
0877: * @throws ProjectionException if a point can't be transformed. This method try
0878: * to transform every points even if some of them can't be transformed.
0879: * Non-transformable points will have value {@link Float#NaN}. If more
0880: * than one point can't be transformed, then this exception may be about
0881: * an arbitrary point.
0882: */
0883: public final void transform(final float[] src, int srcOffset,
0884: final float[] dest, int dstOffset, int numPts)
0885: throws ProjectionException {
0886: final boolean reverse = (src == dest && srcOffset < dstOffset && srcOffset
0887: + (2 * numPts) > dstOffset);
0888: if (reverse) {
0889: srcOffset += 2 * numPts;
0890: dstOffset += 2 * numPts;
0891: }
0892: final Point2D.Double point = new Point2D.Double();
0893: ProjectionException firstException = null;
0894: while (--numPts >= 0) {
0895: try {
0896: point.x = src[srcOffset++];
0897: point.y = src[srcOffset++];
0898: transform(point, point);
0899: dest[dstOffset++] = (float) point.x;
0900: dest[dstOffset++] = (float) point.y;
0901: } catch (ProjectionException exception) {
0902: dest[dstOffset++] = Float.NaN;
0903: dest[dstOffset++] = Float.NaN;
0904: if (firstException == null) {
0905: firstException = exception;
0906: }
0907: }
0908: if (reverse) {
0909: srcOffset -= 4;
0910: dstOffset -= 4;
0911: }
0912: }
0913: if (firstException != null) {
0914: throw firstException;
0915: }
0916: }
0917:
0918: /**
0919: * Inverse of a map projection. Will be created by {@link MapProjection#inverse()} only when
0920: * first required. Implementation of {@code transform(...)} methods are mostly identical
0921: * to {@code MapProjection.transform(...)}, except that they will invokes
0922: * {@link MapProjection#inverseTransformNormalized} instead of
0923: * {@link MapProjection#transformNormalized}.
0924: *
0925: * @version $Id: MapProjection.java 26518 2007-08-10 17:05:01Z desruisseaux $
0926: * @author Martin Desruisseaux
0927: */
0928: private final class Inverse extends AbstractMathTransform.Inverse
0929: implements MathTransform2D {
0930: /**
0931: * Default constructor.
0932: */
0933: public Inverse() {
0934: MapProjection.this .super ();
0935: }
0936:
0937: /**
0938: * Inverse transforms the specified {@code ptSrc} and stores the result in {@code ptDst}.
0939: * <p>
0940: *
0941: * This method standardizes the {@code ptSrc} by removing the
0942: * {@link #falseEasting} and {@link #falseNorthing} and dividing by
0943: * {@link #globalScale} before invoking
0944: * <code>{@link #inverseTransformNormalized inverseTransformNormalized}(x, y, ptDst)</code>.
0945: * It then adds the {@link #centralMeridian} to the {@code x} of the
0946: * point returned by the {@code inverseTransformNormalized} call.
0947: *
0948: * @param ptSrc the specified coordinate point to be transformed.
0949: * Ordinates must be in metres.
0950: * @param ptDst the specified coordinate point that stores the result of transforming
0951: * {@code ptSrc}, or {@code null}. Ordinates will be in decimal degrees.
0952: * @return the coordinate point after transforming {@code ptSrc}
0953: * and stroring the result in {@code ptDst}.
0954: * @throws ProjectionException if the point can't be transformed.
0955: */
0956: public final Point2D transform(final Point2D ptSrc,
0957: Point2D ptDst) throws ProjectionException {
0958: final double x0 = ptSrc.getX();
0959: final double y0 = ptSrc.getY();
0960: ptDst = inverseTransformNormalized((x0 - falseEasting)
0961: / globalScale, (y0 - falseNorthing) / globalScale,
0962: ptDst);
0963: /*
0964: * Makes sure that the longitude after conversion stay within +/- PI radians. As a
0965: * special case, we do not check the range if no rotation were applied on the longitude.
0966: * This is because the user may have a big area ranging from -180° to +180°. With the
0967: * slight rounding errors related to map projections, the 180° longitude may be slightly
0968: * over the limit. Rolling the longitude would changes its sign. For example a bounding
0969: * box from 30° to +180° would become 30° to -180°, which is probably not what the user
0970: * wanted.
0971: */
0972: final double x = Math
0973: .toDegrees(centralMeridian != 0 ? rollLongitude(ptDst
0974: .getX()
0975: + centralMeridian)
0976: : ptDst.getX());
0977: final double y = Math.toDegrees(ptDst.getY());
0978: ptDst.setLocation(x, y);
0979:
0980: // Note: the following tests should not fails for NaN values.
0981: if (x < (Longitude.MIN_VALUE - ANGLE_TOLERANCE)
0982: || x > (Longitude.MAX_VALUE + ANGLE_TOLERANCE)) {
0983: throw new PointOutsideEnvelopeException(Errors.format(
0984: ErrorKeys.LONGITUDE_OUT_OF_RANGE_$1,
0985: new Longitude(x)));
0986: }
0987: if (y < (Latitude.MIN_VALUE - ANGLE_TOLERANCE)
0988: || y > (Latitude.MAX_VALUE + ANGLE_TOLERANCE)) {
0989: throw new PointOutsideEnvelopeException(Errors.format(
0990: ErrorKeys.LATITUDE_OUT_OF_RANGE_$1,
0991: new Latitude(y)));
0992: }
0993: assert checkReciprocal(ptDst, (ptSrc != ptDst) ? ptSrc
0994: : new Point2D.Double(x0, y0), false);
0995: return ptDst;
0996: }
0997:
0998: /**
0999: * Inverse transforms a list of coordinate point ordinal values.
1000: * Ordinates must be (<var>x</var>,<var>y</var>) pairs in metres.
1001: *
1002: * @throws ProjectionException if a point can't be transformed. This method try
1003: * to transform every points even if some of them can't be transformed.
1004: * Non-transformable points will have value {@link Double#NaN}. If more
1005: * than one point can't be transformed, then this exception may be about
1006: * an arbitrary point.
1007: */
1008: public final void transform(final double[] src, int srcOffset,
1009: final double[] dest, int dstOffset, int numPts)
1010: throws TransformException {
1011: /*
1012: * Vérifie s'il faudra parcourir le tableau en sens inverse.
1013: * Ce sera le cas si les tableaux source et destination se
1014: * chevauchent et que la destination est après la source.
1015: */
1016: final boolean reverse = (src == dest
1017: && srcOffset < dstOffset && srcOffset
1018: + (2 * numPts) > dstOffset);
1019: if (reverse) {
1020: srcOffset += 2 * numPts;
1021: dstOffset += 2 * numPts;
1022: }
1023: final Point2D.Double point = new Point2D.Double();
1024: ProjectionException firstException = null;
1025: while (--numPts >= 0) {
1026: try {
1027: point.x = src[srcOffset++];
1028: point.y = src[srcOffset++];
1029: transform(point, point);
1030: dest[dstOffset++] = point.x;
1031: dest[dstOffset++] = point.y;
1032: } catch (ProjectionException exception) {
1033: dest[dstOffset++] = Double.NaN;
1034: dest[dstOffset++] = Double.NaN;
1035: if (firstException == null) {
1036: firstException = exception;
1037: }
1038: }
1039: if (reverse) {
1040: srcOffset -= 4;
1041: dstOffset -= 4;
1042: }
1043: }
1044: if (firstException != null) {
1045: throw firstException;
1046: }
1047: }
1048:
1049: /**
1050: * Inverse transforms a list of coordinate point ordinal values.
1051: * Ordinates must be (<var>x</var>,<var>y</var>) pairs in metres.
1052: *
1053: * @throws ProjectionException if a point can't be transformed. This method try
1054: * to transform every points even if some of them can't be transformed.
1055: * Non-transformable points will have value {@link Float#NaN}. If more
1056: * than one point can't be transformed, then this exception may be about
1057: * an arbitrary point.
1058: */
1059: public final void transform(final float[] src, int srcOffset,
1060: final float[] dest, int dstOffset, int numPts)
1061: throws ProjectionException {
1062: final boolean reverse = (src == dest
1063: && srcOffset < dstOffset && srcOffset
1064: + (2 * numPts) > dstOffset);
1065: if (reverse) {
1066: srcOffset += 2 * numPts;
1067: dstOffset += 2 * numPts;
1068: }
1069: final Point2D.Double point = new Point2D.Double();
1070: ProjectionException firstException = null;
1071: while (--numPts >= 0) {
1072: try {
1073: point.x = src[srcOffset++];
1074: point.y = src[srcOffset++];
1075: transform(point, point);
1076: dest[dstOffset++] = (float) point.x;
1077: dest[dstOffset++] = (float) point.y;
1078: } catch (ProjectionException exception) {
1079: dest[dstOffset++] = Float.NaN;
1080: dest[dstOffset++] = Float.NaN;
1081: if (firstException == null) {
1082: firstException = exception;
1083: }
1084: }
1085: if (reverse) {
1086: srcOffset -= 4;
1087: dstOffset -= 4;
1088: }
1089: }
1090: if (firstException != null) {
1091: throw firstException;
1092: }
1093: }
1094: }
1095:
1096: /**
1097: * Returns the inverse of this map projection.
1098: */
1099: public final MathTransform inverse() {
1100: // No synchronization. Not a big deal if this method is invoked in
1101: // the same time by two threads resulting in two instances created.
1102: if (inverse == null) {
1103: inverse = new Inverse();
1104: }
1105: return inverse;
1106: }
1107:
1108: /**
1109: * Maximal error (in metres) tolerated for assertions, if enabled. When assertions are enabled,
1110: * every direct projection is followed by an inverse projection, and the result is compared to
1111: * the original coordinate. If a distance greater than the tolerance level is found, then an
1112: * {@link ProjectionException} will be thrown. Subclasses should override this method if they
1113: * need to relax the tolerance level.
1114: *
1115: * @param longitude The longitude in decimal degrees.
1116: * @param latitude The latitude in decimal degrees.
1117: * @return The tolerance level for assertions, in meters.
1118: */
1119: protected double getToleranceForAssertions(final double longitude,
1120: final double latitude) {
1121: final double delta = Math.abs(longitude - centralMeridian) / 2
1122: + Math.abs(latitude - latitudeOfOrigin);
1123: if (delta > 40) {
1124: // When far from the valid area, use a larger tolerance.
1125: return 1;
1126: }
1127: // Be less strict when the point is near an edge.
1128: return (Math.abs(longitude) > 179) || (Math.abs(latitude) > 89) ? 1E-1
1129: : EPSILON;
1130: }
1131:
1132: //////////////////////////////////////////////////////////////////////////////////////////
1133: //////// ////////
1134: //////// IMPLEMENTATION OF Object AND MathTransform2D STANDARD METHODS ////////
1135: //////// ////////
1136: //////////////////////////////////////////////////////////////////////////////////////////
1137:
1138: /**
1139: * Returns a hash value for this map projection.
1140: */
1141: public int hashCode() {
1142: long code = Double.doubleToLongBits(semiMajor);
1143: code = code * 37 + Double.doubleToLongBits(semiMinor);
1144: code = code * 37 + Double.doubleToLongBits(centralMeridian);
1145: code = code * 37 + Double.doubleToLongBits(latitudeOfOrigin);
1146: return (int) code ^ (int) (code >>> 32);
1147: }
1148:
1149: /**
1150: * Compares the specified object with this map projection for equality.
1151: */
1152: public boolean equals(final Object object) {
1153: // Do not check 'object==this' here, since this
1154: // optimization is usually done in subclasses.
1155: if (super .equals(object)) {
1156: final MapProjection that = (MapProjection) object;
1157: return equals(this .semiMajor, that.semiMajor)
1158: && equals(this .semiMinor, that.semiMinor)
1159: && equals(this .centralMeridian,
1160: that.centralMeridian)
1161: && equals(this .latitudeOfOrigin,
1162: that.latitudeOfOrigin)
1163: && equals(this .scaleFactor, that.scaleFactor)
1164: && equals(this .falseEasting, that.falseEasting)
1165: && equals(this .falseNorthing, that.falseNorthing);
1166: }
1167: return false;
1168: }
1169:
1170: /**
1171: * Returns {@code true} if the two specified value are equals.
1172: * Two {@link Double#NaN NaN} values are considered equals.
1173: */
1174: static boolean equals(final double value1, final double value2) {
1175: return Double.doubleToLongBits(value1) == Double
1176: .doubleToLongBits(value2);
1177: }
1178:
1179: //////////////////////////////////////////////////////////////////////////////////////////
1180: //////// ////////
1181: //////// FORMULAS FROM SNYDER ////////
1182: //////// ////////
1183: //////////////////////////////////////////////////////////////////////////////////////////
1184:
1185: /**
1186: * Iteratively solve equation (7-9) from Snyder.
1187: */
1188: final double cphi2(final double ts) throws ProjectionException {
1189: final double eccnth = 0.5 * excentricity;
1190: double phi = (Math.PI / 2) - 2.0 * Math.atan(ts);
1191: for (int i = 0; i < MAXIMUM_ITERATIONS; i++) {
1192: final double con = excentricity * Math.sin(phi);
1193: final double dphi = (Math.PI / 2)
1194: - 2.0
1195: * Math.atan(ts
1196: * Math.pow((1 - con) / (1 + con), eccnth))
1197: - phi;
1198: phi += dphi;
1199: if (Math.abs(dphi) <= ITERATION_TOLERANCE) {
1200: return phi;
1201: }
1202: }
1203: throw new ProjectionException(Errors
1204: .format(ErrorKeys.NO_CONVERGENCE));
1205: }
1206:
1207: /**
1208: * Compute function <code>f(s,c,e²) = c/sqrt(1 - s²×e²)</code> needed for the true scale
1209: * latitude (Snyder 14-15), where <var>s</var> and <var>c</var> are the sine and cosine of
1210: * the true scale latitude, and <var>e²</var> is the {@linkplain #excentricitySquared
1211: * eccentricity squared}.
1212: */
1213: final double msfn(final double s, final double c) {
1214: return c / Math.sqrt(1.0 - (s * s) * excentricitySquared);
1215: }
1216:
1217: /**
1218: * Compute function (15-9) and (9-13) from Snyder.
1219: * Equivalent to negative of function (7-7).
1220: */
1221: final double tsfn(final double phi, double sinphi) {
1222: sinphi *= excentricity;
1223: /*
1224: * NOTE: change sign to get the equivalent of Snyder (7-7).
1225: */
1226: return Math.tan(0.5 * ((Math.PI / 2) - phi))
1227: / Math.pow((1 - sinphi) / (1 + sinphi),
1228: 0.5 * excentricity);
1229: }
1230:
1231: //////////////////////////////////////////////////////////////////////////////////////////
1232: //////// ////////
1233: //////// PROVIDER ////////
1234: //////// ////////
1235: //////////////////////////////////////////////////////////////////////////////////////////
1236:
1237: /**
1238: * The base provider for {@link MapProjection}s.
1239: *
1240: * @version $Id: MapProjection.java 26518 2007-08-10 17:05:01Z desruisseaux $
1241: * @author Martin Desruisseaux
1242: */
1243: public static abstract class AbstractProvider extends
1244: MathTransformProvider {
1245: /**
1246: * Serial number for interoperability with different versions.
1247: */
1248: private static final long serialVersionUID = 6280666068007678702L;
1249:
1250: /**
1251: * The operation parameter descriptor for the {@linkplain #semiMajor semi major} parameter
1252: * value. Valid values range is from 0 to infinity. This parameter is mandatory.
1253: *
1254: * @todo Would like to start range from 0 <u>exclusive</u>.
1255: */
1256: public static final ParameterDescriptor SEMI_MAJOR = createDescriptor(
1257: new NamedIdentifier[] {
1258: new NamedIdentifier(Citations.OGC, "semi_major"),
1259: new NamedIdentifier(Citations.EPSG,
1260: "semi-major axis")
1261: // EPSG does not specifically define the above parameter
1262: }, Double.NaN, 0, Double.POSITIVE_INFINITY, SI.METER);
1263:
1264: /**
1265: * The operation parameter descriptor for the {@linkplain #semiMinor semi minor} parameter
1266: * value. Valid values range is from 0 to infinity. This parameter is mandatory.
1267: *
1268: * @todo Would like to start range from 0 <u>exclusive</u>.
1269: */
1270: public static final ParameterDescriptor SEMI_MINOR = createDescriptor(
1271: new NamedIdentifier[] {
1272: new NamedIdentifier(Citations.OGC, "semi_minor"),
1273: new NamedIdentifier(Citations.EPSG,
1274: "semi-minor axis")
1275: // EPSG does not specifically define the above parameter
1276: }, Double.NaN, 0, Double.POSITIVE_INFINITY, SI.METER);
1277:
1278: /**
1279: * The operation parameter descriptor for the {@linkplain #centralMeridian central meridian}
1280: * parameter value. Valid values range is from -180 to 180°. Default value is 0.
1281: */
1282: public static final ParameterDescriptor CENTRAL_MERIDIAN = createDescriptor(
1283: new NamedIdentifier[] {
1284: new NamedIdentifier(Citations.OGC,
1285: "central_meridian"),
1286: new NamedIdentifier(Citations.EPSG,
1287: "Longitude of natural origin"),
1288: new NamedIdentifier(Citations.EPSG,
1289: "Longitude of false origin"),
1290: new NamedIdentifier(Citations.EPSG,
1291: "Longitude of origin"),
1292: new NamedIdentifier(Citations.ESRI,
1293: "Longitude_Of_Origin"),
1294: new NamedIdentifier(Citations.ESRI,
1295: "Longitude_Of_Center"),
1296: new NamedIdentifier(Citations.GEOTIFF,
1297: "NatOriginLong")
1298: // ESRI uses "Longitude_Of_Origin" in orthographic (not to
1299: // be confused with "Longitude_Of_Center" in oblique mercator).
1300: }, 0, -180, 180, NonSI.DEGREE_ANGLE);
1301:
1302: /**
1303: * The operation parameter descriptor for the {@linkplain #latitudeOfOrigin latitude of origin}
1304: * parameter value. Valid values range is from -90 to 90°. Default value is 0.
1305: */
1306: public static final ParameterDescriptor LATITUDE_OF_ORIGIN = createDescriptor(
1307: new NamedIdentifier[] {
1308: new NamedIdentifier(Citations.OGC,
1309: "latitude_of_origin"),
1310: new NamedIdentifier(Citations.EPSG,
1311: "Latitude of false origin"),
1312: new NamedIdentifier(Citations.EPSG,
1313: "Latitude of natural origin"),
1314: new NamedIdentifier(Citations.ESRI,
1315: "Latitude_Of_Center"),
1316: new NamedIdentifier(Citations.GEOTIFF,
1317: "NatOriginLat")
1318: // ESRI uses "Latitude_Of_Center" in orthographic.
1319: }, 0, -90, 90, NonSI.DEGREE_ANGLE);
1320:
1321: /**
1322: * The operation parameter descriptor for the {@linkplain Mercator#standardParallel standard
1323: * parallel} parameter value. Valid values range is from -90 to 90°. Default value is 0.
1324: *
1325: * @deprecated Use {@link #STANDARD_PARALLEL_1} instead.
1326: * Not to be confused with {@link PolarStereographic.ProviderB#STANDARD_PARALLEL}.
1327: */
1328: public static final ParameterDescriptor STANDARD_PARALLEL = createDescriptor(
1329: new NamedIdentifier[] {
1330: new NamedIdentifier(Citations.OGC,
1331: "standard_parallel_1"),
1332: new NamedIdentifier(Citations.EPSG,
1333: "Latitude of 1st standard parallel"),
1334: new NamedIdentifier(Citations.GEOTIFF,
1335: "StdParallel1") }, 0, -90, 90,
1336: NonSI.DEGREE_ANGLE);
1337:
1338: /**
1339: * The operation parameter descriptor for the standard parallel 1 parameter value.
1340: * Valid values range is from -90 to 90°. Default value is 0.
1341: */
1342: public static final ParameterDescriptor STANDARD_PARALLEL_1 = createDescriptor(
1343: new NamedIdentifier[] {
1344: new NamedIdentifier(Citations.OGC,
1345: "standard_parallel_1"),
1346: new NamedIdentifier(Citations.EPSG,
1347: "Latitude of 1st standard parallel"),
1348: new NamedIdentifier(Citations.GEOTIFF,
1349: "StdParallel1") }, 0, -90, 90,
1350: NonSI.DEGREE_ANGLE);
1351:
1352: /**
1353: * The operation parameter descriptor for the standard parallel 2 parameter value.
1354: * Valid values range is from -90 to 90°. Default value is 0.
1355: */
1356: public static final ParameterDescriptor STANDARD_PARALLEL_2 = createOptionalDescriptor(
1357: new NamedIdentifier[] {
1358: new NamedIdentifier(Citations.OGC,
1359: "standard_parallel_2"),
1360: new NamedIdentifier(Citations.EPSG,
1361: "Latitude of 2nd standard parallel"),
1362: new NamedIdentifier(Citations.GEOTIFF,
1363: "StdParallel2") }, -90, 90,
1364: NonSI.DEGREE_ANGLE);
1365:
1366: /**
1367: * The operation parameter descriptor for the {@link #scaleFactor scaleFactor}
1368: * parameter value. Valid values range is from 0 to infinity. Default value is 1.
1369: *
1370: * @todo Would like to start range from 0 <u>exclusive</u>.
1371: */
1372: public static final ParameterDescriptor SCALE_FACTOR = createDescriptor(
1373: new NamedIdentifier[] {
1374: new NamedIdentifier(Citations.OGC,
1375: "scale_factor"),
1376: new NamedIdentifier(Citations.EPSG,
1377: "Scale factor at natural origin"),
1378: new NamedIdentifier(Citations.EPSG,
1379: "Scale factor on initial line"),
1380: new NamedIdentifier(Citations.GEOTIFF,
1381: "ScaleAtNatOrigin"),
1382: new NamedIdentifier(Citations.GEOTIFF,
1383: "ScaleAtCenter") }, 1, 0,
1384: Double.POSITIVE_INFINITY, Unit.ONE);
1385:
1386: /**
1387: * The operation parameter descriptor for the {@link #falseEasting falseEasting}
1388: * parameter value. Valid values range is unrestricted. Default value is 0.
1389: */
1390: public static final ParameterDescriptor FALSE_EASTING = createDescriptor(
1391: new NamedIdentifier[] {
1392: new NamedIdentifier(Citations.OGC,
1393: "false_easting"),
1394: new NamedIdentifier(Citations.EPSG,
1395: "False easting"),
1396: new NamedIdentifier(Citations.EPSG,
1397: "Easting at false origin"),
1398: new NamedIdentifier(Citations.EPSG,
1399: "Easting at projection centre"),
1400: new NamedIdentifier(Citations.GEOTIFF,
1401: "FalseEasting") }, 0,
1402: Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY,
1403: SI.METER);
1404:
1405: /**
1406: * The operation parameter descriptor for the {@link #falseNorthing falseNorthing}
1407: * parameter value. Valid values range is unrestricted. Default value is 0.
1408: */
1409: public static final ParameterDescriptor FALSE_NORTHING = createDescriptor(
1410: new NamedIdentifier[] {
1411: new NamedIdentifier(Citations.OGC,
1412: "false_northing"),
1413: new NamedIdentifier(Citations.EPSG,
1414: "False northing"),
1415: new NamedIdentifier(Citations.EPSG,
1416: "Northing at false origin"),
1417: new NamedIdentifier(Citations.EPSG,
1418: "Northing at projection centre"),
1419: new NamedIdentifier(Citations.GEOTIFF,
1420: "FalseNorthing") }, 0,
1421: Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY,
1422: SI.METER);
1423:
1424: /**
1425: * Constructs a math transform provider from a set of parameters. The provider
1426: * {@linkplain #getIdentifiers identifiers} will be the same than the parameter
1427: * ones.
1428: *
1429: * @param parameters The set of parameters (never {@code null}).
1430: */
1431: public AbstractProvider(
1432: final ParameterDescriptorGroup parameters) {
1433: super (2, 2, parameters);
1434: }
1435:
1436: /**
1437: * Returns the operation type for this map projection.
1438: */
1439: public Class getOperationType() {
1440: return Projection.class;
1441: }
1442:
1443: /**
1444: * Returns {@code true} is the parameters use a spherical datum.
1445: */
1446: static boolean isSpherical(final ParameterValueGroup values) {
1447: try {
1448: return doubleValue(SEMI_MAJOR, values) == doubleValue(
1449: SEMI_MINOR, values);
1450: } catch (IllegalStateException exception) {
1451: // Probably could not find the requested values -- gobble error and be forgiving.
1452: // The error will probably be thrown at MapProjection construction time, which is
1453: // less surprising to some users.
1454: return false;
1455: }
1456: }
1457:
1458: /**
1459: * Returns the parameter value for the specified operation parameter in standard units.
1460: * Values are automatically converted into the standard units specified by the supplied
1461: * {@code param} argument, except {@link NonSI#DEGREE_ANGLE degrees} which are converted
1462: * to {@link SI#RADIAN radians}. This conversion is performed because the radians units
1463: * are standard for all internal computations in the map projection package. For example
1464: * they are the standard units for {@link MapProjection#latitudeOfOrigin latitudeOfOrigin}
1465: * and {@link MapProjection#centralMeridian centralMeridian} fields in the
1466: * {@link MapProjection} class.
1467: *
1468: * @param param The parameter to look for.
1469: * @param group The parameter value group to search into.
1470: * @return The requested parameter value.
1471: * @throws ParameterNotFoundException if the parameter is not found.
1472: */
1473: protected static double doubleValue(
1474: final ParameterDescriptor param,
1475: final ParameterValueGroup group)
1476: throws ParameterNotFoundException {
1477: double v = MathTransformProvider.doubleValue(param, group);
1478: if (NonSI.DEGREE_ANGLE.equals(param.getUnit())) {
1479: v = Math.toRadians(v);
1480: }
1481: return v;
1482: }
1483: }
1484: }
|