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 2019/04/23 17:45:59 UTC
[sis] branch geoapi-4.0 updated: Implement MeridionalDistance.meridianArc(…).
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 8a492bc Implement MeridionalDistance.meridianArc(…).
8a492bc is described below
commit 8a492bc47a272e7ccefbe690aa44dd9377e9149d
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Tue Apr 23 19:45:00 2019 +0200
Implement MeridionalDistance.meridianArc(…).
---
.../projection/MeridionalDistanceBased.java | 102 ++++++++++++++++++++-
.../projection/MeridionalDistanceTest.java | 9 +-
2 files changed, 105 insertions(+), 6 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
index 0d2b2ed..192a097 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
@@ -62,6 +62,16 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
/**
+ * Coefficients for the formula implemented by the {@link #meridianArc(double, double, double)} method.
+ * Values are computed by {@link #computeCoefficients()} at construction time or after deserialization.
+ * The values depends on the form of equation implemented by {@code meridianArc(…)}. We do not use the
+ * form published commonly found in publication since a few algebraic operations allow to replace the
+ * sin(2φ), sin(4φ), sin(6φ) and sin(8φ) terms by only sin²(φ), which is faster to calculate.
+ * See {@link #computeCoefficients()} comments in method body for more information.
+ */
+ private transient double c0, c1, c2, c3, c4;
+
+ /**
* Creates a new normalized projection from the parameters computed by the given initializer.
*/
MeridionalDistanceBased(final Initializer initializer) {
@@ -74,6 +84,11 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
*/
MeridionalDistanceBased(final MeridionalDistanceBased other) {
super(other);
+ c0 = other.c0;
+ c1 = other.c1;
+ c2 = other.c2;
+ c3 = other.c3;
+ c4 = other.c4;
}
/**
@@ -89,7 +104,90 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
* This method shall be invoked after {@code MeridionalDistanceBased} construction or deserialization.
*/
private void computeCoefficients() {
- final double e2 = eccentricitySquared;
- // TODO
+ final double e2 = eccentricitySquared;
+ final double e4 = e2*e2;
+ final double e6 = e2*e4;
+ final double e8 = e4*e4;
+ final double e10 = e2*e8;
+ /*
+ * Snyder 3-21 and EPSG guidance notes #7 part 2 (February 2019) give us the following equation:
+ *
+ * M = a[(1 – e²/4 – 3e⁴/64 – 5e⁶/256 – …)⋅φ
+ * – (3e²/8 + 3e⁴/32 + 45e⁶/1024 + …)⋅sin2φ
+ * + (15e⁴/256 + 45e⁶/1024 + …)⋅sin4φ
+ * – (35e⁶/3072 + …)⋅sin6φ
+ * + …]
+ *
+ * But we do not use this equation as-is. Instead, we start from a series providing two more terms (e⁸ and e¹⁰)
+ * from the following source (also available with less terms on https://en.wikipedia.org/wiki/Meridian_arc):
+ *
+ * Kawase, Kazushige (2011). A General Formula for Calculating Meridian Arc Length and its Application to Coordinate
+ * Conversion in the Gauss-Krüger Projection. Bulletin of the Geospatial Information Authority of Japan, Vol.59.
+ *
+ * That more accurate formula is implemented in MeridionalDistanceTest for comparison purposes.
+ * Then we transform that formula as below:
+ *
+ * 1) Multiply by b²/a = (1 - e²). This is done by combining some terms. For example (1 - e²)⋅(1 + ¾e²) =
+ * (1 + ¾e²) - (1e² + ¾e⁴) = 1 - ¼e² - ¾e⁴. Note that we get the first two terms of EPSG formula, which
+ * already include the multiplication by (1 - e²).
+ *
+ * 2) Application of trigonometric identities:
+ *
+ * replace: f(φ) = A⋅sin(2φ) + B⋅sin(4φ) + C⋅sin(6φ) + D⋅sin(8φ) Snyder 3-34
+ * by: f(φ) = sin(2φ)⋅(A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′ + D′⋅cos(2φ)))) Snyder 3-35
+ * with: A′ = A - C
+ * B′ = 2B - 4D
+ * C′ = 4C
+ * D′ = 8D
+ *
+ * 3) Replacement of sin2φ by 2⋅sinφ⋅cosφ and cos2φ by 1 - 2sin²φ:
+ *
+ * f(φ) = sinφ⋅cosφ⋅(C₁ + sin²φ⋅(C₂ + sin²φ⋅(C₃ + C₄⋅sin²φ)))
+ * with: C₁ = 2⋅(A′ + B′ + C′ + D′)
+ * C₂ = 2⋅(−2B′ − 4C′ − 6D′)
+ * C₃ = 2⋅(4C′ + 12D′)
+ * C₄ = 2⋅(-8D′)
+ *
+ * They are the coefficients computed below (the C₀ coefficient is not affected by trigonometric identities).
+ * We add the terms from smallest values to largest values for accuracy reasons. Most denominators are power
+ * of 2, which result in exact representations in IEEE 754 double type. Derivation from Kawase formula can be
+ * viewed at: https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods
+ */
+ c0 = -441./0x10000 * e10 + -175./0x4000 * e8 + -5./0x100 * e6 + -3./0x40 * e4 + -1./4 * e2 + 1;
+ c1 = 1575./0x20000 * e10 + 175./0x4000 * e8 + 5./0x100 * e6 + 3./0x40 * e4 + -3./4 * e2;
+ c2 = -2625./0x8000 * e10 + 175./0x6000 * e8 + 5120./0x60000 * e6 + -15./0x20 * e4;
+ c3 = 735./0x800 * e10 + 2240./0x60000 * e8 + -35./0x60 * e6;
+ c4 = -2205./0x1000 * e10 + -315./0x400 * e8;
+ }
+
+ /**
+ * Calculates the meridian distance (M). This is the distance along the central meridian from the equator to {@code φ}.
+ * Special cases:
+ * <ul>
+ * <li>If <var>φ</var> is 0°, then this method returns 0.</li>
+ * </ul>
+ *
+ * @param φ latitude to calculate meridian distance for, in radians.
+ * @param sinφ value of sin(φ).
+ * @param cosφ value of cos(φ).
+ * @return meridian distance for the given latitude, on an ellipsoid of semi-major axis of 1.
+ */
+ final double meridianArc(final double φ, final double sinφ, final double cosφ) {
+ final double sinφcosφ = cosφ * sinφ;
+ final double sinφ2 = sinφ * sinφ;
+ return c0*φ + sinφcosφ*(c1 + sinφ2*(c2 + sinφ2*(c3 + sinφ2*c4))); // TODO: use Math.fma with JDK9.
+ }
+
+ /**
+ * Gets the derivative of this {@link #meridianArc(double, double, double)} method.
+ *
+ * @return the derivative at the specified latitude.
+ */
+ final double dM_dφ(final double sinφ2, final double cosφ2) {
+ return c0 +
+ c1 * (sinφ2 - cosφ2) + sinφ2*(
+ c2 * (sinφ2 - 3*cosφ2) + sinφ2*(
+ c3 * (sinφ2 - 5*cosφ2) + sinφ2*
+ c4 * ( 1 - 7*cosφ2)));
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
index 35ee297..2db8320 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
@@ -78,7 +78,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
/**
* Series expansion with more terms. We use this formulas as a reference for testing accuracy of the formula
- * implemented by {@link MeridionalDistanceBased} after making sure that
+ * implemented by {@link MeridionalDistanceBased#meridianArc(double, double, double)} after making sure that
* this value is in agreement with {@link #reference(double)}.
*
* <p>References:</p>
@@ -108,7 +108,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
}
/**
- * Compares {@link MeridionalDistanceBased} calculation with formulas taken as references.
+ * Compares {@link MeridionalDistanceBased#meridianArc(double, double, double)} with formulas taken as references.
*/
@Test
public void compareWithReference() {
@@ -118,8 +118,9 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
final double φ = random.nextDouble() * PI - PI/2;
final double reference = reference(φ);
final double accurate = referenceMoreAccurate(φ);
- assertEquals("Accurate formula disagree with reference.", reference, accurate, 2E-10);
- // TODO: invoke MeridionalDistanceBased method.
+ final double actual = projection.meridianArc(φ, sin(φ), cos(φ));
+ assertEquals("Accurate formula disagrees with reference.", reference, accurate, 2E-10);
+ assertEquals("Implementation disagrees with reference.", accurate, actual, 1E-13);
}
}
}