This is an automated email from the ASF dualhosted git repository.
desruisseaux pushed a commit to branch geoapi4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
The following commit(s) were added to refs/heads/geoapi4.0 by this push:
new b0c3203 After having verified that `TransverseMercator` now behave like a monotonic function up to 90° of longitude, remove the limit previously set at 70°. Add notes saying that results for longitude larger than 60° are inexact. We nevertheless allow computations between 60° and 90° because they are required for handling some raster data.
b0c3203 is described below
commit b0c32031038121918c8d3d1df090b5df978033e1
Author: Martin Desruisseaux
AuthorDate: Thu Feb 25 00:48:23 2021 +0100
After having verified that `TransverseMercator` now behave like a monotonic function up to 90° of longitude,
remove the limit previously set at 70°. Add notes saying that results for longitude larger than 60° are inexact.
We nevertheless allow computations between 60° and 90° because they are required for handling some raster data.

.../operation/projection/TransverseMercator.java  99 +++++++++
.../projection/TransverseMercatorTest.java  55 ++++++++++
2 files changed, 85 insertions(+), 69 deletions()
diff git a/core/sisreferencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java b/core/sisreferencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
index fb00694..94d03f6 100644
 a/core/sisreferencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
+++ b/core/sisreferencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
@@ 57,6 +57,13 @@ import static org.apache.sis.internal.referencing.provider.TransverseMercator.*;
* and the latitude of origin is the equator. A scale factor of 0.9996 and false easting of 500000 metres is used for
* all zones and a false northing of 10000000 metres is used for zones in the southern hemisphere.
*
+ * Domain of validity
+ * The difference between longitude values λ and the central meridian λ₀ should be less than 60°.
+ * Differences larger than 90° of longitude cause a {@link ProjectionException} to be thrown.
+ * Differences between 60° and 90° are not rejected by Apache SIS but should be avoided.
+ * See the {@linkplain #transform(double[], int, double[], int, boolean) projection method}
+ * for more information.
+ *
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @version 1.1
@@ 74,41 +81,6 @@ public class TransverseMercator extends NormalizedProjection {
private static final long serialVersionUID = 627685138188387835L;
/**
 * Distance from central meridian, in degrees, at which errors are considered too important.
 * This threshold is determined by comparisons of computed values against values provided by
 * Karney (2009) Test data for the transverse Mercator projection data file.
 * On the WGS84 ellipsoid we observed at equator:
 *
 *
 *  For ∆λ below 60°, errors below 1 centimetre.
 *  For ∆λ between 60° and 66°, errors up to 0.1 metre.
 *  For ∆λ between 66° and 70°, errors up to 1 metre.
 *  For ∆λ between 70° and 72°, errors up to 2 metres.
 *  For ∆λ between 72° and 74°, errors up to 10 metres.
 *  For ∆λ between 74° and 76°, errors up to 30 metres.
 *  For ∆λ between 76° and 78°, errors up to 1 kilometre.
 *  For ∆λ greater than 85°, results are chaotic.
 *
 *
 * For latitudes greater than 20°, errors are less than 0.5 meters for all ∆λ < 83°.
 * Errors become suddenly much higher at ∆λ > 82.6° no matter the latitude.
 */
 static final double DOMAIN_OF_VALIDITY = 82.5;

 /**
 * The domain of validity at equator. We keep using the limit at all latitudes up to
 * {@value #LATITUDE_OF_REDUCED_DOMAIN} degrees, even if actually the limit could be
 * progressively relaxed until it reaches {@value #DOMAIN_OF_VALIDITY} degrees.
 */
 static final double DOMAIN_OF_VALIDITY_AT_EQUATOR = 70;

 /**
 * The maximal latitude (exclusive) where to use {@link #DOMAIN_OF_VALIDITY_AT_EQUATOR}
 * instead of {@link #DOMAIN_OF_VALIDITY}.
 */
 static final double LATITUDE_OF_REDUCED_DOMAIN = 20;

 /**
* {@code false} for using the original formulas as published by EPSG, or {@code true} for using formulas
* modified using trigonometric identities. The use of trigonometric identities is for reducing the amount
* of calls to the {@link Math#sin(double)} and similar methods. Some identities used are:
@@ 382,6 +354,31 @@ public class TransverseMercator extends NormalizedProjection {
* Converts the specified (λ,φ) coordinate (units in radians) and stores the result in {@code dstPts}.
* In addition, opportunistically computes the projection derivative if {@code derivate} is {@code true}.
*
+ * Accuracy and domain of validity
+ * Projection errors depend on the difference ∆λ between longitude λ and the central meridian λ₀.
+ * All Universal Transverse Mercator (UTM) projections aim for ∆λ ≤ 3°, but this implementation
+ * can nevertheless handle larger values. Results have been compared with values provided by
+ * Karney, C.F.F. (2009).
+ * Test data for the transverse Mercator projection [Data set]. Zenodo.
+ * On the WGS84 ellipsoid we observed the following errors compared to Karney's data:
+ *
+ *
+ *  Errors less than 1 centimetre for ∆λ < 60° at all latitudes.
+ *  At latitudes far enough from equator (φ ≥ 20°), the domain can be extended up to ∆λ <
+ * (1 − ℯ)⋅90° (≈ 82.63627282416406551° on WGS84) with errors less than 70 centimetres.
+ *
+ *
+ * Case of 82.6…° < ∆λ ≤ 90°
+ * Karney (2009) uses an “extended” domain of transverse Mercator projection for ∆λ ≥ (1 − ℯ)⋅90°,
+ * but Apache SIS does not support such extension. Consequently ∆λ values between (1 − ℯ)⋅90° and 90°
+ * should be considered invalid but are not rejected by Apache SIS. Note that those invalid values are
+ * consistent with the {@linkplain #inverseTransform(double[], int, double[], int) inverse projection}
+ * (i.e. applying a projection followed by an inverse projection gives approximately the original values).
+ *
+ * Case of ∆λ > 90°
+ * Longitude values at a distance greater than 90° from the central meridian are rejected.
+ * A {@link ProjectionException} is thrown in that case.
+ *
* @return the matrix of the projection derivative at the given source position,
* or {@code null} if the {@code derivate} argument is {@code false}.
* @throws ProjectionException if the coordinate can not be converted.
@@ 393,7 +390,7 @@ public class TransverseMercator extends NormalizedProjection {
{
final double λ = srcPts[srcOff];
final double φ = srcPts[srcOff+1];
 if (abs(λ) >= DOMAIN_OF_VALIDITY_AT_EQUATOR * (PI/180)) {
+ if (abs(λ) > 90 * (PI/180)) {
/*
* The Transverse Mercator projection is conceptually a Mercator projection rotated by 90°.
* In Mercator projection, y values tend toward infinity for latitudes close to ±90°.
@@ 407,15 +404,11 @@ public class TransverseMercator extends NormalizedProjection {
* Since a distance of 90° from central meridian is far outside the Transverse Mercator
* domain of validity anyway, we do not let the user go further.
*
 * In the particular case of ellipsoidal formulas, we put a limit of 70° instead of 90°
 * because experience shows that results close to equator become chaotic after 85° when
 * using WGS84 ellipsoid. We do not need to reduce the limit for the spherical formulas,
 * because the mathematic are simpler and the function still smooth until 90°.
+ * Historical note: in a previous version, we used a limit of 70° instead of 90° because
+ * esults became chaotic after 85°. That limit has been removed in later version because
+ * this method now behaves like a monotonic function.
*/
 final double limit = (abs(φ) < LATITUDE_OF_REDUCED_DOMAIN * (PI/180))
 ? (PI/180) * DOMAIN_OF_VALIDITY_AT_EQUATOR
 : (PI/180) * DOMAIN_OF_VALIDITY;
 if (Math.abs(IEEEremainder(λ, 2*PI)) >= limit) {
+ if (Math.abs(IEEEremainder(λ, 2*PI)) > 90 * (PI/180)) { // More costly check.
throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
}
}
@@ 705,22 +698,8 @@ public class TransverseMercator extends NormalizedProjection {
final double cosλ = cos(λ);
/*
* The Transverse Mercator projection is conceptually a Mercator projection rotated by 90°.
 * In Mercator projection, y values tend toward infinity for latitudes close to ±90° while
 * in Transverse Mercator, x values tend toward infinity for longitudes close to ±90° from
 * central meridian (and at equator). After we pass the 90° limit, the Transverse Mercator
 * results at (90° + Δ) are the same as for (90°  Δ).
 *
 * Problem is that 90° is an ordinary longitude value, not even close to the limit of longitude
 * values range (±180°). So having f(π/2+Δ, φ) = f(π/2Δ, φ) results in wrong behavior in some
 * algorithms like the one used by Envelopes.transform(CoordinateOperation, Envelope). Since a
 * distance of 90° from central meridian is outside projection domain of validity anyway, we do
 * not let the user go further.
 *
 * Note that in the case of ellipsoidal formulas (implemented by parent class), we put the
 * limit to a lower value (DOMAIN_OF_VALIDITY) because experience shows that results close
 * to equator become chaotic after 85° when using WGS84 ellipsoid. We do not need to reduce
 * the limit for the spherical formulas, because the mathematic are simpler and the function
 * still smooth until 90°.
+ * See comment in the `super.transform(…)` implementation for more information about why we
+ * need to reject ∆λ > 90°.
*/
if (cosλ < 0) { // Implies Math.abs(IEEEremainder(λ, 2*PI)) > PI/2
throw new ProjectionException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
diff git a/core/sisreferencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java b/core/sisreferencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
index e526e3f..d512b4f 100644
 a/core/sisreferencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
+++ b/core/sisreferencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
@@ 43,13 +43,52 @@ import org.opengis.test.CalculationType;
* Tests the {@link TransverseMercator} class.
*
* @author Martin Desruisseaux (Geomatys)
 * @version 1.0
+ * @version 1.1
* @since 0.6
* @module
*/
@DependsOn(NormalizedProjectionTest.class)
public final strictfp class TransverseMercatorTest extends MapProjectionTestCase {
/**
+ * Distance from central meridian, in degrees, at which errors are considered too important.
+ * This threshold is determined by comparisons of computed values against values provided by
+ * Karney (2009) Test data for the transverse Mercator projection data file.
+ * On the WGS84 ellipsoid we observed close to equator:
+ *
+ *
+ *  For ∆λ below 60°, errors below 1 centimetre.
+ *  For ∆λ between 60° and 66°, errors up to 0.1 metre.
+ *  For ∆λ between 66° and 70°, errors up to 1 metre.
+ *  For ∆λ between 70° and 72°, errors up to 2 metres.
+ *  For ∆λ between 72° and 74°, errors up to 10 metres.
+ *  For ∆λ between 74° and 76°, errors up to 30 metres.
+ *  For ∆λ between 76° and 78°, errors up to 1 kilometre.
+ *  For ∆λ greater than 85°, errors grow exponentially.
+ *
+ *
+ * On the WGS84 ellipsoid at latitudes greater than 20°, we found errors less than 1 meter
+ * for all ∆λ < (1 − ℯ)⋅90° (82.63627282416406551 in WGS84 case). For larger ∆λ values
+ * Karney (2009) uses an “extended” domain of transverse Mercator projection, but Apache SIS
+ * does not support such extension. Consequently ∆λ values between (1 − ℯ)⋅90° and 90° should
+ * be considered invalid but are not rejected by Apache SIS. Note that even for those invalid
+ * values, the inverse projection continue to gives back the original values.
+ */
+ private static final double DOMAIN_OF_VALIDITY = 82.63627282416406551; // (1 − ℯ)⋅90°
+
+ /**
+ * The domain of validity at equator. We keep using the limit at all latitudes up to
+ * {@value #LATITUDE_OF_REDUCED_DOMAIN} degrees, even if actually the limit could be
+ * progressively relaxed until it reaches {@value #DOMAIN_OF_VALIDITY} degrees.
+ */
+ private static final double DOMAIN_OF_VALIDITY_AT_EQUATOR = 70;
+
+ /**
+ * The maximal latitude (exclusive) where to use {@link #DOMAIN_OF_VALIDITY_AT_EQUATOR}
+ * instead of {@link #DOMAIN_OF_VALIDITY}.
+ */
+ private static final double LATITUDE_OF_REDUCED_DOMAIN = 20;
+
+ /**
* Creates a new instance of {@link TransverseMercator}.
*
* @param ellipsoidal {@code false} for a sphere, or {@code true} for WGS84 ellipsoid.
@@ 192,6 +231,8 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
* @throws IOException if an error occurred while reading the test file.
* @throws FactoryException if an error occurred while creating the map projection.
* @throws TransformException if an error occurred while transforming coordinates.
+ *
+ * @see OptionalTestData#TRANSVERSE_MERCATOR
*/
@Test
public void compareAgainstDataset() throws IOException, FactoryException, TransformException {
@@ 219,18 +260,14 @@ public final strictfp class TransverseMercatorTest extends MapProjectionTestCase
// Relax tolerance for longitudes very far from central meridian.
final double longitude = abs(source[0]);
final double latitude = abs(source[1]);
 final double limit = (latitude < TransverseMercator.LATITUDE_OF_REDUCED_DOMAIN
 ? TransverseMercator.DOMAIN_OF_VALIDITY_AT_EQUATOR
 : TransverseMercator.DOMAIN_OF_VALIDITY);
+ final double limit = (latitude < LATITUDE_OF_REDUCED_DOMAIN
+ ? DOMAIN_OF_VALIDITY_AT_EQUATOR
+ : DOMAIN_OF_VALIDITY);
if (longitude < limit) {
if (latitude >= 89.9) tolerance = 0.1;
else if (longitude <= 60) tolerance = Formulas.LINEAR_TOLERANCE;
else if (longitude <= 66) tolerance = 0.1;
 else if (longitude <= 70) tolerance = 1;
 else if (longitude <= 72) tolerance = 2;
 else if (longitude <= 74) tolerance = 10;
 else if (longitude <= 76) tolerance = 30;
 else tolerance = 1000;
+ else tolerance = 0.7;
transform.transform(source, 0, source, 0, 1);
assertCoordinateEquals(line, target, source, reader.getLineNumber(), CalculationType.DIRECT_TRANSFORM);
}