sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/03: Merge the implementation of Cassini-Soldner projection from "geoapi-4.0" development branch. JIRA: https://issues.apache.org/jira/browse/SIS-218 IP-review: https://svn.apache.org/repos/asf/sis/ip-review/CassiniSoldner.html
Date Wed, 29 Apr 2020 16:07:35 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sis.git

commit f2e8a22c7c7198d8195d7dc5a9afdc11bedf0258
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Apr 29 18:05:28 2020 +0200

    Merge the implementation of Cassini-Soldner projection from "geoapi-4.0" development branch.
    JIRA: https://issues.apache.org/jira/browse/SIS-218
    IP-review: https://svn.apache.org/repos/asf/sis/ip-review/CassiniSoldner.html
---
 .../referencing/provider/CassiniSoldner.java       | 169 ++++++++++++
 .../operation/projection/CassiniSoldner.java       | 287 +++++++++++++++++++++
 ...g.opengis.referencing.operation.OperationMethod |   1 +
 .../referencing/provider/ProvidersTest.java        |   1 +
 .../operation/projection/CassiniSoldnerTest.java   | 149 +++++++++++
 .../transform/MathTransformFactoryMock.java        |   2 +-
 .../sis/test/suite/ReferencingTestSuite.java       |   1 +
 7 files changed, 609 insertions(+), 1 deletion(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/CassiniSoldner.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/CassiniSoldner.java
new file mode 100644
index 0000000..e40850f
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/CassiniSoldner.java
@@ -0,0 +1,169 @@
+/*
+ * 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 org.opengis.parameter.ParameterDescriptor;
+import org.opengis.parameter.ParameterDescriptorGroup;
+import org.apache.sis.metadata.iso.citation.Citations;
+import org.apache.sis.referencing.operation.projection.NormalizedProjection;
+import org.apache.sis.parameter.Parameters;
+
+
+/**
+ * The provider for "<cite>Cassini-Soldner</cite>" projection (EPSG:9806).
+ * This projection is similar to {@link TransverseMercator}.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ *
+ * @see <a href="http://geotiff.maptools.org/proj_list/cassini_soldner.html">GeoTIFF
parameters for Cassini-Soldner</a>
+ *
+ * @since 1.1
+ * @module
+ */
+@XmlTransient
+public final class CassiniSoldner extends MapProjection {
+    /**
+     * For cross-version compatibility.
+     */
+    private static final long serialVersionUID = 7280273456465057368L;
+
+    /**
+     * The operation parameter descriptor for the <cite>Latitude of natural origin</cite>
(φ₀) parameter value.
+     * Valid values range is [-90 … 90]° and default value is 0°.
+     *
+     * <!-- Generated by ParameterNameTableGenerator -->
+     * <table class="sis">
+     *   <caption>Parameter names</caption>
+     *   <tr><td> EPSG:    </td><td> Latitude of natural origin </td></tr>
+     *   <tr><td> OGC:     </td><td> latitude_of_origin </td></tr>
+     *   <tr><td> ESRI:    </td><td> Latitude_Of_Origin </td></tr>
+     *   <tr><td> NetCDF:  </td><td> latitude_of_projection_origin
</td></tr>
+     *   <tr><td> GeoTIFF: </td><td> NatOriginLat </td></tr>
+     *   <tr><td> Proj4:   </td><td> lat_0 </td></tr>
+     * </table>
+     */
+    public static final ParameterDescriptor<Double> LATITUDE_OF_ORIGIN = TransverseMercator.LATITUDE_OF_ORIGIN;
+
+    /**
+     * The operation parameter descriptor for the <cite>Longitude of natural origin</cite>
(λ₀) parameter value.
+     * Valid values range is [-180 … 180]° and default value is 0°.
+     *
+     * <!-- Generated by ParameterNameTableGenerator -->
+     * <table class="sis">
+     *   <caption>Parameter names</caption>
+     *   <tr><td> EPSG:    </td><td> Longitude of natural origin
</td></tr>
+     *   <tr><td> OGC:     </td><td> central_meridian </td></tr>
+     *   <tr><td> ESRI:    </td><td> Central_Meridian </td></tr>
+     *   <tr><td> NetCDF:  </td><td> longitude_of_central_meridian
</td></tr>
+     *   <tr><td> GeoTIFF: </td><td> NatOriginLong </td></tr>
+     *   <tr><td> Proj4:   </td><td> lon_0 </td></tr>
+     * </table>
+     */
+    public static final ParameterDescriptor<Double> LONGITUDE_OF_ORIGIN = TransverseMercator.LONGITUDE_OF_ORIGIN;
+
+    /**
+     * The operation parameter descriptor for the <cite>Scale factor at natural origin</cite>
(k₀) parameter value.
+     * Valid values range is (0 … ∞) and default value is 1. This is not formally a parameter
of this projection.
+     *
+     * <!-- Generated by ParameterNameTableGenerator -->
+     * <table class="sis">
+     *   <caption>Parameter names</caption>
+     *   <tr><td> OGC:     </td><td> scale_factor </td></tr>
+     *   <tr><td> ESRI:    </td><td> Scale_Factor </td></tr>
+     *   <tr><td> Proj4:   </td><td> k </td></tr>
+     * </table>
+     * <b>Notes:</b>
+     * <ul>
+     *   <li>Optional</li>
+     * </ul>
+     */
+    public static final ParameterDescriptor<Double> SCALE_FACTOR = Mercator2SP.SCALE_FACTOR;
+
+    /**
+     * The operation parameter descriptor for the <cite>False easting</cite>
(FE) parameter value.
+     * Valid values range is unrestricted and default value is 0 metre.
+     *
+     * <!-- Generated by ParameterNameTableGenerator -->
+     * <table class="sis">
+     *   <caption>Parameter names</caption>
+     *   <tr><td> EPSG:    </td><td> False easting </td></tr>
+     *   <tr><td> OGC:     </td><td> false_easting </td></tr>
+     *   <tr><td> ESRI:    </td><td> False_Easting </td></tr>
+     *   <tr><td> NetCDF:  </td><td> false_easting </td></tr>
+     *   <tr><td> GeoTIFF: </td><td> FalseEasting </td></tr>
+     *   <tr><td> Proj4:   </td><td> x_0 </td></tr>
+     * </table>
+     */
+    public static final ParameterDescriptor<Double> FALSE_EASTING = TransverseMercator.FALSE_EASTING;
+
+    /**
+     * The operation parameter descriptor for the <cite>False northing</cite>
(FN) parameter value.
+     * Valid values range is unrestricted and default value is 0 metre.
+     *
+     * <!-- Generated by ParameterNameTableGenerator -->
+     * <table class="sis">
+     *   <caption>Parameter names</caption>
+     *   <tr><td> EPSG:    </td><td> False northing </td></tr>
+     *   <tr><td> OGC:     </td><td> false_northing </td></tr>
+     *   <tr><td> ESRI:    </td><td> False_Northing </td></tr>
+     *   <tr><td> NetCDF:  </td><td> false_northing </td></tr>
+     *   <tr><td> GeoTIFF: </td><td> FalseNorthing </td></tr>
+     *   <tr><td> Proj4:   </td><td> y_0 </td></tr>
+     * </table>
+     */
+    public static final ParameterDescriptor<Double> FALSE_NORTHING = TransverseMercator.FALSE_NORTHING;
+
+    /**
+     * The group of all parameters expected by this coordinate operation.
+     */
+    static final ParameterDescriptorGroup PARAMETERS;
+    static {
+        PARAMETERS = builder()
+                .addIdentifier(              "9806")
+                .addName(                    "Cassini-Soldner")
+                .addName(Citations.OGC,      "Cassini_Soldner")
+                .addName(Citations.ESRI,     "Cassini_Soldner")
+                .addName(Citations.GEOTIFF,  "CT_CassiniSoldner")
+                .addName(Citations.PROJ4,    "cass")
+                .addIdentifier(Citations.GEOTIFF,  "18")
+                .createGroupForMapProjection(
+                        LATITUDE_OF_ORIGIN,
+                        LONGITUDE_OF_ORIGIN,
+                        SCALE_FACTOR,
+                        FALSE_EASTING,
+                        FALSE_NORTHING);
+    }
+
+    /**
+     * Constructs a new provider.
+     */
+    public CassiniSoldner() {
+        super(PARAMETERS);
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * @return the map projection created from the given parameter values.
+     */
+    @Override
+    protected NormalizedProjection createProjection(final Parameters parameters) {
+        return new org.apache.sis.referencing.operation.projection.CassiniSoldner(this, parameters);
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java
new file mode 100644
index 0000000..e857a29
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java
@@ -0,0 +1,287 @@
+/*
+ * 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 java.util.EnumMap;
+import org.opengis.util.FactoryException;
+import org.opengis.parameter.ParameterDescriptor;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.MathTransformFactory;
+import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.OperationMethod;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
+import org.apache.sis.referencing.operation.transform.ContextualParameters;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.util.Workaround;
+
+import static java.lang.Math.*;
+import static org.apache.sis.internal.referencing.provider.CassiniSoldner.*;
+
+
+/**
+ * <cite>Cassini-Soldner</cite> projection (EPSG code 9806).
+ * See the following references for an overview:
+ * <ul>
+ *   <li><a href="https://en.wikipedia.org/wiki/Cassini_projection">Cassini projection
on Wikipedia</a></li>
+ *   <li><a href="https://mathworld.wolfram.com/CassiniProjection.html">Cassini
projection on MathWorld</a></li>
+ * </ul>
+ *
+ * The ellipsoidal form of this projection is suitable only within a few degrees (3° or
4° of longitude)
+ * to either side of the central meridian.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @author  Rémi Maréchal (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public class CassiniSoldner extends MeridianArcBased {
+    /**
+     * For cross-version compatibility.
+     */
+    private static final long serialVersionUID = 3007155786839466950L;
+
+    /**
+     * Creates a Cassini-Soldner projection from the given parameters.
+     * The {@code method} argument can be the description of one of the following:
+     *
+     * <ul>
+     *   <li><cite>"Cassini-Soldner"</cite>.</li>
+     * </ul>
+     *
+     * @param  method      description of the projection parameters.
+     * @param  parameters  the parameter values of the projection to create.
+     */
+    public CassiniSoldner(final OperationMethod method, final Parameters parameters) {
+        this(initializer(method, parameters));
+    }
+
+    /**
+     * 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.7")
+    private static Initializer initializer(final OperationMethod method, final Parameters
parameters) {
+        final EnumMap<ParameterRole, ParameterDescriptor<Double>> roles = new
EnumMap<>(ParameterRole.class);
+        roles.put(ParameterRole.CENTRAL_MERIDIAN, LONGITUDE_OF_ORIGIN);
+        roles.put(ParameterRole.SCALE_FACTOR,     SCALE_FACTOR);
+        roles.put(ParameterRole.FALSE_EASTING,    FALSE_EASTING);
+        roles.put(ParameterRole.FALSE_NORTHING,   FALSE_NORTHING);
+        return new Initializer(method, parameters, roles, (byte) 0);
+    }
+
+    /**
+     * Creates a new Cassini-Soldner projection from the given initializer.
+     */
+    CassiniSoldner(final Initializer initializer) {
+        super(initializer);
+        final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN));
+        final MatrixSIS denormalize = getContextualParameters().getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
+        denormalize.convertBefore(1, null, -distance(φ0, sin(φ0), cos(φ0)));
+    }
+
+    /**
+     * Creates a new projection initialized to the same parameters than the given one.
+     */
+    CassiniSoldner(final CassiniSoldner other) {
+        super(other);
+    }
+
+    /**
+     * Returns the sequence of <cite>normalization</cite> → {@code this} →
<cite>denormalization</cite> transforms
+     * as a whole. The transform returned by this method expects (<var>longitude</var>,
<var>latitude</var>)
+     * coordinates in <em>degrees</em> and returns (<var>x</var>,<var>y</var>)
coordinates in <em>metres</em>.
+     * The non-linear part of the returned transform will be {@code this} transform, except
if the ellipsoid
+     * is spherical. In the later case, {@code this} transform will be replaced by a simplified
implementation.
+     *
+     * @param  factory  the factory to use for creating the transform.
+     * @return the map projection from (λ,φ) to (<var>x</var>,<var>y</var>)
coordinates.
+     * @throws FactoryException if an error occurred while creating a transform.
+     */
+    @Override
+    public MathTransform createMapProjection(final MathTransformFactory factory) throws FactoryException
{
+        CassiniSoldner kernel = this;
+        if (eccentricity == 0) {
+            kernel = new Spherical(this);
+        }
+        return context.completeTransform(factory, kernel);
+    }
+
+    /**
+     * Converts the specified (λ,φ) coordinate (units in radians) and stores the result
in {@code dstPts}.
+     * In addition, opportunistically computes the projection derivative if {@code derivate}
is {@code true}.
+     *
+     * @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 coordinate can not be converted.
+     */
+    @Override
+    public Matrix transform(final double[] srcPts, final int srcOff,
+                            final double[] dstPts, final int dstOff,
+                            final boolean derivate) throws ProjectionException
+    {
+        /*
+         * Formula from IOGP Publication 373-7-2 — Geomatics Guidance Note
+         * — Coordinate Conversions and Transformations including Formulas.
+         *
+         * Following symbols are from EPSG formulas with a=1, λ₀=0 and M₀=0:
+         * λ, φ, T, A, C and rν=1/ν.
+         *
+         * Following symbol are specific to this implementation: Q and S.
+         * Those symbols are defined because of terms reused in derivatives.
+         */
+        final double λ     = srcPts[srcOff  ];
+        final double φ     = srcPts[srcOff+1];
+        final double cosφ  = cos(φ);
+        final double sinφ  = sin(φ);
+        final double sinφ2 = sinφ * sinφ;
+        final double cosφ2 = cosφ * cosφ;
+        final double tanφ  = sinφ / cosφ;
+        final double T     = tanφ * tanφ;
+        final double A     = λ  * cosφ;
+        final double A2    = A  * A;
+        final double A3    = A2 * A;
+        final double rν2   = 1 - sinφ2*eccentricitySquared;
+        final double C     = cosφ2*eccentricitySquared / (1 - eccentricitySquared);
+        final double Q     = T - 8*(C + 1);
+        final double S     = ((5-T)/6 + C)*A2;
+        final double rν    = sqrt(rν2);
+        if (dstPts != null) {
+            dstPts[dstOff  ] = (A - T*A3/6 + Q*T*(A3*A2) / 120) / rν;
+            dstPts[dstOff+1] = distance(φ, sinφ, cosφ) + tanφ*A2*(0.5 + S/4) / rν;
+        }
+        if (!derivate) {
+            return null;
+        }
+        /*
+         * Following formulas have been derived with WxMaxima, then simplified by hand.
+         * See: https://svn.apache.org/repos/asf/sis/analysis/README.html
+         */
+        final double λ2 = λ * λ;
+        final double B  = λ * sinφ;                             // Like A but with sin(φ)
instead of cos(φ).
+        final double B2 = B*B;                                  // = T*A2
+        final double D  = cosφ2*eccentricitySquared / rν2;      // Like C but with a sin²φ
in denominator.
+        final double V  = A2*Q/60 - 1./3;
+        final double W  = B2*(Q*A2/24 - 0.5);
+        return new Matrix2(
+                  (cosφ/rν) * (W + 1),
+                        (B) * (λ2*V - W + B2*(λ2 + 8*C*A2)/60 +(D*(B2*V/2 + 1) - 1)/rν),
+                (B*cosφ/rν) * (S + 1),
+                    (B2/rν) * ((S/4 + 0.5)*(D + 1/sinφ2) - (S+1 + A2*C/2 + λ2/12)) + dM_dφ(sinφ2));
+    }
+
+    /**
+     * Converts the specified (<var>x</var>,<var>y</var>) coordinates
+     * and stores the result in {@code dstPts} (angles in radians).
+     *
+     * @throws ProjectionException if the point can not be converted.
+     */
+    @Override
+    protected void inverseTransform(final double[] srcPts, final int srcOff,
+                                    final double[] dstPts, final int dstOff)
+            throws ProjectionException
+    {
+        final double x     = srcPts[srcOff  ];
+        final double y     = srcPts[srcOff+1];
+        final double φ1    = latitude(y);
+        final double sinφ1 = sin(φ1);
+        final double cosφ1 = cos(φ1);
+        final double tanφ1 = sinφ1 / cosφ1;
+        final double rν2   =  1 - eccentricitySquared*(sinφ1*sinφ1);    // ν₁⁻²
+        final double ρ1_ν1 = (1 - eccentricitySquared) / rν2;           // ρ₁/ν₁
+        final double D     = x  * sqrt(rν2);                            // x/ν₁
+        final double D2    = D  * D;
+        final double D4    = D2 * D2;
+        final double T1    = tanφ1 * tanφ1;
+        dstPts[dstOff  ]   = (D - T1*(D2*D)/3 + (1 + 3*T1)*T1*(D4*D)/15) / cosφ1;
+        dstPts[dstOff+1]   = φ1 - (tanφ1 / ρ1_ν1) * (D2/2 - (1 + 3*T1)*D4 / 24);
+    }
+
+
+    /**
+     * Provides the transform equations for the spherical case of the Cassini-Soldner projection.
+     * Formulas are available on <a href="https://en.wikipedia.org/wiki/Cassini_projection">Wikipedia</a>.
+     *
+     * @author  Martin Desruisseaux (Geomatys)
+     * @author  Rémi Maréchal (Geomatys)
+     * @version 1.1
+     * @since   1.1
+     * @module
+     */
+    static final class Spherical extends CassiniSoldner {
+        /**
+         * For cross-version compatibility.
+         */
+        private static final long serialVersionUID = -8131887916538379858L;
+
+        /**
+         * Constructs a new map projection from the parameters of the given projection.
+         *
+         * @param  other  the other projection (usually ellipsoidal) from which to copy the
parameters.
+         */
+        protected Spherical(final CassiniSoldner other) {
+            super(other);
+        }
+
+        /**
+         * {@inheritDoc}
+         */
+        @Override
+        public Matrix transform(final double[] srcPts, final int srcOff,
+                                final double[] dstPts, final int dstOff,
+                                final boolean derivate)
+        {
+            final double λ    = srcPts[srcOff  ];
+            final double φ    = srcPts[srcOff+1];
+            final double sinλ = sin(λ);
+            final double cosλ = cos(λ);
+            final double sinφ = sin(φ);
+            final double cosφ = cos(φ);
+            final double tanφ = sinφ / cosφ;
+            if (dstPts != null) {
+                dstPts[dstOff  ] = asin (cosφ * sinλ);
+                dstPts[dstOff+1] = atan2(tanφ,  cosλ);
+            }
+            if (!derivate) {
+                return null;
+            }
+            final double cosφ2 = cosφ * cosφ;
+            final double sinλ2 = sinλ * sinλ;
+            final double dxden = sqrt(1 - sinλ2*cosφ2);
+            final double dyden = tanφ*tanφ + cosλ*cosλ;
+            return new Matrix2(
+                    cosλ*cosφ  / dxden,
+                   -sinλ*sinφ  / dxden,
+                    sinλ*tanφ  / dyden,
+                    cosλ/cosφ2 / dyden);
+        }
+
+        /**
+         * {@inheritDoc}
+         */
+        @Override
+        protected void inverseTransform(final double[] srcPts, final int srcOff,
+                                        final double[] dstPts, final int dstOff)
+        {
+            final double x = srcPts[srcOff  ];
+            final double y = srcPts[srcOff+1];
+            dstPts[dstOff  ] = atan2(tan(x),  cos(y));
+            dstPts[dstOff+1] = asin (sin(y) * cos(x));
+        }
+    }
+}
diff --git a/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
b/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
index 11937c1..49cb955 100644
--- a/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
+++ b/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
@@ -44,6 +44,7 @@ org.apache.sis.internal.referencing.provider.LambertCylindricalEqualAreaSpherica
 org.apache.sis.internal.referencing.provider.AlbersEqualArea
 org.apache.sis.internal.referencing.provider.TransverseMercator
 org.apache.sis.internal.referencing.provider.TransverseMercatorSouth
+org.apache.sis.internal.referencing.provider.CassiniSoldner
 org.apache.sis.internal.referencing.provider.PolarStereographicA
 org.apache.sis.internal.referencing.provider.PolarStereographicB
 org.apache.sis.internal.referencing.provider.PolarStereographicC
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
index fd483a4..824c123 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
@@ -94,6 +94,7 @@ public final strictfp class ProvidersTest extends TestCase {
             AlbersEqualArea.class,
             TransverseMercator.class,
             TransverseMercatorSouth.class,
+            CassiniSoldner.class,
             PolarStereographicA.class,
             PolarStereographicB.class,
             PolarStereographicC.class,
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CassiniSoldnerTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CassiniSoldnerTest.java
new file mode 100644
index 0000000..d8992b2
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CassiniSoldnerTest.java
@@ -0,0 +1,149 @@
+/*
+ * 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.util.FactoryException;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.referencing.provider.MapProjection;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.geometry.DirectPosition2D;
+import org.apache.sis.parameter.Parameters;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+import static java.lang.Double.NaN;
+import static java.lang.StrictMath.*;
+
+
+/**
+ * Tests the {@link CassiniSoldner} class.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @author  Rémi Maréchal (Geomatys)
+ * @version 1.1
+ * @since   0.8
+ * @module
+ */
+public final strictfp class CassiniSoldnerTest extends MapProjectionTestCase {
+    /**
+     * Returns the operation method for the projection tested in this class.
+     */
+    private static MapProjection method() {
+        return new org.apache.sis.internal.referencing.provider.CassiniSoldner();
+    }
+
+    /**
+     * Creates a new instance of {@link CassiniSoldner}. This instance expects
+     * angles in radians and a φ₀ value of zero; this is not a complete transform.
+     *
+     * @param  ellipse  {@code false} for a sphere, or {@code true} for WGS84 ellipsoid.
+     * @return newly created projection.
+     */
+    private static CassiniSoldner create(final boolean ellipse) {
+        final MapProjection provider = method();
+        final Parameters pg = Parameters.castOrWrap(provider.getParameters().createValue());
+        if (ellipse) {
+            pg.parameter("semi-major").setValue(WGS84_A);
+            pg.parameter("semi-minor").setValue(WGS84_B);
+            return new CassiniSoldner(provider, pg);
+        } else {
+            pg.parameter("semi-major").setValue(RADIUS);
+            pg.parameter("semi-minor").setValue(RADIUS);
+            return new CassiniSoldner.Spherical(new CassiniSoldner(provider, pg));
+        }
+    }
+
+    /**
+     * Tests some identities related to the {@link CassiniSoldner#distance(double, double,
double)} method.
+     *
+     * @throws TransformException if an error occurred while projecting a coordinate.
+     */
+    @Test
+    public void testSimplePoint() throws TransformException {
+        transform = create(false);
+        /*
+         * Fix φ=45°, which implies tan(φ)=1.
+         * Test using the CassiniSoldner spherical equation.
+         */
+        final DirectPosition2D point = new DirectPosition2D();
+        final double domain = toRadians(5);
+        final double step   = domain / 20;
+        for (double λ = -domain; λ <= domain; λ += step) {
+            final double yFromSimplified = PI/2 - atan(cos(λ));
+            point.x = λ;
+            point.y = PI/4;
+            assertSame(point, transform.transform(point, point));
+            assertEquals("Given excentricity=0 and φ=45°, the spherical equation should
simplify to a "
+                    + "very simple expression, which we are testing here.", yFromSimplified,
point.y, 1E-9);
+        }
+    }
+
+    /**
+     * Tests the point given in EPSG example. This is the same test than {@link #runGeoapiTest()}
+     * but is repeated here for easier debugging.
+     *
+     * @throws FactoryException if an error occurred while creating the map projection.
+     * @throws TransformException if an error occurred while projecting a coordinate.
+     */
+    @Test
+    public void testExampleEPSG() throws FactoryException, TransformException {
+        createCompleteProjection(method(),
+                31706587.88,                            // Semi-major axis (Clarke's links)
+                31706587.88 * (20855233./20926348),     // Semi-minor axis (Clarke's links)
+                -(61 + 20./60),                         // Longitude of natural origin
+                10 + (26 + 30./60)/60,                  // Latitude of natural origin
+                NaN,                                    // Standard parallel 1 (none)
+                NaN,                                    // Standard parallel 2 (none)
+                NaN,                                    // Scale factor (none)
+                430000,                                 // False easting  (Clarke's links)
+                325000);                                // False northing (Clarke's links)
+
+        final DirectPosition2D p = new DirectPosition2D(-62, 10);
+        assertSame(p, transform.transform(p, p));
+        assertEquals(66644.94, p.x, 0.005);
+        assertEquals(82536.22, p.y, 0.005);
+
+        assertSame(p, transform.inverse().transform(p, p));
+        assertEquals(-62, p.x, Formulas.ANGULAR_TOLERANCE);
+        assertEquals(+10, p.y, Formulas.ANGULAR_TOLERANCE);
+    }
+
+    /**
+     * Creates a projection and tests the derivatives at a few points.
+     *
+     * @throws TransformException if an error occurred while projecting a coordinate.
+     */
+    @Test
+    public void testDerivative() throws TransformException {
+        final double delta = toRadians(1. / 60) / 1852;         // Approximately 1 meter.
+        derivativeDeltas = new double[] {delta, delta};
+
+        // Tests spherical formulas
+        tolerance = 1E-9;
+        transform = create(false);
+        validate();
+        verifyDerivative(toRadians(+3), toRadians(-6));
+        verifyDerivative(toRadians(-4), toRadians(40));
+
+        // Tests ellipsoidal formulas
+        transform = create(true);
+        validate();
+        verifyDerivative(toRadians(+3), toRadians(-6));
+        verifyDerivative(toRadians(+3), toRadians(-10));
+        verifyDerivative(toRadians(-4), toRadians(+10));
+    }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformFactoryMock.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformFactoryMock.java
index 0be16de..dd148fd 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformFactoryMock.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformFactoryMock.java
@@ -50,7 +50,7 @@ public final strictfp class MathTransformFactoryMock implements MathTransformFac
     /**
      * Creates a new mock for the given operation method.
      *
-     * @param method The operation method to put in this factory.
+     * @param  method  the operation method to put in this factory.
      */
     public MathTransformFactoryMock(final DefaultOperationMethod method) {
         this.method = method;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 3cafaa9..ca6e15c 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -181,6 +181,7 @@ import org.junit.BeforeClass;
     org.apache.sis.referencing.operation.projection.LambertConicConformalTest.class,
     org.apache.sis.referencing.operation.projection.TransverseMercatorTest.class,
     org.apache.sis.referencing.operation.projection.ZonedGridSystemTest.class,
+    org.apache.sis.referencing.operation.projection.CassiniSoldnerTest.class,
     org.apache.sis.referencing.operation.projection.PolarStereographicTest.class,
     org.apache.sis.referencing.operation.projection.ObliqueStereographicTest.class,
     org.apache.sis.referencing.operation.projection.ObliqueMercatorTest.class,


Mime
View raw message