sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [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.
Date Tue, 01 Oct 2019 17:25:15 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

commit a8a3962f58f1a4d401448a6e7758b99e31357c99
Author: Matthieu_Bastianelli <matthieu.bastianelli@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);
+    }
+
 
 }


Mime
View raw message