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.bastiane...@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); + } + }