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