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
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new 42eed3c Change strategy regarding
https://issues.apache.org/jira/browse/SIS-486 fix by shifting longitude only if
needed for the anti-meridian problem and not in other cases. The intent is to
reduce discontinuities in the mathematical function.
42eed3c is described below
commit 42eed3c1d1c94451ff78aaa81655849bd93b1df9
Author: Martin Desruisseaux <[email protected]>
AuthorDate: Sun Dec 22 22:15:33 2019 +0100
Change strategy regarding https://issues.apache.org/jira/browse/SIS-486 fix
by shifting longitude only if needed for the anti-meridian problem and not in
other cases. The intent is to reduce discontinuities in the mathematical
function.
---
.../operation/projection/AlbersEqualArea.java | 21 +++++------
.../operation/projection/Initializer.java | 41 +++++++++++++++-------
.../projection/LambertConicConformal.java | 19 +++++-----
.../operation/projection/NormalizedProjection.java | 36 ++++++++++++++++++-
.../operation/projection/ObliqueStereographic.java | 29 +++++++--------
.../operation/projection/SatelliteTracking.java | 16 ++++-----
.../operation/projection/AlbersEqualAreaTest.java | 1 -
.../org/apache/sis/internal/util/DoubleDouble.java | 8 ++---
.../apache/sis/internal/util/DoubleDoubleTest.java | 2 +-
9 files changed, 113 insertions(+), 60 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 88476ee..92f1538 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
@@ -82,12 +82,13 @@ public class AlbersEqualArea extends EqualAreaProjection {
final double C;
/**
- * Size of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ
values.
- * We need to ensure that θ values are inside that range before to use it
in trigonometric functions.
+ * 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#spanOfScaledLongitude(DoubleDouble)
+ * @see Initializer#boundOfScaledLongitude(DoubleDouble)
*/
- final double θ_span;
+ final double θ_bound;
/**
* Creates an Albers Equal Area projection from the given parameters.
@@ -172,7 +173,7 @@ public class AlbersEqualArea extends EqualAreaProjection {
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`.
- θ_span = Initializer.spanOfScaledLongitude(rn);
+ θ_bound = initializer.boundOfScaledLongitude(rn);
}
/**
@@ -180,9 +181,9 @@ public class AlbersEqualArea extends EqualAreaProjection {
*/
AlbersEqualArea(final AlbersEqualArea other) {
super(other);
- nm = other.nm;
- C = other.C;
- θ_span = other.θ_span;
+ nm = other.nm;
+ C = other.C;
+ θ_bound = other.θ_bound;
}
/**
@@ -238,7 +239,7 @@ public class AlbersEqualArea extends EqualAreaProjection {
final boolean derivate) throws ProjectionException
{
// θ = n⋅λ reduced to [−n⋅π … n⋅π] range.
- final double θ = IEEEremainder(srcPts[srcOff], θ_span);
+ final double θ = wraparoundScaledLongitude(srcPts[srcOff], θ_bound);
final double φ = srcPts[srcOff + 1];
final double cosθ = cos(θ);
final double sinθ = sin(θ);
@@ -327,7 +328,7 @@ public class AlbersEqualArea extends EqualAreaProjection {
final boolean derivate)
{
// θ = n⋅λ reduced to [−n⋅π … n⋅π] range.
- final double θ = IEEEremainder(srcPts[srcOff], θ_span);
+ final double θ = wraparoundScaledLongitude(srcPts[srcOff],
θ_bound);
final double φ = srcPts[srcOff + 1];
final double cosθ = cos(θ);
final double sinθ = sin(θ);
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 9e88f90..ff0ac5c 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
@@ -88,6 +88,11 @@ final class Initializer {
final DoubleDouble eccentricitySquared;
/**
+ * Sign of central meridian: -1 if negative, 0 if zero, +1 if positive.
+ */
+ private final byte signum_λ0;
+
+ /**
* Map projection variant.
* Values from 0 to 127 inclusive are convenience values at the discretion
of {@link NormalizedProjection} subclasses.
* Values from 128 to 255 inclusive are values handled in a special way by
{@link Initializer} constructor.
@@ -143,6 +148,8 @@ final class Initializer {
final double fn = getAndStore(roles.get(ParameterRole.FALSE_NORTHING))
- getAndStore(roles.get(ParameterRole.FALSE_SOUTHING));
+ signum_λ0 = (λ0 > 0) ? (byte) +1 :
+ (λ0 < 0) ? (byte) -1 : 0;
eccentricitySquared = new DoubleDouble();
DoubleDouble k = DoubleDouble.createAndGuessError(a); // The value by
which to multiply all results of normalized projection.
if (a != b) {
@@ -381,27 +388,37 @@ final class Initializer {
}
/**
- * Returns the size of the [−n⋅π … n⋅π] range, which is the valid range of
θ = n⋅λ values.
+ * 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, them that
implicit wraparound is not the
+ * 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.
- * Example:
- *
- * {@preformat java
- * double spanθ = spanOfScaledLongitude(n); // Should be computed
only once.
- * double θ = Math.IEEEremainder(λn, spanθ); // λ without n is
typically unknown.
- * }
*
* @param n the factor by which longitude values are multiplied before
use in trigonometry.
- * @return size of the [−n⋅π … n⋅π] range, for use in {@link
Math#IEEEremainder(double, double)}.
+ * @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>
*/
- static double spanOfScaledLongitude(final DoubleDouble n) {
- final DoubleDouble r = DoubleDouble.createTwicePi();
+ 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);
- return abs(r.doubleValue());
+ 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 961eea2..91209c3 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
@@ -141,12 +141,13 @@ public class LambertConicConformal extends
ConformalProjection {
final double n;
/**
- * Size of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ
values.
- * We need to ensure that θ values are inside that range before to use it
in trigonometric functions.
+ * 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#spanOfScaledLongitude(DoubleDouble)
+ * @see Initializer#boundOfScaledLongitude(DoubleDouble)
*/
- final double θ_span;
+ final double θ_bound;
/**
* Creates a Lambert projection from the given parameters.
@@ -347,7 +348,7 @@ public class LambertConicConformal extends
ConformalProjection {
normalize .convertAfter(1, sφ, null);
denormalize.convertBefore(0, F, null); F.negate();
denormalize.convertBefore(1, F, rF);
- θ_span = Initializer.spanOfScaledLongitude(sλ);
+ θ_bound = initializer.boundOfScaledLongitude(sλ);
}
/**
@@ -355,8 +356,8 @@ public class LambertConicConformal extends
ConformalProjection {
*/
LambertConicConformal(final LambertConicConformal other) {
super(other);
- n = other.n;
- θ_span = other.θ_span;
+ n = other.n;
+ θ_bound = other.θ_bound;
}
/**
@@ -416,7 +417,7 @@ 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 θ = IEEEremainder(srcPts[srcOff], θ_span); // θ =
λ⋅n (ignoring longitude of origin)
+ final double θ = wraparoundScaledLongitude(srcPts[srcOff],
θ_bound); // θ = Δλ⋅n
final double φ = srcPts[srcOff + 1]; //
Sign may be reversed
final double absφ = abs(φ);
final double sinθ = sin(θ);
@@ -527,7 +528,7 @@ public class LambertConicConformal extends
ConformalProjection {
final double[] dstPts, final int dstOff,
final boolean derivate)
{
- final double θ = IEEEremainder(srcPts[srcOff], θ_span); //
θ = λ⋅n (ignoring longitude of origin)
+ final double θ = wraparoundScaledLongitude(srcPts[srcOff],
θ_bound); // θ = Δλ⋅n
final double φ = srcPts[srcOff + 1]; //
Sign may be reversed
final double absφ = abs(φ);
final double sinθ = sin(θ);
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 b6e2a0d..013e20c 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
@@ -126,7 +126,7 @@ import static java.lang.Math.*;
* @author André Gosselin (MPO)
* @author Rueben Schulz (UBC)
* @author Rémi Maréchal (Geomatys)
- * @version 0.8
+ * @version 1.1
*
* @see ContextualParameters
* @see <a href="http://mathworld.wolfram.com/MapProjection.html";>Map
projections on MathWorld</a>
@@ -647,6 +647,40 @@ 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 spanning 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)}.
+ *
+ * @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 θ.
+ *
+ * @see Initializer#boundOfScaledLongitude(double)
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-486";>SIS-486</a>
+ */
+ static double wraparoundScaledLongitude(double θ, final double θ_bound) {
+ if (θ_bound < 0 ? θ < θ_bound : θ > θ_bound) {
+ θ -= 2*θ_bound;
+ }
+ return θ;
+ }
+
+ /**
* Converts a single coordinate in {@code srcPts} at the given offset and
stores the result
* in {@code dstPts} at the given offset. In addition, opportunistically
computes the
* transform derivative if requested.
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 be81ee1..64adcd0 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
@@ -93,12 +93,13 @@ public class ObliqueStereographic extends
NormalizedProjection {
private final double g, h;
/**
- * Size of the [−n⋅π … n⋅π] range, which is the valid range of Λ = n⋅λ
values.
- * We need to ensure that Λ values are inside that range before to use it
in trigonometric functions.
+ * 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#spanOfScaledLongitude(DoubleDouble)
+ * @see Initializer#boundOfScaledLongitude(DoubleDouble)
*/
- private final double θ_span;
+ private final double θ_bound;
/**
* Creates an Oblique Stereographic projection from the given parameters.
@@ -181,7 +182,7 @@ public class ObliqueStereographic extends
NormalizedProjection {
final double R2 = 2 * initializer.radiusOfConformalSphere(sinφ0);
denormalize.convertBefore(0, R2, null);
denormalize.convertBefore(1, R2, null);
- θ_span = n * (2*PI);
+ θ_bound = initializer.boundOfScaledLongitude(n);
}
/**
@@ -189,14 +190,14 @@ public class ObliqueStereographic extends
NormalizedProjection {
*/
ObliqueStereographic(final ObliqueStereographic other) {
super(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;
- θ_span = other.θ_span;
+ χ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;
}
/**
@@ -267,7 +268,7 @@ public class ObliqueStereographic extends
NormalizedProjection {
final boolean derivate) throws ProjectionException
{
// Λ = λ⋅n (see below), ignoring longitude of origin.
- final double Λ = IEEEremainder(srcPts[srcOff], θ_span);
+ final double Λ = wraparoundScaledLongitude(srcPts[srcOff],
θ_bound);
final double φ = srcPts[srcOff + 1];
final double sinφ = sin(φ);
final double ℯsinφ = eccentricity * 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 1b9d68b..430e92f 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
@@ -99,12 +99,13 @@ public class SatelliteTracking extends NormalizedProjection
{
private final boolean isConic;
/**
- * Size of the [−n⋅π … n⋅π] range, which is the valid range of θ = n⋅λ
values.
- * We need to ensure that θ values are inside that range before to use it
in trigonometric functions.
+ * 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#spanOfScaledLongitude(DoubleDouble)
+ * @see Initializer#boundOfScaledLongitude(DoubleDouble)
*/
- private final double θ_span;
+ private final double θ_bound;
/**
* Work around for RFE #4093999 in Sun's bug database ("Relax constraint on
@@ -202,7 +203,7 @@ public class SatelliteTracking extends NormalizedProjection
{
normalize .convertAfter (0, n, null);
denormalize.convertBefore(0, +ρf, null);
denormalize.convertBefore(1, -ρf, ρ0);
- θ_span = n * (2*PI);
+ θ_bound = initializer.boundOfScaledLongitude(n);
} else {
/*
* Cylindrical projection case. The equations are (ignoring R and
λ₀):
@@ -214,8 +215,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.
*/
- θ_span = 2*PI;
- n = s0 = Double.NaN;
+ n = s0 = θ_bound = 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);
@@ -308,7 +308,7 @@ public class SatelliteTracking extends NormalizedProjection
{
double x = srcPts[srcOff];
double y = λt - p2_on_p1 * λpm;
if (isConic) {
- x = IEEEremainder(x, θ_span);
+ 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/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
index a60e290..5ca8775 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
@@ -24,7 +24,6 @@ import static java.lang.StrictMath.*;
import static java.lang.Double.NaN;
// Test dependencies
-
import org.opengis.test.ToleranceModifier;
import org.apache.sis.test.DependsOnMethod;
import org.apache.sis.test.DependsOn;
diff --git
a/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
b/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
index b5a9b02..b30d0cd 100644
---
a/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
+++
b/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
@@ -309,12 +309,12 @@ public final class DoubleDouble extends Number {
}
/**
- * Returns a new {@code DoubleDouble} instance initialized to the 2π value.
+ * Returns a new {@code DoubleDouble} instance initialized to the π value.
*
- * @return an instance initialized to the
6.28318530717958647692528676655901 value.
+ * @return an instance initialized to the
3.14159265358979323846264338327950 value.
*/
- public static DoubleDouble createTwicePi() {
- return new DoubleDouble(6.28318530717958647692528676655901,
2.4492935982947064E-16);
+ public static DoubleDouble createPi() {
+ return new DoubleDouble(3.14159265358979323846264338327950,
1.2246467991473532E-16);
}
/**
diff --git
a/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
b/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
index ebdf365..c08fc24 100644
---
a/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
+++
b/core/sis-utility/src/test/java/org/apache/sis/internal/util/DoubleDoubleTest.java
@@ -444,7 +444,7 @@ public final strictfp class DoubleDoubleTest extends
TestCase {
for (int i=0; ; i++) {
final DoubleDouble dd;
switch (i) {
- case 0: dd = DoubleDouble.createTwicePi(); break;
+ case 0: dd = DoubleDouble.createPi(); break;
case 1: dd = DoubleDouble.createRadiansToDegrees(); break;
case 2: dd = DoubleDouble.createDegreesToRadians(); break;
case 3: dd = DoubleDouble.createSecondsToRadians(); break;
