sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 01/09: feat (sis-referencing & internal) : implement cylindrical and conic satellite-tracking projections from Snyder's 'Map Projections - a working Manual'. Don't compute derivatives yet. Very few tests.
Date Tue, 01 Oct 2019 17:25:12 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 d6502364fe43513d836c23c2016844d2ca315531
Author: Matthieu_Bastianelli <matthieu.bastianelli@geomatys.com>
AuthorDate: Tue Jul 23 09:31:22 2019 +0200

    feat (sis-referencing & internal) : implement cylindrical and conic satellite-tracking
projections from Snyder's 'Map Projections - a working Manual'. Don't compute derivatives
yet. Very few tests.
---
 .../referencing/provider/SatelliteTracking.java    | 200 +++++++++++
 .../projection/ConicSatelliteTracking.java         | 398 +++++++++++++++++++++
 .../projection/CylindricalSatelliteTracking.java   | 197 ++++++++++
 .../projection/ConicSatelliteTrackingTest.java     | 127 +++++++
 .../CylindricalSatelliteTrackingTest.java          | 117 ++++++
 5 files changed, 1039 insertions(+)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/SatelliteTracking.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/SatelliteTracking.java
new file mode 100644
index 0000000..3e9382d
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/SatelliteTracking.java
@@ -0,0 +1,200 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+package org.apache.sis.internal.referencing.provider;
+
+import javax.xml.bind.annotation.XmlTransient;
+import static org.apache.sis.internal.referencing.provider.AbstractProvider.createZeroConstant;
+import static org.apache.sis.internal.referencing.provider.Mercator1SP.LATITUDE_OF_ORIGIN;
+import org.apache.sis.measure.Units;
+import org.apache.sis.metadata.iso.citation.Citations;
+import org.apache.sis.parameter.ParameterBuilder;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.referencing.operation.projection.ConicSatelliteTracking;
+import org.apache.sis.referencing.operation.projection.CylindricalSatelliteTracking;
+import org.apache.sis.referencing.operation.projection.NormalizedProjection;
+import org.apache.sis.util.ArgumentChecks;
+import org.opengis.parameter.ParameterDescriptor;
+import org.opengis.parameter.ParameterDescriptorGroup;
+import org.opengis.parameter.ParameterNotFoundException;
+
+/**
+ * The provider for <cite>Satellite-Tracking projections"</cite>.
+ *
+ * <cite>
+ * - All groundtracks for satellites orbiting the Earth with the same orbital
+ * parameters are shown as straight lines on the map.
+ *
+ * - Cylindrical {@link CylindricalSatelliteTracking} or conical
+ * {@link ConicSatelliteTracking} form available.
+ *
+ * - Neither conformal nor equal-area.
+ *
+ * - All meridians are equally spaced straight lines, parallel on cylindrical
+ * form and converging to a common point on conical form.
+ *
+ * - All parallels are straight and parallel on cylindrical form and are
+ * concentric circular arcs on conical form. Parallels are unequally spaced.
+ *
+ * - Conformality occurs along two chosen parallels. Scale is correct along one
+ * of these parameters on the conical form and along both on the cylindrical
+ * form.
+ *
+ * Developed 1977 by Snyder
+ * </cite>
+ *
+ * <cite> These formulas are confined to circular orbits and the SPHERICAL
+ * Earth.</cite>
+ *
+ * <cite>The ascending and descending groundtracks meet at the northern an
+ * southern tracking limits, lats. 80.9°N and S for landsat 1, 2 and 3. The map
+ * Projection does not extend closer to the poles.</cite>
+ *
+ * This projection method has no associated EPSG code.
+ *
+ * Earth radius is normalized. Its value is 1 and is'nt an input parameter.
+ *
+ * =============================================================================
+ * REMARK : The parameters associated with the satellite (and its orbit) could
+ * be aggregate in class of the kind : Satellite or SatelliteOrbit.
+ * =============================================================================
+ *
+ * @see <cite>Map Projections - A Working Manual</cite> By John P. Snyder
+ * @author Matthieu Bastianelli (Geomatys)
+ * @version 1.0
+ */
+@XmlTransient
+public class SatelliteTracking extends MapProjection {
+
+    /**
+     * Name of this projection.
+     */
+    public static final String NAME = "Satellite-Tracking";
+
+    /**
+     * The operation parameter descriptor for the <cite>Longitude of projection
+     * center</cite> (λ₀) parameter value. Valid values range is [-180 … 180]°
+     * and default value is 0°.
+     */
+    public static final ParameterDescriptor<Double> CENTRAL_MERIDIAN = ESRI.CENTRAL_MERIDIAN;
+
+    /**
+     * The operation parameter descriptor for the <cite>Latitude of 1st standard
+     * parallel</cite>. For conical satellite-tracking projection : first
+     * parallel of conformality with true scale.
+     *
+     * Valid values range is [-90 … 90]° and default value is the value given to
+     * the {@link #LATITUDE_OF_FALSE_ORIGIN} parameter.
+     */
+    public static final ParameterDescriptor<Double> STANDARD_PARALLEL_1 = LambertConformal2SP.STANDARD_PARALLEL_1;
+
+    /**
+     * The operation parameter descriptor for the <cite>second parallel of
+     * conformality but without true scale</cite> parameter value for conic
+     * projection. Valid values range is [-90 … 90]° and default value is the
+     * opposite value given to the {@link #STANDARD_PARALLEL_1} parameter.
+     */
+    public static final ParameterDescriptor<Double> STANDARD_PARALLEL_2 = LambertConformal2SP.STANDARD_PARALLEL_2;
+
+    /**
+     * Latitude Crossing the central meridian at the desired origin of
+     * rectangular coordinates.
+     */
+    public static final ParameterDescriptor<Double> LATITUDE_OF_ORIGIN = TransverseMercator.LATITUDE_OF_ORIGIN;;
+
+    /**
+     * The operation parameter descriptor for the <cite> Angle of inclination
+     * between the plane of the Earth's Equator and the plane of the satellite
+     * orbit</cite> parameter value. <cite> It's measured counterclockwise from
+     * the Equatorto the orbit plane at the ascending node (99.092° for Landsat
+     * 1, 2, 3</cite>
+     */
+    public static final ParameterDescriptor<Double> SATELLITE_ORBIT_INCLINATION;
+
+    /**
+     * The operation parameter descriptor for the <cite> time required for
+     * revolution of the satellite</cite> parameter value.
+     */
+    public static final ParameterDescriptor<Double> SATELLITE_ORBITAL_PERIOD;
+
+    /**
+     * The operation parameter descriptor for the <cite> length of earth's
+     * rotation with respect to the precessed ascending node</cite> parameter
+     * value. The ascending node is the point on the satellite orbit at which
+     * the satellite crosses the Earth's equatorial plane in a northerly
+     * direction.
+     */
+    public static final ParameterDescriptor<Double> ASCENDING_NODE_PERIOD; //Constant
for landsat but is it the case for every satellite with circular orbital?
+
+    /**
+     * The group of all parameters expected by this coordinate operation.
+     */
+    static final ParameterDescriptorGroup PARAMETERS;
+
+    static {
+        final ParameterBuilder builder = new ParameterBuilder();
+        builder.setCodeSpace(null, null)
+                .setRequired(true);
+
+        SATELLITE_ORBIT_INCLINATION = builder
+                .addName("satellite_orbit_inclination")
+                .setDescription("Angle of inclination between the plane of the Earth's Equator
and the plane of the satellite orbit")
+                .create(0, Units.RADIAN);
+
+        SATELLITE_ORBITAL_PERIOD = builder
+                .addName("satellite_orbital_period")
+                .setDescription("Time required for revolution of the satellite")
+                .createStrictlyPositive(103.267, Units.MINUTE); //Or hours? Seconds?
+
+        ASCENDING_NODE_PERIOD = builder
+                .addName("ascending_node_period")
+                .setDescription("Length of Earth's rotation with respect to the precessed
ascending node")
+                .createStrictlyPositive(98.884, Units.MINUTE);
+                
+        PARAMETERS = builder.addName(NAME)
+                .createGroupForMapProjection(CENTRAL_MERIDIAN,
+                        STANDARD_PARALLEL_1, STANDARD_PARALLEL_2,
+                        SATELLITE_ORBIT_INCLINATION, SATELLITE_ORBITAL_PERIOD,
+                        ASCENDING_NODE_PERIOD, LATITUDE_OF_ORIGIN);
+
+    }
+
+    /**
+     * Constructs a new provider.
+     */
+    public SatelliteTracking() {
+        super(PARAMETERS);
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * @return the map projection created from the given parameter values.
+     */
+    @Override
+    protected NormalizedProjection createProjection(final Parameters parameters) throws ParameterNotFoundException
{
+        ArgumentChecks.ensureNonNull("Parameters", parameters);
+
+        if (parameters.getValue(STANDARD_PARALLEL_2) == -parameters.getValue(STANDARD_PARALLEL_1))
{ 
+            return new org.apache.sis.referencing.operation.projection.CylindricalSatelliteTracking(this,
parameters);
+        } else {
+            return new org.apache.sis.referencing.operation.projection.ConicSatelliteTracking(this,
parameters);
+        }
+    }
+
+}
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
new file mode 100644
index 0000000..3c8d556
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTracking.java
@@ -0,0 +1,398 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import static java.lang.Double.NaN;
+import java.util.EnumMap;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.util.Workaround;
+import org.opengis.parameter.ParameterDescriptor;
+import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.OperationMethod;
+
+import static java.lang.Math.*;
+import java.util.function.DoubleUnaryOperator;
+import org.apache.sis.internal.referencing.Resources;
+import static org.apache.sis.internal.referencing.provider.SatelliteTracking.*;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
+import org.apache.sis.referencing.operation.transform.ContextualParameters;
+
+/**
+ * <cite>Cylindrical Satellite-Tracking projections</cite>.
+ * 
+ * <cite>
+ * - All groundtracks for satellites orbiting the Earth with the same orbital
+ * parameters are shown as straight lines on the map.
+ *
+ * - Cylindrical {@link CylindricalSatelliteTracking}
+ * or conical {@link ConicSatelliteTracking} form available. 
+ *
+ * - Neither conformal nor equal-area.
+ *
+ * - All meridians are equally spaced straight lines, parallel on cylindrical
+ * form and converging to a common point on conical form.
+ *
+ * - All parallels are straight and parallel on cylindrical form and are
+ * concentric circular arcs on conical form. Parallels are unequally spaced.
+ *
+ * - Conformality occurs along two chosen parallels. Scale is correct along one
+ * of these parameters on the conical form and along both on the cylindrical
+ * form.
+ *
+ * Developed 1977 by Snyder
+ * </cite>
+ *
+ * <cite> These formulas are confined to circular orbits and the SPHERICAL
+ * Earth.</cite>
+ *
+ * <cite>The ascending and descending groundtracks meet at the northern an
+ * southern tracking limits, lats. 80.9°N and S for landsat 1, 2 and 3. The map
+ * Projection does not extend closer to the poles.</cite>
+ *
+ * This projection method has no associated EPSG code.
+ *
+ * Earth radius is normalized. Its value is 1 and is'nt an input parameter.
+ *
+ * =============================================================================
+ * REMARK : The parameters associated with the satellite (and its orbit) could
+ * be aggregate in class of the kind : Satellite or SatelliteOrbit.
+ * =============================================================================
+ *
+ * @see <cite>Map Projections - A Working Manual</cite> By John P. Snyder
+ * @author Matthieu Bastianelli (Geomatys)
+ * @version 1.0
+ */
+public class ConicSatelliteTracking extends NormalizedProjection{
+    
+    /**
+     * {@code SATELLITE_ORBIT_INCLINATION}
+     */
+    final double i;
+    
+    /**
+     * Coefficients for both cylindrical and conic Satellite-Tracking Projection.
+     */
+    final double cos_i, sin_i, cos2_i, cos_φ1, p2_on_p1;
+    
+//    /**
+//     * Radius of the circle radius of the circle to which groundtracks
+//     * are tangent on the map.
+//     */
+//    private final double ρs;
+    
+    /**
+     * Projection Cone's constant.
+     */
+    private final double n;
+    /**
+     * Approximation of the limiting lattitude to which the {@code L coefficient}
+     * is associated with a particular case {@code L equals -s0/n}.
+     */
+    private final double lattitudeLimit;
+    /**
+     * Boolean attribute indicating if the projection cone's constant is positive.
+     * */
+    private final boolean positiveN;
+    /**
+     * Coefficients for the Conic Satellite-Tracking Projection.
+     */
+    private final double cos_φ1xsin_F1, s0, ρ0;
+
+    /**
+     * Work around for RFE #4093999 in Sun's bug database ("Relax constraint on
+     * placement of this()/super() call in constructors").
+     */
+    @Workaround(library = "JDK", version = "1.8")
+    static Initializer initializer(final OperationMethod method, final Parameters parameters)
{
+        final EnumMap<NormalizedProjection.ParameterRole, ParameterDescriptor<Double>>
roles = new EnumMap<>(NormalizedProjection.ParameterRole.class);
+        roles.put(NormalizedProjection.ParameterRole.CENTRAL_MERIDIAN, CENTRAL_MERIDIAN);
+        return new Initializer(method, parameters, roles, (byte) 0);
+    }
+
+    /**
+     * Create a Cylindrical Satellite Tracking Projection from the given
+     * parameters.
+     *
+     * The parameters are described in <cite>Map Projections - A Working
+     * Manual</cite> By John P. Snyder.
+     *
+     * @param method : description of the projection parameters.
+     * @param parameters : the parameter values of the projection to create.
+     */
+    public ConicSatelliteTracking(final OperationMethod method, final Parameters parameters)
{
+        this(initializer(method, parameters));
+    }
+
+    /**
+     * Constructor for ConicSatelliteTracking.
+     * 
+     * Calculation are based on <cite>28 .SATTELITE-TRACKING PROJECTIONS </cite>

+     * in <cite> Map Projections - A Working Manual</cite> By John P. Snyder.
+     * 
+     * @param initializer 
+     */
+    ConicSatelliteTracking(final Initializer initializer) {
+        super(initializer);
+
+        //======================================================================
+        // Common for both cylindrical and conic sattelite tracking projections :
+        //======================================================================
+        i        = toRadians(initializer.getAndStore(SATELLITE_ORBIT_INCLINATION)); // Radian
input value.
+        cos_i    = cos(i);
+        sin_i    = sin(i);
+        cos2_i   = cos_i * cos_i;
+        p2_on_p1 = initializer.getAndStore(SATELLITE_ORBITAL_PERIOD) / initializer.getAndStore(ASCENDING_NODE_PERIOD);
+
+        final double φ1      = toRadians(initializer.getAndStore(STANDARD_PARALLEL_1));
 //appropriated use of toRadians??
+        cos_φ1               = cos(φ1);
+        final double cos2_φ1 = cos_φ1 * cos_φ1;
+        
+        //======================================================================
+        // For conic projection :
+        //=======================
+        if(!(this instanceof CylindricalSatelliteTracking)) {
+            final double sin_φ1 = sin(φ1);
+            final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN));  
//appropriated use of toRadians??
+            final double cos_φ0 = cos(φ0);
+            final double cos2_φ0 = cos_φ0 * cos_φ0;
+            final double sin_φ0 = sin(φ0);
+            final double φ2 = toRadians(initializer.getAndStore(STANDARD_PARALLEL_2)); 
 //appropriated use of toRadians??
+
+//        final double F0 = atan( (p2_on_p1*cos2_φ0 - cos_i)/sqrt(cos2_φ0 - cos2_i) );
+//        final double F1 = atan( (p2_on_p1*cos2_φ1 - cos_i)/sqrt(cos2_φ1 - cos2_i) );
+//        final double F2 = atan( (p2_on_p1*cos2_φ2 - cos_i)/sqrt(cos2_φ2 - cos2_i) );
+//        final double dλ0 = -asin(sin_φ0/sin_i);
+//        final double dλ1 = -asin(sin_φ1/sin_i);
+//        final double dλ2 = -asin(sin_φ2/sin_i);   
+//        final double λt0 = atan( tan(dλ0) * cos_i );
+//        final double λt1 = atan( tan(dλ1) * cos_i );
+//        final double λt2 = atan( tan(dλ2) * cos_i );
+            final DoubleUnaryOperator computeFn = (cos2_φn) -> atan((p2_on_p1 * cos2_φn
- cos_i) / sqrt(cos2_φn - cos2_i)); // eq.28-9 in Snyder 
+            final double F0 = computeFn.applyAsDouble(cos2_φ0);
+            /*
+         *Inclination of the groundtrack to the meridian at latitude φ1
+             */
+            final double F1 = computeFn.applyAsDouble(cos2_φ1);
+
+            cos_φ1xsin_F1 = cos_φ1 * sin(F1);
+
+            final DoubleUnaryOperator computedλn = (sin_φn) -> -asin(sin_φn / sin_i);
// eq.28-2a in Snyder 
+            final double dλ0 = computedλn.applyAsDouble(sin_φ0);
+            final double dλ1 = computedλn.applyAsDouble(sin_φ1);
+
+            final DoubleUnaryOperator computeλtn = (dλn) -> atan(tan(dλn) * cos_i);
// eq.28-3a in Snyder 
+            final double λt0 = computeλtn.applyAsDouble(dλ0);
+            final double λt1 = computeλtn.applyAsDouble(dλ1);
+
+            final double L0 = λt0 - p2_on_p1 * dλ0;
+            final double L1 = λt1 - p2_on_p1 * dλ1;
+
+            if (φ1 == PI / 2 - i) { //tracking limit
+                final double factor = (p2_on_p1 * cos_i - 1);
+                final double factor2 = factor * factor;
+                n = sin_i / (factor2); //eq. 28-18 in Snyder 
+            } else if (φ2 != φ1) {
+                final double cos_φ2 = cos(φ2);
+                final double cos2_φ2 = cos_φ2 * cos_φ2;
+                final double sin_φ2 = sin(φ2);
+                final double dλ2 = computedλn.applyAsDouble(sin_φ2);
+                final double λt2 = computeλtn.applyAsDouble(dλ2);
+                final double L2 = λt2 - p2_on_p1 * dλ2;
+                final double F2 = computeFn.applyAsDouble(cos2_φ2);
+                n = (F2 - F1) / (L2 - L1);
+            } 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 
+            }
+            s0 = F1 - n * L1;
+            ρ0 = cos_φ1xsin_F1 / (n * sin(n * L0 + s0)); // *R in eq.28-12 in Snyder 
+
+            //======================== Unsure ======================================
+            // Aim to assess the limit lattitude associated with -s0/n L-value.
+            //
+            lattitudeLimit = lattitudeFromNewtonMethod(-s0 / n);
+            positiveN = (n >= 0);
+            //======================================================================
+
+//            //Additionally we can compute the radius of the circle to which groundtracks
+//            //are tangent on the map :
+//            ρs = cos_φ1xsin_F1 / n; //*R
+
+            //======================================================================    
   
+            /*
+         * At this point, all parameters have been processed. Now process to their
+         * validation and the initialization of (de)normalize affine transforms.
+             */
+        final MatrixSIS normalize   = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
+//        final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
+        normalize  .convertAfter (0, n, null);  //For conic tracking
+//        denormalize.convertBefore(0, 1/PI,     null);
+//        denormalize.convertBefore(1, , null);
+        }else{
+            n = lattitudeLimit = cos_φ1xsin_F1 = s0 = ρ0 = NaN;
+            positiveN = false;
+        }
+    }
+
+    /**
+     * Converts the specified (λ,φ) coordinate (units in radians) and stores the result
in {@code dstPts}
+     * (linear distance on a unit sphere). In addition, opportunistically computes the projection
derivative
+     * if {@code derivate}
+     * is {@code true}.
+     *
+     * <cite> The Yaxis lies along the central meridian λ0, y increasing
+     * northerly, and X axis intersects perpendicularly at LATITUDE_OF_ORIGIN
+     * φ0, x increasing easterly.
+     * </cite>
+     *
+     * @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 coordinates can not be converted.
+     */
+    @Override
+    public Matrix transform(double[] srcPts, int srcOff, 
+                            double[] dstPts, int dstOff,
+                            boolean derivate) throws ProjectionException {
+        
+        final double λ      = srcPts[srcOff];
+        final double φ      = srcPts[srcOff + 1];
+        
+        /*
+         * According to the Snyder (page 236) in those cases cannot or should not be plotted.
+         */
+        if( ((positiveN)&&(φ<=lattitudeLimit)) || ((!positiveN)&&(φ>=lattitudeLimit))
){
+            throw new ProjectionException(Resources.format(Resources.Keys.CanNotTransformCoordinates_2));
+        }
+        
+        final double sin_φ  = sin(φ);
+        final double dλ     = -asin(sin_φ / sin_i);
+        final double λt     = atan(tan(dλ) * cos_i);
+        final double L      = λt - p2_on_p1 * dλ;
+        
+        /*
+         * As {@code lattitudeLimit} is an approximation we repeat the test here.
+         */
+        if( ((positiveN)&&(L<=-s0/n)) || ((!positiveN)&&(L>=-s0/n))
){
+            //TODO if usefull, could we add :
+            // lattitudeLimit = φ;
+            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 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;
+        }
+        
+         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φ = ?;
+        //======================================================================
+        
+        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);
+    }
+
+    /**
+     * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
+     * and stores the result in {@code dstPts} (angles in radians).
+     *
+     * @throws ProjectionException if the coordinates can not be converted.
+     */
+    @Override
+    protected void inverseTransform(double[] srcPts, int srcOff,
+                                    double[] dstPts, int dstOff)
+                                    throws ProjectionException {
+        
+        final double x   = srcPts[srcOff];
+        final double y   = srcPts[srcOff + 1] - ρ0;
+        
+        //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.
+        }
+               
+        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
+        
+        //TODO ensure that λ0 will be added. ;  In eq. Snyder 28-23 : λ0 + θ/n
+        dstPts[dstOff  ] =  θ ;//λ 
+        dstPts[dstOff+1] = lattitudeFromNewtonMethod(L); //φ
+    }
+    
+    protected final double lattitudeFromNewtonMethod(final double l){
+        double dλ = -PI/2;
+        double dλn = Double.MIN_VALUE;
+        double A, A2, Δdλ;
+        
+        int count=0;
+        int maxIteration = 100; //TODO : check the value.
+                
+        while((count==0) || dλ-dλn >= (0.1*PI/180)){ //TODO : check the condition. It
is considered here that convergence is reached when improvement is lower than 0.1°.
+            if(count >= maxIteration){
+                throw new RuntimeException(Resources.format(Resources.Keys.NoConvergence));
+            }
+            dλn = dλ;
+
+//            λt = l + p2_on_p1 * dλ_n; 
+//            dλ = atan(tan(λt) / 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); 
+    }
+    
+    
+    public double getLattitudeLimit(){
+        return lattitudeLimit;
+    }
+    
+//    /**
+//     * 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
new file mode 100644
index 0000000..70ea6fe
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTracking.java
@@ -0,0 +1,197 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import org.opengis.referencing.operation.Matrix;
+
+import org.apache.sis.parameter.Parameters;
+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.MatrixSIS;
+import org.apache.sis.referencing.operation.transform.ContextualParameters;
+
+/**
+ * <cite>Cylindrical Satellite-Tracking projections</cite>. 
+ * 
+ * <cite>
+ * - All groundtracks
+ * for satellites orbiting the Earth with the same orbital parameters are shown
+ * as straight lines on the map.
+ *
+ * - Cylindrical {@link CylindricalSatelliteTracking}
+ * or conical {@link ConicSatelliteTracking} form available. 
+ *
+ * - Neither conformal nor equal-area.
+ *
+ * - All meridians are equally spaced straight lines, parallel on cylindrical
+ * form and converging to a common point on conical form.
+ *
+ * - All parallels are straight and parallel on cylindrical form and are
+ * concentric circular arcs on conical form. Parallels are unequally spaced.
+ *
+ * - Conformality occurs along two chosen parallels. Scale is correct along one
+ * of these parameters on the conical form and along both on the cylindrical
+ * form.
+ *
+ * Developed 1977 by Snyder
+ * </cite>
+ *
+ * <cite> These formulas are confined to circular orbits and the SPHERICAL
+ * Earth.</cite>
+ *
+ * <cite>The ascending and descending groundtracks meet at the northern an
+ * southern tracking limits, lats. 80.9°N and S for landsat 1, 2 and 3. The map
+ * Projection does not extend closer to the poles.</cite>
+ *
+ * This projection method has no associated EPSG code.
+ *
+ * Earth radius is normalized. Its value is 1 and is'nt an input parameter.
+ *
+ * =============================================================================
+ * REMARK : The parameters associated with the satellite (and its orbit) could
+ * be aggregate in class of the kind : Satellite or SatelliteOrbit.
+ * =============================================================================
+ *
+ * @see <cite>Map Projections - A Working Manual</cite> By John P. Snyder
+ * @author Matthieu Bastianelli (Geomatys)
+ * @version 1.0
+ */
+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.
+     *
+     * The parameters are described in <cite>Map Projections - A Working
+     * Manual</cite> By John P. Snyder.
+     *
+     * @param method : description of the projection parameters.
+     * @param parameters : the parameter values of the projection to create.
+     */
+    public CylindricalSatelliteTracking(final OperationMethod method, final Parameters parameters)
{
+        this(initializer(method, parameters));
+    }
+
+    private CylindricalSatelliteTracking(final Initializer initializer) {
+        super(initializer);
+        
+        final double cos2_φ1 = cos_φ1 * cos_φ1;
+        dF1 = (p2_on_p1 * cos2_φ1 - cos_i) / sqrt(cos2_φ1 - cos2_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);
+
+    }
+
+    /**
+     * Converts the specified (λ,φ) coordinate (units in radians) and stores the result
in {@code dstPts}
+     * (linear distance on a unit sphere). In addition, opportunistically computes the projection
derivative
+     * if {@code derivate} is {@code true}.
+     * 
+     * <cite> The Yaxis lies along the central meridian λ0, y increasing northerly,
+     * and X axis intersects perpendicularly at O_PARALLEL φ0, x increasing easterly.
+     * </cite> 
+     *
+     * @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 coordinates can not be converted.
+     */
+    @Override
+    public Matrix transform(double[] srcPts, int srcOff,
+            double[] dstPts, int dstOff,
+            boolean derivate) throws ProjectionException {
+        
+        final double λ = srcPts[srcOff];
+        final double φ = srcPts[srcOff + 1];
+
+        // TODO : check the condition and the thrown exception.
+        if (abs(φ) > PI / 2 - abs(PI / 2 - i)) {  // Exceed tracking limit 
+            throw new ProjectionException(Resources.format(Resources.Keys.CanNotTransformCoordinates_2));
+        }
+
+//        final double cos_φ  = cos(φ);
+//        final double cos2_φ = cos_φ * cos_φ;
+        final double sin_φ  = sin(φ);
+
+        final double dλ     = -asin(sin_φ / sin_i);
+        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
+        }
+
+        /* =====================================================================
+        * Uncomputed scale factors :
+        *===========================
+        * F' : tangente of the angle on the globe between the groundtrack and 
+        * 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
+        ===================================================================== */
+        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φ =
+                
+        throw new UnsupportedOperationException("Not supported yet."); //To change body of
generated methods, choose Tools | Templates.
+    }
+
+    /**
+     * Transforms the specified (<var>x</var>,<var>y</var>) coordinates
+     * and stores the result in {@code dstPts} (angles in radians).
+     *
+     * @throws ProjectionException if the coordinates can not be converted.
+     */
+    @Override
+    protected void inverseTransform(double[] srcPts, int srcOff, 
+                                    double[] dstPts, int dstOff) 
+            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
+        
+        dstPts[dstOff  ] =  x;
+        dstPts[dstOff+1] = lattitudeFromNewtonMethod(L); //φ
+    }
+
+}
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
new file mode 100644
index 0000000..28d83df
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConicSatelliteTrackingTest.java
@@ -0,0 +1,127 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import static java.lang.Double.NaN;
+import java.util.Collections;
+import org.apache.sis.internal.referencing.NilReferencingObject;
+import org.apache.sis.internal.referencing.provider.SatelliteTracking;
+import org.apache.sis.measure.Units;
+import org.apache.sis.referencing.datum.DefaultEllipsoid;
+import org.apache.sis.referencing.operation.transform.MathTransformFactoryMock;
+import org.apache.sis.test.DependsOn;
+import static org.junit.Assert.assertTrue;
+import org.junit.Test;
+import org.opengis.parameter.ParameterValueGroup;
+import org.opengis.referencing.operation.TransformException;
+import org.opengis.util.FactoryException;
+
+/**
+ * Tests coordiantes computed by applying a conic satellite-tracking projection
+ * with {@link ConicSatelliteTracking}.
+ *
+ * @author Matthieu Bastianelli (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+@DependsOn(ConformalProjectionTest.class)
+public class ConicSatelliteTrackingTest extends MapProjectionTestCase {
+
+    /**
+     * Creates a new instance of {@link ConicSatelliteTracking} concatenated
+     * with the (de)normalization matrices. The new instance is stored in the
+     * inherited {@link #transform} field.
+     *
+     * @param i : Angle of inclination between the plane of the Earth's Equator
+     * and the plane of the satellite orbit.
+     * @param orbitalT : Time required for revolution of the satellite.
+     * @param ascendingNodeT : Length of Earth's rotation with respect to the
+     * precessed ascending node.
+     * @param λ0 : central meridian.
+     * @param φ1 : first parallel of conformality, with true scale.
+     * @param φ2 : second parallel of conformality, without true scale.
+     * @param φ0 : lattitude_of_origin : latitude Crossing the central meridian
+     * at the desired origin of rectangular coordinates (null or NaN for
+     * cylindrical satellite tracking projection.)
+     * @return
+     */
+    void createProjection(final double i,
+            final double orbitalT, final double ascendingNodeT,
+            final double λ0,       final double φ1, 
+            final double φ2,       final double φ0) 
+            throws FactoryException {
+
+        final SatelliteTracking provider = new SatelliteTracking();
+        final ParameterValueGroup values = provider.getParameters().createValue();
+        final DefaultEllipsoid sphere = DefaultEllipsoid.createEllipsoid(
+                Collections.singletonMap(DefaultEllipsoid.NAME_KEY, NilReferencingObject.UNNAMED),
+                1, 1, Units.METRE);
+        
+        values.parameter("semi_major").setValue(sphere.getSemiMajorAxis());
+        values.parameter("semi_minor").setValue(sphere.getSemiMinorAxis());
+        values.parameter("satellite_orbit_inclination").setValue(i);
+        values.parameter("satellite_orbital_period").setValue(orbitalT);
+        values.parameter("ascending_node_period").setValue(ascendingNodeT);
+        values.parameter("central_meridian").setValue(λ0);
+        values.parameter("standard_parallel_1").setValue(φ1);
+        
+        if (!Double.isNaN(φ2)) {
+            values.parameter("standard_parallel_2").setValue(φ2);
+        }else{
+            values.parameter("standard_parallel_2").setValue(-φ1); //Cylindrical case
+        }
+        if (!Double.isNaN(φ0)) {
+            values.parameter("latitude_of_origin").setValue(φ0);
+        }
+        
+        transform = new MathTransformFactoryMock(provider).createParameterizedTransform(values);
+        validate();
+    }
+
+    /**
+     * Tests the projection of a few points on a sphere.
+     *
+     * @throws FactoryException if an error occurred while creating the map
+     * projection.
+     * @throws TransformException if an error occurred while projecting a point.
+     */
+    @Test
+    public void testTransform() throws FactoryException, TransformException {
+        tolerance = 1E-5;
+        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       //lattitude_of_origin
+        );
+        assertTrue(isInverseTransformSupported);
+        verifyTransform(
+                new double[]{ // (λ,φ) coordinates in degrees to project.
+                    -75, 40
+                },
+                new double[]{ // Expected (x,y) results in metres.
+                    0.2001910, 0.2121685
+                });
+    }
+
+}
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
new file mode 100644
index 0000000..0a86ba6
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalSatelliteTrackingTest.java
@@ -0,0 +1,117 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements.  See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership.  The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied.  See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import static java.lang.Double.NaN;
+import org.apache.sis.internal.referencing.Formulas;
+import static org.junit.Assert.assertTrue;
+import org.junit.Test;
+import org.opengis.referencing.operation.TransformException;
+import org.opengis.util.FactoryException;
+
+/**
+ * Tests coordiantes computed by applying a cylindrical satellite-tracking projection
+ * with {@link CylindricalSatelliteTracking}.
+ *
+ * @author Matthieu Bastianelli (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+public class CylindricalSatelliteTrackingTest extends ConicSatelliteTrackingTest {
+    
+    
+    
+    /**
+     * Creates a new instance of {@link ConicSatelliteTracking} concatenated
+     * with the (de)normalization matrices. The new instance is stored in the
+     * inherited {@link #transform} field.
+     *
+     * @param i : Angle of inclination between the plane of the Earth's Equator
+     * and the plane of the satellite orbit.
+     * @param orbitalT : Time required for revolution of the satellite.
+     * @param ascendingNodeT : Length of Earth's rotation with respect to the
+     * precessed ascending node.
+     * @param λ0 : central meridian.
+     * @param φ1 : first parallel of conformality, with true scale.
+     * @param φ2 : second parallel of conformality, without true scale.
+     * @param φ0 : lattitude_of_origin : latitude Crossing the central meridian
+     * at the desired origin of rectangular coordinates (null or NaN for
+     * cylindrical satellite tracking projection.)
+     * @return
+     */
+    private void createProjection(final double i,
+            final double orbitalT, final double ascendingNodeT,
+            final double λ0,       final double φ1) 
+            throws FactoryException {
+        super.createProjection(i, orbitalT, ascendingNodeT, λ0, φ1, NaN, NaN);
+    }
+    
+    /**
+     * Tests the projection of a few points on a sphere.
+     *
+     * @throws FactoryException if an error occurred while creating the map
+     * projection.
+     * @throws TransformException if an error occurred while projecting a point.
+     */
+    @Test
+    @Override
+    public void testTransform() throws FactoryException, TransformException {
+        tolerance = 1E-5;
+        createProjection(
+                99.092,  //satellite_orbit_inclination
+                103.267, //satellite_orbital_period
+                1440.0,  //ascending_node_period
+                -90,     //central_meridian
+                30       //standard_parallel_1
+        );
+        assertTrue(isInverseTransformSupported);
+        verifyTransform(
+                new double[]{ // (λ,φ) coordinates in degrees to project.
+                    -75, 40
+                },
+                new double[]{ // Expected (x,y) results in metres.
+                    0.2267249, 0.6459071
+                });
+        
+        createProjection(
+                99.092,  //satellite_orbit_inclination
+                103.267, //satellite_orbital_period
+                1440.0,  //ascending_node_period
+                0,       //central_meridian
+                0        //standard_parallel_1
+        );
+        
+        tolerance = Formulas.LINEAR_TOLERANCE; //Don't pass with the former tolerance.
+        
+        final double xConverterFactor=0.017453;
+        verifyTransform(
+                new double[]{ // (λ,φ) coordinates in degrees to project.
+                      0,  0,
+                     10,  0, 
+                    -10, 10
+                },
+                new double[]{ // Expected (x,y) results in metres.
+                    0, 0,
+                    xConverterFactor *  10,  0,
+                    xConverterFactor * -10,  0.17579,
+                });
+    }
+    
+}


Mime
View raw message