| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.sis.referencing.operation.projection; |
| |
| import java.util.EnumMap; |
| import org.opengis.util.FactoryException; |
| import org.opengis.parameter.ParameterDescriptor; |
| import org.opengis.referencing.operation.MathTransform; |
| import org.opengis.referencing.operation.MathTransformFactory; |
| import org.opengis.referencing.operation.Matrix; |
| import org.opengis.referencing.operation.OperationMethod; |
| import org.apache.sis.referencing.operation.matrix.Matrix2; |
| import org.apache.sis.referencing.operation.matrix.MatrixSIS; |
| import org.apache.sis.referencing.operation.transform.ContextualParameters; |
| import org.apache.sis.internal.referencing.provider.TransverseMercatorSouth; |
| import org.apache.sis.internal.referencing.Resources; |
| import org.apache.sis.internal.util.DoubleDouble; |
| import org.apache.sis.parameter.Parameters; |
| import org.apache.sis.util.resources.Errors; |
| import org.apache.sis.util.Workaround; |
| |
| import static java.lang.Math.*; |
| import static org.apache.sis.math.MathFunctions.asinh; |
| import static org.apache.sis.math.MathFunctions.atanh; |
| import static org.apache.sis.internal.referencing.provider.TransverseMercator.*; |
| |
| |
| /** |
| * <cite>Transverse Mercator</cite> projection (EPSG codes 9807). |
| * This class implements the "JHS formulas" reproduced in |
| * IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015. |
| * |
| * <div class="section">Description</div> |
| * This is a cylindrical projection, in which the cylinder has been rotated 90°. |
| * Instead of being tangent to the equator (or to an other standard latitude), it is tangent to a central meridian. |
| * Deformation are more important as we are going further from the central meridian. |
| * The Transverse Mercator projection is appropriate for region which have a greater extent north-south than east-west. |
| * |
| * <p>There are a number of versions of the Transverse Mercator projection including the Universal (UTM) |
| * and Modified (MTM) Transverses Mercator projections. In these cases the earth is divided into zones. |
| * For the UTM the zones are 6 degrees wide, numbered from 1 to 60 proceeding east from 180 degrees longitude, |
| * and between latitude 84 degrees North and 80 degrees South. The central meridian is taken as the center of the zone |
| * and the latitude of origin is the equator. A scale factor of 0.9996 and false easting of 500000 metres is used for |
| * all zones and a false northing of 10000000 metres is used for zones in the southern hemisphere. |
| * |
| * @author Martin Desruisseaux (Geomatys) |
| * @author Rémi Maréchal (Geomatys) |
| * @version 1.0 |
| * |
| * @see Mercator |
| * @see ObliqueMercator |
| * |
| * @since 0.6 |
| * @module |
| */ |
| public class TransverseMercator extends NormalizedProjection { |
| /** |
| * For cross-version compatibility. |
| */ |
| private static final long serialVersionUID = -627685138188387835L; |
| |
| /** |
| * Distance from central meridian, in degrees, at which errors are considered too important. |
| * This threshold is determined by comparisons of computed values against values provided by |
| * <cite>Karney (2009) Test data for the transverse Mercator projection</cite> data file. |
| * On the WGS84 ellipsoid we observed: |
| * |
| * <ul> |
| * <li>For ∆λ below 60°, errors below 1 centimetre.</li> |
| * <li>For ∆λ between 60° and 66°, errors up to 0.1 metre.</li> |
| * <li>For ∆λ between 66° and 70°, errors up to 1 metre.</li> |
| * <li>For ∆λ between 70° and 72°, errors up to 2 metres.</li> |
| * <li>For ∆λ between 72° and 74°, errors up to 10 metres.</li> |
| * <li>For ∆λ between 74° and 76°, errors up to 30 metres.</li> |
| * <li>For ∆λ between 76° and 78°, errors up to 1 kilometre.</li> |
| * <li>For ∆λ greater than 85°, results are chaotic.</li> |
| * </ul> |
| */ |
| static final double DOMAIN_OF_VALIDITY = 70; |
| |
| /** |
| * {@code false} for using the original formulas as published by EPSG, or {@code true} for using formulas |
| * modified using trigonometric identities. The use of trigonometric identities is for reducing the amount |
| * of calls to the {@link Math#sin(double)} and similar methods. Some identities used are: |
| * |
| * <ul> |
| * <li>sin(2θ) = 2⋅sinθ⋅cosθ</li> |
| * <li>cos(2θ) = cos²θ - sin²θ</li> |
| * <li>sin(3θ) = (3 - 4⋅sin²θ)⋅sinθ</li> |
| * <li>cos(3θ) = (4⋅cos³θ) - 3⋅cosθ</li> |
| * <li>sin(4θ) = (4 - 8⋅sin²θ)⋅sinθ⋅cosθ</li> |
| * <li>cos(4θ) = (8⋅cos⁴θ) - (8⋅cos²θ) + 1</li> |
| * </ul> |
| * |
| * Hyperbolic formulas: |
| * |
| * <ul> |
| * <li>sinh(2θ) = 2⋅sinhθ⋅coshθ</li> |
| * <li>cosh(2θ) = cosh²θ + sinh²θ = 2⋅cosh²θ - 1 = 1 + 2⋅sinh²θ</li> |
| * <li>sinh(3θ) = (3 + 4⋅sinh²θ)⋅sinhθ</li> |
| * <li>cosh(3θ) = ((4⋅cosh²θ) - 3)⋅coshθ</li> |
| * <li>sinh(4θ) = (1 + 2⋅sinh²θ)⋅4.sinhθ⋅coshθ |
| * = 4.cosh(2θ).sinhθ⋅coshθ</li> |
| * <li>cosh(4θ) = (8⋅cosh⁴θ) - (8⋅cosh²θ) + 1 |
| * = 8⋅cosh²(θ) ⋅ (cosh²θ - 1) + 1 |
| * = 8⋅cosh²(θ) ⋅ sinh²(θ) + 1 |
| * = 2⋅sinh²(2θ) + 1</li> |
| * </ul> |
| * |
| * Note that since this boolean is static final, the compiler should exclude the code in the branch that is never |
| * executed (no need to comment-out that code). |
| * |
| * @see #identityEquals(double, double) |
| */ |
| private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true; |
| |
| /** |
| * Verifies if a trigonometric identity produced the expected value. This method is used in assertions only, |
| * for values close to the [-1 … +1] range. The tolerance threshold is approximately 1.5E-12 (note that it |
| * still about 7000 time greater than {@code Math.ulp(1.0)}). |
| * |
| * @see #ALLOW_TRIGONOMETRIC_IDENTITIES |
| */ |
| private static boolean identityEquals(final double actual, final double expected) { |
| return !(abs(actual - expected) > // Use !(a > b) instead of (a <= b) in order to tolerate NaN. |
| (ANGULAR_TOLERANCE / 1000) * max(1, abs(expected))); // Increase tolerance for values outside the [-1 … +1] range. |
| } |
| |
| /** |
| * Coefficients in the series expansion of the forward projection, |
| * depending only on {@linkplain #eccentricity eccentricity} value. |
| * The series expansion is of the following form: |
| * |
| * <blockquote>cf₂⋅ｆ(2θ) + cf₄⋅ｆ(4θ) + cf₆⋅ｆ(6θ) + cf₈⋅ｆ(8θ)</blockquote> |
| * |
| * Those coefficients are named h₁, h₂, h₃ and h₄ in §1.3.5.1 of |
| * IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015. |
| * |
| * <div class="note"><b>Serialization note:</b> |
| * we do not strictly need to serialize those fields since they could be computed after deserialization. |
| * Bu we serialize them anyway in order to simplify a little bit this class (it allows us to keep those |
| * fields final) and because values computed after deserialization could be slightly different than the |
| * ones computed after construction since the constructor uses the double-double values provided by |
| * {@link Initializer}.</div> |
| */ |
| private final double cf2, cf4, cf6, cf8; |
| |
| /** |
| * Coefficients in the series expansion of the inverse projection, |
| * depending only on {@linkplain #eccentricity eccentricity} value. |
| */ |
| private final double ci2, ci4, ci6, ci8; |
| |
| /** |
| * Creates a Transverse Mercator projection from the given parameters. |
| * The {@code method} argument can be the description of one of the following: |
| * |
| * <ul> |
| * <li><cite>"Transverse Mercator"</cite>.</li> |
| * <li><cite>"Transverse Mercator (South Orientated)"</cite>.</li> |
| * </ul> |
| * |
| * @param method description of the projection parameters. |
| * @param parameters the parameter values of the projection to create. |
| */ |
| public TransverseMercator(final OperationMethod method, final Parameters parameters) { |
| this(initializer(method, parameters)); |
| } |
| |
| /** |
| * Work around for RFE #4093999 in Sun's bug database |
| * ("Relax constraint on placement of this()/super() call in constructors"). |
| */ |
| @Workaround(library="JDK", version="1.7") |
| private static Initializer initializer(final OperationMethod method, final Parameters parameters) { |
| final boolean isSouth = identMatch(method, "(?i).*\\bSouth\\b.*", TransverseMercatorSouth.IDENTIFIER); |
| final EnumMap<ParameterRole, ParameterDescriptor<Double>> roles = new EnumMap<>(ParameterRole.class); |
| ParameterRole xOffset = ParameterRole.FALSE_EASTING; |
| ParameterRole yOffset = ParameterRole.FALSE_NORTHING; |
| if (isSouth) { |
| xOffset = ParameterRole.FALSE_WESTING; |
| yOffset = ParameterRole.FALSE_SOUTHING; |
| } |
| roles.put(ParameterRole.CENTRAL_MERIDIAN, LONGITUDE_OF_ORIGIN); |
| roles.put(ParameterRole.SCALE_FACTOR, SCALE_FACTOR); |
| roles.put(xOffset, FALSE_EASTING); |
| roles.put(yOffset, FALSE_NORTHING); |
| return new Initializer(method, parameters, roles, isSouth ? (byte) 1 : (byte) 0); |
| } |
| |
| /** |
| * Creates a new Transverse Mercator projection from the given initializer. |
| * This constructor is used also by {@link ZonedGridSystem}. |
| */ |
| TransverseMercator(final Initializer initializer) { |
| super(initializer); |
| final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN)); |
| /* |
| * Opportunistically use double-double arithmetic for computation of B since we will store |
| * it in the denormalization matrix, and there is no sine/cosine functions involved here. |
| */ |
| double cf4, cf6, cf8, ci4, ci6, ci8; |
| final DoubleDouble B; |
| { // For keeping the `t` and `n` variables locale. |
| /* |
| * The n parameter is defined by the EPSG guide as: |
| * |
| * n = f / (2-f) |
| * |
| * Where f is the flattening factor (1 - b/a). This equation can be rewritten as: |
| * |
| * n = (1 - b/a) / (1 + b/a) |
| * |
| * As much as possible, b/a should be computed from the map projection parameters. |
| * However if those parameters are not available anymore, then they can be computed |
| * from the eccentricity as: |
| * |
| * b/a = √(1 - ℯ²) |
| */ |
| final DoubleDouble t = initializer.axisLengthRatio(); // t = b/a |
| t.ratio_1m_1p(); // t = (1 - t) / (1 + t) |
| final double n = t.doubleValue(); // n = f / (2-f) |
| final double n2 = n * n; |
| final double n3 = n2 * n; |
| final double n4 = n2 * n2; |
| /* |
| * Coefficients for the forward projections. |
| * Add the smallest values first in order to reduce rounding errors. |
| */ |
| cf2 = ( 41. / 180)*n4 + ( 5. / 16)*n3 + (-2. / 3)*n2 + n/2; |
| cf4 = ( 557. / 1440)*n4 + (-3. / 5)*n3 + (13. / 48)*n2; |
| cf6 = ( -103. / 140)*n4 + (61. / 240)*n3; |
| cf8 = (49561. / 161280)*n4; |
| /* |
| * Coefficients for the inverse projections. |
| * Add the smallest values first in order to reduce rounding errors. |
| */ |
| ci2 = ( -1. / 360)*n4 + (37. / 96)*n3 + (-2. / 3)*n2 + n/2; |
| ci4 = (-437. / 1440)*n4 + ( 1. / 15)*n3 + ( 1. / 48)*n2; |
| ci6 = ( -37. / 840)*n4 + (17. / 480)*n3; |
| ci8 = (4397. / 161280)*n4; |
| /* |
| * Compute B = (1 + n²/4 + n⁴/64) / (1 + n) |
| */ |
| B = new DoubleDouble(t); // B = n |
| B.square(); |
| B.series(1, 0.25, 1./64); // B = (1 + n²/4 + n⁴/64) |
| t.add(1); |
| B.divide(t); // B = (1 + n²/4 + n⁴/64) / (1 + n) |
| } |
| /* |
| * Compute M₀ = B⋅(ξ₁ + ξ₂ + ξ₃ + ξ₄) and negate in anticipation for what will be needed |
| * in the denormalization matrix. We opportunistically use double-double arithmetic but |
| * only for the final multiplication by B, for consistency with the translation term to |
| * be stored in the denormalization matrix. It is not worth to use double-double in the |
| * sum of sine functions because the extra digits would be meaningless. |
| * |
| * NOTE: the EPSG documentation makes special cases for φ₀ = 0 or ±π/2. This is not |
| * needed here; we verified that the code below produces naturally the expected values. |
| */ |
| final double Q = asinh(tan(φ0)) - eccentricity * atanh(eccentricity * sin(φ0)); |
| final double β = atan(sinh(Q)); |
| final DoubleDouble M0 = new DoubleDouble(); |
| M0.value = cf8 * sin(8*β) |
| + cf6 * sin(6*β) |
| + cf4 * sin(4*β) |
| + cf2 * sin(2*β) |
| + β; |
| M0.multiply(B); |
| M0.negate(); |
| /* |
| * At this point, all parameters have been processed. Now store |
| * the linear operations in the (de)normalize affine transforms: |
| * |
| * Normalization: |
| * - Subtract central meridian to longitudes (done by the super-class constructor). |
| * - Convert longitudes and latitudes from degrees to radians (done by the super-class constructor) |
| * |
| * Denormalization |
| * - Scale x and y by B. |
| * - Subtract M₀ to the northing. |
| * - Multiply by the scale factor (done by the super-class constructor). |
| * - Add false easting and false northing (done by the super-class constructor). |
| */ |
| final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION); |
| denormalize.convertBefore(0, B, null); |
| denormalize.convertBefore(1, B, M0); |
| /* |
| * When rewriting equations using trigonometric identities, some constants appear. |
| * For example sin(2θ) = 2⋅sinθ⋅cosθ, so we can factor out the 2 constant into the |
| * corresponding 'c' field. Note: this factorization can only be performed after |
| * the constructor finished to compute other constants. |
| */ |
| if (ALLOW_TRIGONOMETRIC_IDENTITIES) { |
| // Multiplication by powers of 2 does not bring any additional rounding error. |
| cf4 *= 4; ci4 *= 4; |
| cf6 *= 16; ci6 *= 16; |
| cf8 *= 64; ci8 *= 64; |
| } |
| this.cf4 = cf4; |
| this.cf6 = cf6; |
| this.cf8 = cf8; |
| this.ci4 = ci4; |
| this.ci6 = ci6; |
| this.ci8 = ci8; |
| } |
| |
| /** |
| * Creates a new projection initialized to the same parameters than the given one. |
| */ |
| TransverseMercator(final TransverseMercator other) { |
| super(other); |
| cf2 = other.cf2; |
| cf4 = other.cf4; |
| cf6 = other.cf6; |
| cf8 = other.cf8; |
| ci2 = other.ci2; |
| ci4 = other.ci4; |
| ci6 = other.ci6; |
| ci8 = other.ci8; |
| } |
| |
| /** |
| * Returns the sequence of <cite>normalization</cite> → {@code this} → <cite>denormalization</cite> transforms |
| * as a whole. The transform returned by this method expects (<var>longitude</var>, <var>latitude</var>) |
| * coordinates in <em>degrees</em> and returns (<var>x</var>,<var>y</var>) coordinates in <em>metres</em>. |
| * |
| * <p>The non-linear part of the returned transform will be {@code this} transform, except if the ellipsoid |
| * is spherical. In the later case, {@code this} transform will be replaced by a simplified implementation.</p> |
| * |
| * @param factory the factory to use for creating the transform. |
| * @return the map projection from (λ,φ) to (<var>x</var>,<var>y</var>) coordinates. |
| * @throws FactoryException if an error occurred while creating a transform. |
| */ |
| @Override |
| public MathTransform createMapProjection(final MathTransformFactory factory) throws FactoryException { |
| TransverseMercator kernel = this; |
| if (eccentricity == 0) { |
| kernel = new Spherical(this); |
| } |
| return context.completeTransform(factory, kernel); |
| } |
| |
| /** |
| * Converts the specified (λ,φ) coordinate (units in radians) and stores the result in {@code dstPts}. |
| * In addition, opportunistically computes the projection derivative if {@code derivate} is {@code true}. |
| * |
| * @return the matrix of the projection derivative at the given source position, |
| * or {@code null} if the {@code derivate} argument is {@code false}. |
| * @throws ProjectionException if the coordinate can not be converted. |
| */ |
| @Override |
| public Matrix transform(final double[] srcPts, final int srcOff, |
| final double[] dstPts, final int dstOff, |
| final boolean derivate) throws ProjectionException |
| { |
| final double λ = srcPts[srcOff]; |
| if (abs(λ) >= DOMAIN_OF_VALIDITY * (PI/180)) { |
| /* |
| * The Transverse Mercator projection is conceptually a Mercator projection rotated by 90°. |
| * In Mercator projection, y values tend toward infinity for latitudes close to ±90°. |
| * In Transverse Mercator, x values tend toward infinity for longitudes close to ±90° |
| * at equator and after subtraction of central meridian. After we pass the 90° limit, |
| * the Transverse Mercator results at (90° + Δ) are the same as for (90° - Δ). |
| * |
| * Problem is that 90° is an ordinary longitude value, not even close to the limit of longitude |
| * values range (±180°). So having f(π/2+Δ, φ) = f(π/2-Δ, φ) results in wrong behavior in some |
| * algorithms like the one used by Envelopes.transform(CoordinateOperation, Envelope). |
| * Since a distance of 90° from central meridian is far outside the Transverse Mercator |
| * domain of validity anyway, we do not let the user go further. |
| * |
| * In the particular case of ellipsoidal formulas, we put a limit of 81° instead of 90° |
| * because experience shows that results close to equator become chaotic after 85° when |
| * using WGS84 ellipsoid. We do not need to reduce the limit for the spherical formulas, |
| * because the mathematic are simpler and the function still smooth until 90°. |
| */ |
| throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity)); |
| } |
| final double φ = srcPts[srcOff+1]; |
| final double sinλ = sin(λ); |
| final double ℯsinφ = sin(φ) * eccentricity; |
| final double Q = asinh(tan(φ)) - atanh(ℯsinφ) * eccentricity; |
| final double coshQ = cosh(Q); // Can not be smaller than 1. |
| final double η0 = atanh(sinλ / coshQ); // Tends toward ±∞ if λ → ±90°. |
| /* |
| * Original formula: η0 = atanh(sin(λ) * cos(β)) where |
| * cos(β) = cos(atan(sinh(Q))) |
| * = 1 / sqrt(1 + sinh²(Q)) |
| * = 1 / (sqrt(cosh²(Q) - sinh²(Q) + sinh²(Q))) |
| * = 1 / sqrt(cosh²(Q)) |
| * = 1 / cosh(Q) |
| * |
| * So η0 = atanh(sin(λ) / cosh(Q)) |
| */ |
| final double coshη0 = cosh(η0); |
| final double ξ0 = asin(tanh(Q) * coshη0); |
| /* |
| * Compute sin(2⋅ξ₀), sin(4⋅ξ₀), sin(6⋅ξ₀), sin(8⋅ξ₀) and same for cos, but using the following |
| * trigonometric identities in order to reduce the number of calls to Math.sin and cos methods. |
| */ |
| final double sin_2ξ0 = sin(2*ξ0); |
| final double cos_2ξ0 = cos(2*ξ0); |
| final double sin_4ξ0, sin_6ξ0, sin_8ξ0, |
| cos_4ξ0, cos_6ξ0, cos_8ξ0; |
| if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { |
| sin_4ξ0 = sin(4*ξ0); |
| cos_4ξ0 = cos(4*ξ0); |
| sin_6ξ0 = sin(6*ξ0); |
| cos_6ξ0 = cos(6*ξ0); |
| sin_8ξ0 = sin(8*ξ0); |
| cos_8ξ0 = cos(8*ξ0); |
| } else { |
| final double sin2 = sin_2ξ0 * sin_2ξ0; |
| final double cos2 = cos_2ξ0 * cos_2ξ0; |
| sin_4ξ0 = sin_2ξ0 * cos_2ξ0; assert identityEquals(sin_4ξ0, sin(4*ξ0) / 2) : ξ0; |
| cos_4ξ0 = (cos2 - sin2) * 0.5; assert identityEquals(cos_4ξ0, cos(4*ξ0) / 2) : ξ0; |
| sin_6ξ0 = (0.75 - sin2) * sin_2ξ0; assert identityEquals(sin_6ξ0, sin(6*ξ0) / 4) : ξ0; |
| cos_6ξ0 = (cos2 - 0.75) * cos_2ξ0; assert identityEquals(cos_6ξ0, cos(6*ξ0) / 4) : ξ0; |
| sin_8ξ0 = sin_4ξ0 * cos_4ξ0; assert identityEquals(sin_8ξ0, sin(8*ξ0) / 8) : ξ0; |
| cos_8ξ0 = 0.125 - sin_4ξ0 * sin_4ξ0; assert identityEquals(cos_8ξ0, cos(8*ξ0) / 8) : ξ0; |
| } |
| /* |
| * Compute sinh(2⋅ξ₀), sinh(4⋅ξ₀), sinh(6⋅ξ₀), sinh(8⋅ξ₀) and same for cosh, but using the following |
| * hyperbolic identities in order to reduce the number of calls to Math.sinh and cosh methods. |
| * Note that the formulas are very similar to the above ones, with only some signs reversed. |
| */ |
| final double sinh_2η0 = sinh(2*η0); |
| final double cosh_2η0 = cosh(2*η0); |
| final double sinh_4η0, sinh_6η0, sinh_8η0, |
| cosh_4η0, cosh_6η0, cosh_8η0; |
| if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { |
| sinh_4η0 = sinh(4*η0); |
| cosh_4η0 = cosh(4*η0); |
| sinh_6η0 = sinh(6*η0); |
| cosh_6η0 = cosh(6*η0); |
| sinh_8η0 = sinh(8*η0); |
| cosh_8η0 = cosh(8*η0); |
| } else { |
| final double sinh2 = sinh_2η0 * sinh_2η0; |
| final double cosh2 = cosh_2η0 * cosh_2η0; |
| cosh_4η0 = (cosh2 + sinh2) * 0.5; assert identityEquals(cosh_4η0, cosh(4*η0) / 2) : η0; |
| sinh_4η0 = cosh_2η0 * sinh_2η0; assert identityEquals(sinh_4η0, sinh(4*η0) / 2) : η0; |
| cosh_6η0 = cosh_2η0 * (cosh2 - 0.75); assert identityEquals(cosh_6η0, cosh(6*η0) / 4) : η0; |
| sinh_6η0 = sinh_2η0 * (sinh2 + 0.75); assert identityEquals(sinh_6η0, sinh(6*η0) / 4) : η0; |
| cosh_8η0 = sinh_4η0 * sinh_4η0 + 0.125; assert identityEquals(cosh_8η0, cosh(8*η0) / 8) : η0; |
| sinh_8η0 = sinh_4η0 * cosh_4η0; assert identityEquals(sinh_8η0, sinh(8*η0) / 8) : η0; |
| } |
| /* |
| * The projection of (λ,φ) is given by (η⋅B, ξ⋅B+M₀) — ignoring scale factors and false easting/northing. |
| * But the B and M₀ parameters have been merged by the constructor with other linear operations in the |
| * "denormalization" matrix. Consequently we only need to compute (η,ξ) below. |
| */ |
| if (dstPts != null) { |
| // η(λ,φ) |
| dstPts[dstOff ] = cf8 * cos_8ξ0 * sinh_8η0 |
| + cf6 * cos_6ξ0 * sinh_6η0 |
| + cf4 * cos_4ξ0 * sinh_4η0 |
| + cf2 * cos_2ξ0 * sinh_2η0 |
| + η0; |
| // ξ(λ,φ) |
| dstPts[dstOff+1] = cf8 * sin_8ξ0 * cosh_8η0 |
| + cf6 * sin_6ξ0 * cosh_6η0 |
| + cf4 * sin_4ξ0 * cosh_4η0 |
| + cf2 * sin_2ξ0 * cosh_2η0 |
| + ξ0; |
| } |
| if (!derivate) { |
| return null; |
| } |
| /* |
| * Now compute the derivative, if the user asked for it. |
| */ |
| final double cosλ = cos(λ); //-- λ |
| final double cosφ = cos(φ); //-- φ |
| final double cosh2Q = coshQ * coshQ; //-- Q |
| final double sinhQ = sinh(Q); |
| final double tanhQ = tanh(Q); |
| final double cosh2Q_sin2λ = cosh2Q - (sinλ * sinλ); //-- Qλ |
| final double sinhη0 = sinh(η0); //-- η0 |
| final double sqrt1_thQchη0 = sqrt(1 - (tanhQ * tanhQ) * (coshη0 * coshη0)); //-- Qη0 |
| |
| //-- dQ_dλ = 0; |
| final double dQ_dφ = 1 / cosφ - eccentricitySquared * cosφ / (1 - ℯsinφ * ℯsinφ); |
| |
| final double dη0_dλ = cosλ * coshQ / cosh2Q_sin2λ; |
| final double dη0_dφ = -dQ_dφ * sinλ * sinhQ / cosh2Q_sin2λ; |
| |
| final double dξ0_dλ = sinhQ * sinhη0 * cosλ / (cosh2Q_sin2λ * sqrt1_thQchη0); |
| final double dξ0_dφ = (dQ_dφ * coshη0 / cosh2Q + dη0_dφ * sinhη0 * tanhQ) / sqrt1_thQchη0; |
| /* |
| * Jac(Proj(λ,φ)) is the Jacobian matrix of Proj(λ,φ) function. |
| * So the derivative of Proj(λ,φ) is defined by: |
| * |
| * ┌ ┐ |
| * │ ∂η(λ,φ)/∂λ, ∂η(λ,φ)/∂φ │ |
| * Jac = │ │ |
| * (Proj(λ,φ)) │ ∂ξ(λ,φ)/∂λ, ∂ξ(λ,φ)/∂φ │ |
| * └ ┘ |
| */ |
| //-- dξ(λ, φ) / dλ |
| final double dξ_dλ = dξ0_dλ |
| + 2 * (cf2 * (dξ0_dλ * cos_2ξ0 * cosh_2η0 + dη0_dλ * sinh_2η0 * sin_2ξ0) |
| + 3 * cf6 * (dξ0_dλ * cos_6ξ0 * cosh_6η0 + dη0_dλ * sinh_6η0 * sin_6ξ0) |
| + 2 * (cf4 * (dξ0_dλ * cos_4ξ0 * cosh_4η0 + dη0_dλ * sinh_4η0 * sin_4ξ0) |
| + 2 * cf8 * (dξ0_dλ * cos_8ξ0 * cosh_8η0 + dη0_dλ * sinh_8η0 * sin_8ξ0))); |
| |
| //-- dξ(λ, φ) / dφ |
| final double dξ_dφ = dξ0_dφ |
| + 2 * (cf2 * (dξ0_dφ * cos_2ξ0 * cosh_2η0 + dη0_dφ * sinh_2η0 * sin_2ξ0) |
| + 3 * cf6 * (dξ0_dφ * cos_6ξ0 * cosh_6η0 + dη0_dφ * sinh_6η0 * sin_6ξ0) |
| + 2 * (cf4 * (dξ0_dφ * cos_4ξ0 * cosh_4η0 + dη0_dφ * sinh_4η0 * sin_4ξ0) |
| + 2 * cf8 * (dξ0_dφ * cos_8ξ0 * cosh_8η0 + dη0_dφ * sinh_8η0 * sin_8ξ0))); |
| |
| //-- dη(λ, φ) / dλ |
| final double dη_dλ = dη0_dλ |
| + 2 * (cf2 * (dη0_dλ * cosh_2η0 * cos_2ξ0 - dξ0_dλ * sin_2ξ0 * sinh_2η0) |
| + 3 * cf6 * (dη0_dλ * cosh_6η0 * cos_6ξ0 - dξ0_dλ * sin_6ξ0 * sinh_6η0) |
| + 2 * (cf4 * (dη0_dλ * cosh_4η0 * cos_4ξ0 - dξ0_dλ * sin_4ξ0 * sinh_4η0) |
| + 2 * cf8 * (dη0_dλ * cosh_8η0 * cos_8ξ0 - dξ0_dλ * sin_8ξ0 * sinh_8η0))); |
| |
| //-- dη(λ, φ) / dφ |
| final double dη_dφ = dη0_dφ |
| + 2 * (cf2 * (dη0_dφ * cosh_2η0 * cos_2ξ0 - dξ0_dφ * sin_2ξ0 * sinh_2η0) |
| + 3 * cf6 * (dη0_dφ * cosh_6η0 * cos_6ξ0 - dξ0_dφ * sin_6ξ0 * sinh_6η0) |
| + 2 * (cf4 * (dη0_dφ * cosh_4η0 * cos_4ξ0 - dξ0_dφ * sin_4ξ0 * sinh_4η0) |
| + 2 * cf8 * (dη0_dφ * cosh_8η0 * cos_8ξ0 - dξ0_dφ * sin_8ξ0 * sinh_8η0))); |
| |
| return new Matrix2(dη_dλ, dη_dφ, |
| dξ_dλ, dξ_dφ); |
| } |
| |
| /** |
| * Transforms the specified (η, ξ) coordinates and stores the result in {@code dstPts} (angles in radians). |
| * |
| * @throws ProjectionException if the point can not be converted. |
| */ |
| @Override |
| protected void inverseTransform(final double[] srcPts, final int srcOff, |
| final double[] dstPts, final int dstOff) |
| throws ProjectionException |
| { |
| final double η = srcPts[srcOff ]; |
| final double ξ = srcPts[srcOff+1]; |
| /* |
| * Following calculation of sin_2ξ, sin_4ξ, etc. is basically a copy-and-paste of the code in transform(…). |
| * Its purpose is the same than for transform(…): reduce the amount of calls to Math.sin(double) and other |
| * methods. |
| */ |
| final double sin_2ξ = sin (2*ξ); |
| final double cos_2ξ = cos (2*ξ); |
| final double sinh_2η = sinh(2*η); |
| final double cosh_2η = cosh(2*η); |
| final double sin_4ξ, sin_6ξ, sin_8ξ, |
| cos_4ξ, cos_6ξ, cos_8ξ, |
| sinh_4η, sinh_6η, sinh_8η, |
| cosh_4η, cosh_6η, cosh_8η; |
| if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { |
| sin_4ξ = sin(4*ξ); |
| cos_4ξ = cos(4*ξ); |
| sin_6ξ = sin(6*ξ); |
| cos_6ξ = cos(6*ξ); |
| sin_8ξ = sin(8*ξ); |
| cos_8ξ = cos(8*ξ); |
| |
| sinh_4η = sinh(4*η); |
| cosh_4η = cosh(4*η); |
| sinh_6η = sinh(6*η); |
| cosh_6η = cosh(6*η); |
| sinh_8η = sinh(8*η); |
| cosh_8η = cosh(8*η); |
| } else { |
| final double sin2 = sin_2ξ * sin_2ξ; |
| final double cos2 = cos_2ξ * cos_2ξ; |
| sin_4ξ = sin_2ξ * cos_2ξ; assert identityEquals(sin_4ξ, sin(4*ξ) / 2) : ξ; |
| cos_4ξ = (cos2 - sin2) * 0.5; assert identityEquals(cos_4ξ, cos(4*ξ) / 2) : ξ; |
| sin_6ξ = (0.75 - sin2) * sin_2ξ; assert identityEquals(sin_6ξ, sin(6*ξ) / 4) : ξ; |
| cos_6ξ = (cos2 - 0.75) * cos_2ξ; assert identityEquals(cos_6ξ, cos(6*ξ) / 4) : ξ; |
| sin_8ξ = sin_4ξ * cos_4ξ; assert identityEquals(sin_8ξ, sin(8*ξ) / 8) : ξ; |
| cos_8ξ = 0.125 - sin_4ξ * sin_4ξ; assert identityEquals(cos_8ξ, cos(8*ξ) / 8) : ξ; |
| |
| final double sinh2 = sinh_2η * sinh_2η; |
| final double cosh2 = cosh_2η * cosh_2η; |
| cosh_4η = (cosh2 + sinh2) * 0.5; assert identityEquals(cosh_4η, cosh(4*η) / 2) : η; |
| sinh_4η = cosh_2η * sinh_2η; assert identityEquals(sinh_4η, sinh(4*η) / 2) : η; |
| cosh_6η = cosh_2η * (cosh2 - 0.75); assert identityEquals(cosh_6η, cosh(6*η) / 4) : η; |
| sinh_6η = sinh_2η * (sinh2 + 0.75); assert identityEquals(sinh_6η, sinh(6*η) / 4) : η; |
| cosh_8η = sinh_4η * sinh_4η + 0.125; assert identityEquals(cosh_8η, cosh(8*η) / 8) : η; |
| sinh_8η = sinh_4η * cosh_4η; assert identityEquals(sinh_8η, sinh(8*η) / 8) : η; |
| } |
| /* |
| * The actual inverse transform. |
| */ |
| final double ξ0 = ξ - (ci8 * sin_8ξ * cosh_8η |
| + ci6 * sin_6ξ * cosh_6η |
| + ci4 * sin_4ξ * cosh_4η |
| + ci2 * sin_2ξ * cosh_2η); |
| |
| final double η0 = η - (ci8 * cos_8ξ * sinh_8η |
| + ci6 * cos_6ξ * sinh_6η |
| + ci4 * cos_4ξ * sinh_4η |
| + ci2 * cos_2ξ * sinh_2η); |
| |
| final double β = asin(sin(ξ0) / cosh(η0)); |
| final double Q = asinh(tan(β)); |
| /* |
| * Following usually converges in 4 iterations. |
| * The first iteration is unrolled. |
| */ |
| double p = eccentricity * atanh(eccentricity * tanh(Q)); |
| double Qp = Q + p; |
| for (int it=0; it<MAXIMUM_ITERATIONS; it++) { |
| final double c = eccentricity * atanh(eccentricity * tanh(Qp)); |
| Qp = Q + c; |
| if (abs(c - p) <= ITERATION_TOLERANCE) { |
| dstPts[dstOff ] = asin(tanh(η0) / cos(β)); |
| dstPts[dstOff+1] = atan(sinh(Qp)); |
| return; |
| } |
| p = c; |
| } |
| throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence)); |
| } |
| |
| |
| |
| |
| /** |
| * Provides the transform equations for the spherical case of the Transverse Mercator projection. |
| * |
| * @author André Gosselin (MPO) |
| * @author Martin Desruisseaux (IRD, Geomatys) |
| * @author Rueben Schulz (UBC) |
| * @version 0.6 |
| * @since 0.6 |
| * @module |
| */ |
| private static final class Spherical extends TransverseMercator { |
| /** |
| * For cross-version compatibility. |
| */ |
| private static final long serialVersionUID = 8903592710452235162L; |
| |
| /** |
| * Constructs a new map projection from the parameters of the given projection. |
| * |
| * @param other the other projection (usually ellipsoidal) from which to copy the parameters. |
| */ |
| protected Spherical(final TransverseMercator other) { |
| super(other); |
| } |
| |
| /** |
| * {@inheritDoc} |
| */ |
| @Override |
| public Matrix transform(final double[] srcPts, final int srcOff, |
| final double[] dstPts, final int dstOff, |
| final boolean derivate) throws ProjectionException |
| { |
| final double λ = srcPts[srcOff]; |
| /* |
| * The Transverse Mercator projection is conceptually a Mercator projection rotated by 90°. |
| * In Mercator projection, y values tend toward infinity for latitudes close to ±90° while |
| * in Transverse Mercator, x values tend toward infinity for longitudes close to ±90° from |
| * central meridian (and at equator). After we pass the 90° limit, the Transverse Mercator |
| * results at (90° + Δ) are the same as for (90° - Δ). |
| * |
| * Problem is that 90° is an ordinary longitude value, not even close to the limit of longitude |
| * values range (±180°). So having f(π/2+Δ, φ) = f(π/2-Δ, φ) results in wrong behavior in some |
| * algorithms like the one used by Envelopes.transform(CoordinateOperation, Envelope). Since a |
| * distance of 90° from central meridian is outside projection domain of validity anyway, we do |
| * not let the user go further. |
| * |
| * Note that in the case of ellipsoidal formulas (implemented by parent class), we put the |
| * limit to a lower value (DOMAIN_OF_VALIDITY) because experience shows that results close |
| * to equator become chaotic after 85° when using WGS84 ellipsoid. We do not need to reduce |
| * the limit for the spherical formulas, because the mathematic are simpler and the function |
| * still smooth until 90°. |
| */ |
| if (abs(λ) > PI/2) { |
| throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity)); |
| } |
| final double φ = srcPts[srcOff+1]; |
| final double sinλ = sin(λ); |
| final double cosλ = cos(λ); |
| final double sinφ = sin(φ); |
| final double cosφ = cos(φ); |
| final double tanφ = sinφ / cosφ; |
| final double B = cosφ * sinλ; |
| /* |
| * Using Snyder's equation for calculating y, instead of the one used in Proj.4. |
| * Potential problems when y and x = 90 degrees, but behaves ok in tests. |
| */ |
| if (dstPts != null) { |
| dstPts[dstOff ] = atanh(B); // Snyder 8-1; |
| dstPts[dstOff+1] = atan2(tanφ, cosλ); // Snyder 8-3; |
| } |
| if (!derivate) { |
| return null; |
| } |
| final double Bm = B*B - 1; |
| final double sct = cosλ*cosλ + tanφ*tanφ; |
| return new Matrix2(-(cosφ * cosλ) / Bm, // ∂x/∂λ |
| (sinφ * sinλ) / Bm, // ∂x/∂φ |
| (tanφ * sinλ) / sct, // ∂y/∂λ |
| cosλ / (cosφ * cosφ * sct)); // ∂y/∂φ |
| } |
| |
| /** |
| * {@inheritDoc} |
| */ |
| @Override |
| protected void inverseTransform(final double[] srcPts, final int srcOff, |
| final double[] dstPts, final int dstOff) |
| { |
| final double x = srcPts[srcOff ]; |
| final double y = srcPts[srcOff+1]; |
| final double sinhx = sinh(x); |
| final double cosy = cos(y); |
| // 'copySign' corrects for the fact that we made everything positive using sqrt(…) |
| dstPts[dstOff ] = atan2(sinhx, cosy); |
| dstPts[dstOff+1] = copySign(asin(sqrt((1 - cosy*cosy) / (1 + sinhx*sinhx))), y); |
| } |
| } |
| } |
