| /* |
| * 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 org.apache.sis.internal.referencing.Resources; |
| |
| import static java.lang.Math.*; |
| |
| |
| /** |
| * Base class of {@link LambertConicConformal}, {@link Mercator} and {@link PolarStereographic} projections. |
| * Those projections have in common the property of being <cite>conformal</cite>, i.e. they preserve angles locally. |
| * However we do not put this base class in public API because not all conformal projections extend this base class. |
| * For example the {@link TransverseMercator} projection, despite being conformal, uses very different formulas. |
| * |
| * <p>Note that no projection can be both conformal and equal-area. This restriction is implemented in class |
| * hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct classes.</p> |
| * |
| * <p>This base class can been seen as a generalization of <cite>Lambert Conic Conformal</cite> projection, |
| * which includes some other projections like Mercator and Polar Stereographic as special cases. |
| * For this base class, the Mercator projection is considered as <cite>"a special limiting case of the |
| * Lambert Conic Conformal map projection with the equator as the single standard parallel."</cite> |
| * (Source: §1.3.3 in IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015). |
| * Indeed, those two projections have some equation in commons which are provided in this base class.</p> |
| * |
| * <p>The Polar Stereographic projection is not documented as a special case of Lambert Conic Conformal, |
| * but the equations in the {@code PolarStereographic.transform(…)} and {@code inverseTransform(…)} methods |
| * appear to be the same with the <var>n</var> factor fixed to 1 or -1, so we leverage the code provided by |
| * this base class. This class hierarchy is only an implementation convenience and not part of public API.</p> |
| * |
| * <p>The Transverse Mercator projection is also conformal, but does not use the formulas provided in this class. |
| * It will instead compute the coefficients itself and use its own, more complex, formulas with those coefficients. |
| * However the formulas provided in this {@code ConformalProjection} class can be seen as a special case of |
| * Transverse Mercator formulas for <var>x</var> = 0.</p> |
| * |
| * <div class="note"><b>Reference:</b> |
| * “Lambert developed the regular Conformal Conic as the oblique aspect of a family containing the previously |
| * known polar Stereographic and regular Mercator projections. (…) If the standard parallels are symmetrical |
| * about the Equator, the regular Mercator results (although formulas must be revised). If the only standard |
| * parallel is a pole, the polar Stereographic results.” (Snyder, page 105)</div> |
| * |
| * @author Martin Desruisseaux (Geomatys) |
| * @version 1.0 |
| * @since 0.6 |
| * @module |
| */ |
| abstract class ConformalProjection extends NormalizedProjection { |
| /** |
| * For cross-version compatibility. |
| */ |
| private static final long serialVersionUID = -292755890184138414L; |
| |
| /** |
| * The threshold value of {@link #eccentricity} at which we consider the accuracy of the |
| * series expansion insufficient. This threshold is determined empirically with the help |
| * of the {@code MercatorMethodComparison} class in the test directory. |
| * We choose the value where: |
| * |
| * <ul> |
| * <li>the average error of series expansion become greater than {@link NormalizedProjection#ITERATION_TOLERANCE},</li> |
| * <li>the maximal error of series expansion become greater than {@link NormalizedProjection#ANGULAR_TOLERANCE}.</li> |
| * </ul> |
| */ |
| static final double ECCENTRICITY_THRESHOLD = 0.16; |
| |
| /** |
| * {@code true} if the {@link #eccentricity} value is greater than or equals to {@link #ECCENTRICITY_THRESHOLD}, |
| * in which case the {@link #φ(double)} method will need to use an iterative method. |
| */ |
| private final boolean useIterations; |
| |
| /** |
| * Coefficients of the first terms in the series expansion of the inverse projection. |
| * Values of those coefficients depend only on {@linkplain #eccentricity eccentricity} value. |
| * The series expansion is published under the following form, where χ is the <cite>conformal latitude</cite>: |
| * |
| * <blockquote>c₂⋅sin(2χ) + c₄⋅sin(4χ) + c₆⋅sin(6χ) + c₈⋅sin(8χ)</blockquote> |
| * |
| * However we rewrite above series expansion for taking advantage of trigonometric identities. |
| * The equation become (with different <var>c</var> values): |
| * |
| * <blockquote>sin(2χ)⋅(c₂ + cos(2χ)⋅(c₄ + cos(2χ)⋅(c₆ + cos(2χ)⋅c₈)))</blockquote> |
| * |
| * <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 if a future version of the constructor uses the double-double values |
| * provided by {@link Initializer}.</div> |
| */ |
| private final double c2χ, c4χ, c6χ, c8χ; |
| |
| /** |
| * Creates a new normalized projection from the parameters computed by the given initializer. |
| * This constructor computes the coefficients in the series expansion from the {@link #eccentricitySquared} value. |
| * |
| * @param initializer the initializer for computing map projection internal parameters. |
| */ |
| ConformalProjection(final Initializer initializer) { |
| super(initializer); |
| useIterations = (eccentricity >= ECCENTRICITY_THRESHOLD); |
| final double e2 = eccentricitySquared; |
| final double e4 = e2 * e2; |
| final double e6 = e2 * e4; |
| final double e8 = e4 * e4; |
| /* |
| * Published equation is of the following form: |
| * |
| * φ = χ + A⋅sin(2χ) + B⋅sin(4χ) + C⋅sin(6χ) + D⋅sin(8χ) (Snyder 3-34) |
| * |
| * We rewrite as: |
| * |
| * φ = χ + sin(2χ)⋅(A′ + cos(2χ)⋅(B′ + cos(2χ)⋅(C′ + cos(2χ)⋅D′))) (Snyder 3-35) |
| * A′ = A - C |
| * B′ = 2B - 4D |
| * C′ = 4C |
| * D′ = 8D |
| * |
| * Published coefficients are: |
| * |
| * A = 13/360 ⋅ℯ⁸ + 1/12 ⋅ℯ⁶ + 5/24⋅ℯ⁴ + ℯ²/2 |
| * B = 811/11520 ⋅ℯ⁸ + 29/240⋅ℯ⁶ + 7/48⋅ℯ⁴ |
| * C = 81/1120 ⋅ℯ⁸ + 7/120⋅ℯ⁶ |
| * D = 4279/161280⋅ℯ⁸ |
| * |
| * We replace A B C D by A′ B′ C′ D′ using Snyder 3-35. The result are coefficients below. |
| * For each line, we add the smallest values first in order to reduce rounding errors. |
| * The smallest values are the ones using the eccentricity raised to the highest power. |
| */ |
| c2χ = -73./ 2016 * e8 + 1./40 * e6 + 5./24 * e4 + 1./2 * e2; |
| c4χ = 233./ 6720 * e8 + 29./120 * e6 + 7./24 * e4; |
| c6χ = 81./ 280 * e8 + 7./30 * e6; |
| c8χ = 4279./20160 * e8; |
| } |
| |
| /** |
| * Creates a new projection initialized to the values of the given one. This constructor may be invoked after |
| * we determined that the default implementation can be replaced by another one, for example using spherical |
| * formulas instead than the ellipsoidal ones. This constructor allows to transfer all parameters to the new |
| * instance without recomputing them. |
| */ |
| ConformalProjection(final ConformalProjection other) { |
| super(other); |
| useIterations = other.useIterations; |
| c2χ = other.c2χ; |
| c4χ = other.c4χ; |
| c6χ = other.c6χ; |
| c8χ = other.c8χ; |
| } |
| |
| /** |
| * Computes the latitude for a value closely related to the <var>y</var> value of a Mercator projection. |
| * This formula is also part of other projections, since Mercator can be considered as a special case of |
| * Lambert Conic Conformal for instance. |
| * |
| * <p>This function is <em>almost</em> the converse of the {@link #expΨ(double, double)} function. |
| * In a Mercator inverse projection, the value of the {@code rexpΨ} argument is {@code exp(-Ψ)}.</p> |
| * |
| * <p>The input should be a positive number, otherwise the result will be either outside |
| * the [-π/2 … π/2] range, or will be NaN. Its behavior at some particular points is:</p> |
| * |
| * <ul> |
| * <li>φ(0) = π/2</li> |
| * <li>φ(1) = 0</li> |
| * <li>φ(∞) = -π/2.</li> |
| * </ul> |
| * |
| * <b>Note:</b> §1.3.3 in Geomatics Guidance Note number 7 part 2 (April 2015) uses a series expansion |
| * while USGS used an iterative method. The series expansion is twice faster than the iterative method |
| * for the same precision, but this precision is achieved "only" for relatively small eccentricity like |
| * the Earth's one. See the {@code MercatorMethodComparison} class in the test package for more discussion. |
| * |
| * @param rexpΨ the <em>reciprocal</em> of the value returned by {@link #expΨ}. |
| * @return the latitude in radians. |
| * @throws ProjectionException if the iteration does not converge. |
| * |
| * @see #expΨ(double, double) |
| * @see #dy_dφ(double, double) |
| */ |
| final double φ(final double rexpΨ) throws ProjectionException { |
| /* |
| * Get a first approximation of φ from Snyder (7-11). The result below would be exact if the |
| * ellipsoid was actually a sphere. But if the eccentricity is different than 0, then we will |
| * need to add a correction. |
| * |
| * Note that the φ value computed by the line below is called χ in EPSG guide. |
| * We name it φ in our code because we will modify that value in-place in order to get φ. |
| */ |
| double φ = (PI/2) - 2*atan(rexpΨ); // at this point == χ (conformal latitude). |
| /* |
| * Add a correction for the flattened shape of the Earth. The correction can be represented by an |
| * infinite series. Here, we apply only the first 4 terms. Those terms are given by §1.3.3 in the |
| * EPSG guidance note and modified by application of Snyder 3-35. We add those terms in reverse order, |
| * beginning with the smallest values, for reducing rounding errors due to IEEE 754 arithmetic. |
| * |
| * For the original formula as published by EPSG, see MercatorMethodComparison.bySeriesExpansion() |
| * in the test suite. |
| */ |
| final double sin_2φ = sin(2*φ); |
| final double cos_2φ = cos(2*φ); |
| φ += sin_2φ * (c2χ + cos_2φ * (c4χ + cos_2φ * (c6χ + cos_2φ * c8χ))); |
| /* |
| * Note: a previous version checked if the value of the smallest term c8χ⋅sin(8φ) was smaller than |
| * the iteration tolerance. But this was not reliable enough. We use now a hard coded threshold |
| * determined empirically by MercatorMethodComparison. |
| */ |
| if (!useIterations) { |
| return φ; |
| } |
| /* |
| * We should never reach this point for map projections on Earth. But if the ellipsoid is for some |
| * other planet having a high eccentricity, then the above series expansion may not be sufficient. |
| * Try to improve by iteratively solving equation (7-9) from Snyder. However instead than using |
| * Snyder (7-11) as the starting point, we take the result of above calculation as the initial φ. |
| * Assuming that it is closer to the real φ value, this save us some iteration loops and usually |
| * gives us more accurate results (according MercatorMethodComparison tests). |
| */ |
| final double hℯ = 0.5 * eccentricity; |
| for (int it=0; it<MAXIMUM_ITERATIONS; it++) { |
| final double ℯsinφ = eccentricity * sin(φ); |
| final double Δφ = φ - (φ = PI/2 - 2*atan(rexpΨ * pow((1 - ℯsinφ)/(1 + ℯsinφ), hℯ))); |
| if (!(abs(Δφ) > ITERATION_TOLERANCE)) { // Use '!' for accepting NaN. |
| return φ; |
| } |
| } |
| throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence)); |
| } |
| |
| /** |
| * Computes <code>{@linkplain Math#exp(double) exp}(Ψ)</code> where Ψ is the isometric latitude. |
| * This is part of the Mercator projection for the given latitude. This formula is also part of |
| * Lambert Conic Conformal projection, since Mercator can be considered as a special case of that |
| * Lambert projection with the equator as the single standard parallel. |
| * |
| * <p>The isometric latitude is given by the {@linkplain Math#log(double) natural logarithm} of the value returned |
| * by this method. The <em>reciprocal</em> of this function {@code 1/expΨ(φ)} is the converse of {@link #φ(double)}. |
| * Isometric latitude Ψ is related to conformal latitude χ by {@literal χ(φ) = gd(Ψ(φ))} where {@literal gd(x)} is the |
| * Gudermannian function. There is many representations of that function, e.g. {@literal gd(Ψ) = 2⋅atan(expΨ) - π/2}.</p> |
| * |
| * <p>In IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015, |
| * a function closely related to this method has the letter <var>t</var>.</p> |
| * |
| * |
| * <h4>Properties</h4> |
| * This function is used with φ values in the [-π/2 … π/2] range and has a periodicity of 2π. |
| * The result is always a positive number when the φ argument is inside the above-cited range. |
| * If, after removal of any 2π periodicity, φ is still outside the [-π/2 … π/2] range, then the |
| * result is a negative number. In a Mercator projection, such negative number will result in NaN. |
| * |
| * <p>Some values are:</p> |
| * <ul> |
| * <li>expΨ(NaN) = NaN</li> |
| * <li>expΨ(±∞) = NaN</li> |
| * <li>expΨ(-π/2) = 0</li> |
| * <li>expΨ( 0 ) = 1</li> |
| * <li>expΨ(+π/2) → ∞ (actually some large value like 1.633E+16)</li> |
| * <li>expΨ(-φ) = 1 / expΨ(φ)</li> |
| * </ul> |
| * |
| * |
| * <h4>The π/2 special case</h4> |
| * The value at {@code Math.PI/2} is not exactly infinity because there is no exact representation of π/2. |
| * However since the conversion of 90° to radians gives {@code Math.PI/2}, we can presume that the user was |
| * expecting infinity. The caller should check for the PI/2 special case himself if desired, as this method |
| * does nothing special about it. |
| * |
| * <p>Note that the result for the φ value after {@code Math.PI/2} (as given by {@link Math#nextUp(double)}) |
| * is still positive, maybe because {@literal PI/2 < π/2 < nextUp(PI/2)}. Only the {@code nextUp(nextUp(PI/2))} |
| * value become negative. Callers may need to take this behavior in account: special check for {@code Math.PI/2} |
| * is not sufficient, the check needs to include at least the {@code nextUp(Math.PI/2)} case.</p> |
| * |
| * |
| * <h4>Relationship with Snyder</h4> |
| * This function is related to the following functions from Snyder: |
| * |
| * <ul> |
| * <li>(7-7) in the <cite>Mercator projection</cite> chapter.</li> |
| * <li>Reciprocal of (9-13) in the <cite>Oblique Mercator projection</cite> chapter.</li> |
| * <li>Reciprocal of (15-9) in the <cite>Lambert Conformal Conic projection</cite> chapter.</li> |
| * </ul> |
| * |
| * @param φ the latitude in radians. |
| * @param ℯsinφ the sine of the φ argument multiplied by {@link #eccentricity}. |
| * @return {@code Math.exp(Ψ)} where Ψ is the isometric latitude. |
| * |
| * @see #φ(double) |
| * @see #dy_dφ(double, double) |
| * @see <a href="https://en.wikipedia.org/wiki/Latitude#Isometric_latitude">Isometric latitude on Wikipedia</a> |
| */ |
| final double expΨ(final double φ, final double ℯsinφ) { |
| /* |
| * Note: tan(π/4 - φ/2) = 1 / tan(π/4 + φ/2) |
| * |
| * A + sign in the equation favorises slightly the accuracy in South hemisphere, while a - sign |
| * favorises slightly the North hemisphere (but the differences are very small). In Apache SIS, |
| * we handle that by changing the sign of some terms in the (de)normalisation matrices. |
| */ |
| return tan(PI/4 + 0.5*φ) * pow((1 - ℯsinφ) / (1 + ℯsinφ), 0.5*eccentricity); |
| } |
| |
| /** |
| * Computes the partial derivative of a Mercator projection at the given latitude. This formula is also part of |
| * other projections, since Mercator can be considered as a special case of Lambert Conic Conformal for instance. |
| * |
| * <p>In order to get the derivative of the {@link #expΨ(double, double)} function, caller can multiply |
| * the returned value by {@code expΨ}.</p> |
| * |
| * @param sinφ the sine of latitude. |
| * @param cosφ the cosine of latitude. |
| * @return the partial derivative of a Mercator projection at the given latitude. |
| * |
| * @see #expΨ(double, double) |
| * @see #φ(double) |
| */ |
| final double dy_dφ(final double sinφ, final double cosφ) { |
| return (1 / cosφ) - eccentricitySquared * cosφ / (1 - eccentricitySquared * (sinφ*sinφ)); |
| } |
| } |
