[28000] | 1 | /*
|
---|
| 2 | * GeoTools - The Open Source Java GIS Toolkit
|
---|
| 3 | * http://geotools.org
|
---|
| 4 | *
|
---|
| 5 | * (C) 1999-2008, Open Source Geospatial Foundation (OSGeo)
|
---|
| 6 | *
|
---|
| 7 | * This library is free software; you can redistribute it and/or
|
---|
| 8 | * modify it under the terms of the GNU Lesser General Public
|
---|
| 9 | * License as published by the Free Software Foundation;
|
---|
| 10 | * version 2.1 of the License.
|
---|
| 11 | *
|
---|
| 12 | * This library is distributed in the hope that it will be useful,
|
---|
| 13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
---|
| 15 | * Lesser General Public License for more details.
|
---|
| 16 | *
|
---|
| 17 | * This package contains formulas from the PROJ package of USGS.
|
---|
| 18 | * USGS's work is fully acknowledged here. This derived work has
|
---|
| 19 | * been relicensed under LGPL with Frank Warmerdam's permission.
|
---|
| 20 | */
|
---|
| 21 | package org.geotools.referencing.operation.projection;
|
---|
| 22 |
|
---|
| 23 | import java.awt.geom.Point2D;
|
---|
| 24 | import java.io.Serializable;
|
---|
| 25 | import java.util.Collection;
|
---|
| 26 | import java.util.logging.Level;
|
---|
| 27 | import java.util.logging.Logger;
|
---|
| 28 | import java.util.logging.LogRecord;
|
---|
| 29 |
|
---|
| 30 | import javax.measure.unit.NonSI;
|
---|
| 31 | import javax.measure.unit.SI;
|
---|
| 32 | import javax.measure.unit.Unit;
|
---|
| 33 |
|
---|
| 34 | import org.opengis.parameter.InvalidParameterValueException;
|
---|
| 35 | import org.opengis.parameter.GeneralParameterDescriptor;
|
---|
| 36 | import org.opengis.parameter.ParameterDescriptor;
|
---|
| 37 | import org.opengis.parameter.ParameterDescriptorGroup;
|
---|
| 38 | import org.opengis.parameter.ParameterNotFoundException;
|
---|
| 39 | import org.opengis.parameter.ParameterValueGroup;
|
---|
| 40 | import org.opengis.referencing.operation.MathTransform2D;
|
---|
| 41 | import org.opengis.referencing.operation.NoninvertibleTransformException;
|
---|
| 42 | import org.opengis.referencing.operation.Projection;
|
---|
| 43 | import org.opengis.referencing.operation.TransformException;
|
---|
| 44 |
|
---|
| 45 | import org.geotools.math.XMath;
|
---|
| 46 | import org.geotools.measure.Latitude;
|
---|
| 47 | import org.geotools.measure.Longitude;
|
---|
| 48 | import org.geotools.metadata.iso.citation.Citations;
|
---|
| 49 | import org.geotools.referencing.NamedIdentifier;
|
---|
| 50 | import org.geotools.referencing.operation.MathTransformProvider;
|
---|
| 51 | import org.geotools.referencing.operation.transform.AbstractMathTransform;
|
---|
| 52 | import org.geotools.resources.i18n.Errors;
|
---|
| 53 | import org.geotools.resources.i18n.ErrorKeys;
|
---|
| 54 | import org.geotools.util.Utilities;
|
---|
| 55 | import org.geotools.util.logging.Logging;
|
---|
| 56 |
|
---|
| 57 | import static java.lang.Math.*;
|
---|
| 58 |
|
---|
| 59 |
|
---|
| 60 | /**
|
---|
| 61 | * Base class for transformation services between ellipsoidal and cartographic projections.
|
---|
| 62 | * This base class provides the basic feature needed for all methods (no need to overrides
|
---|
| 63 | * methods). Subclasses must "only" implements the following methods:
|
---|
| 64 | * <ul>
|
---|
| 65 | * <li>{@link #getParameterValues}</li>
|
---|
| 66 | * <li>{@link #transformNormalized}</li>
|
---|
| 67 | * <li>{@link #inverseTransformNormalized}</li>
|
---|
| 68 | * </ul>
|
---|
| 69 | * <p>
|
---|
| 70 | * <strong>NOTE:</strong>Serialization of this class is appropriate for short-term storage
|
---|
| 71 | * or RMI use, but will probably not be compatible with future version. For long term storage,
|
---|
| 72 | * WKT (Well Know Text) or XML (not yet implemented) are more appropriate.
|
---|
| 73 | *
|
---|
| 74 | * @since 2.0
|
---|
| 75 | * @version $Id: MapProjection.java 37299 2011-05-25 05:21:24Z mbedward $
|
---|
| 76 | *
|
---|
| 77 | * @source $URL: http://svn.osgeo.org/geotools/branches/2.7.x/modules/library/referencing/src/main/java/org/geotools/referencing/operation/projection/MapProjection.java $
|
---|
| 78 | * @author André Gosselin
|
---|
| 79 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
| 80 | * @author Rueben Schulz
|
---|
| 81 | *
|
---|
| 82 | * @see <A HREF="http://mathworld.wolfram.com/MapProjection.html">Map projections on MathWorld</A>
|
---|
| 83 | * @see <A HREF="http://atlas.gc.ca/site/english/learningresources/carto_corner/map_projections.html">Map projections on the atlas of Canada</A>
|
---|
| 84 | * @tutorial http://www.geotools.org/display/GEOTOOLS/How+to+add+new+projections
|
---|
| 85 | */
|
---|
| 86 | public abstract class MapProjection extends AbstractMathTransform
|
---|
| 87 | implements MathTransform2D, Serializable
|
---|
| 88 | {
|
---|
| 89 | /**
|
---|
| 90 | * For cross-version compatibility.
|
---|
| 91 | */
|
---|
| 92 | private static final long serialVersionUID = -406751619777246914L;
|
---|
| 93 |
|
---|
| 94 | /**
|
---|
| 95 | * The projection package logger
|
---|
| 96 | */
|
---|
| 97 | protected static final Logger LOGGER = Logging.getLogger(MapProjection.class);
|
---|
| 98 |
|
---|
| 99 | /**
|
---|
| 100 | * Maximum difference allowed when comparing real numbers.
|
---|
| 101 | * This field is private because subclasses may use different threshold value.
|
---|
| 102 | */
|
---|
| 103 | private static final double EPSILON = 1E-6;
|
---|
| 104 |
|
---|
| 105 | /**
|
---|
| 106 | * Maximum difference allowed when comparing longitudes or latitudes in degrees.
|
---|
| 107 | * This tolerance do not apply to angle in radians. A tolerance of 1E-4 is about
|
---|
| 108 | * 10 kilometers.
|
---|
| 109 | */
|
---|
| 110 | private static final double ANGLE_TOLERANCE = 1E-4;
|
---|
| 111 |
|
---|
| 112 | /**
|
---|
| 113 | * Difference allowed in iterative computations.
|
---|
| 114 | */
|
---|
| 115 | private static final double ITERATION_TOLERANCE = 1E-10;
|
---|
| 116 |
|
---|
| 117 | /**
|
---|
| 118 | * Relative iteration precision used in the <code>mlfn<code> method
|
---|
| 119 | */
|
---|
| 120 | private static final double MLFN_TOL = 1E-11;
|
---|
| 121 |
|
---|
| 122 |
|
---|
| 123 | /**
|
---|
| 124 | * Maximum number of iterations for iterative computations.
|
---|
| 125 | */
|
---|
| 126 | private static final int MAXIMUM_ITERATIONS = 15;
|
---|
| 127 |
|
---|
| 128 | /**
|
---|
| 129 | * Constants used to calculate {@link #en0}, {@link #en1},
|
---|
| 130 | * {@link #en2}, {@link #en3}, {@link #en4}.
|
---|
| 131 | */
|
---|
| 132 | private static final double C00= 1.0,
|
---|
| 133 | C02= 0.25,
|
---|
| 134 | C04= 0.046875,
|
---|
| 135 | C06= 0.01953125,
|
---|
| 136 | C08= 0.01068115234375,
|
---|
| 137 | C22= 0.75,
|
---|
| 138 | C44= 0.46875,
|
---|
| 139 | C46= 0.01302083333333333333,
|
---|
| 140 | C48= 0.00712076822916666666,
|
---|
| 141 | C66= 0.36458333333333333333,
|
---|
| 142 | C68= 0.00569661458333333333,
|
---|
| 143 | C88= 0.3076171875;
|
---|
| 144 |
|
---|
| 145 | /**
|
---|
| 146 | * Ellipsoid excentricity, equals to <code>sqrt({@link #excentricitySquared})</code>.
|
---|
| 147 | * Value 0 means that the ellipsoid is spherical.
|
---|
| 148 | *
|
---|
| 149 | * @see #excentricitySquared
|
---|
| 150 | * @see #isSpherical
|
---|
| 151 | */
|
---|
| 152 | protected final double excentricity;
|
---|
| 153 |
|
---|
| 154 | /**
|
---|
| 155 | * The square of excentricity: e² = (a²-b²)/a² where
|
---|
| 156 | * <var>e</var> is the {@linkplain #excentricity excentricity},
|
---|
| 157 | * <var>a</var> is the {@linkplain #semiMajor semi major} axis length and
|
---|
| 158 | * <var>b</var> is the {@linkplain #semiMinor semi minor} axis length.
|
---|
| 159 | *
|
---|
| 160 | * @see #excentricity
|
---|
| 161 | * @see #semiMajor
|
---|
| 162 | * @see #semiMinor
|
---|
| 163 | * @see #isSpherical
|
---|
| 164 | */
|
---|
| 165 | protected final double excentricitySquared;
|
---|
| 166 |
|
---|
| 167 | /**
|
---|
| 168 | * {@code true} if this projection is spherical. Spherical model has identical
|
---|
| 169 | * {@linkplain #semiMajor semi major} and {@linkplain #semiMinor semi minor} axis
|
---|
| 170 | * length, and an {@linkplain #excentricity excentricity} zero.
|
---|
| 171 | *
|
---|
| 172 | * @see #excentricity
|
---|
| 173 | * @see #semiMajor
|
---|
| 174 | * @see #semiMinor
|
---|
| 175 | */
|
---|
| 176 | protected final boolean isSpherical;
|
---|
| 177 |
|
---|
| 178 | /**
|
---|
| 179 | * Length of semi-major axis, in metres. This is named '<var>a</var>' or '<var>R</var>'
|
---|
| 180 | * (Radius in spherical cases) in Snyder.
|
---|
| 181 | *
|
---|
| 182 | * @see #excentricity
|
---|
| 183 | * @see #semiMinor
|
---|
| 184 | */
|
---|
| 185 | protected final double semiMajor;
|
---|
| 186 |
|
---|
| 187 | /**
|
---|
| 188 | * Length of semi-minor axis, in metres. This is named '<var>b</var>' in Snyder.
|
---|
| 189 | *
|
---|
| 190 | * @see #excentricity
|
---|
| 191 | * @see #semiMajor
|
---|
| 192 | */
|
---|
| 193 | protected final double semiMinor;
|
---|
| 194 |
|
---|
| 195 | /**
|
---|
| 196 | * Central longitude in <u>radians</u>. Default value is 0, the Greenwich meridian.
|
---|
| 197 | * This is called '<var>lambda0</var>' in Snyder.
|
---|
| 198 | * <p>
|
---|
| 199 | * <strong>Consider this field as final</strong>. It is not final only
|
---|
| 200 | * because some classes need to modify it at construction time.
|
---|
| 201 | */
|
---|
| 202 | protected double centralMeridian;
|
---|
| 203 |
|
---|
| 204 | /**
|
---|
| 205 | * Latitude of origin in <u>radians</u>. Default value is 0, the equator.
|
---|
| 206 | * This is called '<var>phi0</var>' in Snyder.
|
---|
| 207 | * <p>
|
---|
| 208 | * <strong>Consider this field as final</strong>. It is not final only
|
---|
| 209 | * because some classes need to modify it at construction time.
|
---|
| 210 | */
|
---|
| 211 | protected double latitudeOfOrigin;
|
---|
| 212 |
|
---|
| 213 | /**
|
---|
| 214 | * The scale factor. Default value is 1. Named '<var>k</var>' in Snyder.
|
---|
| 215 | * <p>
|
---|
| 216 | * <strong>Consider this field as final</strong>. It is not final only
|
---|
| 217 | * because some classes need to modify it at construction time.
|
---|
| 218 | */
|
---|
| 219 | protected double scaleFactor;
|
---|
| 220 |
|
---|
| 221 | /**
|
---|
| 222 | * False easting, in metres. Default value is 0.
|
---|
| 223 | */
|
---|
| 224 | protected final double falseEasting;
|
---|
| 225 |
|
---|
| 226 | /**
|
---|
| 227 | * False northing, in metres. Default value is 0.
|
---|
| 228 | */
|
---|
| 229 | protected final double falseNorthing;
|
---|
| 230 |
|
---|
| 231 | /**
|
---|
| 232 | * Global scale factor. Default value {@code globalScale} is equal
|
---|
| 233 | * to {@link #semiMajor}×{@link #scaleFactor}.
|
---|
| 234 | * <p>
|
---|
| 235 | * <strong>Consider this field as final</strong>. It is not final only
|
---|
| 236 | * because some classes need to modify it at construction time.
|
---|
| 237 | */
|
---|
| 238 | protected double globalScale;
|
---|
| 239 |
|
---|
| 240 | /**
|
---|
| 241 | * The inverse of this map projection. Will be created only when needed.
|
---|
| 242 | */
|
---|
| 243 | private transient MathTransform2D inverse;
|
---|
| 244 |
|
---|
| 245 | /**
|
---|
| 246 | * Constant needed for the <code>mlfn<code> method.
|
---|
| 247 | * Setup at construction time.
|
---|
| 248 | */
|
---|
| 249 | protected double en0,en1,en2,en3,en4;
|
---|
| 250 |
|
---|
| 251 | /**
|
---|
| 252 | * When different than {@link #globalRangeCheckSemaphore}, coordinate ranges will be
|
---|
| 253 | * checked and a {@code WARNING} log will be issued if they are out of their natural
|
---|
| 254 | * ranges (-180/180° for longitude, -90/90° for latitude).
|
---|
| 255 | *
|
---|
| 256 | * @see #verifyCoordinateRanges()
|
---|
| 257 | * @see #warningLogged()
|
---|
| 258 | */
|
---|
| 259 | private transient int rangeCheckSemaphore;
|
---|
| 260 |
|
---|
| 261 | /**
|
---|
| 262 | * The value to be checked against {@link #rangeCheckSemaphore} in order to determine
|
---|
| 263 | * if coordinates ranges should be checked.
|
---|
| 264 | */
|
---|
| 265 | private static int globalRangeCheckSemaphore = 1;
|
---|
| 266 |
|
---|
| 267 | /**
|
---|
| 268 | * Marks if the projection is invertible. The vast majority is, subclasses can override.
|
---|
| 269 | */
|
---|
| 270 | protected boolean invertible = true;
|
---|
| 271 |
|
---|
| 272 | /**
|
---|
| 273 | * Constructs a new map projection from the suplied parameters.
|
---|
| 274 | *
|
---|
| 275 | * @param values The parameter values in standard units.
|
---|
| 276 | * The following parameter are recognized:
|
---|
| 277 | * <ul>
|
---|
| 278 | * <li>"semi_major" (mandatory: no default)</li>
|
---|
| 279 | * <li>"semi_minor" (mandatory: no default)</li>
|
---|
| 280 | * <li>"central_meridian" (default to 0°)</li>
|
---|
| 281 | * <li>"latitude_of_origin" (default to 0°)</li>
|
---|
| 282 | * <li>"scale_factor" (default to 1 )</li>
|
---|
| 283 | * <li>"false_easting" (default to 0 )</li>
|
---|
| 284 | * <li>"false_northing" (default to 0 )</li>
|
---|
| 285 | * </ul>
|
---|
| 286 | * @throws ParameterNotFoundException if a mandatory parameter is missing.
|
---|
| 287 | */
|
---|
| 288 | protected MapProjection(final ParameterValueGroup values) throws ParameterNotFoundException {
|
---|
| 289 | this(values, null);
|
---|
| 290 | }
|
---|
| 291 |
|
---|
| 292 | /**
|
---|
| 293 | * Constructor invoked by sub-classes when we can't rely on {@link #getParameterDescriptors}
|
---|
| 294 | * before the construction is completed. This is the case when the later method depends on
|
---|
| 295 | * the value of some class's attribute, which has not yet been set. An example is
|
---|
| 296 | * {@link ObliqueMercator#getParameterDescriptors}.
|
---|
| 297 | * <p>
|
---|
| 298 | * This method is not public because it is not a very elegant hack, and a work around exists.
|
---|
| 299 | * For example {@code ObliqueMercator} two-points case could be implemented by as a separated
|
---|
| 300 | * classes, in which case {@link #getParameterDescriptors} returns a constant and can be safely
|
---|
| 301 | * invoked in a constructor.
|
---|
| 302 | */
|
---|
| 303 | MapProjection(final ParameterValueGroup values, Collection<GeneralParameterDescriptor> expected)
|
---|
| 304 | throws ParameterNotFoundException
|
---|
| 305 | {
|
---|
| 306 | if (expected == null) {
|
---|
| 307 | expected = getParameterDescriptors().descriptors();
|
---|
| 308 | }
|
---|
| 309 | semiMajor = doubleValue(expected, AbstractProvider.SEMI_MAJOR, values);
|
---|
| 310 | semiMinor = doubleValue(expected, AbstractProvider.SEMI_MINOR, values);
|
---|
| 311 | centralMeridian = doubleValue(expected, AbstractProvider.CENTRAL_MERIDIAN, values);
|
---|
| 312 | latitudeOfOrigin = doubleValue(expected, AbstractProvider.LATITUDE_OF_ORIGIN, values);
|
---|
| 313 | scaleFactor = doubleValue(expected, AbstractProvider.SCALE_FACTOR, values);
|
---|
| 314 | falseEasting = doubleValue(expected, AbstractProvider.FALSE_EASTING, values);
|
---|
| 315 | falseNorthing = doubleValue(expected, AbstractProvider.FALSE_NORTHING, values);
|
---|
| 316 | isSpherical = (semiMajor == semiMinor);
|
---|
| 317 | excentricitySquared = 1.0 - (semiMinor * semiMinor) / (semiMajor * semiMajor);
|
---|
| 318 | excentricity = sqrt(excentricitySquared);
|
---|
| 319 | globalScale = scaleFactor * semiMajor;
|
---|
| 320 | ensureLongitudeInRange(AbstractProvider.CENTRAL_MERIDIAN, centralMeridian, true);
|
---|
| 321 | ensureLatitudeInRange (AbstractProvider.LATITUDE_OF_ORIGIN, latitudeOfOrigin, true);
|
---|
| 322 |
|
---|
| 323 | // Compute constants for the mlfn
|
---|
| 324 | double t;
|
---|
| 325 | en0 = C00 - excentricitySquared * (C02 + excentricitySquared *
|
---|
| 326 | (C04 + excentricitySquared * (C06 + excentricitySquared * C08)));
|
---|
| 327 | en1 = excentricitySquared * (C22 - excentricitySquared *
|
---|
| 328 | (C04 + excentricitySquared * (C06 + excentricitySquared * C08)));
|
---|
| 329 | en2 = (t = excentricitySquared * excentricitySquared) *
|
---|
| 330 | (C44 - excentricitySquared * (C46 + excentricitySquared * C48));
|
---|
| 331 | en3 = (t *= excentricitySquared) * (C66 - excentricitySquared * C68);
|
---|
| 332 | en4 = t * excentricitySquared * C88;
|
---|
| 333 | }
|
---|
| 334 |
|
---|
| 335 | /**
|
---|
| 336 | * Returns {@code true} if the specified parameter can apply to this map projection.
|
---|
| 337 | * The set of expected parameters must be supplied. The default implementation just
|
---|
| 338 | * invokes {@code expected.contains(param)}. Some subclasses will override this method
|
---|
| 339 | * in order to handle {@link ModifiedParameterDescriptor} in a special way.
|
---|
| 340 | *
|
---|
| 341 | * @see #doubleValue
|
---|
| 342 | * @see #set
|
---|
| 343 | */
|
---|
| 344 | boolean isExpectedParameter(final Collection<GeneralParameterDescriptor> expected,
|
---|
| 345 | final ParameterDescriptor param)
|
---|
| 346 | {
|
---|
| 347 | return expected.contains(param);
|
---|
| 348 | }
|
---|
| 349 |
|
---|
| 350 | /**
|
---|
| 351 | * Returns the parameter value for the specified operation parameter. Values are
|
---|
| 352 | * automatically converted into the standard units specified by the supplied
|
---|
| 353 | * {@code param} argument, except {@link NonSI#DEGREE_ANGLE degrees} which
|
---|
| 354 | * are converted to {@link SI#RADIAN radians}.
|
---|
| 355 | *
|
---|
| 356 | * @param expected The value returned by {@code getParameterDescriptors().descriptors()}.
|
---|
| 357 | * @param param The parameter to look for.
|
---|
| 358 | * @param group The parameter value group to search into.
|
---|
| 359 | * @return The requested parameter value, or {@code NaN} if {@code param} is
|
---|
| 360 | * {@linkplain MathTransformProvider#createOptionalDescriptor optional}
|
---|
| 361 | * and the user didn't provided any value.
|
---|
| 362 | * @throws ParameterNotFoundException if the parameter is not found.
|
---|
| 363 | *
|
---|
| 364 | * @see MathTransformProvider#doubleValue
|
---|
| 365 | */
|
---|
| 366 | final double doubleValue(final Collection<GeneralParameterDescriptor> expected,
|
---|
| 367 | final ParameterDescriptor param,
|
---|
| 368 | final ParameterValueGroup group)
|
---|
| 369 | throws ParameterNotFoundException
|
---|
| 370 | {
|
---|
| 371 | if (isExpectedParameter(expected, param)) {
|
---|
| 372 | /*
|
---|
| 373 | * Gets the value supplied by the user. The conversion from decimal
|
---|
| 374 | * degrees to radians (if needed) is performed by AbstractProvider.
|
---|
| 375 | */
|
---|
| 376 | return AbstractProvider.doubleValue(param, group);
|
---|
| 377 | }
|
---|
| 378 | /*
|
---|
| 379 | * The constructor asked for a parameter value that do not apply to the type of the
|
---|
| 380 | * projection to be created. Returns a default value common to all projection types,
|
---|
| 381 | * but this value should not be used in projection computations.
|
---|
| 382 | */
|
---|
| 383 | double v;
|
---|
| 384 | final Object value = param.getDefaultValue();
|
---|
| 385 | if (value instanceof Number) {
|
---|
| 386 | v = ((Number) value).doubleValue();
|
---|
| 387 | if (NonSI.DEGREE_ANGLE.equals(param.getUnit())) {
|
---|
| 388 | v = toRadians(v);
|
---|
| 389 | }
|
---|
| 390 | } else {
|
---|
| 391 | v = Double.NaN;
|
---|
| 392 | }
|
---|
| 393 | return v;
|
---|
| 394 | }
|
---|
| 395 |
|
---|
| 396 | /**
|
---|
| 397 | * Ensures that this projection has equals semi-major and semi-minor axis. This method is
|
---|
| 398 | * invoked by constructors of classes implementing only spherical formulas.
|
---|
| 399 | */
|
---|
| 400 | final void ensureSpherical() throws IllegalArgumentException {
|
---|
| 401 | if (!isSpherical) {
|
---|
| 402 | throw new IllegalArgumentException(Errors.format(ErrorKeys.ELLIPTICAL_NOT_SUPPORTED));
|
---|
| 403 | }
|
---|
| 404 | }
|
---|
| 405 |
|
---|
| 406 | /**
|
---|
| 407 | * Ensures that the latitude is within allowed limits (±π/2).
|
---|
| 408 | * This method is useful to check the validity of projection parameters,
|
---|
| 409 | * like {@link #latitudeOfOrigin}.
|
---|
| 410 | *
|
---|
| 411 | * @param y Latitude to check, in radians.
|
---|
| 412 | * @param edge {@code true} to accept latitudes of ±π/2.
|
---|
| 413 | * @throws IllegalArgumentException if the latitude is out of range.
|
---|
| 414 | */
|
---|
| 415 | static void ensureLatitudeInRange(final ParameterDescriptor name, double y, final boolean edge)
|
---|
| 416 | throws IllegalArgumentException
|
---|
| 417 | {
|
---|
| 418 | if (edge ? (y >= Latitude.MIN_VALUE * PI/180 && y <= Latitude.MAX_VALUE * PI/180) :
|
---|
| 419 | (y > Latitude.MIN_VALUE * PI/180 && y < Latitude.MAX_VALUE * PI/180))
|
---|
| 420 | {
|
---|
| 421 | return;
|
---|
| 422 | }
|
---|
| 423 | y = toDegrees(y);
|
---|
| 424 | throw new InvalidParameterValueException(Errors.format(ErrorKeys.LATITUDE_OUT_OF_RANGE_$1,
|
---|
| 425 | new Latitude(y)), name.getName().getCode(), y);
|
---|
| 426 | }
|
---|
| 427 |
|
---|
| 428 | /**
|
---|
| 429 | * Ensures that the longitue is within allowed limits (±π).
|
---|
| 430 | * This method is used to check the validity of projection parameters,
|
---|
| 431 | * like {@link #centralMeridian}.
|
---|
| 432 | *
|
---|
| 433 | * @param x Longitude to verify, in radians.
|
---|
| 434 | * @param edge {@code true} for accepting longitudes of ±π.
|
---|
| 435 | * @throws IllegalArgumentException if the longitude is out of range.
|
---|
| 436 | */
|
---|
| 437 | static void ensureLongitudeInRange(final ParameterDescriptor name, double x, final boolean edge)
|
---|
| 438 | throws IllegalArgumentException
|
---|
| 439 | {
|
---|
| 440 | if (edge ? (x >= Longitude.MIN_VALUE * PI/180 && x <= Longitude.MAX_VALUE * PI/180) :
|
---|
| 441 | (x > Longitude.MIN_VALUE * PI/180 && x < Longitude.MAX_VALUE * PI/180))
|
---|
| 442 | {
|
---|
| 443 | return;
|
---|
| 444 | }
|
---|
| 445 | x = toDegrees(x);
|
---|
| 446 | throw new InvalidParameterValueException(Errors.format(ErrorKeys.LONGITUDE_OUT_OF_RANGE_$1,
|
---|
| 447 | new Longitude(x)), name.getName().getCode(), x);
|
---|
| 448 | }
|
---|
| 449 |
|
---|
| 450 | /**
|
---|
| 451 | * Verifies if the given coordinates are in the range of geographic coordinates. If they are
|
---|
| 452 | * not, then this method logs a warning and returns {@code true}. Otherwise this method does
|
---|
| 453 | * nothing and returns {@code false}.
|
---|
| 454 | *
|
---|
| 455 | * @param tr The caller.
|
---|
| 456 | * @param x The longitude in decimal degrees.
|
---|
| 457 | * @param y The latitude in decimal degrees.
|
---|
| 458 | * @return {@code true} if the coordinates are not in the geographic range, in which case
|
---|
| 459 | * a warning has been logged.
|
---|
| 460 | */
|
---|
| 461 | private static boolean verifyGeographicRanges(final AbstractMathTransform tr,
|
---|
| 462 | final double x, final double y)
|
---|
| 463 | {
|
---|
| 464 | // Note: the following tests should not fails for NaN values.
|
---|
| 465 | final boolean xOut, yOut;
|
---|
| 466 | xOut = (x < (Longitude.MIN_VALUE - ANGLE_TOLERANCE) || x > (Longitude.MAX_VALUE + ANGLE_TOLERANCE));
|
---|
| 467 | yOut = (y < (Latitude .MIN_VALUE - ANGLE_TOLERANCE) || y > (Latitude .MAX_VALUE + ANGLE_TOLERANCE));
|
---|
| 468 | if (!xOut && !yOut) {
|
---|
| 469 | return false;
|
---|
| 470 | }
|
---|
| 471 | final String lineSeparator = System.getProperty("line.separator", "\n");
|
---|
| 472 | final StringBuilder buffer = new StringBuilder();
|
---|
| 473 | buffer.append(Errors.format(ErrorKeys.OUT_OF_PROJECTION_VALID_AREA_$1, tr.getName()));
|
---|
| 474 | if (xOut) {
|
---|
| 475 | buffer.append(lineSeparator);
|
---|
| 476 | buffer.append(Errors.format(ErrorKeys.LONGITUDE_OUT_OF_RANGE_$1, new Longitude(x)));
|
---|
| 477 | }
|
---|
| 478 | if (yOut) {
|
---|
| 479 | buffer.append(lineSeparator);
|
---|
| 480 | buffer.append(Errors.format(ErrorKeys.LATITUDE_OUT_OF_RANGE_$1, new Latitude(y)));
|
---|
| 481 | }
|
---|
| 482 | final LogRecord record = new LogRecord(Level.WARNING, buffer.toString());
|
---|
| 483 | final String classe;
|
---|
| 484 | if (tr instanceof Inverse) {
|
---|
| 485 | classe = ((Inverse) tr).inverse().getClass().getName() + ".Inverse";
|
---|
| 486 | } else {
|
---|
| 487 | classe = tr.getClass().getName();
|
---|
| 488 | }
|
---|
| 489 | record.setSourceClassName(classe);
|
---|
| 490 | record.setSourceMethodName("transform");
|
---|
| 491 | record.setLoggerName(LOGGER.getName());
|
---|
| 492 | LOGGER.log(record);
|
---|
| 493 | return true;
|
---|
| 494 | }
|
---|
| 495 |
|
---|
| 496 | /**
|
---|
| 497 | * Sets the value in a parameter group. This convenience method is used
|
---|
| 498 | * by subclasses for {@link #getParameterValues} implementation. Values
|
---|
| 499 | * are automatically converted from radians to decimal degrees if needed.
|
---|
| 500 | *
|
---|
| 501 | * @param expected The value returned by {@code getParameterDescriptors().descriptors()}.
|
---|
| 502 | * @param param One of the {@link AbstractProvider} constants.
|
---|
| 503 | * @param group The group in which to set the value.
|
---|
| 504 | * @param value The value to set.
|
---|
| 505 | */
|
---|
| 506 | final void set(final Collection<GeneralParameterDescriptor> expected,
|
---|
| 507 | final ParameterDescriptor<?> param,
|
---|
| 508 | final ParameterValueGroup group,
|
---|
| 509 | double value)
|
---|
| 510 | {
|
---|
| 511 | if (isExpectedParameter(expected, param)) {
|
---|
| 512 | if (NonSI.DEGREE_ANGLE.equals(param.getUnit())) {
|
---|
| 513 | /*
|
---|
| 514 | * Converts radians to degrees and try to fix rounding error
|
---|
| 515 | * (e.g. -61.500000000000014 --> -61.5). This is necessary
|
---|
| 516 | * in order to avoid a bias when formatting a transform and
|
---|
| 517 | * parsing it again.
|
---|
| 518 | */
|
---|
| 519 | value = toDegrees(value);
|
---|
| 520 | double old = value;
|
---|
| 521 | value = XMath.trimDecimalFractionDigits(value, 4, 12);
|
---|
| 522 | if (value == old) {
|
---|
| 523 | /*
|
---|
| 524 | * The attempt to fix rounding error failed. Try again with the
|
---|
| 525 | * assumption that the true value is a multiple of 1/3 of angle
|
---|
| 526 | * (e.g. 51.166666666666664 --> 51.166666666666666), which is
|
---|
| 527 | * common in the EPSG database.
|
---|
| 528 | */
|
---|
| 529 | old *= 3;
|
---|
| 530 | final double test = XMath.trimDecimalFractionDigits(old, 4, 12);
|
---|
| 531 | if (test != old) {
|
---|
| 532 | value = test/3;
|
---|
| 533 | }
|
---|
| 534 | }
|
---|
| 535 | }
|
---|
| 536 | group.parameter(param.getName().getCode()).setValue(value);
|
---|
| 537 | }
|
---|
| 538 | }
|
---|
| 539 |
|
---|
| 540 | /**
|
---|
| 541 | * Returns the parameter descriptors for this map projection.
|
---|
| 542 | * This is used for a providing a default implementation of
|
---|
| 543 | * {@link #getParameterValues}, as well as arguments checking.
|
---|
| 544 | */
|
---|
| 545 | @Override
|
---|
| 546 | public abstract ParameterDescriptorGroup getParameterDescriptors();
|
---|
| 547 |
|
---|
| 548 | /**
|
---|
| 549 | * Returns the parameter values for this map projection.
|
---|
| 550 | *
|
---|
| 551 | * @return A copy of the parameter values for this map projection.
|
---|
| 552 | */
|
---|
| 553 | @Override
|
---|
| 554 | public ParameterValueGroup getParameterValues() {
|
---|
| 555 | final ParameterDescriptorGroup descriptor = getParameterDescriptors();
|
---|
| 556 | final Collection<GeneralParameterDescriptor> expected = descriptor.descriptors();
|
---|
| 557 | final ParameterValueGroup values = descriptor.createValue();
|
---|
| 558 | set(expected, AbstractProvider.SEMI_MAJOR, values, semiMajor );
|
---|
| 559 | set(expected, AbstractProvider.SEMI_MINOR, values, semiMinor );
|
---|
| 560 | set(expected, AbstractProvider.CENTRAL_MERIDIAN, values, centralMeridian );
|
---|
| 561 | set(expected, AbstractProvider.LATITUDE_OF_ORIGIN, values, latitudeOfOrigin);
|
---|
| 562 | set(expected, AbstractProvider.SCALE_FACTOR, values, scaleFactor );
|
---|
| 563 | set(expected, AbstractProvider.FALSE_EASTING, values, falseEasting );
|
---|
| 564 | set(expected, AbstractProvider.FALSE_NORTHING, values, falseNorthing );
|
---|
| 565 | return values;
|
---|
| 566 | }
|
---|
| 567 |
|
---|
| 568 | /**
|
---|
| 569 | * Returns the dimension of input points.
|
---|
| 570 | */
|
---|
| 571 | public final int getSourceDimensions() {
|
---|
| 572 | return 2;
|
---|
| 573 | }
|
---|
| 574 |
|
---|
| 575 | /**
|
---|
| 576 | * Returns the dimension of output points.
|
---|
| 577 | */
|
---|
| 578 | public final int getTargetDimensions() {
|
---|
| 579 | return 2;
|
---|
| 580 | }
|
---|
| 581 |
|
---|
| 582 |
|
---|
| 583 |
|
---|
| 584 |
|
---|
| 585 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 586 | //////// ////////
|
---|
| 587 | //////// TRANSFORMATION METHODS ////////
|
---|
| 588 | //////// Includes an inner class for inverse projections. ////////
|
---|
| 589 | //////// ////////
|
---|
| 590 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 591 |
|
---|
| 592 | /**
|
---|
| 593 | * Returns the orthodromic distance between the two specified points using a spherical
|
---|
| 594 | * approximation. This is used for assertions only.
|
---|
| 595 | */
|
---|
| 596 | private double orthodromicDistance(final Point2D source, final Point2D target) {
|
---|
| 597 | final double y1 = toRadians(source.getY());
|
---|
| 598 | final double y2 = toRadians(target.getY());
|
---|
| 599 | final double dx = toRadians(abs(target.getX() - source.getX()) % 360);
|
---|
| 600 | double rho = sin(y1)*sin(y2) + cos(y1)*cos(y2)*cos(dx);
|
---|
| 601 | if (rho > +1) {assert rho <= +(1+EPSILON) : rho; rho = +1;}
|
---|
| 602 | if (rho < -1) {assert rho >= -(1+EPSILON) : rho; rho = -1;}
|
---|
| 603 | return acos(rho) * semiMajor;
|
---|
| 604 | }
|
---|
| 605 |
|
---|
| 606 | /**
|
---|
| 607 | * Check point for private use by {@link #checkReciprocal}. This class is necessary in order
|
---|
| 608 | * to avoid never-ending loop in {@code assert} statements (when an {@code assert}
|
---|
| 609 | * calls {@code transform(...)}, which calls {@code inverse.transform(...)}, which
|
---|
| 610 | * calls {@code transform(...)}, etc.).
|
---|
| 611 | */
|
---|
| 612 | @SuppressWarnings("serial")
|
---|
| 613 | private static final class CheckPoint extends Point2D.Double {
|
---|
| 614 | public CheckPoint(final Point2D point) {
|
---|
| 615 | super(point.getX(), point.getY());
|
---|
| 616 | }
|
---|
| 617 | }
|
---|
| 618 |
|
---|
| 619 | /**
|
---|
| 620 | * Check if the transform of {@code point} is close enough to {@code target}.
|
---|
| 621 | * "Close enough" means that the two points are separated by a distance shorter than
|
---|
| 622 | * {@link #getToleranceForAssertions}. This method is used for assertions with J2SE 1.4.
|
---|
| 623 | *
|
---|
| 624 | * @param point Point to transform, in decimal degrees if {@code inverse} is {@code false}.
|
---|
| 625 | * @param target Point to compare to, in metres if {@code inverse} is {@code false}.
|
---|
| 626 | * @param inverse {@code true} for an inverse transform instead of a direct one.
|
---|
| 627 | * @return {@code true} if the two points are close enough.
|
---|
| 628 | */
|
---|
| 629 | private boolean checkReciprocal(Point2D point, final Point2D target, final boolean inverse)
|
---|
| 630 | throws ProjectionException
|
---|
| 631 | {
|
---|
| 632 | if (!(point instanceof CheckPoint)) try {
|
---|
| 633 | point = new CheckPoint(point);
|
---|
| 634 | final double longitude;
|
---|
| 635 | final double latitude;
|
---|
| 636 | final double distance;
|
---|
| 637 | if (inverse) {
|
---|
| 638 | // Computes orthodromic distance (spherical model) in metres.
|
---|
| 639 | point = inverse().transform(point, point);
|
---|
| 640 | distance = orthodromicDistance(point, target);
|
---|
| 641 | longitude = point.getX();
|
---|
| 642 | latitude = point.getY();
|
---|
| 643 | } else {
|
---|
| 644 | // Computes cartesian distance in metres.
|
---|
| 645 | longitude = point.getX();
|
---|
| 646 | latitude = point.getY();
|
---|
| 647 | point = transform(point, point);
|
---|
| 648 | distance = point.distance(target);
|
---|
| 649 | }
|
---|
| 650 | if (distance > getToleranceForAssertions(longitude, latitude)) {
|
---|
| 651 | /*
|
---|
| 652 | * Do not fail for NaN values. For other cases we must throw a ProjectionException,
|
---|
| 653 | * not an AssertionError, because some code like CRS.transform(CoordinateOperation,
|
---|
| 654 | * ...) will project points that are know to be suspicious by surrounding them in
|
---|
| 655 | * "try ... catch" statements. Failure are normal in their case and we want to let
|
---|
| 656 | * them handle the exception the way they are used to.
|
---|
| 657 | */
|
---|
| 658 | throw new ProjectionException(Errors.format(ErrorKeys.PROJECTION_CHECK_FAILED_$4,
|
---|
| 659 | distance,
|
---|
| 660 | new Longitude(longitude - toDegrees(centralMeridian )),
|
---|
| 661 | new Latitude (latitude - toDegrees(latitudeOfOrigin)),
|
---|
| 662 | getName()));
|
---|
| 663 | }
|
---|
| 664 | } catch (ProjectionException exception) {
|
---|
| 665 | throw exception;
|
---|
| 666 | } catch (TransformException exception) {
|
---|
| 667 | throw new ProjectionException(exception);
|
---|
| 668 | }
|
---|
| 669 | return true;
|
---|
| 670 | }
|
---|
| 671 |
|
---|
| 672 | /**
|
---|
| 673 | * Checks if transform using spherical formulas produces the same result
|
---|
| 674 | * than ellipsoidal formulas. This method is invoked during assertions only.
|
---|
| 675 | *
|
---|
| 676 | * @param x The easting computed by spherical formulas, in metres.
|
---|
| 677 | * @param y The northing computed by spherical formulas, in metres.
|
---|
| 678 | * @param expected The (easting,northing) computed by ellipsoidal formulas.
|
---|
| 679 | * @param tolerance The tolerance (optional).
|
---|
| 680 | */
|
---|
| 681 | static boolean checkTransform(final double x, final double y,
|
---|
| 682 | final Point2D expected, final double tolerance)
|
---|
| 683 | {
|
---|
| 684 | compare("x", expected.getX(), x, tolerance);
|
---|
| 685 | compare("y", expected.getY(), y, tolerance);
|
---|
| 686 | return tolerance < Double.POSITIVE_INFINITY;
|
---|
| 687 | }
|
---|
| 688 |
|
---|
| 689 | /**
|
---|
[28055] | 690 | * Default version of {@link #checkTransform(double,double,Point2D,double)}.
|
---|
| 691 | */
|
---|
| 692 | static boolean checkTransform(final double x, final double y, final Point2D expected) {
|
---|
| 693 | return checkTransform(x, y, expected, EPSILON);
|
---|
| 694 | }
|
---|
| 695 |
|
---|
| 696 | /**
|
---|
[28000] | 697 | * Checks if inverse transform using spherical formulas produces the same result
|
---|
| 698 | * than ellipsoidal formulas. This method is invoked during assertions only.
|
---|
| 699 | * <p>
|
---|
| 700 | * <strong>Note:</strong> this method ignores the longitude if the latitude is
|
---|
| 701 | * at a pole, because in such case the longitude is meanless.
|
---|
| 702 | *
|
---|
| 703 | * @param longitude The longitude computed by spherical formulas, in radians.
|
---|
| 704 | * @param latitude The latitude computed by spherical formulas, in radians.
|
---|
| 705 | * @param expected The (longitude,latitude) computed by ellipsoidal formulas.
|
---|
| 706 | * @param tolerance The tolerance (optional).
|
---|
| 707 | */
|
---|
| 708 | static boolean checkInverseTransform(final double longitude, final double latitude,
|
---|
| 709 | final Point2D expected, final double tolerance)
|
---|
| 710 | {
|
---|
| 711 | compare("latitude", expected.getY(), latitude, tolerance);
|
---|
| 712 | if (abs(PI/2 - abs(latitude)) > EPSILON) {
|
---|
| 713 | compare("longitude", expected.getX(), longitude, tolerance);
|
---|
| 714 | }
|
---|
| 715 | return tolerance < Double.POSITIVE_INFINITY;
|
---|
| 716 | }
|
---|
| 717 |
|
---|
| 718 | /**
|
---|
[28055] | 719 | * Default version of {@link #checkInverseTransform(double,double,Point2D,double)}.
|
---|
| 720 | */
|
---|
| 721 | static boolean checkInverseTransform(double longitude, double latitude, Point2D expected) {
|
---|
| 722 | return checkInverseTransform(longitude, latitude, expected, EPSILON);
|
---|
| 723 | }
|
---|
| 724 |
|
---|
| 725 | /**
|
---|
[28000] | 726 | * Compares two value for equality up to some tolerance threshold. This is used during
|
---|
| 727 | * assertions only. The comparaison do not fails if at least one value to compare is
|
---|
| 728 | * {@link Double#NaN}.
|
---|
| 729 | * <p>
|
---|
| 730 | * <strong>Hack:</strong> if the {@code variable} name starts by lower-case {@code L}
|
---|
| 731 | * (as in "longitude" and "latitude"), then the value is assumed to be an angle in
|
---|
| 732 | * radians. This is used for formatting an error message, if needed.
|
---|
| 733 | */
|
---|
| 734 | private static void compare(String variable, double expected, double actual, double tolerance) {
|
---|
| 735 | if (abs(expected - actual) > tolerance) {
|
---|
| 736 | if (variable.charAt(0) == 'l') {
|
---|
| 737 | actual = toDegrees(actual);
|
---|
| 738 | expected = toDegrees(expected);
|
---|
| 739 | }
|
---|
| 740 | throw new AssertionError(Errors.format(ErrorKeys.TEST_FAILURE_$3, variable, expected, actual));
|
---|
| 741 | }
|
---|
| 742 | }
|
---|
| 743 |
|
---|
| 744 | /**
|
---|
| 745 | * Transforms the specified coordinate and stores the result in {@code ptDst}. This method
|
---|
| 746 | * returns longitude as <var>x</var> values in the range {@code [-PI..PI]} and latitude as
|
---|
| 747 | * <var>y</var> values in the range {@code [-PI/2..PI/2]}. It will be checked by the caller,
|
---|
| 748 | * so this method doesn't need to performs this check.
|
---|
| 749 | * <p>
|
---|
| 750 | * Input coordinates have the {@link #falseEasting} and {@link #falseNorthing} removed and are
|
---|
| 751 | * divided by {@link #globalScale} before this method is invoked. After this method is invoked,
|
---|
| 752 | * the {@link #centralMeridian} is added to the {@code x} results in {@code ptDst}. This means
|
---|
| 753 | * that projections that implement this method are performed on an ellipse (or sphere) with a
|
---|
| 754 | * semi-major axis of 1.
|
---|
| 755 | * <p>
|
---|
| 756 | * In <A HREF="http://www.remotesensing.org/proj/">PROJ.4</A>, the same standardization,
|
---|
| 757 | * described above, is handled by {@code pj_inv.c}. Therefore when porting projections
|
---|
| 758 | * from PROJ.4, the inverse transform equations can be used directly here with minimal
|
---|
| 759 | * change. In the equations of Snyder, {@link #falseEasting}, {@link #falseNorthing} and
|
---|
| 760 | * {@link #scaleFactor} are usually not given. When implementing these equations here, you
|
---|
| 761 | * will not need to add the {@link #centralMeridian} to the output longitude or remove the
|
---|
| 762 | * {@link #semiMajor} (<var>a</var> or <var>R</var>).
|
---|
| 763 | *
|
---|
| 764 | * @param x The easting of the coordinate, linear distance on a unit sphere or ellipse.
|
---|
| 765 | * @param y The northing of the coordinate, linear distance on a unit sphere or ellipse.
|
---|
| 766 | * @param ptDst the specified coordinate point that stores the result of transforming
|
---|
| 767 | * {@code ptSrc}, or {@code null}. Ordinates will be in <strong>radians</strong>.
|
---|
| 768 | * @return the coordinate point after transforming {@code x}, {@code y}
|
---|
| 769 | * and storing the result in {@code ptDst}.
|
---|
| 770 | * @throws ProjectionException if the point can't be transformed.
|
---|
| 771 | */
|
---|
| 772 | protected abstract Point2D inverseTransformNormalized(double x, double y, final Point2D ptDst)
|
---|
| 773 | throws ProjectionException;
|
---|
| 774 |
|
---|
| 775 | /**
|
---|
| 776 | * Transforms the specified coordinate and stores the result in {@code ptDst}. This method is
|
---|
| 777 | * usually (but <strong>not</strong> guaranteed) to be invoked with values of <var>x</var> in
|
---|
| 778 | * the range {@code [-PI..PI]} and values of <var>y</var> in the range {@code [-PI/2..PI/2]}.
|
---|
| 779 | * Values outside those ranges are accepted (sometime with a warning logged) on the assumption
|
---|
| 780 | * that most implementations use those values only in trigonometric functions like
|
---|
| 781 | * {@linkplain Math#sin sin} and {@linkplain Math#cos cos}.
|
---|
| 782 | * <p>
|
---|
| 783 | * Coordinates have the {@link #centralMeridian} removed from <var>lambda</var> before this
|
---|
| 784 | * method is invoked. After this method is invoked, the results in {@code ptDst} are multiplied
|
---|
| 785 | * by {@link #globalScale}, and the {@link #falseEasting} and {@link #falseNorthing} are added.
|
---|
| 786 | * This means that projections that implement this method are performed on an ellipse (or sphere)
|
---|
| 787 | * with a semi-major axis of 1.
|
---|
| 788 | * <p>
|
---|
| 789 | * In <A HREF="http://www.remotesensing.org/proj/">PROJ.4</A>, the same standardization,
|
---|
| 790 | * described above, is handled by {@code pj_fwd.c}. Therefore when porting projections
|
---|
| 791 | * from PROJ.4, the forward transform equations can be used directly here with minimal
|
---|
| 792 | * change. In the equations of Snyder, {@link #falseEasting}, {@link #falseNorthing} and
|
---|
| 793 | * {@link #scaleFactor} are usually not given. When implementing these equations here,
|
---|
| 794 | * you will not need to remove the {@link #centralMeridian} from <var>lambda</var> or apply
|
---|
| 795 | * the {@link #semiMajor} (<var>a</var> or <var>R</var>).
|
---|
| 796 | *
|
---|
| 797 | * @param lambda The longitude of the coordinate, in <strong>radians</strong>.
|
---|
| 798 | * @param phi The latitude of the coordinate, in <strong>radians</strong>.
|
---|
| 799 | * @param ptDst the specified coordinate point that stores the result of transforming
|
---|
| 800 | * {@code ptSrc}, or {@code null}. Ordinates will be in a
|
---|
| 801 | * dimensionless unit, as a linear distance on a unit sphere or ellipse.
|
---|
| 802 | * @return the coordinate point after transforming ({@code lambda}, {@code phi})
|
---|
| 803 | * and storing the result in {@code ptDst}.
|
---|
| 804 | * @throws ProjectionException if the point can't be transformed.
|
---|
| 805 | */
|
---|
| 806 | protected abstract Point2D transformNormalized(double lambda, double phi, final Point2D ptDst)
|
---|
| 807 | throws ProjectionException;
|
---|
| 808 |
|
---|
| 809 | /**
|
---|
| 810 | * Transforms the specified {@code ptSrc} and stores the result in {@code ptDst}.
|
---|
| 811 | * <p>
|
---|
| 812 | * This method standardizes the source {@code x} coordinate
|
---|
| 813 | * by removing the {@link #centralMeridian}, before invoking
|
---|
| 814 | * <code>{@link #transformNormalized transformNormalized}(x, y, ptDst)</code>.
|
---|
| 815 | * It also multiplies by {@link #globalScale} and adds the {@link #falseEasting} and
|
---|
| 816 | * {@link #falseNorthing} to the point returned by the {@code transformNormalized(...)} call.
|
---|
| 817 | *
|
---|
| 818 | * @param ptSrc the specified coordinate point to be transformed.
|
---|
| 819 | * Ordinates must be in decimal degrees.
|
---|
| 820 | * @param ptDst the specified coordinate point that stores the result of transforming
|
---|
| 821 | * {@code ptSrc}, or {@code null}. Ordinates will be in metres.
|
---|
| 822 | * @return the coordinate point after transforming {@code ptSrc} and storing
|
---|
| 823 | * the result in {@code ptDst}.
|
---|
| 824 | * @throws ProjectionException if the point can't be transformed.
|
---|
| 825 | */
|
---|
| 826 | @Override
|
---|
| 827 | public final Point2D transform(final Point2D ptSrc, Point2D ptDst) throws ProjectionException {
|
---|
| 828 | final double x = ptSrc.getX();
|
---|
| 829 | final double y = ptSrc.getY();
|
---|
| 830 | if (verifyCoordinateRanges()) {
|
---|
| 831 | if (verifyGeographicRanges(this, x, y)) {
|
---|
| 832 | warningLogged();
|
---|
| 833 | }
|
---|
| 834 | }
|
---|
| 835 | /*
|
---|
| 836 | * Makes sure that the longitude before conversion stay within +/- PI radians. As a
|
---|
| 837 | * special case, we do not check the range if no rotation were applied on the longitude.
|
---|
| 838 | * This is because the user may have a big area ranging from -180° to +180°. With the
|
---|
| 839 | * slight rounding errors related to map projections, the 180° longitude may be slightly
|
---|
| 840 | * over the limit. Rolling the longitude would changes its sign. For example a bounding
|
---|
| 841 | * box from 30° to +180° would become 30° to -180°, which is probably not what the user
|
---|
| 842 | * wanted.
|
---|
| 843 | */
|
---|
| 844 | ptDst = transformNormalized(centralMeridian != 0 ?
|
---|
| 845 | rollLongitude(toRadians(x) - centralMeridian) : toRadians(x), toRadians(y), ptDst);
|
---|
| 846 | ptDst.setLocation(globalScale*ptDst.getX() + falseEasting,
|
---|
| 847 | globalScale*ptDst.getY() + falseNorthing);
|
---|
| 848 |
|
---|
| 849 | if(invertible) {
|
---|
| 850 | assert checkReciprocal(ptDst, (ptSrc!=ptDst) ? ptSrc : new Point2D.Double(x,y), true);
|
---|
| 851 | }
|
---|
| 852 | return ptDst;
|
---|
| 853 | }
|
---|
| 854 |
|
---|
| 855 | /**
|
---|
| 856 | * Transforms a list of coordinate point ordinal values. Ordinates must be
|
---|
| 857 | * (<var>longitude</var>,<var>latitude</var>) pairs in decimal degrees.
|
---|
| 858 | *
|
---|
| 859 | * @throws ProjectionException if a point can't be transformed. This method tries to transform
|
---|
| 860 | * every points even if some of them can't be transformed. Non-transformable points will
|
---|
| 861 | * have value {@link Double#NaN}. If more than one point can't be transformed, then this
|
---|
| 862 | * exception may be about an arbitrary point.
|
---|
| 863 | */
|
---|
| 864 | public final void transform(final double[] srcPts, int srcOff,
|
---|
| 865 | final double[] dstPts, int dstOff, int numPts)
|
---|
| 866 | throws ProjectionException
|
---|
| 867 | {
|
---|
| 868 | /*
|
---|
| 869 | * Vérifie s'il faudra parcourir le tableau en sens inverse.
|
---|
| 870 | * Ce sera le cas si les tableaux source et destination se
|
---|
| 871 | * chevauchent et que la destination est après la source.
|
---|
| 872 | */
|
---|
| 873 | final boolean reverse = (srcPts == dstPts && srcOff < dstOff &&
|
---|
| 874 | srcOff + (2*numPts) > dstOff);
|
---|
| 875 | if (reverse) {
|
---|
| 876 | srcOff += 2*numPts;
|
---|
| 877 | dstOff += 2*numPts;
|
---|
| 878 | }
|
---|
| 879 | final Point2D.Double point = new Point2D.Double();
|
---|
| 880 | ProjectionException firstException = null;
|
---|
| 881 | while (--numPts >= 0) {
|
---|
| 882 | try {
|
---|
| 883 | point.x = srcPts[srcOff++];
|
---|
| 884 | point.y = srcPts[srcOff++];
|
---|
| 885 | transform(point, point);
|
---|
| 886 | dstPts[dstOff++] = point.x;
|
---|
| 887 | dstPts[dstOff++] = point.y;
|
---|
| 888 | } catch (ProjectionException exception) {
|
---|
| 889 | dstPts[dstOff++] = Double.NaN;
|
---|
| 890 | dstPts[dstOff++] = Double.NaN;
|
---|
| 891 | if (firstException == null) {
|
---|
| 892 | firstException = exception;
|
---|
| 893 | }
|
---|
| 894 | }
|
---|
| 895 | if (reverse) {
|
---|
| 896 | srcOff -= 4;
|
---|
| 897 | dstOff -= 4;
|
---|
| 898 | }
|
---|
| 899 | }
|
---|
| 900 | if (firstException != null) {
|
---|
| 901 | throw firstException;
|
---|
| 902 | }
|
---|
| 903 | }
|
---|
| 904 |
|
---|
| 905 | /**
|
---|
| 906 | * Transforms a list of coordinate point ordinal values. Ordinates must be
|
---|
| 907 | * (<var>longitude</var>,<var>latitude</var>) pairs in decimal degrees.
|
---|
| 908 | *
|
---|
| 909 | * @throws ProjectionException if a point can't be transformed. This method tries to transform
|
---|
| 910 | * every points even if some of them can't be transformed. Non-transformable points will
|
---|
| 911 | * have value {@link Float#NaN}. If more than one point can't be transformed, then this
|
---|
| 912 | * exception may be about an arbitrary point.
|
---|
| 913 | */
|
---|
| 914 | @Override
|
---|
| 915 | public final void transform(final float[] srcPts, int srcOff,
|
---|
| 916 | final float[] dstPts, int dstOff, int numPts)
|
---|
| 917 | throws ProjectionException
|
---|
| 918 | {
|
---|
| 919 | final boolean reverse = (srcPts == dstPts && srcOff < dstOff &&
|
---|
| 920 | srcOff + (2*numPts) > dstOff);
|
---|
| 921 | if (reverse) {
|
---|
| 922 | srcOff += 2*numPts;
|
---|
| 923 | dstOff += 2*numPts;
|
---|
| 924 | }
|
---|
| 925 | final Point2D.Double point = new Point2D.Double();
|
---|
| 926 | ProjectionException firstException = null;
|
---|
| 927 | while (--numPts >= 0) {
|
---|
| 928 | try {
|
---|
| 929 | point.x = srcPts[srcOff++];
|
---|
| 930 | point.y = srcPts[srcOff++];
|
---|
| 931 | transform(point, point);
|
---|
| 932 | dstPts[dstOff++] = (float) point.x;
|
---|
| 933 | dstPts[dstOff++] = (float) point.y;
|
---|
| 934 | } catch (ProjectionException exception) {
|
---|
| 935 | dstPts[dstOff++] = Float.NaN;
|
---|
| 936 | dstPts[dstOff++] = Float.NaN;
|
---|
| 937 | if (firstException == null) {
|
---|
| 938 | firstException = exception;
|
---|
| 939 | }
|
---|
| 940 | }
|
---|
| 941 | if (reverse) {
|
---|
| 942 | srcOff -= 4;
|
---|
| 943 | dstOff -= 4;
|
---|
| 944 | }
|
---|
| 945 | }
|
---|
| 946 | if (firstException != null) {
|
---|
| 947 | throw firstException;
|
---|
| 948 | }
|
---|
| 949 | }
|
---|
| 950 |
|
---|
| 951 | /**
|
---|
| 952 | * Inverse of a map projection. Will be created by {@link MapProjection#inverse()} only when
|
---|
| 953 | * first required. Implementation of {@code transform(...)} methods are mostly identical
|
---|
| 954 | * to {@code MapProjection.transform(...)}, except that they will invokes
|
---|
| 955 | * {@link MapProjection#inverseTransformNormalized} instead of
|
---|
| 956 | * {@link MapProjection#transformNormalized}.
|
---|
| 957 | *
|
---|
| 958 | * @version $Id: MapProjection.java 37299 2011-05-25 05:21:24Z mbedward $
|
---|
| 959 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
| 960 | */
|
---|
| 961 | private final class Inverse extends AbstractMathTransform.Inverse implements MathTransform2D {
|
---|
| 962 | /**
|
---|
| 963 | * For cross-version compatibility.
|
---|
| 964 | */
|
---|
| 965 | private static final long serialVersionUID = -9138242780765956870L;
|
---|
| 966 |
|
---|
| 967 | /**
|
---|
| 968 | * Default constructor.
|
---|
| 969 | */
|
---|
| 970 | public Inverse() {
|
---|
| 971 | MapProjection.this.super();
|
---|
| 972 | }
|
---|
| 973 |
|
---|
| 974 | /**
|
---|
| 975 | * Inverse transforms the specified {@code ptSrc} and stores the result in {@code ptDst}.
|
---|
| 976 | * <p>
|
---|
| 977 | * This method standardizes the {@code ptSrc} by removing the {@link #falseEasting}
|
---|
| 978 | * and {@link #falseNorthing} and dividing by {@link #globalScale} before invoking
|
---|
| 979 | * <code>{@link #inverseTransformNormalized inverseTransformNormalized}(x, y, ptDst)</code>.
|
---|
| 980 | * It then adds the {@link #centralMeridian} to the {@code x} of the point returned by the
|
---|
| 981 | * {@code inverseTransformNormalized} call.
|
---|
| 982 | *
|
---|
| 983 | * @param ptSrc the specified coordinate point to be transformed.
|
---|
| 984 | * Ordinates must be in metres.
|
---|
| 985 | * @param ptDst the specified coordinate point that stores the result of transforming
|
---|
| 986 | * {@code ptSrc}, or {@code null}. Ordinates will be in decimal degrees.
|
---|
| 987 | * @return the coordinate point after transforming {@code ptSrc}
|
---|
| 988 | * and stroring the result in {@code ptDst}.
|
---|
| 989 | * @throws ProjectionException if the point can't be transformed.
|
---|
| 990 | */
|
---|
| 991 | @Override
|
---|
| 992 | public final Point2D transform(final Point2D ptSrc, Point2D ptDst)
|
---|
| 993 | throws ProjectionException
|
---|
| 994 | {
|
---|
| 995 | final double x0 = ptSrc.getX();
|
---|
| 996 | final double y0 = ptSrc.getY();
|
---|
| 997 | ptDst = inverseTransformNormalized((x0 - falseEasting ) / globalScale,
|
---|
| 998 | (y0 - falseNorthing) / globalScale, ptDst);
|
---|
| 999 | /*
|
---|
| 1000 | * Makes sure that the longitude after conversion stay within +/- PI radians. As a
|
---|
| 1001 | * special case, we do not check the range if no rotation were applied on the longitude.
|
---|
| 1002 | * This is because the user may have a big area ranging from -180° to +180°. With the
|
---|
| 1003 | * slight rounding errors related to map projections, the 180° longitude may be slightly
|
---|
| 1004 | * over the limit. Rolling the longitude would changes its sign. For example a bounding
|
---|
| 1005 | * box from 30° to +180° would become 30° to -180°, which is probably not what the user
|
---|
| 1006 | * wanted.
|
---|
| 1007 | */
|
---|
| 1008 | final double x = toDegrees(centralMeridian != 0 ?
|
---|
| 1009 | rollLongitude(ptDst.getX() + centralMeridian) : ptDst.getX());
|
---|
| 1010 | final double y = toDegrees(ptDst.getY());
|
---|
| 1011 | ptDst.setLocation(x,y);
|
---|
| 1012 | if (verifyCoordinateRanges()) {
|
---|
| 1013 | if (verifyGeographicRanges(this, x, y)) {
|
---|
| 1014 | warningLogged();
|
---|
| 1015 | }
|
---|
| 1016 | }
|
---|
| 1017 | assert checkReciprocal(ptDst, (ptSrc!=ptDst) ? ptSrc : new Point2D.Double(x0, y0), false);
|
---|
| 1018 | return ptDst;
|
---|
| 1019 | }
|
---|
| 1020 |
|
---|
| 1021 | /**
|
---|
| 1022 | * Inverse transforms a list of coordinate point ordinal values.
|
---|
| 1023 | * Ordinates must be (<var>x</var>,<var>y</var>) pairs in metres.
|
---|
| 1024 | *
|
---|
| 1025 | * @throws ProjectionException if a point can't be transformed. This method tries
|
---|
| 1026 | * to transform every points even if some of them can't be transformed.
|
---|
| 1027 | * Non-transformable points will have value {@link Double#NaN}. If more
|
---|
| 1028 | * than one point can't be transformed, then this exception may be about
|
---|
| 1029 | * an arbitrary point.
|
---|
| 1030 | */
|
---|
| 1031 | public final void transform(final double[] src, int srcOffset,
|
---|
| 1032 | final double[] dest, int dstOffset, int numPts)
|
---|
| 1033 | throws TransformException
|
---|
| 1034 | {
|
---|
| 1035 | /*
|
---|
| 1036 | * Vérifie s'il faudra parcourir le tableau en sens inverse.
|
---|
| 1037 | * Ce sera le cas si les tableaux source et destination se
|
---|
| 1038 | * chevauchent et que la destination est après la source.
|
---|
| 1039 | */
|
---|
| 1040 | final boolean reverse = (src==dest && srcOffset<dstOffset &&
|
---|
| 1041 | srcOffset+(2*numPts) > dstOffset);
|
---|
| 1042 | if (reverse) {
|
---|
| 1043 | srcOffset += 2*numPts;
|
---|
| 1044 | dstOffset += 2*numPts;
|
---|
| 1045 | }
|
---|
| 1046 | final Point2D.Double point = new Point2D.Double();
|
---|
| 1047 | ProjectionException firstException = null;
|
---|
| 1048 | while (--numPts>=0) {
|
---|
| 1049 | try {
|
---|
| 1050 | point.x = src[srcOffset++];
|
---|
| 1051 | point.y = src[srcOffset++];
|
---|
| 1052 | transform(point, point);
|
---|
| 1053 | dest[dstOffset++] = point.x;
|
---|
| 1054 | dest[dstOffset++] = point.y;
|
---|
| 1055 | } catch (ProjectionException exception) {
|
---|
| 1056 | dest[dstOffset++] = Double.NaN;
|
---|
| 1057 | dest[dstOffset++] = Double.NaN;
|
---|
| 1058 | if (firstException == null) {
|
---|
| 1059 | firstException = exception;
|
---|
| 1060 | }
|
---|
| 1061 | }
|
---|
| 1062 | if (reverse) {
|
---|
| 1063 | srcOffset -= 4;
|
---|
| 1064 | dstOffset -= 4;
|
---|
| 1065 | }
|
---|
| 1066 | }
|
---|
| 1067 | if (firstException != null) {
|
---|
| 1068 | throw firstException;
|
---|
| 1069 | }
|
---|
| 1070 | }
|
---|
| 1071 |
|
---|
| 1072 | /**
|
---|
| 1073 | * Inverse transforms a list of coordinate point ordinal values.
|
---|
| 1074 | * Ordinates must be (<var>x</var>,<var>y</var>) pairs in metres.
|
---|
| 1075 | *
|
---|
| 1076 | * @throws ProjectionException if a point can't be transformed. This method tries
|
---|
| 1077 | * to transform every points even if some of them can't be transformed.
|
---|
| 1078 | * Non-transformable points will have value {@link Float#NaN}. If more
|
---|
| 1079 | * than one point can't be transformed, then this exception may be about
|
---|
| 1080 | * an arbitrary point.
|
---|
| 1081 | */
|
---|
| 1082 | @Override
|
---|
| 1083 | public final void transform(final float[] src, int srcOffset,
|
---|
| 1084 | final float[] dest, int dstOffset, int numPts)
|
---|
| 1085 | throws ProjectionException
|
---|
| 1086 | {
|
---|
| 1087 | final boolean reverse = (src==dest && srcOffset<dstOffset &&
|
---|
| 1088 | srcOffset+(2*numPts) > dstOffset);
|
---|
| 1089 | if (reverse) {
|
---|
| 1090 | srcOffset += 2*numPts;
|
---|
| 1091 | dstOffset += 2*numPts;
|
---|
| 1092 | }
|
---|
| 1093 | final Point2D.Double point = new Point2D.Double();
|
---|
| 1094 | ProjectionException firstException = null;
|
---|
| 1095 | while (--numPts>=0) {
|
---|
| 1096 | try {
|
---|
| 1097 | point.x = src[srcOffset++];
|
---|
| 1098 | point.y = src[srcOffset++];
|
---|
| 1099 | transform(point, point);
|
---|
| 1100 | dest[dstOffset++] = (float) point.x;
|
---|
| 1101 | dest[dstOffset++] = (float) point.y;
|
---|
| 1102 | } catch (ProjectionException exception) {
|
---|
| 1103 | dest[dstOffset++] = Float.NaN;
|
---|
| 1104 | dest[dstOffset++] = Float.NaN;
|
---|
| 1105 | if (firstException == null) {
|
---|
| 1106 | firstException = exception;
|
---|
| 1107 | }
|
---|
| 1108 | }
|
---|
| 1109 | if (reverse) {
|
---|
| 1110 | srcOffset -= 4;
|
---|
| 1111 | dstOffset -= 4;
|
---|
| 1112 | }
|
---|
| 1113 | }
|
---|
| 1114 | if (firstException!=null) {
|
---|
| 1115 | throw firstException;
|
---|
| 1116 | }
|
---|
| 1117 | }
|
---|
| 1118 |
|
---|
| 1119 | /**
|
---|
| 1120 | * Returns the original map projection.
|
---|
| 1121 | */
|
---|
| 1122 | @Override
|
---|
| 1123 | public MathTransform2D inverse() {
|
---|
| 1124 | return (MathTransform2D) super.inverse();
|
---|
| 1125 | }
|
---|
| 1126 | }
|
---|
| 1127 |
|
---|
| 1128 | /**
|
---|
| 1129 | * Returns the inverse of this map projection.
|
---|
| 1130 | */
|
---|
| 1131 | @Override
|
---|
| 1132 | public final MathTransform2D inverse() throws NoninvertibleTransformException {
|
---|
| 1133 | if(!invertible) {
|
---|
| 1134 | throw new NoninvertibleTransformException(Errors.format(ErrorKeys.NONINVERTIBLE_TRANSFORM));
|
---|
| 1135 | }
|
---|
| 1136 |
|
---|
| 1137 |
|
---|
| 1138 | // No synchronization. Not a big deal if this method is invoked in
|
---|
| 1139 | // the same time by two threads resulting in two instances created.
|
---|
| 1140 | if (inverse == null) {
|
---|
| 1141 | inverse = new Inverse();
|
---|
| 1142 | }
|
---|
| 1143 | return inverse;
|
---|
| 1144 | }
|
---|
| 1145 |
|
---|
| 1146 | /**
|
---|
| 1147 | * Maximal error (in metres) tolerated for assertions, if enabled. When assertions are enabled,
|
---|
| 1148 | * every direct projection is followed by an inverse projection, and the result is compared to
|
---|
| 1149 | * the original coordinate. If a distance greater than the tolerance level is found, then an
|
---|
| 1150 | * {@link ProjectionException} will be thrown. Subclasses should override this method if they
|
---|
| 1151 | * need to relax the tolerance level.
|
---|
| 1152 | *
|
---|
| 1153 | * @param longitude The longitude in decimal degrees.
|
---|
| 1154 | * @param latitude The latitude in decimal degrees.
|
---|
| 1155 | * @return The tolerance level for assertions, in meters.
|
---|
| 1156 | */
|
---|
| 1157 | protected double getToleranceForAssertions(final double longitude, final double latitude) {
|
---|
| 1158 | final double delta = abs(longitude - centralMeridian)/2 +
|
---|
| 1159 | abs(latitude - latitudeOfOrigin);
|
---|
| 1160 | if (delta > 40) {
|
---|
| 1161 | // When far from the valid area, use a larger tolerance.
|
---|
| 1162 | return 1;
|
---|
| 1163 | }
|
---|
| 1164 | // Be less strict when the point is near an edge.
|
---|
| 1165 | return (abs(longitude) > 179) || (abs(latitude) > 89) ? 1E-1 : 1E-5;
|
---|
| 1166 | }
|
---|
| 1167 |
|
---|
| 1168 | /**
|
---|
| 1169 | * When {@code true}, coordinate ranges will be checked, and a {@link Level#WARNING WARNING}
|
---|
| 1170 | * log will be issued if they are out of their natural ranges (-180/180° for longitude,
|
---|
| 1171 | * -90/90° for latitude).
|
---|
| 1172 | * <p>
|
---|
| 1173 | * To avoid excessive logging, this flag will be set to {@code false} after the first
|
---|
| 1174 | * coordinate failing the checks is found.
|
---|
| 1175 | */
|
---|
| 1176 | final boolean verifyCoordinateRanges() {
|
---|
| 1177 | /*
|
---|
| 1178 | * Do not synchronize - doing so would be a major bottleneck since this method will be
|
---|
| 1179 | * invoked thousands of time. The consequence is that a call to {@link #resetWarnings}
|
---|
| 1180 | * may not be reflected immediately in other threads, but the later is defined only on
|
---|
| 1181 | * a "best effort" basis.
|
---|
| 1182 | */
|
---|
| 1183 | return rangeCheckSemaphore != globalRangeCheckSemaphore;
|
---|
| 1184 | }
|
---|
| 1185 |
|
---|
| 1186 | /**
|
---|
| 1187 | * To be invoked after a warning in order to disable subsequent warnings.
|
---|
| 1188 | */
|
---|
| 1189 | final void warningLogged() {
|
---|
| 1190 | synchronized (MapProjection.class) {
|
---|
| 1191 | rangeCheckSemaphore = globalRangeCheckSemaphore;
|
---|
| 1192 | }
|
---|
| 1193 | }
|
---|
| 1194 |
|
---|
| 1195 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 1196 | //////// ////////
|
---|
| 1197 | //////// IMPLEMENTATION OF Object AND MathTransform2D STANDARD METHODS ////////
|
---|
| 1198 | //////// ////////
|
---|
| 1199 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 1200 |
|
---|
| 1201 | /**
|
---|
| 1202 | * Returns a hash value for this map projection.
|
---|
| 1203 | */
|
---|
| 1204 | @Override
|
---|
| 1205 | public int hashCode() {
|
---|
| 1206 | long code = Double.doubleToLongBits(semiMajor);
|
---|
| 1207 | code = code*37 + Double.doubleToLongBits(semiMinor);
|
---|
| 1208 | code = code*37 + Double.doubleToLongBits(centralMeridian);
|
---|
| 1209 | code = code*37 + Double.doubleToLongBits(latitudeOfOrigin);
|
---|
| 1210 | return (int) code ^ (int) (code >>> 32);
|
---|
| 1211 | }
|
---|
| 1212 |
|
---|
| 1213 | /**
|
---|
| 1214 | * Compares the specified object with this map projection for equality.
|
---|
| 1215 | */
|
---|
| 1216 | @Override
|
---|
| 1217 | public boolean equals(final Object object) {
|
---|
| 1218 | // Do not check 'object==this' here, since this
|
---|
| 1219 | // optimization is usually done in subclasses.
|
---|
| 1220 | if (super.equals(object)) {
|
---|
| 1221 | final MapProjection that = (MapProjection) object;
|
---|
| 1222 | return equals(this.semiMajor, that.semiMajor) &&
|
---|
| 1223 | equals(this.semiMinor, that.semiMinor) &&
|
---|
| 1224 | equals(this.centralMeridian, that.centralMeridian) &&
|
---|
| 1225 | equals(this.latitudeOfOrigin, that.latitudeOfOrigin) &&
|
---|
| 1226 | equals(this.scaleFactor, that.scaleFactor) &&
|
---|
| 1227 | equals(this.falseEasting, that.falseEasting) &&
|
---|
| 1228 | equals(this.falseNorthing, that.falseNorthing);
|
---|
| 1229 | }
|
---|
| 1230 | return false;
|
---|
| 1231 | }
|
---|
| 1232 |
|
---|
| 1233 | /**
|
---|
| 1234 | * Returns {@code true} if the two specified value are equals.
|
---|
| 1235 | * Two {@link Double#NaN NaN} values are considered equals.
|
---|
| 1236 | */
|
---|
| 1237 | static boolean equals(final double value1, final double value2) {
|
---|
| 1238 | return Utilities.equals(value1, value2);
|
---|
| 1239 | }
|
---|
| 1240 |
|
---|
| 1241 |
|
---|
| 1242 |
|
---|
| 1243 |
|
---|
| 1244 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 1245 | //////// ////////
|
---|
| 1246 | //////// FORMULAS FROM SNYDER ////////
|
---|
| 1247 | //////// ////////
|
---|
| 1248 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 1249 |
|
---|
| 1250 | /**
|
---|
| 1251 | * Iteratively solve equation (7-9) from Snyder.
|
---|
| 1252 | */
|
---|
| 1253 | final double cphi2(final double ts) throws ProjectionException {
|
---|
| 1254 | final double eccnth = 0.5 * excentricity;
|
---|
| 1255 | double phi = (PI/2) - 2.0 * atan(ts);
|
---|
| 1256 | for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
|
---|
| 1257 | final double con = excentricity * sin(phi);
|
---|
| 1258 | final double dphi = (PI/2) - 2.0*atan(ts * pow((1-con)/(1+con), eccnth)) - phi;
|
---|
| 1259 | phi += dphi;
|
---|
| 1260 | if (abs(dphi) <= ITERATION_TOLERANCE) {
|
---|
| 1261 | return phi;
|
---|
| 1262 | }
|
---|
| 1263 | }
|
---|
| 1264 | throw new ProjectionException(ErrorKeys.NO_CONVERGENCE);
|
---|
| 1265 | }
|
---|
| 1266 |
|
---|
| 1267 | /**
|
---|
| 1268 | * Computes function <code>f(s,c,e²) = c/sqrt(1 - s²×e²)</code> needed for the true scale
|
---|
| 1269 | * latitude (Snyder 14-15), where <var>s</var> and <var>c</var> are the sine and cosine of
|
---|
| 1270 | * the true scale latitude, and <var>e²</var> is the {@linkplain #excentricitySquared
|
---|
| 1271 | * eccentricity squared}.
|
---|
| 1272 | */
|
---|
| 1273 | final double msfn(final double s, final double c) {
|
---|
| 1274 | return c / sqrt(1.0 - (s*s) * excentricitySquared);
|
---|
| 1275 | }
|
---|
| 1276 |
|
---|
| 1277 | /**
|
---|
| 1278 | * Computes function (15-9) and (9-13) from Snyder.
|
---|
| 1279 | * Equivalent to negative of function (7-7).
|
---|
| 1280 | */
|
---|
| 1281 | final double tsfn(final double phi, double sinphi) {
|
---|
| 1282 | sinphi *= excentricity;
|
---|
| 1283 | /*
|
---|
| 1284 | * NOTE: change sign to get the equivalent of Snyder (7-7).
|
---|
| 1285 | */
|
---|
| 1286 | return tan(0.5 * (PI/2 - phi)) / pow((1 - sinphi) / (1 + sinphi), 0.5*excentricity);
|
---|
| 1287 | }
|
---|
| 1288 |
|
---|
| 1289 | /**
|
---|
| 1290 | * Calculates the meridian distance. This is the distance along the central
|
---|
| 1291 | * meridian from the equator to {@code phi}. Accurate to < 1e-5 meters
|
---|
| 1292 | * when used in conjuction with typical major axis values.
|
---|
| 1293 | *
|
---|
| 1294 | * @param phi latitude to calculate meridian distance for.
|
---|
| 1295 | * @param sphi sin(phi).
|
---|
| 1296 | * @param cphi cos(phi).
|
---|
| 1297 | * @return meridian distance for the given latitude.
|
---|
| 1298 | */
|
---|
| 1299 | protected final double mlfn(final double phi, double sphi, double cphi) {
|
---|
| 1300 | cphi *= sphi;
|
---|
| 1301 | sphi *= sphi;
|
---|
| 1302 | return en0 * phi - cphi *
|
---|
| 1303 | (en1 + sphi *
|
---|
| 1304 | (en2 + sphi *
|
---|
| 1305 | (en3 + sphi *
|
---|
| 1306 | (en4))));
|
---|
| 1307 | }
|
---|
| 1308 |
|
---|
| 1309 | /**
|
---|
| 1310 | * Calculates the latitude ({@code phi}) from a meridian distance.
|
---|
| 1311 | * Determines phi to TOL (1e-11) radians, about 1e-6 seconds.
|
---|
| 1312 | *
|
---|
| 1313 | * @param arg meridian distance to calulate latitude for.
|
---|
| 1314 | * @return the latitude of the meridian distance.
|
---|
| 1315 | * @throws ProjectionException if the itteration does not converge.
|
---|
| 1316 | */
|
---|
| 1317 | protected final double inv_mlfn(double arg) throws ProjectionException {
|
---|
| 1318 | double s, t, phi, k = 1.0/(1.0 - excentricitySquared);
|
---|
| 1319 | int i;
|
---|
| 1320 | phi = arg;
|
---|
| 1321 | for (i=MAXIMUM_ITERATIONS; true;) { // rarely goes over 5 iterations
|
---|
| 1322 | if (--i < 0) {
|
---|
| 1323 | throw new ProjectionException(Errors.format(ErrorKeys.NO_CONVERGENCE));
|
---|
| 1324 | }
|
---|
| 1325 | s = Math.sin(phi);
|
---|
| 1326 | t = 1.0 - excentricitySquared * s * s;
|
---|
| 1327 | t = (mlfn(phi, s, Math.cos(phi)) - arg) * (t * Math.sqrt(t)) * k;
|
---|
| 1328 | phi -= t;
|
---|
| 1329 | if (Math.abs(t) < MLFN_TOL) {
|
---|
| 1330 | return phi;
|
---|
| 1331 | }
|
---|
| 1332 | }
|
---|
| 1333 | }
|
---|
| 1334 |
|
---|
| 1335 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 1336 | //////// ////////
|
---|
| 1337 | //////// PROVIDER ////////
|
---|
| 1338 | //////// ////////
|
---|
| 1339 | //////////////////////////////////////////////////////////////////////////////////////////
|
---|
| 1340 |
|
---|
| 1341 | /**
|
---|
| 1342 | * The base provider for {@link MapProjection}s.
|
---|
| 1343 | *
|
---|
| 1344 | * @version $Id: MapProjection.java 37299 2011-05-25 05:21:24Z mbedward $
|
---|
| 1345 | * @author Martin Desruisseaux (PMO, IRD)
|
---|
| 1346 | */
|
---|
| 1347 | public static abstract class AbstractProvider extends MathTransformProvider {
|
---|
| 1348 | /**
|
---|
| 1349 | * Serial number for interoperability with different versions.
|
---|
| 1350 | */
|
---|
| 1351 | private static final long serialVersionUID = 6280666068007678702L;
|
---|
| 1352 |
|
---|
| 1353 | /**
|
---|
| 1354 | * The operation parameter descriptor for the {@linkplain #semiMajor semi major} parameter
|
---|
| 1355 | * value. Valid values range is from 0 to infinity. This parameter is mandatory.
|
---|
| 1356 | *
|
---|
| 1357 | * @todo Would like to start range from 0 <u>exclusive</u>.
|
---|
| 1358 | */
|
---|
| 1359 | public static final ParameterDescriptor SEMI_MAJOR = createDescriptor(
|
---|
| 1360 | new NamedIdentifier[] {
|
---|
| 1361 | new NamedIdentifier(Citations.OGC, "semi_major"),
|
---|
| 1362 | new NamedIdentifier(Citations.EPSG, "semi-major axis")
|
---|
| 1363 | // EPSG does not specifically define the above parameter
|
---|
| 1364 | },
|
---|
| 1365 | Double.NaN, 0, Double.POSITIVE_INFINITY, SI.METER);
|
---|
| 1366 |
|
---|
| 1367 | /**
|
---|
| 1368 | * The operation parameter descriptor for the {@linkplain #semiMinor semi minor} parameter
|
---|
| 1369 | * value. Valid values range is from 0 to infinity. This parameter is mandatory.
|
---|
| 1370 | *
|
---|
| 1371 | * @todo Would like to start range from 0 <u>exclusive</u>.
|
---|
| 1372 | */
|
---|
| 1373 | public static final ParameterDescriptor SEMI_MINOR = createDescriptor(
|
---|
| 1374 | new NamedIdentifier[] {
|
---|
| 1375 | new NamedIdentifier(Citations.OGC, "semi_minor"),
|
---|
| 1376 | new NamedIdentifier(Citations.EPSG, "semi-minor axis")
|
---|
| 1377 | // EPSG does not specifically define the above parameter
|
---|
| 1378 | },
|
---|
| 1379 | Double.NaN, 0, Double.POSITIVE_INFINITY, SI.METER);
|
---|
| 1380 |
|
---|
| 1381 | /**
|
---|
| 1382 | * The operation parameter descriptor for the {@linkplain #centralMeridian central meridian}
|
---|
| 1383 | * parameter value. Valid values range is from -180 to 180°. Default value is 0.
|
---|
| 1384 | */
|
---|
| 1385 | public static final ParameterDescriptor CENTRAL_MERIDIAN = createDescriptor(
|
---|
| 1386 | new NamedIdentifier[] {
|
---|
| 1387 | new NamedIdentifier(Citations.OGC, "central_meridian"),
|
---|
| 1388 | new NamedIdentifier(Citations.EPSG, "Longitude of natural origin"),
|
---|
| 1389 | new NamedIdentifier(Citations.EPSG, "Longitude of false origin"),
|
---|
| 1390 | new NamedIdentifier(Citations.EPSG, "Longitude of origin"),
|
---|
| 1391 | new NamedIdentifier(Citations.ESRI, "Longitude_Of_Center"),
|
---|
| 1392 | new NamedIdentifier(Citations.ESRI, "longitude_of_center"),
|
---|
| 1393 | new NamedIdentifier(Citations.ESRI, "Longitude_Of_Origin"),
|
---|
| 1394 | new NamedIdentifier(Citations.ESRI, "longitude_of_origin"),
|
---|
| 1395 | new NamedIdentifier(Citations.GEOTIFF, "NatOriginLong")
|
---|
| 1396 | // ESRI uses "Longitude_Of_Origin" in orthographic (not to
|
---|
| 1397 | // be confused with "Longitude_Of_Center" in oblique mercator).
|
---|
| 1398 | },
|
---|
| 1399 | 0, -180, 180, NonSI.DEGREE_ANGLE);
|
---|
| 1400 |
|
---|
| 1401 | /**
|
---|
| 1402 | * The operation parameter descriptor for the {@linkplain #latitudeOfOrigin latitude of origin}
|
---|
| 1403 | * parameter value. Valid values range is from -90 to 90°. Default value is 0.
|
---|
| 1404 | */
|
---|
| 1405 | public static final ParameterDescriptor LATITUDE_OF_ORIGIN = createDescriptor(
|
---|
| 1406 | new NamedIdentifier[] {
|
---|
| 1407 | new NamedIdentifier(Citations.OGC, "latitude_of_origin"),
|
---|
| 1408 | new NamedIdentifier(Citations.EPSG, "Latitude of false origin"),
|
---|
| 1409 | new NamedIdentifier(Citations.EPSG, "Latitude of natural origin"),
|
---|
| 1410 | new NamedIdentifier(Citations.ESRI, "Latitude_Of_Origin"),
|
---|
| 1411 | new NamedIdentifier(Citations.ESRI, "latitude_of_origin"),
|
---|
| 1412 | new NamedIdentifier(Citations.ESRI, "Latitude_Of_Center"),
|
---|
| 1413 | new NamedIdentifier(Citations.ESRI, "latitude_of_center"),
|
---|
| 1414 | new NamedIdentifier(Citations.GEOTIFF, "NatOriginLat")
|
---|
| 1415 | // ESRI uses "Latitude_Of_Center" in orthographic.
|
---|
| 1416 | },
|
---|
| 1417 | 0, -90, 90, NonSI.DEGREE_ANGLE);
|
---|
| 1418 |
|
---|
| 1419 | /**
|
---|
| 1420 | * The operation parameter descriptor for the standard parallel 1 parameter value.
|
---|
| 1421 | * Valid values range is from -90 to 90°. Default value is 0.
|
---|
| 1422 | */
|
---|
| 1423 | public static final ParameterDescriptor STANDARD_PARALLEL_1 = createDescriptor(
|
---|
| 1424 | new NamedIdentifier[] {
|
---|
| 1425 | new NamedIdentifier(Citations.OGC, "standard_parallel_1"),
|
---|
| 1426 | new NamedIdentifier(Citations.EPSG, "Latitude of 1st standard parallel"),
|
---|
| 1427 | new NamedIdentifier(Citations.ESRI, "Standard_Parallel_1"),
|
---|
| 1428 | new NamedIdentifier(Citations.ESRI, "standard_parallel_1"),
|
---|
| 1429 | new NamedIdentifier(Citations.GEOTIFF, "StdParallel1")
|
---|
| 1430 | },
|
---|
| 1431 | 0, -90, 90, NonSI.DEGREE_ANGLE);
|
---|
| 1432 |
|
---|
| 1433 | /**
|
---|
| 1434 | * The operation parameter descriptor for the standard parallel 2 parameter value.
|
---|
| 1435 | * Valid values range is from -90 to 90°. Default value is 0.
|
---|
| 1436 | */
|
---|
| 1437 | public static final ParameterDescriptor STANDARD_PARALLEL_2 = createOptionalDescriptor(
|
---|
| 1438 | new NamedIdentifier[] {
|
---|
| 1439 | new NamedIdentifier(Citations.OGC, "standard_parallel_2"),
|
---|
| 1440 | new NamedIdentifier(Citations.EPSG, "Latitude of 2nd standard parallel"),
|
---|
| 1441 | new NamedIdentifier(Citations.ESRI, "Standard_Parallel_2"),
|
---|
| 1442 | new NamedIdentifier(Citations.ESRI, "standard_parallel_2"),
|
---|
| 1443 | new NamedIdentifier(Citations.GEOTIFF, "StdParallel2")
|
---|
| 1444 | },
|
---|
| 1445 | -90, 90, NonSI.DEGREE_ANGLE);
|
---|
| 1446 |
|
---|
| 1447 | /**
|
---|
| 1448 | * The operation parameter descriptor for the {@link #scaleFactor scaleFactor}
|
---|
| 1449 | * parameter value. Valid values range is from 0 to infinity. Default value is 1.
|
---|
| 1450 | *
|
---|
| 1451 | * @todo Would like to start range from 0 <u>exclusive</u>.
|
---|
| 1452 | */
|
---|
| 1453 | public static final ParameterDescriptor SCALE_FACTOR = createDescriptor(
|
---|
| 1454 | new NamedIdentifier[] {
|
---|
| 1455 | new NamedIdentifier(Citations.OGC, "scale_factor"),
|
---|
| 1456 | new NamedIdentifier(Citations.EPSG, "Scale factor at natural origin"),
|
---|
| 1457 | new NamedIdentifier(Citations.EPSG, "Scale factor on initial line"),
|
---|
| 1458 | new NamedIdentifier(Citations.GEOTIFF, "ScaleAtNatOrigin"),
|
---|
| 1459 | new NamedIdentifier(Citations.GEOTIFF, "ScaleAtCenter"),
|
---|
| 1460 | new NamedIdentifier(Citations.ESRI, "Scale_Factor"),
|
---|
| 1461 | new NamedIdentifier(Citations.ESRI, "scale_factor"),
|
---|
| 1462 | },
|
---|
| 1463 | 1, 0, Double.POSITIVE_INFINITY, Unit.ONE);
|
---|
| 1464 |
|
---|
| 1465 | /**
|
---|
| 1466 | * The operation parameter descriptor for the {@link #falseEasting falseEasting}
|
---|
| 1467 | * parameter value. Valid values range is unrestricted. Default value is 0.
|
---|
| 1468 | */
|
---|
| 1469 | public static final ParameterDescriptor FALSE_EASTING = createDescriptor(
|
---|
| 1470 | new NamedIdentifier[] {
|
---|
| 1471 | new NamedIdentifier(Citations.OGC, "false_easting"),
|
---|
| 1472 | new NamedIdentifier(Citations.EPSG, "False easting"),
|
---|
| 1473 | new NamedIdentifier(Citations.EPSG, "Easting at false origin"),
|
---|
| 1474 | new NamedIdentifier(Citations.EPSG, "Easting at projection centre"),
|
---|
| 1475 | new NamedIdentifier(Citations.GEOTIFF, "FalseEasting"),
|
---|
| 1476 | new NamedIdentifier(Citations.ESRI, "false_easting")
|
---|
| 1477 | },
|
---|
| 1478 | 0, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, SI.METER);
|
---|
| 1479 |
|
---|
| 1480 | /**
|
---|
| 1481 | * The operation parameter descriptor for the {@link #falseNorthing falseNorthing}
|
---|
| 1482 | * parameter value. Valid values range is unrestricted. Default value is 0.
|
---|
| 1483 | */
|
---|
| 1484 | public static final ParameterDescriptor FALSE_NORTHING = createDescriptor(
|
---|
| 1485 | new NamedIdentifier[] {
|
---|
| 1486 | new NamedIdentifier(Citations.OGC, "false_northing"),
|
---|
| 1487 | new NamedIdentifier(Citations.EPSG, "False northing"),
|
---|
| 1488 | new NamedIdentifier(Citations.EPSG, "Northing at false origin"),
|
---|
| 1489 | new NamedIdentifier(Citations.EPSG, "Northing at projection centre"),
|
---|
| 1490 | new NamedIdentifier(Citations.GEOTIFF, "FalseNorthing"),
|
---|
| 1491 | new NamedIdentifier(Citations.ESRI, "False_Northing"),
|
---|
| 1492 | new NamedIdentifier(Citations.ESRI, "false_northing")
|
---|
| 1493 | },
|
---|
| 1494 | 0, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, SI.METER);
|
---|
| 1495 |
|
---|
| 1496 | /**
|
---|
| 1497 | * Constructs a math transform provider from a set of parameters. The provider
|
---|
| 1498 | * {@linkplain #getIdentifiers identifiers} will be the same than the parameter
|
---|
| 1499 | * ones.
|
---|
| 1500 | *
|
---|
| 1501 | * @param parameters The set of parameters (never {@code null}).
|
---|
| 1502 | */
|
---|
| 1503 | public AbstractProvider(final ParameterDescriptorGroup parameters) {
|
---|
| 1504 | super(2, 2, parameters);
|
---|
| 1505 | }
|
---|
| 1506 |
|
---|
| 1507 | /**
|
---|
| 1508 | * Returns the operation type for this map projection.
|
---|
| 1509 | */
|
---|
| 1510 | @Override
|
---|
| 1511 | public Class<? extends Projection> getOperationType() {
|
---|
| 1512 | return Projection.class;
|
---|
| 1513 | }
|
---|
| 1514 |
|
---|
| 1515 | /**
|
---|
| 1516 | * Returns {@code true} is the parameters use a spherical datum.
|
---|
| 1517 | */
|
---|
| 1518 | static boolean isSpherical(final ParameterValueGroup values) {
|
---|
| 1519 | try {
|
---|
| 1520 | return doubleValue(SEMI_MAJOR, values) == doubleValue(SEMI_MINOR, values);
|
---|
| 1521 | } catch (IllegalStateException exception) {
|
---|
| 1522 | // Probably could not find the requested values -- gobble error and be forgiving.
|
---|
| 1523 | // The error will probably be thrown at MapProjection construction time, which is
|
---|
| 1524 | // less surprising to some users.
|
---|
| 1525 | return false;
|
---|
| 1526 | }
|
---|
| 1527 | }
|
---|
| 1528 |
|
---|
| 1529 | /**
|
---|
| 1530 | * Returns the parameter value for the specified operation parameter in standard units.
|
---|
| 1531 | * Values are automatically converted into the standard units specified by the supplied
|
---|
| 1532 | * {@code param} argument, except {@link NonSI#DEGREE_ANGLE degrees} which are converted
|
---|
| 1533 | * to {@link SI#RADIAN radians}. This conversion is performed because the radians units
|
---|
| 1534 | * are standard for all internal computations in the map projection package. For example
|
---|
| 1535 | * they are the standard units for {@link MapProjection#latitudeOfOrigin latitudeOfOrigin}
|
---|
| 1536 | * and {@link MapProjection#centralMeridian centralMeridian} fields in the
|
---|
| 1537 | * {@link MapProjection} class.
|
---|
| 1538 | *
|
---|
| 1539 | * @param param The parameter to look for.
|
---|
| 1540 | * @param group The parameter value group to search into.
|
---|
| 1541 | * @return The requested parameter value.
|
---|
| 1542 | * @throws ParameterNotFoundException if the parameter is not found.
|
---|
| 1543 | */
|
---|
| 1544 | protected static double doubleValue(final ParameterDescriptor param,
|
---|
| 1545 | final ParameterValueGroup group)
|
---|
| 1546 | throws ParameterNotFoundException
|
---|
| 1547 | {
|
---|
| 1548 | double v = MathTransformProvider.doubleValue(param, group);
|
---|
| 1549 | if (NonSI.DEGREE_ANGLE.equals(param.getUnit())) {
|
---|
| 1550 | v = toRadians(v);
|
---|
| 1551 | }
|
---|
| 1552 | return v;
|
---|
| 1553 | }
|
---|
| 1554 | }
|
---|
| 1555 | }
|
---|