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 c65774a Document better `GeodeticCalculator` accuracy considerations
and replace the special cases in JUnit tests by a `KnownProblem` enumeration in
an attempt to make clearer what the problem can be and how the tests handle
them.
c65774a is described below
commit c65774ab7961b2bc0c9322083bb5a492de3118d7
Author: Martin Desruisseaux <[email protected]>
AuthorDate: Mon Jan 18 14:36:20 2021 +0100
Document better `GeodeticCalculator` accuracy considerations and replace
the special cases in JUnit tests by a `KnownProblem` enumeration in an attempt
to make clearer what the problem can be and how the tests handle them.
---
.../apache/sis/internal/referencing/Formulas.java | 6 ++
.../sis/referencing/GeodesicsOnEllipsoid.java | 15 ++-
.../apache/sis/referencing/GeodeticCalculator.java | 4 +-
.../sis/referencing/GeodesicsOnEllipsoidTest.java | 47 +++------
.../sis/referencing/GeodeticCalculatorTest.java | 105 ++++++++++++---------
5 files changed, 94 insertions(+), 83 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
index ee258e4..80c250a 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
@@ -43,6 +43,12 @@ public final class Formulas extends Static {
* assuming that the unit of measurement is metre. This constant
determines also
* (indirectly) the minimum accuracy of iterative methods in map
projections.
*
+ * <h4>Maintenance</h4>
+ * If this value is modified, then all usages of this constant should be
verified.
+ * Some usages may need to be compensated. For example {@code
GeodesicsOnEllipsoid}
+ * uses a millimetric precision by dividing the tolerance by 10 or more.
We way want
+ * to keep the same precision there even if {@code LINEAR_TOLERANCE} was
made smaller.
+ *
* @see #ANGULAR_TOLERANCE
* @see org.apache.sis.internal.util.Numerics#COMPARISON_THRESHOLD
*/
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
index 52fd130..2fb94d0 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
@@ -88,16 +88,21 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
/**
* Accuracy threshold for iterative computations, in radians.
- * We take a finer accuracy than default SIS configuration in order to met
the accuracy of numbers
- * published in Karney (2013). If this value is modified, the effect can
be verified by executing
- * the {@code GeodesicsOnEllipsoidTest} methods that compare computed
values against Karney's tables.
- * Remember to update {@link GeodeticCalculator} class javadoc if this
value is changed.
+ * This accuracy must be at least {@value Formulas#ANGULAR_TOLERANCE}
degrees (converted to radians) for
+ * conformance with the accuracy reported in {@link GeodeticCalculator}
class javadoc. Actually we take
+ * a finer accuracy than above value in order to met the accuracy of
numbers published in Karney (2013),
+ * but this extra accuracy is not guaranteed because it is hard to achieve
in all cases.
*
* <p><b>Note:</b> when the iteration loop detects that it reached this
requested accuracy, the loop
* completes the iteration step which was in progress. Consequently the
final accuracy is one iteration
* better than the accuracy computed from this value.</p>
+ *
+ * <h4>Maintenance</h4>
+ * If this value is modified, the effect can be verified by executing the
{@code GeodesicsOnEllipsoidTest}
+ * methods that compare computed values against Karney's tables. The
{@link GeodeticCalculator} javadoc may
+ * need to be edited accordingly.
*/
- static final double ITERATION_TOLERANCE = Formulas.ANGULAR_TOLERANCE *
(PI/180) / 20;
+ static final double ITERATION_TOLERANCE = (Formulas.ANGULAR_TOLERANCE /
20) * (PI/180);
/**
* Difference between ending point and antipode of starting point for
considering them as nearly antipodal.
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
index 6161778..06d29a7 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
@@ -84,7 +84,7 @@ import static
org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEqui
* </li>
* </ul>
*
- * <h2>Algorithms and accuracy</h2>
+ * <h2>Algorithms</h2>
* {@code GeodeticCalculator} uses two set of formulas, depending if the
figure of the Earth
* {@linkplain Ellipsoid#isSphere() is a sphere} or an ellipsoid.
* Publications relevant to this class are:
@@ -104,7 +104,9 @@ import static
org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEqui
* for the reference implementation.</li>
* </ul>
*
+ * <h2>Accuracy</h2>
* {@code GeodeticCalculator} aims for a positional accuracy of one centimetre.
+ * The accuracy is often better (about one millimetre), but not everywhere.
* Azimuthal accuracy corresponds to an error of one centimetre at a distance
of one kilometer,
* except for nearly antipodal points (less than 1° of longitude and latitude
from antipode)
* and points close to the poles where the azimuthal errors are larger.
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
index c39d1ec..75ccf4e 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
@@ -126,7 +126,7 @@ public final strictfp class GeodesicsOnEllipsoidTest
extends GeodeticCalculatorT
* {@code true} if iteration stopped before to reach the desired
accuracy because of limitation
* in {@code double} precision. This field must be reset to {@code
false} before any new point.
*
- * @see #iterationReachedPrecisionLimit()
+ * @see #relaxIfConfirmed(KnownProblem)
* @see #clear()
*/
private boolean iterationReachedPrecisionLimit;
@@ -143,8 +143,7 @@ public final strictfp class GeodesicsOnEllipsoidTest
extends GeodeticCalculatorT
@Override void store(final String name, final double value) {
if (name.equals("dα₁ ≪ α₁")) {
iterationReachedPrecisionLimit = true;
- }
- if (name.length() == 1) {
+ } else if (name.length() == 1) {
switch (name.charAt(0)) {
case 'x': x = value; break;
case 'y': y = value; break;
@@ -523,40 +522,24 @@ public final strictfp class GeodesicsOnEllipsoidTest
extends GeodeticCalculatorT
}
/**
- * Returns {@code true} if iteration stopped before to reach the desired
accuracy because of limitation in
- * {@code double} precision. This problem may happen in the {@link
GeodesicsOnEllipsoid#computeDistance()}
+ * Returns a value {@literal > 1} if iteration stopped before to reach the
desired accuracy because of limitation
+ * in {@code double} precision. This problem may happen in the {@link
GeodesicsOnEllipsoid#computeDistance()}
* method when {@literal dα₁ ≪ α₁}. If locale variable storage is enabled,
this situation is flagged by the
* {@code "dα₁ ≪ α₁"} key. Otherwise we conservatively assume that this
situation occurred.
*/
@Override
- boolean iterationReachedPrecisionLimit() {
- if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) {
- return testedEarth.iterationReachedPrecisionLimit;
- } else {
- // Conservative value in absence of information.
- return true;
- }
- }
-
- /**
- * Tells whether failure to compute geodesic for the given data should
cause the test case to fail.
- * This is invoked by {@link #compareAgainstDataset()}.
- *
- * @param expected a row from the {@code $SIS_DATA/Tests/GeodTest.dat}
file.
- * @return whether the JUnit test should fail.
- */
- @Override
- boolean isFailure(final double[] expected) {
- final double φ1 = expected[COLUMN_φ1];
- final double φ2 = expected[COLUMN_φ2];
- final double Δλ = expected[COLUMN_λ2] - expected[COLUMN_λ1];
- if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) {
- return true; // Ignore equatorial case
in previous version (not anymore).
- }
- if (Δλ > 179 && abs(φ1 + φ2) < 0.002) {
- return false; // Ignore antipodal case.
+ double relaxIfConfirmed(final KnownProblem potentialProblem) {
+ if (potentialProblem ==
KnownProblem.ITERATION_REACHED_PRECISION_LIMIT) {
+ if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) {
+ if (testedEarth.iterationReachedPrecisionLimit) {
+ return 2;
+ }
+ } else {
+ // No information about whether the problem really occurred.
+ return 2;
+ }
}
- return super.isFailure(expected);
+ return super.relaxIfConfirmed(potentialProblem);
}
/**
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
index 1df1bb6..e7b62eb 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
@@ -463,24 +463,25 @@ public strictfp class GeodeticCalculatorTest extends
TestCase {
* We execute only one test for each row instead than
executing both tests,
* for making sure that `GeodeticCalculator` never see the
expected values.
*/
+ KnownProblem potentialProblem = null;
final boolean isTestingInverse = random.nextBoolean();
- final double cosφ1 = abs(cos(toRadians(expected[COLUMN_φ1])));
// For adjusting longitude tolerance.
- final double cosφ2 = abs(cos(toRadians(expected[COLUMN_φ2])));
+ final double φ1 = expected[COLUMN_φ1];
+ final double φ2 = expected[COLUMN_φ2];
+ final double cosφ1 = abs(cos(toRadians(φ1))); // For
adjusting longitude tolerance.
+ final double cosφ2 = abs(cos(toRadians(φ2)));
+ final double Δλ = abs(expected[COLUMN_λ2] -
expected[COLUMN_λ1]);
double linearTolerance, latitudeTolerance, longitudeTolerance,
azimuthTolerance;
- final boolean isLongArcOnEquator;
if (isSphere) {
/*
* When spherical formulas are used instead than
ellipsoidal formulas, an error up to 1% is expected
* in distance calculations (source: Wikipedia). A yet
larger error is observed for azimuth values,
* especially near poles and between antipodal points.
Following are empirical thresholds.
*/
- isLongArcOnEquator = false;
// Not a problem in spherical case.
linearTolerance = expected[COLUMN_Δs] * 0.01;
latitudeTolerance = toDegrees(linearTolerance /
c.semiMajorAxis);
- longitudeTolerance = expected[COLUMN_φ2] > 89.5 ? 180 :
latitudeTolerance / cosφ2;
+ longitudeTolerance = φ2 > 89.5 ? 180 : latitudeTolerance /
cosφ2;
azimuthTolerance = 0.5;
// About 8.8 metres at distance of 1 km.
if (isTestingInverse) {
- final double Δλ = abs(expected[COLUMN_λ2] -
expected[COLUMN_λ1]);
if (Δλ > 179) azimuthTolerance = 100;
else if (Δλ > 178) azimuthTolerance = 20;
else if (Δλ > 175) azimuthTolerance = 10;
@@ -498,18 +499,16 @@ public strictfp class GeodeticCalculatorTest extends
TestCase {
longitudeTolerance = Formulas.ANGULAR_TOLERANCE / cosφ2;
azimuthTolerance = Formulas.LINEAR_TOLERANCE * (180/PI)
/ 10000;
if (isTestingInverse) {
- final double Δ = max(abs(180 - abs(expected[COLUMN_λ2]
- expected[COLUMN_λ1])),
- abs(expected[COLUMN_φ1]
+ expected[COLUMN_φ2]));
- if (Δ < 1) {
+ if (max(abs(180 - Δλ), abs(φ1 + φ2)) < 1) {
azimuthTolerance = 1 * (180/PI) / 10000;
// 1 meter for 10 km.
}
+ if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) {
+ potentialProblem =
KnownProblem.ITERATION_REACHED_PRECISION_LIMIT;
+ }
+ if (Δλ > 179 && abs(φ1 + φ2) < 0.002) {
+ potentialProblem =
KnownProblem.NO_CONVERGENCE_ON_ANTIPODAL_POINTS;
+ }
}
- /*
- * If the start and end points are both within about 3
meters from equator and
- * their distance is more than 10,000 km, we may need to
relax the tolerance.
- * This is an empirical adjustment based on test failures
observation.
- */
- isLongArcOnEquator = Math.min(cosφ1, cosφ2) >=
0.9999999999999 && expected[COLUMN_Δs] > 1E+7;
assertEquals("Consistency with accuracy reported in
Javadoc.", 0.001, linearTolerance, 0.0005);
}
/*
@@ -536,21 +535,24 @@ public strictfp class GeodeticCalculatorTest extends
TestCase {
assertEquals("λ₂", expected[COLUMN_λ2],
end.getOrdinate(1), longitudeTolerance);
assertEquals("α₂", expected[COLUMN_α2],
c.getEndingAzimuth(), azimuthTolerance / cosφ2);
assertEquals("∆s", expected[COLUMN_Δs],
c.getGeodesicDistance(), linearTolerance *
- (isLongArcOnEquator &&
iterationReachedPrecisionLimit() ? 2 : 1));
+
relaxIfConfirmed(potentialProblem));
clear();
} catch (GeodeticException | AssertionError e) {
- if (!isTestingInverse || e instanceof AssertionError ||
isFailure(expected)) {
- out.printf("Test failure at line %d: %s%n"
- + "The values provided in the test file are:%n"
- + "(φ₁,λ₁) = %16.12f %16.12f%n"
- + "(φ₂,λ₂) = %16.12f %16.12f%n"
- + "The values computed by the geodesic
calculator are:%n",
- reader.getLineNumber(),
e.getLocalizedMessage(),
- expected[0], expected[1], expected[3],
expected[4]);
- out.println(c);
- throw e;
+ if (e instanceof GeodeticException) {
+ if (potentialProblem ==
KnownProblem.NO_CONVERGENCE_ON_ANTIPODAL_POINTS) {
+ noConvergenceCount++;
+ continue;
+ }
}
- noConvergenceCount++;
+ out.printf("Test failure at line %d: %s%n"
+ + "The values provided in the test file are:%n"
+ + "(φ₁,λ₁) = %16.12f %16.12f%n"
+ + "(φ₂,λ₂) = %16.12f %16.12f%n"
+ + "The values computed by the geodesic calculator
are:%n",
+ reader.getLineNumber(), e.getLocalizedMessage(),
+ expected[0], expected[1], expected[3],
expected[4]);
+ out.println(c);
+ throw e;
}
}
}
@@ -558,31 +560,44 @@ public strictfp class GeodeticCalculatorTest extends
TestCase {
}
/**
- * Returns {@code true} if iteration stopped before to reach the desired
accuracy because the correction
- * to apply in iterative steps is smaller than what can be applied using
IEEE 754 double arithmetic.
- * It never happen in the spherical case because there is no iteration in
that case.
+ * Returns the factor by which to relax the linear tolerance, or 1 for no
relaxation.
+ * This method is invoked after {@link
GeodeticCalculator#computeDistance()} has been invoked.
+ * It should check if the potential problem really occurred, and if yes
return a value greater
+ * than 1. If no problem (other than IEEE 754 rounding errors) is expected
to occur, then this
+ * method should return exactly 1.
*
- * <p>This method is invoked only in situations empirically identified as
susceptible to have this problem:</p>
- * <uL>
- * <li>When start point and end point are both within about 3 meters
from equator
- * and their distance is more than 10,000 km.</li>
- * </ul>
+ * @param potentialProblem the problem that may happen, or {@code null}
if none.
+ * @return factor by which to relax the linear tolerance threshold, or 1
if no relaxation.
*/
- boolean iterationReachedPrecisionLimit() {
- return false;
+ double relaxIfConfirmed(KnownProblem potentialProblem) {
+ return 1;
}
/**
- * Tells whether failure to compute geodesic for the given data should
cause the test case to fail.
- * The default implementation always return {@code true}. Subclass can
override if some points are
- * known to fail.
+ * Known problems with our implementation of inverse geodetic with
ellipsoidal formulas.
+ * Those problems may occur in {@link #compareAgainstDataset()} test when
ellipsoidal formulas
+ * are used (those problems do not occur with spherical formulas). This is
an enumeration
+ * of cases where {@link GeodesicsOnEllipsoid#computeDistance()} does not
converge.
*
- * @param expected a row from the {@code $SIS_DATA/Tests/GeodTest.dat}
file.
- * Use {@code COLUMN_*} constant for accessing values by column
indices.
- * @return whether the JUnit test should fail.
+ * @see #compareAgainstDataset()
+ * @see #relaxIfConfirmed(KnownProblem)
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-467";>SIS-467</a>
*/
- boolean isFailure(final double[] expected) {
- return true;
+ enum KnownProblem {
+ /**
+ * Iteration stopped before to reach the desired accuracy because the
correction to apply in iterative
+ * steps is smaller than what can be applied using IEEE 754 double
arithmetic. This problem occurs in
+ * {@literal α₁ -= dα₁} expression when {@literal dα₁ ≪ α₁}, in which
case the subtraction has no effect.
+ * It has been observed on the equator between 2 close points. This is
not considered a failure because
+ * the precision is still in 1 cm precision target. We only need to
relax the 1 mm precision check.
+ */
+ ITERATION_REACHED_PRECISION_LIMIT,
+
+ /**
+ * Iteration failed with a "no convergence error". It sometime happens
during distance calculation between
+ * antipodal points.
+ */
+ NO_CONVERGENCE_ON_ANTIPODAL_POINTS;
}
/**
