sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Add skeleton base class for meridional distance calculation on ellipsoid. Add formulas to be used as reference in the test class (the formula actually used by implementation may differ). https://issues.apache.org/jira/browse/SIS-450
Date Mon, 22 Apr 2019 15:59:51 GMT
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 9566eb7  Add skeleton base class for meridional distance calculation on ellipsoid.
Add formulas to be used as reference in the test class (the formula actually used by implementation
may differ). https://issues.apache.org/jira/browse/SIS-450
9566eb7 is described below

commit 9566eb7e7894483da943682782124e19f9895d1c
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Apr 22 17:54:28 2019 +0200

    Add skeleton base class for meridional distance calculation on ellipsoid.
    Add formulas to be used as reference in the test class (the formula actually used by implementation
may differ).
    https://issues.apache.org/jira/browse/SIS-450
---
 .../projection/MeridionalDistanceBased.java        |  95 ++++++++++++++++
 .../operation/projection/Sinusoidal.java           |   2 +-
 .../projection/MeridionalDistanceTest.java         | 125 +++++++++++++++++++++
 .../operation/projection/SinusoidalTest.java       |   2 +-
 .../sis/test/suite/ReferencingTestSuite.java       |   1 +
 5 files changed, 223 insertions(+), 2 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
new file mode 100644
index 0000000..0d2b2ed
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
@@ -0,0 +1,95 @@
+/*
+ * 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.io.IOException;
+import java.io.ObjectInputStream;
+
+
+/**
+ * Base class of map projections based on distance along the meridian from equator to latitude
φ.
+ *
+ * @author  Martin Desruisseaux (MPO, IRD, Geomatys)
+ * @author  Rémi Maréchal (Geomatys)
+ * @version 1.0
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Meridian_arc">Meridian arc on Wikipedia</a>
+ *
+ * @since 1.0
+ * @module
+ */
+abstract class MeridionalDistanceBased extends NormalizedProjection {
+    /**
+     * {@code false} for using the original formulas as published by EPSG, or {@code true}
for using formulas
+     * modified using trigonometric identities. Use of trigonometric identities reduces the
amount of calls to
+     * {@link Math#sin(double)} and similar methods. Snyder 3-34 to 3-39 give the following
identities:
+     *
+     * <pre>
+     *     If:     f(φ) = A⋅sin(2φ) + B⋅sin(4φ) + C⋅sin(6φ) + D⋅sin(8φ)
+     *     Then:   f(φ) = sin(2φ)⋅(A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′ + D′⋅cos(2φ))))
+     *     Where:  A′ = A - C
+     *             B′ = 2B - 4D
+     *             C′ = 4C
+     *             D′ = 8D
+     * </pre>
+     *
+     * Similar, but with cosine instead than sin and the addition of a constant:
+     *
+     * <pre>
+     *     If:     f(φ) = A + B⋅cos(2φ) + C⋅cos(4φ) + D⋅cos(6φ) + E⋅cos(8φ)
+     *     Then:   f(φ) = A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′ + cos(2φ)⋅(D′
+ E′⋅cos(2φ))))
+     *     Where:  A′ = A - C + E
+     *             B′ = B - 3D
+     *             C′ = 2C - 8E
+     *             D′ = 4D
+     *             E′ = 8E
+     * </pre>
+     */
+    private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
+
+    /**
+     * Creates a new normalized projection from the parameters computed by the given initializer.
+     */
+    MeridionalDistanceBased(final Initializer initializer) {
+        super(initializer);
+        computeCoefficients();
+    }
+
+    /**
+     * Creates a new projection initialized to the same parameters than the given one.
+     */
+    MeridionalDistanceBased(final MeridionalDistanceBased other) {
+        super(other);
+    }
+
+    /**
+     * Restores transient fields after deserialization.
+     */
+    private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException
{
+        in.defaultReadObject();
+        computeCoefficients();
+    }
+
+    /**
+     * Computes the coefficients in the series expansions from the {@link #eccentricitySquared}
value.
+     * This method shall be invoked after {@code MeridionalDistanceBased} construction or
deserialization.
+     */
+    private void computeCoefficients() {
+        final double e2 = eccentricitySquared;
+        // TODO
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
index fc94ef6..64b5fac 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Sinusoidal.java
@@ -42,7 +42,7 @@ import static org.apache.sis.internal.referencing.provider.Sinusoidal.*;
  * @since   1.0
  * @module
  */
-public class Sinusoidal extends NormalizedProjection {
+public class Sinusoidal extends MeridionalDistanceBased {
     /**
      * For cross-version compatibility.
      */
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
new file mode 100644
index 0000000..35ee297
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
@@ -0,0 +1,125 @@
+/*
+ * 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.Random;
+import org.apache.sis.referencing.operation.DefaultOperationMethod;
+import org.apache.sis.test.TestUtilities;
+import org.apache.sis.test.DependsOn;
+import org.junit.Test;
+
+import static java.lang.StrictMath.*;
+import static org.junit.Assert.assertEquals;
+
+
+/**
+ * Tests the meridional distances computed by {@link MeridionalDistanceBased}.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since   1.0
+ * @module
+ */
+@DependsOn(NormalizedProjectionTest.class)
+public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase {
+    /**
+     * Creates a new test case.
+     */
+    public MeridionalDistanceTest() {
+        final DefaultOperationMethod provider = new org.apache.sis.internal.referencing.provider.Sinusoidal();
+        transform = new Sinusoidal(provider, parameters(provider, true));
+        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;     // = linear tolerance on
a sphere of radius 1.
+    }
+
+    /**
+     * Computes the meridional distance using equation given in EPSG guidance notes, which
is also from Snyder book.
+     * The equation is given in Snyder 3-21. We use this equation as a reference for testing
validity of other forms.
+     * The equation is:
+     *
+     * {@preformat math
+     *   M = a[(1 – e²/4 – 3e⁴/64  –  5e⁶/256  – …)⋅φ
+     *          – (3e²/8 + 3e⁴/32  + 45e⁶/1024 + …)⋅sin2φ
+     *                 + (15e⁴/256 + 45e⁶/1024 + …)⋅sin4φ
+     *                            – (35e⁶/3072 + …)⋅sin6φ
+     *                                         + …]
+     * }
+     *
+     * @param  φ  latitude in radians.
+     * @return meridional distance from equator to the given latitude on an ellipsoid with
semi-major axis of 1.
+     */
+    private double reference(final double φ) {
+        final MeridionalDistanceBased projection = (MeridionalDistanceBased) transform;
+        final double e2 = projection.eccentricitySquared;
+        final double e4 = e2*e2;
+        final double e6 = e2*e4;
+        /*
+         * Smallest terms of the series should be added first for better accuracy.
+         * The final multiplication by (1 - e2) is already included in all terms.
+         */
+        return - (35./3072*e6)                        * sin(6*φ)
+               + (45./1024*e6 + 15./256*e4)           * sin(4*φ)
+               - (45./1024*e6 +  3./32 *e4 + 3./8*e2) * sin(2*φ)
+               + (-5./256 *e6 -  3./64 *e4 - 1./4*e2 + 1)  *  φ;
+    }
+
+    /**
+     * 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
+     * this value is in agreement with {@link #reference(double)}.
+     *
+     * <p>References:</p>
+     * <ul>
+     *   <li>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.
+     *       <a href="http://www.gsi.go.jp/common/000062452.pdf">(download)</a></li>
+     *   <li><a href="https://en.wikipedia.org/wiki/Meridian_arc#Series_expansions">Meridian
arc — series expansions</a>
+     *       on Wikipedia.</li>
+     * </ul>
+     */
+    private double referenceMoreAccurate(final double φ) {
+        final MeridionalDistanceBased projection = (MeridionalDistanceBased) transform;
+        final double e2  = projection.eccentricitySquared;
+        final double e4  = e2*e2;
+        final double e6  = e2*e4;
+        final double e8  = e4*e4;
+        final double e10 = e2*e8;
+        final double C1  = 43659./ 65536*e10 + 11025./16384*e8 + 175./256*e6 + 45./64*e4
+ 3./4*e2 + 1;
+        final double C2  = 72765./ 65536*e10 +  2205./ 2048*e8 + 525./512*e6 + 15./16*e4
+ 3./4*e2;
+        final double C3  = 10395./ 16384*e10 +  2205./ 4096*e8 + 105./256*e6 + 15./64*e4;
+        final double C4  = 31185./131072*e10 +   315./ 2048*e8 +  35./512*e6;
+        final double C5  =  3465./ 65536*e10 +   315./16384*e8;
+        final double C6  =   693./131072*e10;
+        final double v   = -C6*sin(10*φ)/10 + C5*sin(8*φ)/8 - C4*sin(6*φ)/6 + C3*sin(4*φ)/4
- C2*sin(2*φ)/2 + C1*φ;
+        return v * (1 - e2);
+    }
+
+    /**
+     * Compares {@link MeridionalDistanceBased} calculation with formulas taken as references.
+     */
+    @Test
+    public void compareWithReference() {
+        final MeridionalDistanceBased projection = (MeridionalDistanceBased) transform;
+        final Random random = TestUtilities.createRandomNumberGenerator();
+        for (int i=0; i<100; i++) {
+            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.
+        }
+    }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
index acd8869..a51ebcb 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
@@ -34,7 +34,7 @@ import static org.junit.Assert.*;
  * @since   1.0
  * @module
  */
-@DependsOn(NormalizedProjectionTest.class)
+@DependsOn(MeridionalDistanceTest.class)
 public final strictfp class SinusoidalTest extends MapProjectionTestCase {
     /**
      * Creates a new instance of {@link Sinusoidal} concatenated with the (de)normalization
matrices.
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 972fcdb..e72ed1f 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -175,6 +175,7 @@ import org.junit.BeforeClass;
     org.apache.sis.referencing.operation.projection.ObliqueMercatorTest.class,
     org.apache.sis.referencing.operation.projection.CylindricalEqualAreaTest.class,
     org.apache.sis.referencing.operation.projection.AlbersEqualAreaTest.class,
+    org.apache.sis.referencing.operation.projection.MeridionalDistanceTest.class,
     org.apache.sis.referencing.operation.projection.SinusoidalTest.class,
     org.apache.sis.referencing.operation.projection.MollweideTest.class,
 


Mime
View raw message