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: Implement derivate for "Azimuthal Equidistant (Spherical)" projection.
Date Fri, 20 Mar 2020 21:57:19 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 0045413  Implement derivate for "Azimuthal Equidistant (Spherical)" projection.
0045413 is described below

commit 0045413191f4d6b9d75470945b087abd870a0407
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri Mar 20 22:55:40 2020 +0100

    Implement derivate for "Azimuthal Equidistant (Spherical)" projection.
---
 .../operation/projection/AzimuthalEquidistant.java | 52 +++++++++++++++++-----
 .../projection/ModifiedAzimuthalEquidistant.java   | 20 ++-------
 .../projection/AzimuthalEquidistantTest.java       | 46 ++++++++++++++++++-
 .../ModifiedAzimuthalEquidistantTest.java          | 10 +++++
 4 files changed, 99 insertions(+), 29 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistant.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistant.java
index c613410..6bd1b4a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistant.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistant.java
@@ -20,6 +20,7 @@ import java.util.EnumMap;
 import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.OperationMethod;
 import org.opengis.parameter.ParameterDescriptor;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.util.Workaround;
 
@@ -103,6 +104,11 @@ public class AzimuthalEquidistant extends NormalizedProjection {
     /**
      * Converts the specified (λ,φ) coordinate and stores the (<var>x</var>,<var>y</var>)
result in {@code dstPts}.
      *
+     * @param  srcPts    source point coordinate, as (<var>longitude</var>, <var>latitude</var>)
in radians.
+     * @param  srcOff    the offset of the single coordinate to be converted in the source
array.
+     * @param  dstPts    the array into which the converted coordinate is returned (may be
the same than {@code srcPts}).
+     * @param  dstOff    the offset of the location of the converted coordinate that is stored
in the destination array.
+     * @param  derivate  {@code true} for computing the derivative, or {@code false} if not
needed.
      * @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.
@@ -112,27 +118,51 @@ public class AzimuthalEquidistant extends NormalizedProjection {
                             final double[] dstPts, final int dstOff,
                             final boolean derivate) throws ProjectionException
     {
-        final double λ    = srcPts[srcOff  ];
-        final double φ    = srcPts[srcOff+1];
-        final double cosλ = cos(λ);
-        final double sinλ = sin(λ);
-        final double cosφ = cos(φ);
-        final double sinφ = sin(φ);
-        final double c    = acos(sinφ0*sinφ + cosφ0*cosφ*cosλ);
-        final double k    = abs(c) >= ANGULAR_TOLERANCE ? c/sin(c) : 1;
+        final double  λ     = srcPts[srcOff  ];
+        final double  φ     = srcPts[srcOff+1];
+        final double  cosλ  = cos(λ);
+        final double  sinλ  = sin(λ);
+        final double  cosφ  = cos(φ);
+        final double  sinφ  = sin(φ);
+        final double  cosc  = sinφ0*sinφ + cosφ0*cosφ*cosλ;
+        final double  c     = acos(cosc);
+        final boolean ind   = abs(c) < ANGULAR_TOLERANCE;
+        final double  k     = ind ? 1 : c/sin(c);
+        final double  cφcλ  = cosφ * cosλ;
+        final double  cφsλ  = cosφ * sinλ;
+        final double  x     = k * cφsλ;
+        final double  y     = k * (cosφ0*sinφ - sinφ0*cφcλ);
         if (dstPts != null) {
-            dstPts[dstOff  ] = k * cosφ*sinλ;
-            dstPts[dstOff+1] = k * (cosφ0*sinφ - sinφ0*cosφ*cosλ);
+            dstPts[dstOff  ] = x;
+            dstPts[dstOff+1] = y;
         }
         if (!derivate) {
             return null;
         }
-        throw new ProjectionException("Derivative not yet implemented.");
+        /*
+         * Formulas below can be verified with Maxima.
+         *
+         * https://svn.apache.org/repos/asf/sis/analysis/Azimuthal%20Equidistant%20(Spherical).wxmx
+         */
+        final double t    = ind ? 1./3 : (1/k - cosc) / (1 - cosc*cosc);
+        final double sφcλ = sinφ * cosλ;
+        final double tλ   = cφsλ * cosφ0*t;
+        final double tφ   = (cosφ0*sφcλ - sinφ0*cosφ) * t;
+        return new Matrix2(x*tλ + k*cφcλ,                           // ∂x/∂λ
+                           x*tφ - k*sinλ*sinφ,                      // ∂x/∂φ
+                           y*tλ + x*sinφ0,                          // ∂y/∂λ
+                           y*tφ + k*(sinφ0*sφcλ + cosφ0*cosφ));     // ∂y/∂φ
     }
 
     /**
      * Converts the specified (<var>x</var>,<var>y</var>) coordinates
      * and stores the result in {@code dstPts} (angles in radians).
+     *
+     * @param  srcPts  the array containing the source point coordinate, as linear distance
on a unit sphere or ellipse.
+     * @param  srcOff  the offset of the point to be converted in the source array.
+     * @param  dstPts  the array into which the converted point coordinate is returned (may
be the same than {@code srcPts}).
+     * @param  dstOff  the offset of the location of the converted point that is stored in
the destination array.
+     * @throws ProjectionException if the point can not be converted.
      */
     @Override
     protected void inverseTransform(final double[] srcPts, final int srcOff,
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
index b957554..5e26b4c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
@@ -23,8 +23,6 @@ import org.opengis.parameter.ParameterDescriptor;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.referencing.operation.transform.ContextualParameters;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
-import org.apache.sis.internal.util.Numerics;
-import org.apache.sis.util.ComparisonMode;
 import org.apache.sis.util.Workaround;
 
 import static java.lang.Math.*;
@@ -118,10 +116,11 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
     @Workaround(library="JDK", version="1.8")
     private ModifiedAzimuthalEquidistant(final Initializer initializer) {
         super(initializer);
-        final double ν0, f;
+        final double axisRatio, ν0, f;
+        axisRatio   = initializer.axisLengthRatio().doubleValue();
         ν0          = initializer.radiusOfCurvature(sinφ0);
         ℯ2_ν0_sinφ0 = eccentricitySquared * ν0 * sinφ0;
-        f           = eccentricity / sqrt(1 - eccentricitySquared);
+        f           = eccentricity / axisRatio;                 // √(1 - ℯ²) = b/a
         G           = f * sinφ0;
         Hp          = f * cosφ0;
         Bp          = 3*eccentricitySquared * (sinφ0*cosφ0) / (1 - eccentricitySquared);
@@ -181,6 +180,7 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
                + (s3/8   *  GH*(1 - 2*H2))
                + (s4/120 * (H2*(4 - 7*H2) - 3*(G*G)*(1 - 7*H2)))
                - (s5/48  * GH);
+
         if (dstPts != null) {
             dstPts[dstOff  ] = c * sinα;
             dstPts[dstOff+1] = c * cosα;
@@ -220,16 +220,4 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
         dstPts[dstOff+1]  = atan((1 - eccentricitySquared*sinφ0*K / sinΨ) * (sinΨ/cosΨ)
                                / (1 - eccentricitySquared));
     }
-
-    /**
-     * Compares the given object with this transform for equivalence.
-     */
-    @Override
-    public boolean equals(final Object object, final ComparisonMode mode) {
-        if (super.equals(object, mode)) {
-            final ModifiedAzimuthalEquidistant that = (ModifiedAzimuthalEquidistant) object;
-            return Numerics.epsilonEqual(ℯ2_ν0_sinφ0, that.ℯ2_ν0_sinφ0, mode);
-        }
-        return false;
-    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
index 717eb29..c716b81 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
@@ -16,11 +16,12 @@
  */
 package org.apache.sis.referencing.operation.projection;
 
+import org.opengis.util.FactoryException;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.internal.referencing.provider.MapProjection;
+import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.test.DependsOn;
 import org.junit.Test;
-import org.opengis.util.FactoryException;
 
 
 /**
@@ -108,7 +109,9 @@ public strictfp class AzimuthalEquidistantTest extends MapProjectionTestCase
{
                 Double.NaN,                     // Scale factor
                 40000,                          // False easting
                 60000);                         // False Northing
-
+        /*
+         * Test point given in EPSG guidance note.
+         */
         verifyTransform(new double[] {
             138 + (11 + 34.908/60)/60,          // 138°11'34.908"E
               9 + (35 + 47.493/60)/60           //   9°35'47.493"N
@@ -116,5 +119,44 @@ public strictfp class AzimuthalEquidistantTest extends MapProjectionTestCase
{
             42665.90,
             65509.82
         });
+        /*
+         * North of map origin, for entering in the special case for c/sin(c)
+         * when c is close to zero. This point is not given by EPSG guidance
+         * notes; this is an anti-regression test.
+         */
+        verifyTransform(new double[] {
+            138 + (10 +  7.48/60)/60 + 0.00000000001,
+              9 + (32 + 48.15/60)/60 + 0.01
+        }, new double[] {
+            40000.00,
+            61105.98
+        });
+    }
+
+    /**
+     * Tests the derivatives at a few points on a sphere. This method compares the derivatives
computed
+     * by the projection with an estimation of derivatives computed by the finite differences
method.
+     *
+     * @throws FactoryException if an error occurred while creating the map projection.
+     * @throws TransformException if an error occurred while projecting a point.
+     */
+    @Test
+    public void testDerivative() throws FactoryException, TransformException {
+        createCompleteProjection(method(),
+                CLARKE_A,
+                CLARKE_B,
+                 40,                // Longitude of natural origin (central-meridian)
+                 25,                // Latitude of natural origin
+                Double.NaN,         // Standard parallel 1
+                Double.NaN,         // Standard parallel 2
+                Double.NaN,         // Scale factor
+                40000,              // False easting
+                60000);             // False Northing
+        final double delta = (1.0 / 60) / 1852;                 // Approximately 1 metre.
+        derivativeDeltas = new double[] {delta, delta};
+        tolerance = Formulas.LINEAR_TOLERANCE / 100;
+        verifyDerivative(30, 27);
+        verifyDerivative(27, 20);
+        verifyDerivative(40, 25);
     }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistantTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistantTest.java
index 1001168..37f2065 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistantTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistantTest.java
@@ -83,4 +83,14 @@ public final strictfp class ModifiedAzimuthalEquidistantTest extends AzimuthalEq
     public void runGeoapiTest() throws FactoryException, TransformException {
         createGeoApiTest(method()).testModifiedAzimuthalEquidistant();
     }
+
+    /**
+     * Not yet supported.
+     */
+    @Test
+    @Override
+    @org.junit.Ignore("Implementation not yet completed")
+    public void testDerivative() throws FactoryException, TransformException {
+        // TODO
+    }
 }


Mime
View raw message