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/10/01 17:25:15 UTC
[sis] 04/09: feat&refactor (referencing- satellite tracking) :
refactor to extract terms from transform/inverseTransform methods AND try
to compute Jacobian Matrix. Tests not passed for the conic satellite
tracking projection yet. Refactor should be lead to compute L coefficient
and its derivative in a method.
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
commit a8a3962f58f1a4d401448a6e7758b99e31357c99
Author: Matthieu_Bastianelli <ma...@geomatys.com>
AuthorDate: Wed Jul 24 18:31:59 2019 +0200
feat&refactor (referencing- satellite tracking) : refactor to extract terms from transform/inverseTransform methods AND try to compute Jacobian Matrix. Tests not passed for the conic satellite tracking projection yet. Refactor should be lead to compute L coefficient and its derivative in a method.
---
.../projection/ConicSatelliteTracking.java | 97 +++++++++++++---------
.../projection/CylindricalSatelliteTracking.java | 48 +++++------
.../projection/ConicSatelliteTrackingTest.java | 30 +++++++
.../CylindricalSatelliteTrackingTest.java | 27 ++++++
4 files changed, 135 insertions(+), 67 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTracking.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTracking.java
index d650fa2..9b173f1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTracking.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTracking.java
@@ -27,8 +27,10 @@ import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.OperationMethod;
import static java.lang.Math.*;
+import java.util.Map;
import org.apache.sis.internal.referencing.Resources;
import static org.apache.sis.internal.referencing.provider.SatelliteTracking.*;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.ContextualParameters;
@@ -111,12 +113,13 @@ public class ConicSatelliteTracking extends NormalizedProjection{
private final double latitudeLimit;
/**
* Boolean attribute indicating if the projection cone's constant is positive.
- * */
+ */
private final boolean positiveN;
/**
* Coefficients for the Conic Satellite-Tracking Projection.
+ * cosφ1xsinF1_n = cos(φ1)*sin(F1)/n
*/
- private final double cos_φ1xsin_F1, s0, ρ0;
+ private final double cosφ1xsinF1_n, s0, ρ0;
/**
* Work around for RFE #4093999 in Sun's bug database ("Relax constraint on
@@ -194,8 +197,6 @@ public class ConicSatelliteTracking extends NormalizedProjection{
*/
final double F1 = computeFn(cos2_φ1);
- cos_φ1xsin_F1 = cos_φ1 * sin(F1);
-
final double dλ0 = computedλn(sin_φ0);
final double dλ1 = computedλn(sin_φ1);
@@ -223,8 +224,10 @@ public class ConicSatelliteTracking extends NormalizedProjection{
} else {
n = sin_φ1 * (p2_on_p1 * (2 * cos2_i - cos2_φ1) - cos_i) / (p2_on_p1 * cos2_φ1 - cos_i); //eq. 28-17 in Snyder
}
+
+ cosφ1xsinF1_n = cos_φ1 * sin(F1)/n;
s0 = F1 - n * L1;
- ρ0 = cos_φ1xsin_F1 / (n * sin(n * L0 + s0)); // *R in eq.28-12 in Snyder
+ ρ0 = cosφ1xsinF1_n / sin(n * L0 + s0); // *R in eq.28-12 in Snyder
//======================== Unsure ======================================
// Aim to assess the limit latitude associated with -s0/n L-value.
@@ -247,7 +250,7 @@ public class ConicSatelliteTracking extends NormalizedProjection{
final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
denormalize.convertBefore(1, 1, ρ0); //For conic tracking
} else {
- n = latitudeLimit = cos_φ1xsin_F1 = s0 = ρ0 = NaN;
+ n = latitudeLimit = cosφ1xsinF1_n = s0 = ρ0 = NaN;
positiveN = false;
}
}
@@ -283,8 +286,10 @@ public class ConicSatelliteTracking extends NormalizedProjection{
}
final double sin_φ = sin(φ);
- final double dλ = -asin(sin_φ / sin_i);
- final double λt = atan(tan(dλ) * cos_i);
+ final double sinφ_sini = sin_φ / sin_i;
+ final double dλ = -asin(sinφ_sini);
+ final double tan_dλ = tan(dλ);
+ final double λt = atan(tan_dλ * cos_i);
final double L = λt - p2_on_p1 * dλ;
/*
@@ -296,34 +301,46 @@ public class ConicSatelliteTracking extends NormalizedProjection{
throw new ProjectionException(Resources.format(Resources.Keys.CanNotTransformCoordinates_2));
}
- final double ρ = cos_φ1xsin_F1/(n*sin(n*L+s0));
- final double θ = λ; // extracted n
+ final double nLandS0 = n*L+s0;
+ final double sin_nLandS0 = sin(nLandS0);
+
+ final double ρ = cosφ1xsinF1_n/sin_nLandS0;
+// final double ρ = 1/sin(n*L+s0);
- final double sinθ = sin(θ);
- final double cosθ = cos(θ);
+ final double sinλ = sin(λ);
+ final double cosλ = cos(λ);
if (dstPts != null) {
- dstPts[dstOff ] = ρ * sinθ; // x
-// dstPts[dstOff + 1] = ρ0 - ρ*cosθ; // y //TODO : extract ρ0 when ensuring : λ = λ - λ0;
- dstPts[dstOff + 1] = - ρ*cosθ;
+ dstPts[dstOff ] = ρ * sinλ; // x
+ dstPts[dstOff + 1] = - ρ*cosλ; // y
}
if (!derivate) {
return null;
}
- //=========================TO Resolve =================================
-// final double dx_dλ = ρ*n*cosθ;
-// final double dx_dφ =?;
-//
-// final double dy_dλ = ρ*n*sinθ;
-// final double dy_dφ = ?;
+ //=========================To check the resolution =====================
+ final double dρ_dφ = cosφ1xsinF1_n * (-1 / (sin_nLandS0*sin_nLandS0)) *cos(nLandS0) * n
+ // dL/dφ :
+ * ((cos(φ) / sin_i) *(1/sqrt(1-sinφ_sini*sinφ_sini))) // derivative of (-dλ/dφ)
+ * ( p2_on_p1 - ((1+tan_dλ*tan_dλ)*cos_i/(1+λt*λt) ) );
+
+ final double dx_dλ = ρ*cosλ;
+ final double dx_dφ = sinλ * dρ_dφ;
+
+ final double dy_dλ = ρ*sinλ;
+ final double dy_dφ = -cosλ * dρ_dφ;
//======================================================================
- throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates.
+ return new Matrix2(dx_dλ, dy_dλ,
+ dx_dφ, dy_dφ);
+// throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates.
- //Additionally we can compute the scale factors :
- // k = ρ*n/cos(φ); // /R
- // h = k*tan(F)/tan(n*L+s0);
+ /* =====================================================================
+ * Uncomputed scale factors :
+ *===========================
+ * k = ρ*n/cos(φ); // /R
+ * h = k*tan(F)/tan(n*L+s0);
+ ===================================================================== */
}
/**
@@ -338,11 +355,8 @@ public class ConicSatelliteTracking extends NormalizedProjection{
throws ProjectionException {
final double x = srcPts[srcOff];
-// final double y = srcPts[srcOff + 1] - ρ0;
final double y = srcPts[srcOff + 1];
- //TODO : extract - ρ0 : MatrixSIS convertBefore and convertAfter??
-
if ((x== 0) && (y == 0)) {
//TODO : which values does it imply?
throw new UnsupportedOperationException("Not supported yet for those coordinates."); //To change body of generated methods, choose Tools | Templates.
@@ -350,13 +364,22 @@ public class ConicSatelliteTracking extends NormalizedProjection{
final double ρ = positiveN ? hypot(x,y) : -hypot(x,y);
final double θ = atan(x/(-y) ); //undefined if x=y=0
- final double L = (asin(cos_φ1xsin_F1/(ρ*n)) -s0)/n; //undefined if x=y=0 //eq.28-26 in Snyder with R=1
+ final double L = (asin(cosφ1xsinF1_n/ρ) -s0)/n; //undefined if x=y=0 //eq.28-26 in Snyder with R=1
- //TODO ensure that λ0 will be added. ; In eq. Snyder 28-23 : λ0 + θ/n
dstPts[dstOff ] = θ ;//λ
dstPts[dstOff+1] = latitudeFromNewtonMethod(L); //φ
}
+ /**
+ * Return an approximation of the latitude associated with L coefficient
+ * by applying Newton-Raphson method.
+ *
+ * Eq. 28-24, 28-25 and then 28-22 in Snyder's Manual.
+ *
+ * @param l the L coefficient used in (cylindrical and conic)
+ * satellite-tracking projections.
+ * @return Approximation of the associated latitude.
+ */
protected final double latitudeFromNewtonMethod(final double l){
double dλ = -PI/2;
double dλn = Double.MIN_VALUE;
@@ -371,17 +394,17 @@ public class ConicSatelliteTracking extends NormalizedProjection{
}
dλn = dλ;
+ // Alternative calculation with Snyder's eq. 28-20 and 28-21
// λt = l + p2_on_p1 * dλ_n;
// dλ = atan(tan(λt) / cos_i);
- A = tan(l + p2_on_p1*dλn) / cos_i;
+ A = tan(l + p2_on_p1 * dλn) / cos_i;
A2=A*A;
Δdλ = -(dλn-atan(A)) / (1- (A2 + 1/cos2_i) * (p2_on_p1*cos_i/(A2+1)) );
dλ = dλn + Δdλ ;
count++;
}
-// λt = L + p2_on_p1 * dλ;
final double sin_dλ = sin(dλ);
return -asin(sin_dλ * sin_i);
}
@@ -421,13 +444,5 @@ public class ConicSatelliteTracking extends NormalizedProjection{
private double computeλtn(final double dλn) {
return atan(tan(dλn) * cos_i); // eq.28-3a in Snyder
}
-// /**
-// * Radius of the circle radius of the circle to which groundtracks
-// * are tangent on the map.
-// *
-// * @return radius ρs.
-// */
-// public double getRadiusOfTangencyCircle(){
-// return ρs;
-// }
+
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTracking.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTracking.java
index 13f8521..d68ec87 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTracking.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTracking.java
@@ -25,6 +25,7 @@ import org.opengis.referencing.operation.OperationMethod;
import static java.lang.Math.*;
import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.ContextualParameters;
@@ -76,13 +77,6 @@ import org.apache.sis.referencing.operation.transform.ContextualParameters;
*/
public class CylindricalSatelliteTracking extends ConicSatelliteTracking {
-
- /**
- * F1' in Snyder : tangente of the angle on both the globe and on the map between
- * the groundtrack and the meridian at latitude φ1.
- */
- final double dF1;
-
/**
* Create a Cylindrical Satellite Tracking Projection from the given
* parameters.
@@ -101,15 +95,13 @@ public class CylindricalSatelliteTracking extends ConicSatelliteTracking {
super(initializer);
final double cos2_φ1 = cos_φ1 * cos_φ1;
- dF1 = (p2_on_p1 * cos2_φ1 - cos_i) / sqrt(cos2_φ1 - cos2_i);
+ final double cosφ1_dF1 = sqrt(cos2_φ1 - cos2_i) *cos_φ1 / (p2_on_p1 * cos2_φ1 - cos_i);
final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
-// final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
-
normalize .convertAfter (0, cos_φ1, null); //For conic tracking
-// denormalize.convertBefore(0, 1/PI, null);
-// denormalize.convertBefore(1, , null);
+ final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
+ denormalize.convertBefore(1, cosφ1_dF1, null);
}
@@ -139,17 +131,17 @@ public class CylindricalSatelliteTracking extends ConicSatelliteTracking {
throw new ProjectionException(Resources.format(Resources.Keys.CanNotTransformCoordinates_2));
}
-// final double cos_φ = cos(φ);
-// final double cos2_φ = cos_φ * cos_φ;
final double sin_φ = sin(φ);
+ final double sinφ_sini = sin_φ / sin_i;
final double dλ = -asin(sin_φ / sin_i);
+ final double tan_dλ = tan(dλ);
final double λt = atan(tan(dλ) * cos_i);
final double L = λt - p2_on_p1 * dλ;
if (dstPts != null) {
dstPts[dstOff ] = λ; // In eq. Snyder 28-5 : R(λ - λ0) cos_φ1
- dstPts[dstOff + 1] = L * cos_φ1 / dF1; // In eq. Snyder 28-6 : R L cos(φ1)/F'1
+ dstPts[dstOff + 1] = L; // In eq. Snyder 28-6 : R L cos(φ1)/F'1
}
/* =====================================================================
@@ -159,19 +151,25 @@ public class CylindricalSatelliteTracking extends ConicSatelliteTracking {
* the meridian at latitude φ
* final double dF = (p2_on_p1 * cos2_φ - cos_i) / sqrt(cos2_φ - cos2_i);
* k = cos_φ1/cos_φ; // Parallel eq. Snyder 28-7
- * h = k* dF / dF1; // Meridian eq. Snyder 28-8
+ * h = k* dF / cosφ1_dF1; // Meridian eq. Snyder 28-8
===================================================================== */
if (!derivate) {
return null;
}
-// final double dx_dλ = cos_φ1; //*R
-// final double dx_dφ = 0;
-//
-// final double dy_dλ = 0;
-// final double dy_dφ =
+ //=========================To check the resolution =====================
+ final double dx_dλ = 1; //*R
+ final double dx_dφ = 0;
+
+ final double dy_dλ = 0;
+ final double dy_dφ = ((cos(φ) / sin_i) *(1/sqrt(1-sinφ_sini*sinφ_sini))) // derivative of (-dλ/dφ)
+ * ( p2_on_p1 - ((1+tan_dλ*tan_dλ)*cos_i/(1+λt*λt) ) );
+ //======================================================================
- throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates.
+ return new Matrix2(dx_dλ, dy_dλ,
+ dx_dφ, dy_dφ);
+
+// throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates.
}
/**
@@ -186,12 +184,10 @@ public class CylindricalSatelliteTracking extends ConicSatelliteTracking {
throws ProjectionException {
final double x = srcPts[srcOff];
- final double y = srcPts[srcOff + 1];
-
- final double L = y * dF1 / cos_φ1; // In eq. Snyder 28-19 : L = y * dF1 / R . cos_φ1
+ final double y = srcPts[srcOff + 1]; // In eq. Snyder 28-19 : y = yinit * cosφ1_dF1 / R . cos_φ1
dstPts[dstOff ] = x;
- dstPts[dstOff+1] = latitudeFromNewtonMethod(L); //φ
+ dstPts[dstOff+1] = latitudeFromNewtonMethod(y); //φ
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTrackingTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTrackingTest.java
index 7152538..242e024 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTrackingTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTrackingTest.java
@@ -33,6 +33,7 @@ import org.opengis.referencing.operation.TransformException;
import org.opengis.util.FactoryException;
import static java.lang.StrictMath.toRadians;
+import org.apache.sis.test.DependsOnMethod;
import static org.junit.Assert.assertTrue;
/**
@@ -243,4 +244,33 @@ public class ConicSatelliteTrackingTest extends MapProjectionTestCase {
0.21642 * sin(n * (toRadians(-120 - -90))), 2.46908 //Tracking limit
});
}
+
+
+ /**
+ * 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
+ @DependsOnMethod("testInverseDerivative")
+ public void testDerivativeOnSphere() throws FactoryException, TransformException {
+ createProjection(
+ 99.092, //satellite_orbit_inclination
+ 103.267, //satellite_orbital_period
+ 1440.0, //ascending_node_period
+ -90, //central_meridian
+ 45, //standard_parallel_1
+ 70, //standard_parallel_2
+ 30 //latitude_of_origin
+ );
+ final double delta = (1.0 / 60) / 1852; // Approximately 1 metre.
+ derivativeDeltas = new double[] {delta, delta};
+ tolerance = Formulas.LINEAR_TOLERANCE / 10;
+ verifyDerivative( 0, -10);
+ verifyDerivative(-100, 3);
+ verifyDerivative( -56, 50);
+ verifyDerivative( -20, 47);
+ }
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTrackingTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTrackingTest.java
index 02c1b09..844ba7f 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTrackingTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTrackingTest.java
@@ -20,6 +20,7 @@ package org.apache.sis.referencing.operation.projection;
import static java.lang.Double.NaN;
import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.test.DependsOnMethod;
import static org.junit.Assert.assertTrue;
import org.junit.Test;
import org.opengis.referencing.operation.TransformException;
@@ -205,5 +206,31 @@ public class CylindricalSatelliteTrackingTest extends ConicSatelliteTrackingTest
}
+ /**
+ * 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
+ @DependsOnMethod("testInverseDerivative")
+ @Override
+ public void testDerivativeOnSphere() throws FactoryException, TransformException {
+ createProjection(
+ 99.092, //satellite_orbit_inclination
+ 103.267, //satellite_orbital_period
+ 1440.0, //ascending_node_period
+ -90, //central_meridian
+ 30 //standard_parallel_1
+ );
+ final double delta = (1.0 / 60) / 1852; // Approximately 1 metre.
+ derivativeDeltas = new double[] {delta, delta};
+ tolerance = Formulas.LINEAR_TOLERANCE / 10;
+ verifyDerivative(-100, 3);
+ verifyDerivative( -56, 50);
+ verifyDerivative( -20, 47);
+ }
+
}