You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2022/08/20 14:27:27 UTC
[sis] 02/02: Partially revert commit 3fa0138 ("Apply a longitude wraparound on the Mercator projection"). Instead, the Mercator projection applies no wraparound by default (same behavior as SIS 1.2) and the wraparound can be requested by invoking `createWithWraparound(…)`.
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
commit fa0971363889cfea057aae4ae258273b56d3da24
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Sat Jul 30 11:31:54 2022 +0200
Partially revert commit 3fa0138 ("Apply a longitude wraparound on the Mercator projection").
Instead, the Mercator projection applies no wraparound by default (same behavior as SIS 1.2)
and the wraparound can be requested by invoking `createWithWraparound(…)`.
https://issues.apache.org/jira/browse/SIS-547
https://issues.apache.org/jira/browse/SIS-486
---
.../operation/projection/AlbersEqualArea.java | 31 +-
.../operation/projection/AuthalicMercator.java | 19 +-
.../operation/projection/Initializer.java | 37 +--
.../projection/LambertConicConformal.java | 27 +-
.../operation/projection/LongitudeWraparound.java | 347 +++++++++++++++++++++
.../referencing/operation/projection/Mercator.java | 90 +++---
.../operation/projection/Mollweide.java | 6 +-
.../operation/projection/NormalizedProjection.java | 114 +++++--
.../operation/projection/ObliqueStereographic.java | 39 +--
.../operation/projection/SatelliteTracking.java | 40 ++-
.../operation/projection/TransverseMercator.java | 2 +-
.../operation/projection/ZonedGridSystem.java | 2 +-
.../projection/LambertConicConformalTest.java | 2 +-
.../operation/projection/MercatorTest.java | 8 +-
.../projection/PolarStereographicTest.java | 2 +-
.../projection/TransverseMercatorTest.java | 2 +-
16 files changed, 553 insertions(+), 215 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
index fc1e246ffb..3267480e3d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -50,7 +50,7 @@ import static org.apache.sis.internal.referencing.provider.AlbersEqualArea.*;
*
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.2
+ * @version 1.3
* @since 0.8
* @module
*/
@@ -58,7 +58,7 @@ public class AlbersEqualArea extends AuthalicConversion {
/**
* For cross-version compatibility.
*/
- private static final long serialVersionUID = -3466040922402982480L;
+ private static final long serialVersionUID = -3024658742514888646L;
/**
* Internal coefficients for computation, depending only on eccentricity and values of standards parallels.
@@ -81,15 +81,6 @@ public class AlbersEqualArea extends AuthalicConversion {
*/
final double C;
- /**
- * A bound of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ values.
- * Some (not all) θ values need to be shifted inside that range before to use them
- * in trigonometric functions.
- *
- * @see Initializer#boundOfScaledLongitude(DoubleDouble)
- */
- final double θ_bound;
-
/**
* Creates an Albers Equal Area projection from the given parameters.
*
@@ -173,7 +164,6 @@ public class AlbersEqualArea extends AuthalicConversion {
denormalize.convertBefore(0, rn, null); rn.negate();
denormalize.convertBefore(1, rn, ρ0); rn.inverseDivide(-1);
normalize.convertAfter(0, rn, null); // On this line, `rn` became `n`.
- θ_bound = initializer.boundOfScaledLongitude(rn);
}
/**
@@ -181,9 +171,8 @@ public class AlbersEqualArea extends AuthalicConversion {
*/
AlbersEqualArea(final AlbersEqualArea other) {
super(other);
- nm = other.nm;
- C = other.C;
- θ_bound = other.θ_bound;
+ nm = other.nm;
+ C = other.C;
}
/**
@@ -222,7 +211,7 @@ public class AlbersEqualArea extends AuthalicConversion {
if (eccentricity == 0 && getClass() == AlbersEqualArea.class) {
kernel = new Spherical(this);
}
- return context.completeTransform(factory, kernel);
+ return kernel.completeWithWraparound(factory);
}
/**
@@ -240,9 +229,8 @@ public class AlbersEqualArea extends AuthalicConversion {
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- // θ = n⋅λ reduced to [−n⋅π … n⋅π] range.
- final double θ = wraparoundScaledLongitude(srcPts[srcOff], θ_bound);
- final double φ = srcPts[srcOff + 1];
+ final double θ = srcPts[srcOff ]; // θ = n⋅λ reduced to [−n⋅π … n⋅π] range.
+ final double φ = srcPts[srcOff+1];
final double cosθ = cos(θ);
final double sinθ = sin(θ);
final double sinφ = sin(φ);
@@ -329,9 +317,8 @@ public class AlbersEqualArea extends AuthalicConversion {
final double[] dstPts, final int dstOff,
final boolean derivate)
{
- // θ = n⋅λ reduced to [−n⋅π … n⋅π] range.
- final double θ = wraparoundScaledLongitude(srcPts[srcOff], θ_bound);
- final double φ = srcPts[srcOff + 1];
+ final double θ = srcPts[srcOff ]; // θ = n⋅λ reduced to [−n⋅π … n⋅π] range.
+ final double φ = srcPts[srcOff+1];
final double cosθ = cos(θ);
final double sinθ = sin(θ);
final double sinφ = sin(φ);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicMercator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicMercator.java
index 656c64a4c2..7eff535735 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicMercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicMercator.java
@@ -16,9 +16,13 @@
*/
package org.apache.sis.referencing.operation.projection;
+import java.util.Optional;
+import org.opengis.geometry.Envelope;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.referencing.operation.matrix.Matrix2;
+import org.apache.sis.referencing.operation.transform.DomainDefinition;
+import org.apache.sis.geometry.Envelope2D;
import static java.lang.Math.*;
import static org.apache.sis.math.MathFunctions.atanh;
@@ -29,7 +33,7 @@ import static org.apache.sis.math.MathFunctions.atanh;
* This is used for implementation of ESRI "Mercator Auxiliary Sphere type 3" projection.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.2
+ * @version 1.3
* @since 1.2
* @module
*/
@@ -53,6 +57,19 @@ final class AuthalicMercator extends AuthalicConversion {
super(null, other);
}
+ /**
+ * Returns the domain of input coordinates.
+ * This method is defined for consistency with {@link Mercator#getDomain(DomainDefinition)}.
+ *
+ * @since 1.3
+ */
+ @Override
+ public Optional<Envelope> getDomain(final DomainDefinition criteria) {
+ return Optional.of(new Envelope2D(null,
+ -LARGE_LONGITUDE_LIMIT, -POLAR_AREA_LIMIT,
+ 2*LARGE_LONGITUDE_LIMIT, 2*POLAR_AREA_LIMIT));
+ }
+
/**
* Projects the specified (λ,φ) coordinates (units in radians) and stores the result in {@code dstPts}.
* In addition, opportunistically computes the projection derivative if {@code derivate} is {@code true}.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
index b3551e0da6..a8107cc6bd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
@@ -56,7 +56,7 @@ import static org.apache.sis.util.ArgumentChecks.ensureNonNull;
*
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.2
+ * @version 1.3
* @since 0.6
* @module
*/
@@ -384,39 +384,4 @@ final class Initializer {
s.inverseDivide(cosφ);
return s.doubleValue();
}
-
- /**
- * Returns a bound of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ values.
- * This method is invoked by map projections that multiply the longitude values by some scale factor before
- * to use them in trigonometric functions. Usually we do not explicitly wraparound the longitude values,
- * because trigonometric functions do that automatically for us. However if the longitude is multiplied
- * by some factor before to be used in trigonometric functions, then that implicit wraparound is not the
- * one we expect. The map projection code needs to perform explicit wraparound in such cases.
- *
- * @param n the factor by which longitude values are multiplied before use in trigonometry.
- * @return a bound of the [−n⋅π … n⋅π] range.
- *
- * @see NormalizedProjection#wraparoundScaledLongitude(double, double)
- * @see <a href="https://issues.apache.org/jira/browse/SIS-486">SIS-486</a>
- */
- final double boundOfScaledLongitude(final double n) {
- return boundOfScaledLongitude(new DoubleDouble(n));
- }
-
- /**
- * Same as {@link #boundOfScaledLongitude(double)} with opportunistic use of double-double precision.
- * This is used when than object is available anyway.
- *
- * @param n the factor by which longitude values are multiplied before use in trigonometry.
- * @return a bound of the [−n⋅π … n⋅π] range.
- */
- final double boundOfScaledLongitude(final DoubleDouble n) {
- if (signum_λ0 == 0 || n.doubleValue() >= 1) {
- return Double.NaN; // Do not apply any wraparound.
- }
- final DoubleDouble r = DoubleDouble.createPi();
- r.multiply(n);
- final double θ_bound = abs(r.doubleValue());
- return (signum_λ0 < 0) ? θ_bound : -θ_bound;
- }
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
index 911a4dd001..581e84ea01 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
@@ -74,7 +74,7 @@ public class LambertConicConformal extends ConformalProjection {
/**
* For cross-version compatibility.
*/
- private static final long serialVersionUID = -8786460743797422415L;
+ private static final long serialVersionUID = 2067358524298002016L;
/**
* Variants of Lambert Conical Conformal projection. Those variants modify the way the projections are constructed
@@ -165,15 +165,6 @@ public class LambertConicConformal extends ConformalProjection {
*/
final double n;
- /**
- * A bound of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ values.
- * Some (not all) θ values need to be shifted inside that range before to use them
- * in trigonometric functions.
- *
- * @see Initializer#boundOfScaledLongitude(DoubleDouble)
- */
- final double θ_bound;
-
/**
* Creates a Lambert projection from the given parameters.
* The {@code method} argument can be the description of one of the following:
@@ -375,7 +366,6 @@ public class LambertConicConformal extends ConformalProjection {
normalize .convertAfter(1, sφ, null);
denormalize.convertBefore(0, F, null); F.negate();
denormalize.convertBefore(1, F, rF);
- θ_bound = initializer.boundOfScaledLongitude(sλ);
}
/**
@@ -383,8 +373,7 @@ public class LambertConicConformal extends ConformalProjection {
*/
LambertConicConformal(final LambertConicConformal other) {
super(other);
- n = other.n;
- θ_bound = other.θ_bound;
+ n = other.n;
}
/**
@@ -423,7 +412,7 @@ public class LambertConicConformal extends ConformalProjection {
if (eccentricity == 0 && getClass() == LambertConicConformal.class) {
kernel = new Spherical(this);
}
- return context.completeTransform(factory, kernel);
+ return kernel.completeWithWraparound(factory);
}
/**
@@ -436,7 +425,7 @@ public class LambertConicConformal extends ConformalProjection {
*/
@Override
public Optional<Envelope> getDomain(final DomainDefinition criteria) {
- final double x = abs(θ_bound);
+ final double x = getWraparoundLongitude();
final double y = copySign(PI/2, -n);
return Optional.of(new Envelope2D(null, -x, Math.min(y, 0), 2*x, Math.abs(y)));
}
@@ -461,8 +450,8 @@ public class LambertConicConformal extends ConformalProjection {
* the first non-linear one moved to the "normalize" affine transform, and the linear operations
* applied after the last non-linear one moved to the "denormalize" affine transform.
*/
- final double θ = wraparoundScaledLongitude(srcPts[srcOff], θ_bound); // θ = Δλ⋅n
- final double φ = srcPts[srcOff + 1]; // Sign may be reversed
+ final double θ = srcPts[srcOff ]; // θ = Δλ⋅n
+ final double φ = srcPts[srcOff+1]; // Sign may be reversed
final double absφ = abs(φ);
final double sinθ = sin(θ);
final double cosθ = cos(θ);
@@ -572,8 +561,8 @@ public class LambertConicConformal extends ConformalProjection {
final double[] dstPts, final int dstOff,
final boolean derivate)
{
- final double θ = wraparoundScaledLongitude(srcPts[srcOff], θ_bound); // θ = Δλ⋅n
- final double φ = srcPts[srcOff + 1]; // Sign may be reversed
+ final double θ = srcPts[srcOff ]; // θ = Δλ⋅n
+ final double φ = srcPts[srcOff+1]; // Sign may be reversed
final double absφ = abs(φ);
final double sinθ = sin(θ);
final double cosθ = cos(θ);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java
new file mode 100644
index 0000000000..98da9b1472
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LongitudeWraparound.java
@@ -0,0 +1,347 @@
+/*
+ * 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.Optional;
+import java.io.Serializable;
+import org.opengis.geometry.Envelope;
+import org.opengis.parameter.ParameterValueGroup;
+import org.opengis.parameter.ParameterDescriptorGroup;
+import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.MathTransform2D;
+import org.opengis.referencing.operation.NoninvertibleTransformException;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.transform.AbstractMathTransform2D;
+import org.apache.sis.referencing.operation.transform.ContextualParameters;
+import org.apache.sis.referencing.operation.transform.DomainDefinition;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
+import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.internal.util.DoubleDouble;
+import org.apache.sis.internal.util.Numerics;
+import org.apache.sis.util.ComparisonMode;
+import org.apache.sis.measure.Longitude;
+
+
+/**
+ * If the scaled longitude θ=n⋅λ is outside the [−n⋅π … n⋅π] range, maybe shifts θ to that range.
+ * This transform intentionally does <strong>not</strong> force θ to be inside that range in all cases.
+ * We avoid explicit wraparounds as much as possible (as opposed to implicit wraparounds performed by
+ * trigonometric functions) because they tend to introduce discontinuities. We perform wraparounds only
+ * when necessary for the problem of area crossing the anti-meridian (±180°).
+ *
+ * <div class="note"><b>Example:</b>
+ * a CRS for Alaska may have the central meridian at λ₀=−154° of longitude. If the point to project is
+ * at λ=177° of longitude, calculations will be performed with Δλ=331° while the correct value that we
+ * need to use is Δλ=−29°.</div>
+ *
+ * In order to avoid wraparound operations as much as possible, we test only the bound where anti-meridian
+ * problem may happen; no wraparound will be applied for the opposite bound. Furthermore we add or subtract
+ * 360° only once. Even if the point did many turns around the Earth, the 360° shift will still be applied
+ * at most once. The desire to apply the minimal amount of shifts is the reason why we do not use
+ * {@link Math#IEEEremainder(double, double)}.
+ *
+ * <h2>When to use</h2>
+ * Many map projections implicitly wraparound longitude values through the use of trigonometric functions
+ * ({@code sin(λ)}, {@code cos(λ)}, <i>etc</i>). For those map projections, the wraparound is unconditional
+ * and this {@code LongitudeWraparound} class is not needed. This class is used only when the wraparound is
+ * not implicitly done and the central meridian is not zero. The latter condition is because subtraction of
+ * central meridian may cause longitude values to go outside the −180° … +180° range.
+ *
+ * <p>This transform is hidden in WKT (it does not appear as a concatenation).</p>
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ *
+ * @see org.apache.sis.referencing.operation.transform.WraparoundTransform
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-486">SIS-486</a>
+ *
+ * @since 1.3
+ * @module
+ */
+final class LongitudeWraparound extends AbstractMathTransform2D implements Serializable {
+ /**
+ * For cross-version compatibility.
+ */
+ private static final long serialVersionUID = -4658152274068444690L;
+
+ /**
+ * The actual map projection to execute after the longitude wraparound.
+ */
+ final NormalizedProjection projection;
+
+ /**
+ * A bound of the [−n⋅π … n⋅π] range which, if exceeded, should cause wraparound.
+ * Some (not all) θ = n⋅λ values need to be shifted inside that range before to
+ * use them in trigonometric functions.
+ *
+ * <p>The sign is significant. A negative value means that the wraparound is applied
+ * only on longitudes less than −180°. A positive value means that the wraparound is
+ * applied only on longitudes greater than +180°.</p>
+ */
+ final double bound;
+
+ /**
+ * Whether the bound is on the side of negative longitudes. This is {@code bound < 0}.
+ * This field is trivial in {@code LongitudeWraparound} case, but more important in
+ * {@link Inverse} case because {@code rotation - bound} may cancel to zero.
+ */
+ final boolean negative;
+
+ /**
+ * The inverse of this operation.
+ *
+ * @see #inverse()
+ */
+ private final Inverse inverse;
+
+ /**
+ * Creates a new transform for wrapping around the longitude values before a map projection.
+ *
+ * @param projection the actual map projection to execute after the longitude wraparound.
+ * @param bound one bound of the [−n⋅π … n⋅π] range, on the side where wraparound needs to be applied.
+ * @param rotation longitude rotation applied by the normalization matrix after conversion to projection domain.
+ */
+ LongitudeWraparound(final NormalizedProjection projection, final double bound, final double rotation) {
+ this.projection = projection;
+ this.bound = bound;
+ negative = bound < 0;
+ inverse = new Inverse(this, rotation - bound);
+ }
+
+ /**
+ * Returns a bound of the [−n⋅π … n⋅π] domain where <var>n</var> is a map projection dependent factor.
+ * The factor is inferred from the {@link NormalizedProjection#context}.
+ *
+ * @param normalize the normalization matrix of the projection for which to get a bound.
+ * @param negative {@code true} for the −180° bound, or {@code false} for the +180° bound.
+ * @return a bound of the [−n⋅π … n⋅π] range.
+ */
+ static double boundOfScaledLongitude(final MatrixSIS normalize, final boolean negative) {
+ final DoubleDouble bound = DoubleDouble.castOrCopy(normalize.getNumber(0, 0));
+ bound.multiply(negative ? Longitude.MIN_VALUE : Longitude.MAX_VALUE);
+ return bound.doubleValue();
+ }
+
+ /**
+ * Returns the ranges of coordinate values which can be used as inputs.
+ * The {@link #projection} domain is used verbatim, without wraparound adjustment.
+ *
+ * @param criteria controls the definition of transform domain.
+ * @return estimation of a domain where this transform is considered numerically applicable.
+ */
+ @Override
+ public Optional<Envelope> getDomain(DomainDefinition criteria) throws TransformException {
+ return projection.getDomain(criteria);
+ }
+
+ /**
+ * Returns the parameter descriptors for this math transform, or {@code null} if unknown.
+ * Delegates to {@link #projection} since this {@code LongitudeWraparound} is hidden.
+ * This is used by WKT formatting.
+ */
+ @Override
+ public ParameterDescriptorGroup getParameterDescriptors() {
+ return projection.getParameterDescriptors();
+ }
+
+ /**
+ * Returns the parameter values for this math transform, or {@code null} if unknown.
+ * Delegates to {@link #projection} since this {@code LongitudeWraparound} is hidden.
+ * This is used by WKT formatting.
+ */
+ @Override
+ public ParameterValueGroup getParameterValues() {
+ return projection.getParameterValues();
+ }
+
+ /**
+ * Returns the parameters for a sequence of <cite>normalize</cite> → {@code this} → <cite>denormalize</cite>.
+ * Delegates to {@link #projection} since this {@code LongitudeWraparound} is hidden.
+ * This is used by WKT formatting.
+ */
+ @Override
+ protected ContextualParameters getContextualParameters() {
+ return projection.getContextualParameters();
+ }
+
+ /**
+ * Transforms a single coordinate tuple in an array, and optionally computes the transform derivative.
+ * The wraparound is applied, if needed, on the longitude value before to delegate to {@link #projection}.
+ */
+ @Override
+ public Matrix transform(final double[] srcPts, final int srcOff,
+ final double[] dstPts, final int dstOff, boolean derivate) throws TransformException
+ {
+ final double λ = srcPts[srcOff];
+ if (negative ? λ < bound : λ > bound) {
+ dstPts[dstOff+1] = srcPts[srcOff+1]; // Must be first.
+ dstPts[dstOff ] = λ - 2*bound;
+ return projection.transform(dstPts, dstOff, dstPts, dstOff, derivate);
+ } else {
+ return projection.transform(srcPts, srcOff, dstPts, dstOff, derivate);
+ }
+ }
+
+ /**
+ * Transforms a list of coordinate tuples. This method is provided for efficiently transforming many points.
+ * Wraparound is applied on all longitude values before to delegate to {@link #projection}.
+ */
+ @Override
+ public void transform(final double[] srcPts, final int srcOff,
+ final double[] dstPts, final int dstOff, final int numPts) throws TransformException
+ {
+ if (srcPts != dstPts || srcOff != dstOff) {
+ System.arraycopy(srcPts, srcOff, dstPts, dstOff, numPts * DIMENSION);
+ }
+ final double period = 2*bound;
+ final int stop = dstOff + numPts * DIMENSION;
+ for (int i=dstOff; i<stop; i+=DIMENSION) {
+ final double λ = dstPts[i];
+ if (negative ? λ < bound : λ > bound) {
+ dstPts[i] = λ - period;
+ }
+ }
+ projection.transform(dstPts, dstOff, dstPts, dstOff, numPts);
+ }
+
+ /**
+ * Returns the inverse transform of this object.
+ */
+ @Override
+ public MathTransform2D inverse() throws NoninvertibleTransformException {
+ return inverse;
+ }
+
+ /**
+ * Longitude wraparound applied on reverse projection.
+ * This is a copy of {@code NormalizedProjection.Inverse} with longitude wraparound added after conversion.
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.3
+ * @since 1.3
+ * @module
+ */
+ private static final class Inverse extends AbstractMathTransform2D.Inverse implements Serializable {
+ /**
+ * For cross-version compatibility.
+ */
+ private static final long serialVersionUID = -543869926271003589L;
+
+ /**
+ * The projection to reverse, which is the enclosing transform.
+ */
+ private final LongitudeWraparound forward;
+
+ /**
+ * {@link LongitudeWraparound#bound} with opposite sign and translated by the longitude rotation.
+ * This is the bound that the reverse projection needs to use before longitude rotation.
+ * This bound <strong>shall not</strong> be used for period computation.
+ */
+ private final double bound;
+
+ /**
+ * Creates a reverse projection for the given forward projection.
+ */
+ Inverse(final LongitudeWraparound forward, final double bound) {
+ this.forward = forward;
+ this.bound = bound;
+ }
+
+ /**
+ * Returns the inverse of this math transform, which is the forward projection.
+ */
+ @Override
+ public MathTransform2D inverse() {
+ return forward;
+ }
+
+ /**
+ * Reverse projects the specified {@code srcPts} and stores the result in {@code dstPts}.
+ * If the derivative has been requested, then this method will delegate the derivative
+ * calculation to the enclosing class and inverts the resulting matrix.
+ */
+ @Override
+ public Matrix transform(final double[] srcPts, final int srcOff,
+ double[] dstPts, int dstOff,
+ final boolean derivate) throws TransformException
+ {
+ if (derivate && dstPts == null) {
+ dstPts = new double[DIMENSION];
+ dstOff = 0;
+ }
+ forward.projection.inverseTransform(srcPts, srcOff, dstPts, dstOff);
+ final Matrix d = derivate ? Matrices.inverse(forward.transform(dstPts, dstOff, null, 0, true)) : null;
+ final double λ = dstPts[dstOff];
+ if (forward.negative ? λ > bound : λ < bound) { // Interpretation of `negative` is inversed.
+ dstPts[dstOff] = λ + 2*forward.bound;
+ }
+ return d;
+ }
+
+ /**
+ * Inverse transforms an arbitrary amount of coordinates. This method optimizes the
+ * case where conversions can be applied by a loop with indices in increasing order.
+ */
+ @Override
+ public void transform(final double[] srcPts, int srcOff,
+ final double[] dstPts, int dstOff, int numPts) throws TransformException
+ {
+ if (srcPts == dstPts && srcOff < dstOff) {
+ super.transform(srcPts, srcOff, dstPts, dstOff, numPts);
+ } else {
+ final double period = 2 * forward.bound;
+ while (--numPts >= 0) {
+ forward.projection.inverseTransform(srcPts, srcOff, dstPts, dstOff);
+ final double λ = dstPts[dstOff];
+ if (forward.negative ? λ > bound : λ < bound) { // Interpretation of `negative` is inversed.
+ dstPts[dstOff] = λ + period;
+ }
+ srcOff += DIMENSION;
+ dstOff += DIMENSION;
+ }
+ }
+ }
+ }
+
+ /*
+ * We do not implement `tryConcatenate` yet because the result of invoking `projection.tryConcatenate(…)`
+ * is either null or a linear transform. In the latter case, the linear transform can not be wrapped by
+ * this `longitudeWraparound` class. Even if we could, it would block concatenation with surrounding
+ * affine transforms. We have no easy solution for now.
+ */
+
+ /**
+ * Computes a hash value for this transform.
+ */
+ @Override
+ protected int computeHashCode() {
+ return projection.hashCode() + Double.hashCode(bound);
+ }
+
+ /**
+ * Compares the specified object with this math transform for equality.
+ */
+ @Override
+ public boolean equals(final Object object, final ComparisonMode mode) {
+ if (object instanceof LongitudeWraparound) {
+ final LongitudeWraparound that = (LongitudeWraparound) object;
+ return Numerics.epsilonEqual(bound, that.bound, mode)
+ && projection.equals(that.projection, mode);
+ }
+ return false;
+ }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
index 1059bd8a6a..d6be730fbb 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
@@ -195,12 +195,11 @@ public class Mercator extends ConformalProjection {
final Variant variant = variant(method, Variant.values(), Variant.TWO_PARALLELS);
final EnumMap<ParameterRole, ParameterDescriptor<Double>> roles = new EnumMap<>(ParameterRole.class);
/*
- * "Longitude of origin" is a parameter of all Mercator projections, but is intentionally omitted from
- * this map because it will be handled in a special way by the Mercator constructor. The "scale factor"
- * is not formally a "Mercator 2SP" argument, but we accept it anyway for all Mercator projections
- * since it may be used in some Well Known Text (WKT).
+ * The "scale factor" is not formally a "Mercator 2SP" argument, but we accept it anyway
+ * for all Mercator projections because it may be used in some Well Known Text (WKT).
*/
- roles.put(ParameterRole.SCALE_FACTOR, Mercator1SP.SCALE_FACTOR);
+ roles.put(ParameterRole.SCALE_FACTOR, Mercator1SP.SCALE_FACTOR);
+ roles.put(ParameterRole.CENTRAL_MERIDIAN, Mercator1SP.LONGITUDE_OF_ORIGIN);
switch (variant) {
case REGIONAL: {
roles.put(ParameterRole.FALSE_EASTING, RegionalMercator.EASTING_AT_FALSE_ORIGIN);
@@ -270,31 +269,10 @@ public class Mercator extends ConformalProjection {
* if they really want, since we sometime see such CRS definitions.
*/
final double φ1 = toRadians(initializer.getAndStore(Mercator2SP.STANDARD_PARALLEL));
- final DoubleDouble k0 = new DoubleDouble(initializer.scaleAtφ(sin(φ1), cos(φ1)));
- final DoubleDouble ku = DoubleDouble.createAndGuessError(PI/0.5);
- ku.multiply(k0);
- /*
- * Note about `ku`: in principle we should convert longitude values to radians, then apply the same
- * scale factor k₀ to both longitude and latitude values. However in this implementation instead of
- * converting the [−180 … 180]° range to the [−π … π] range, we rather convert to a [−½ … ½] range.
- * This is because longitude conversions in Mercator projection do not use trigonometric functions,
- * so we do not get the implicit wraparound performed by trigonometric functions such as sin(λ−λ₀).
- * We will use `Math.rint(…)` instead.
- */
+ final Number k0 = new DoubleDouble(initializer.scaleAtφ(sin(φ1), cos(φ1)));
final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
- final DoubleDouble forUnity = DoubleDouble.createAndGuessError(1.0/360);
- normalize.setNumber(0, 0, forUnity); // Overwrite the "degrees to radians" coefficient.
- if (λ0 != 0) {
- /*
- * Use double-double arithmetic here for consistency with the work done in the normalization matrix.
- * The intent is to have exact value at `double` precision when computing Matrix.invert(). Note that
- * there is no such goal for other parameters computed from sine or consine functions.
- */
- forUnity.multiplyGuessError(-λ0);
- normalize.setNumber(0, 2, forUnity);
- }
- denormalize.convertBefore(0, ku, null);
+ denormalize.convertBefore(0, k0, null);
denormalize.convertBefore(1, k0, null);
if (φ0 != 0) {
denormalize.convertBefore(1, null, new DoubleDouble(-log(expΨ(φ0, eccentricity * sin(φ0)))));
@@ -391,7 +369,7 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
}
kernel = new Spherical(this);
}
- return context.completeTransform(factory, kernel);
+ return kernel.completeWithWraparound(factory);
}
/**
@@ -399,12 +377,17 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
* the limit is arbitrarily set to 84° of latitude North and South. This is consistent with the
* "World Mercator" domain of validity defined by EPSG:3395, which is 80°S to 84°N.
*
+ * <p>The range of longitude values is set to an arbitrary range larger than −180° … +180°,
+ * because the Mercator projection is mathematically capable to handle coordinates beyond that range
+ * even if those coordinates have no real world meaning. This expansion can facilitate the projection
+ * of envelopes, geometries or rasters.</p>
+ *
* @since 1.3
*/
@Override
public Optional<Envelope> getDomain(final DomainDefinition criteria) {
- final double limit = (variant == Variant.MILLER) ? -PI/2 : -PI/2 * (84d/90);
- return Optional.of(new Envelope2D(null, -0.5, limit, 1, -2*limit));
+ final double limit = (variant == Variant.MILLER) ? -PI/2 : -POLAR_AREA_LIMIT;
+ return Optional.of(new Envelope2D(null, -LARGE_LONGITUDE_LIMIT, limit, 2*LARGE_LONGITUDE_LIMIT, -2*limit));
}
/**
@@ -421,7 +404,6 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- final double x = srcPts[srcOff ];
final double φ = srcPts[srcOff+1];
final double sinφ = sin(φ);
if (dstPts != null) {
@@ -447,7 +429,7 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
y = NaN;
}
}
- dstPts[dstOff ] = x - rint(x); // Scale will be applied by the denormalization matrix.
+ dstPts[dstOff ] = srcPts[srcOff]; // Scale will be applied by the denormalization matrix.
dstPts[dstOff+1] = y;
}
/*
@@ -467,21 +449,24 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
final double[] dstPts, int dstOff, int numPts)
throws TransformException
{
- if ((srcPts == dstPts && srcOff < dstOff) || getClass() != Mercator.class) {
+ if (srcPts != dstPts || srcOff != dstOff || getClass() != Mercator.class) {
super.transform(srcPts, srcOff, dstPts, dstOff, numPts);
} else {
+ /*
+ * Override the super-class method only as an optimization in the special case where the target coordinates
+ * are written at the same locations than the source coordinates. In such case, we can take advantage of
+ * the fact that the λ values are not modified by the normalized Mercator projection.
+ */
+ dstOff--;
while (--numPts >= 0) {
- final double x = srcPts[srcOff++];
- final double φ = srcPts[srcOff++];
- final double y;
- if (φ == 0) {
- y = φ;
- } else {
+ final double φ = dstPts[dstOff += DIMENSION]; // Same as srcPts[srcOff + 1].
+ if (φ != 0) {
/*
* See the javadoc of the Spherical inner class for a note
* about why we perform explicit checks for the pole cases.
*/
final double a = abs(φ);
+ final double y;
if (a < PI/2) {
y = log(expΨ(φ, eccentricity * sin(φ)));
} else if (a <= (PI/2 + ANGULAR_TOLERANCE)) {
@@ -489,9 +474,8 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
} else {
y = NaN;
}
+ dstPts[dstOff] = y;
}
- dstPts[dstOff++] = x - rint(x);
- dstPts[dstOff++] = y;
}
}
}
@@ -561,7 +545,6 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
final double[] dstPts, final int dstOff,
final boolean derivate)
{
- final double x = srcPts[srcOff ];
final double φ = srcPts[srcOff+1];
if (dstPts != null) {
/*
@@ -583,7 +566,7 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
y = NaN;
}
}
- dstPts[dstOff ] = x - rint(x);
+ dstPts[dstOff ] = srcPts[srcOff];
dstPts[dstOff+1] = y;
}
return derivate ? new Matrix2(1, 0, 0, 1/cos(φ)) : null;
@@ -601,18 +584,16 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
final double[] dstPts, int dstOff, int numPts)
throws TransformException
{
- if (srcPts == dstPts && srcOff < dstOff) {
+ if (srcPts != dstPts || srcOff != dstOff) {
super.transform(srcPts, srcOff, dstPts, dstOff, numPts);
} else {
+ dstOff--;
while (--numPts >= 0) {
- final double x = srcPts[srcOff++];
- final double φ = srcPts[srcOff++];
- final double y;
- if (φ == 0) {
- y = φ;
- } else {
+ final double φ = dstPts[dstOff += DIMENSION]; // Same as srcPts[srcOff + 1].
+ if (φ != 0) {
// See class javadoc for a note about explicit check for poles.
final double a = abs(φ);
+ final double y;
if (a < PI/2) {
y = log(tan(PI/4 + 0.5*φ)); // Part of Snyder (7-2)
} else if (a <= (PI/2 + ANGULAR_TOLERANCE)) {
@@ -620,9 +601,8 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
} else {
y = NaN;
}
+ dstPts[dstOff] = y;
}
- dstPts[dstOff++] = x - rint(x);
- dstPts[dstOff++] = y;
}
}
}
@@ -644,7 +624,7 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
* Invoked when {@link #tryConcatenate(boolean, MathTransform, MathTransformFactory)} detected a
* (inverse) → (affine) → (this) transforms sequence. If the affine transform in the middle does
* not change the latitude value, then we can take advantage of the fact that longitude conversion
- * is linear (ignoring wraparound).
+ * is linear.
*/
@Override
final MathTransform tryConcatenate(boolean projectedSpace, Matrix affine, MathTransformFactory factory)
@@ -653,7 +633,7 @@ subst: if (variant.spherical || (eccentricity == 0 && getClass() == Mercator.cl
/*
* Verify that the latitude row is an identity conversion except for the sign which is allowed to change
* (but no scale and no translation are allowed). Ignore the longitude row because it just pass through
- * this Mercator projection with no impact on any calculation (if we ignore the wraparound).
+ * this Mercator projection with no impact on any calculation.
*/
if (affine.getElement(1,0) == 0 && affine.getElement(1, DIMENSION) == 0 && Math.abs(affine.getElement(1,1)) == 1) {
if (factory != null) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
index aeb3fa02e5..17a14be96a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
@@ -132,10 +132,10 @@ public class Mollweide extends NormalizedProjection {
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- final double λ = srcPts[srcOff];
+ final double λ = srcPts[srcOff]; // Scaled by 2√2.
final double φ = srcPts[srcOff + 1];
final double sinφ = sin(φ);
- double θp = 2 * asin(φ / (PI/2)); // θ′ in Snyder formulas.
+ double θp = 2 * asin(φ * (2/PI)); // θ′ in Snyder formulas.
/*
* If sin(φ) is 1 or -1 we are on a pole.
* Iteration would produce NaN values.
@@ -182,7 +182,7 @@ public class Mollweide extends NormalizedProjection {
final double y = srcPts[srcOff + 1];
final double θ = asin(y);
final double θp = 2 * θ;
- final double φ = asin((θp + sin(θp)) / PI);
+ final double φ = asin((θp + sin(θp)) * (1/PI));
double λ = x / cos(θ);
if (abs(λ) > 2*SQRT_2*PI) {
λ = Double.NaN;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
index 6717376c01..39121747fd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
@@ -51,6 +51,7 @@ import org.apache.sis.referencing.operation.transform.AbstractMathTransform2D;
import org.apache.sis.referencing.operation.transform.ContextualParameters;
import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory;
import org.apache.sis.referencing.operation.transform.MathTransformProvider;
+import org.apache.sis.referencing.operation.transform.DomainDefinition;
import org.apache.sis.internal.referencing.provider.MapProjection;
import org.apache.sis.internal.referencing.CoordinateOperations;
import org.apache.sis.internal.referencing.Formulas;
@@ -177,6 +178,22 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D imple
*/
static final int MAXIMUM_ITERATIONS = Formulas.MAXIMUM_ITERATIONS;
+ /**
+ * Arbitrary latitude threshold (in radians) for considering that a point is in the polar area.
+ * This is used for implementations of the {@link #getDomain(DomainDefinition)} method.
+ */
+ static final double POLAR_AREA_LIMIT = PI/2 * (84d/90);
+
+ /**
+ * An arbitrarily large longitude value (in radians) for map projections capable to have infinite
+ * extent in the east-west direction. This is the case of {@link Mercator} projection for example.
+ * Longitudes have no real world meaning outside −180° … +180° range (unless wraparound is applied),
+ * but we nevertheless accept large values without wraparound because they make envelope projection easier.
+ * This is used for implementations of the {@link #getDomain(DomainDefinition)} method,
+ * which use arbitrary limits anyway.
+ */
+ static final double LARGE_LONGITUDE_LIMIT = 100*PI;
+
/**
* The internal parameter descriptors. Keys are implementation classes. Values are parameter descriptor groups
* containing at least a parameter for the {@link #eccentricity} value, and optionally other internal parameter
@@ -511,6 +528,48 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D imple
return context.completeTransform(factory, this);
}
+ /**
+ * Returns a transform which may shift scaled longitude θ=n⋅λ inside the [−n⋅π … n⋅π] range.
+ * The transform intentionally does <strong>not</strong> force θ to be inside that range in all cases.
+ * We avoid explicit wraparounds as much as possible (as opposed to implicit wraparounds performed by
+ * trigonometric functions) because they tend to introduce discontinuities. We perform wraparounds only
+ * when necessary for the problem of area crossing the anti-meridian (±180°).
+ *
+ * <div class="note"><b>Example:</b>
+ * a CRS for Alaska may have the central meridian at λ₀=−154° of longitude. If the point to project is
+ * at λ=177° of longitude, calculations will be performed with Δλ=331° while the correct value that we
+ * need to use is Δλ=−29°.</div>
+ *
+ * In order to avoid wraparound operations as much as possible, we test only the bound where anti-meridian
+ * problem may happen; no wraparound will be applied for the opposite bound. Furthermore we add or subtract
+ * 360° only once. Even if the point did many turns around the Earth, the 360° shift will still be applied
+ * at most once. The desire to apply the minimal amount of shifts is the reason why we do not use
+ * {@link Math#IEEEremainder(double, double)}.
+ *
+ * <h4>When to use</h4>
+ * This method is invoked by map projections that multiply the longitude values by some scale factor before
+ * to use them in trigonometric functions. Usually we do not explicitly wraparound the longitude values,
+ * because trigonometric functions do that automatically for us. However if the longitude is multiplied
+ * by some factor before to be used in trigonometric functions, then that implicit wraparound is not the
+ * one we expect. The map projection code needs to perform explicit wraparound in such cases.
+ *
+ * @param factory the factory to use for completing the transform with normalization/denormalization steps.
+ * @return the map projection from (λ,φ) to (<var>x</var>,<var>y</var>) coordinates with wraparound if needed.
+ * @throws FactoryException if an error occurred while creating a transform.
+ *
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-486">SIS-486</a>
+ */
+ final MathTransform completeWithWraparound(final MathTransformFactory factory) throws FactoryException {
+ MathTransform kernel = this;
+ final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
+ final double rotation = normalize.getElement(0, DIMENSION);
+ if (rotation != 0 && Double.isFinite(rotation)) { // Finite check is paranoiac (shall always be true).
+ kernel = new LongitudeWraparound(this,
+ LongitudeWraparound.boundOfScaledLongitude(normalize, rotation < 0), rotation);
+ }
+ return context.completeTransform(factory, kernel);
+ }
+
/**
* If this map projection can not handle the parameters given by the user but an other projection could, delegates
* to the other projection. This method can be invoked by some {@link #createMapProjection(MathTransformFactory)}
@@ -649,37 +708,26 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D imple
}
/**
- * If the given scaled longitude θ=n⋅λ is outside the [−n⋅π … n⋅π] range, maybe shifts θ to that range.
- * This method intentionally does <strong>not</strong> force θ to be inside that range in all cases.
- * We avoid explicit wraparounds as much as possible (as opposed to implicit wraparounds performed by
- * trigonometric functions) because they tend to introduce discontinuities. We perform wraparounds only
- * when necessary for the problem of area crossing the anti-meridian (±180°).
- *
- * <div class="note"><b>Example:</b>
- * a CRS for Alaska may have the central meridian at λ₀=−154° of longitude. If the point to project is
- * at λ=177° of longitude, calculations will be performed with Δλ=331° while the correct value that we
- * need to use is Δλ=−29°.</div>
+ * The longitude value where wraparound is, or would be, applied by this map projection.
+ * This is typically {@link Math#PI} (180° converted to radians) but not necessarily,
+ * because implementations are free to scale the longitude values by an arbitrary factor.
*
- * In order to avoid wraparound operations as much as possible, we test only the bound where anti-meridian
- * problem may happen; no wraparound will be applied for the opposite bound. Furthermore we add or subtract
- * 360° only once. Even if the point did many turns around the Earth, the 360° shift will still be applied
- * at most once. The desire to apply the minimal amount of shifts is the reason why we do not use
- * {@link Math#IEEEremainder(double, double)}.
+ * <p>The wraparound may not be really applied by the {@code transform(…)} methods.
+ * Many map projections implicitly wraparound longitude values through the use of trigonometric functions
+ * ({@code sin(λ)}, {@code cos(λ)}, <i>etc</i>). For those map projections, the wraparound is unconditional.
+ * But some other map projections are capable to handle longitude values beyond the [−180° … +180°] range
+ * as if the world was expanding toward infinity in east and west directions.
+ * The most common example is the {@linkplain Mercator} projection.
+ * In those latter cases, wraparounds are avoided as much as possible in order to facilitate the projection
+ * of envelopes, geometries or rasters, where discontinuities (sudden jumps of 360°) cause artifacts.</p>
*
- * @param θ the scaled longitude value θ=n⋅λ where <var>n</var> is a projection-dependent factor.
- * @param θ_bound minimal (if negative) or maximal (if positive) value of θ before to apply the shift.
- * This is computed by <code>{@linkplain Initializer#boundOfScaledLongitude(double)
- * Initializer.boundOfScaledLongitude}(n)</code>
- * @return θ or shifted θ.
+ * @return the longitude value where wraparound is or would be applied.
*
- * @see Initializer#boundOfScaledLongitude(double)
- * @see <a href="https://issues.apache.org/jira/browse/SIS-486">SIS-486</a>
+ * @see #getDomain(DomainDefinition)
*/
- static double wraparoundScaledLongitude(double θ, final double θ_bound) {
- if (θ_bound < 0 ? θ < θ_bound : θ > θ_bound) {
- θ -= 2*θ_bound;
- }
- return θ;
+ final double getWraparoundLongitude() {
+ return LongitudeWraparound.boundOfScaledLongitude(
+ context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION), false);
}
/*
@@ -775,7 +823,9 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D imple
}
/**
- * Inverse of a normalized map projection.
+ * Reverse of a normalized map projection.
+ * Note that a slightly modified copy of this class is in {@code LongitudeWraparound.Inverse}.
+ * If this is class is modified, consider updating {@code LongitudeWraparound.Inverse} accordingly.
*
* @author Martin Desruisseaux (Geomatys)
* @version 1.0
@@ -789,19 +839,19 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D imple
private static final long serialVersionUID = 6014176098150309651L;
/**
- * The enclosing transform.
+ * The projection to reverse, which is the enclosing transform.
*/
private final NormalizedProjection forward;
/**
- * Default constructor.
+ * Creates a reverse projection for the given forward projection.
*/
Inverse(final NormalizedProjection forward) {
this.forward = forward;
}
/**
- * Returns the inverse of this math transform.
+ * Returns the inverse of this math transform, which is the forward projection.
*/
@Override
public MathTransform2D inverse() {
@@ -809,7 +859,7 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D imple
}
/**
- * Inverse transforms the specified {@code srcPts} and stores the result in {@code dstPts}.
+ * Reverse projects the specified {@code srcPts} and stores the result in {@code dstPts}.
* If the derivative has been requested, then this method will delegate the derivative
* calculation to the enclosing class and inverts the resulting matrix.
*/
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
index 7e71293a4d..9d2f3eb0b2 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
@@ -63,7 +63,7 @@ import static org.apache.sis.internal.referencing.provider.ObliqueStereographic.
* @author Rémi Maréchal (Geomatys)
* @author Martin Desruisseaux (Geomatys)
* @author Emmanuel Giasson (Thales)
- * @version 1.2
+ * @version 1.3
* @since 0.7
* @module
*/
@@ -71,7 +71,7 @@ public class ObliqueStereographic extends NormalizedProjection {
/**
* For cross-version compatibility.
*/
- private static final long serialVersionUID = -1725537881127730658L;
+ private static final long serialVersionUID = -1454098847621943639L;
/**
* Conformal latitude of origin (χ₀), together with its sine and cosine.
@@ -94,15 +94,6 @@ public class ObliqueStereographic extends NormalizedProjection {
*/
private final double g, h;
- /**
- * A bound of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ values.
- * Some (not all) θ values need to be shifted inside that range before to use them
- * in trigonometric functions.
- *
- * @see Initializer#boundOfScaledLongitude(DoubleDouble)
- */
- private final double θ_bound;
-
/**
* Creates an Oblique Stereographic projection from the given parameters.
* The {@code method} argument can be the description of one of the following:
@@ -184,7 +175,6 @@ public class ObliqueStereographic extends NormalizedProjection {
final double R2 = 2 * initializer.radiusOfConformalSphere(sinφ0);
denormalize.convertBefore(0, R2, null);
denormalize.convertBefore(1, R2, null);
- θ_bound = initializer.boundOfScaledLongitude(n);
}
/**
@@ -192,14 +182,13 @@ public class ObliqueStereographic extends NormalizedProjection {
*/
ObliqueStereographic(final ObliqueStereographic other) {
super(null, other);
- χ0 = other.χ0;
- sinχ0 = other.sinχ0;
- cosχ0 = other.cosχ0;
- c = other.c;
- n = other.n;
- g = other.g;
- h = other.h;
- θ_bound = other.θ_bound;
+ χ0 = other.χ0;
+ sinχ0 = other.sinχ0;
+ cosχ0 = other.cosχ0;
+ c = other.c;
+ n = other.n;
+ g = other.g;
+ h = other.h;
}
/**
@@ -249,11 +238,10 @@ public class ObliqueStereographic extends NormalizedProjection {
return delegate(factory, PolarStereographicA.NAME);
}
}
- ObliqueStereographic kernel = this;
if (eccentricity == 0 && getClass() == ObliqueStereographic.class) {
- kernel = new Spherical(this);
+ return context.completeTransform(factory, new Spherical(this));
}
- return context.completeTransform(factory, kernel);
+ return completeWithWraparound(factory);
}
/**
@@ -271,9 +259,8 @@ public class ObliqueStereographic extends NormalizedProjection {
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
- // Λ = λ⋅n (see below), ignoring longitude of origin.
- final double Λ = wraparoundScaledLongitude(srcPts[srcOff], θ_bound);
- final double φ = srcPts[srcOff + 1];
+ final double Λ = srcPts[srcOff ]; // Λ = λ⋅n (see below), ignoring longitude of origin.
+ final double φ = srcPts[srcOff+1];
final double sinφ = sin(φ);
final double ℯsinφ = eccentricity * sinφ;
final double Sa = (1 + sinφ) / (1 - sinφ);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
index 454920c019..039f1ac336 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
@@ -17,10 +17,13 @@
package org.apache.sis.referencing.operation.projection;
import java.util.EnumMap;
+import org.opengis.util.FactoryException;
import org.opengis.parameter.ParameterDescriptor;
+import org.opengis.parameter.InvalidParameterValueException;
import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.MathTransformFactory;
import org.opengis.referencing.operation.OperationMethod;
-import org.opengis.parameter.InvalidParameterValueException;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.referencing.operation.matrix.Matrix2;
@@ -61,7 +64,7 @@ import static org.apache.sis.internal.referencing.provider.SatelliteTracking.*;
*
* @author Matthieu Bastianelli (Geomatys)
* @author Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.3
* @since 1.1
* @module
*/
@@ -69,7 +72,7 @@ public class SatelliteTracking extends NormalizedProjection {
/**
* For cross-version compatibility.
*/
- private static final long serialVersionUID = -209787336760184649L;
+ private static final long serialVersionUID = 859940667477896653L;
/**
* Sines and cosines of inclination between the plane of the Earth's Equator and the plane
@@ -98,15 +101,6 @@ public class SatelliteTracking extends NormalizedProjection {
*/
private final boolean isConic;
- /**
- * A bound of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ values.
- * Some (not all) θ values need to be shifted inside that range before to use them
- * in trigonometric functions.
- *
- * @see Initializer#boundOfScaledLongitude(DoubleDouble)
- */
- private final double θ_bound;
-
/**
* Work around for RFE #4093999 in Sun's bug database ("Relax constraint on
* placement of this()/super() call in constructors").
@@ -203,7 +197,6 @@ public class SatelliteTracking extends NormalizedProjection {
normalize .convertAfter (0, n, null);
denormalize.convertBefore(0, +ρf, null);
denormalize.convertBefore(1, -ρf, ρ0);
- θ_bound = initializer.boundOfScaledLongitude(n);
} else {
/*
* Cylindrical projection case. The equations are (ignoring R and λ₀):
@@ -215,7 +208,7 @@ public class SatelliteTracking extends NormalizedProjection {
* The cosφ₁ (for x at dimension 0) and cosφ₁/F₁′ (for y at dimension 1) factors are computed
* in advance and stored below. The remaining factor to compute in transform(…) method is L.
*/
- n = s0 = θ_bound = Double.NaN;
+ n = s0 = Double.NaN;
final double cotF = sqrt(cos2_φ1 - cos2_i) / (p2_on_p1*cos2_φ1 - cos_i); // Cotangente of F₁.
denormalize.convertBefore(0, cosφ1, null);
denormalize.convertBefore(1, cosφ1*cotF, null);
@@ -225,6 +218,24 @@ public class SatelliteTracking extends NormalizedProjection {
}
}
+ /**
+ * 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>.
+ *
+ * @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 {
+ if (isConic) {
+ return completeWithWraparound(factory);
+ } else {
+ return super.createMapProjection(factory);
+ }
+ }
+
/**
* Returns the names of additional internal parameters which need to be taken in account when
* comparing two {@code SatelliteTracking} projections or formatting them in debug mode.
@@ -330,7 +341,6 @@ public class SatelliteTracking extends NormalizedProjection {
double x = srcPts[srcOff];
double y = λt - p2_on_p1 * λpm;
if (isConic) {
- x = wraparoundScaledLongitude(x, θ_bound);
λpm = n*y + s0; // Use this variable for a new purpose. Needed for derivative.
if ((Double.doubleToRawLongBits(λpm) ^ Double.doubleToRawLongBits(n)) < 0) {
/*
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
index 3f6ebf450e..4696faabb4 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
@@ -386,7 +386,7 @@ public class TransverseMercator extends NormalizedProjection {
public Optional<Envelope> getDomain(final DomainDefinition criteria) {
final Envelope2D domain = new Envelope2D();
domain.x = -PI/2 * (40d/90);
- domain.y = -PI/2 * (84d/90);
+ domain.y = -POLAR_AREA_LIMIT;
domain.width = -2 * domain.x;
domain.height = -2 * domain.y;
return Optional.of(domain);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java
index af043a9996..eff7676b22 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ZonedGridSystem.java
@@ -160,7 +160,7 @@ public class ZonedGridSystem extends AbstractMathTransform2D implements Serializ
*/
@Override
public Optional<Envelope> getDomain(final DomainDefinition criteria) {
- final double y = -PI/2 * (84d/90);
+ final double y = -NormalizedProjection.POLAR_AREA_LIMIT;
return Optional.of(new Envelope2D(null, -PI, y, 2*PI, -2*y));
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
index 0edb4c7586..42d0e846df 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
@@ -335,7 +335,7 @@ public final strictfp class LambertConicConformalTest extends MapProjectionTestC
createCompleteProjection(new LambertConformal1SP(),
6371007, // Semi-major axis length
6371007, // Semi-minor axis length
- 0.5, // Central meridian
+ 0, // Central meridian
40, // Latitude of origin
NaN, // Standard parallel 1
NaN, // Standard parallel 2
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java
index 15d09bdf6e..16bc487a69 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java
@@ -364,7 +364,7 @@ public final strictfp class MercatorTest extends MapProjectionTestCase {
createCompleteProjection(new Mercator1SP(),
6371007, // Semi-major axis length
6371007, // Semi-minor axis length
- 0.5, // Central meridian
+ 0, // Central meridian
NaN, // Latitude of origin (none)
NaN, // Standard parallel 1 (none)
NaN, // Standard parallel 2 (none)
@@ -408,5 +408,11 @@ public final strictfp class MercatorTest extends MapProjectionTestCase {
isInverseTransformSupported = false;
tolerance = Formulas.LINEAR_TOLERANCE;
verifyTransform(coordinates, expected);
+ /*
+ * Replace the 181° longitude by -179° so we can test inverse projection.
+ */
+ coordinates[2] = -179;
+ isInverseTransformSupported = true;
+ verifyTransform(coordinates, expected);
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java
index 4d99553d93..b7d56a93b2 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java
@@ -170,7 +170,7 @@ public final strictfp class PolarStereographicTest extends MapProjectionTestCase
createCompleteProjection(new PolarStereographicA(),
6371007, // Semi-major axis length
6371007, // Semi-minor axis length
- 0.5, // Central meridian
+ 0, // Central meridian
latitudeOfOrigin, // Latitude of origin
NaN, // Standard parallel 1 (none)
NaN, // Standard parallel 2 (none)
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
index d512b4f971..3268687771 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
@@ -146,7 +146,7 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
createCompleteProjection(new org.apache.sis.internal.referencing.provider.TransverseMercator(),
6371007, // Semi-major axis length
6371007, // Semi-minor axis length
- 0.5, // Central meridian
+ 0, // Central meridian
2.5, // Latitude of origin
NaN, // Standard parallel 1 (none)
NaN, // Standard parallel 2 (none)