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 <[email protected]>
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
+ }
}
