sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Initial implementation of Cassini-Soldner projection written from EPSG formulas. This commit provides a direct, non-optimized implementation with almost verbatim copy of formulas published in:
Date Mon, 27 Apr 2020 21:10:58 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 b44873a08fd4c02c8092d97c9ab9d8d9964218c4
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Apr 27 23:01:22 2020 +0200

    Initial implementation of Cassini-Soldner projection written from EPSG formulas.
    This commit provides a direct, non-optimized implementation with almost verbatim
    copy of formulas published in:
    
      IOGP Publication 373-7-2. Geomatics Guidance Note number 7, part 2.
      Coordinate Conversions and Transformations including Formulas.
      Revised September 2019.
      §3.2.2 Cassini-Soldner
    
    Next commits will provide optimizations, accuracy improvements,
    duplicated code removal, tests, etc.
---
 .../referencing/provider/CassiniSoldner.java       | 169 ++++++++++++++++++++
 .../operation/projection/CassiniSoldner.java       | 176 +++++++++++++++++++++
 ...g.opengis.referencing.operation.OperationMethod |   1 +
 .../referencing/provider/ProvidersTest.java        |   1 +
 4 files changed, 347 insertions(+)

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..2748288
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CassiniSoldner.java
@@ -0,0 +1,176 @@
+/*
+ * 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.parameter.ParameterDescriptor;
+import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.OperationMethod;
+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 NormalizedProjection {
+
+    private final double M0;
+
+    /**
+     * 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));
+        M0 = M(φ0);
+    }
+
+    /**
+     * 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
+    {
+        final double λ    = srcPts[srcOff  ];
+        final double φ    = srcPts[srcOff+1];
+        final double cosφ = cos(φ);
+        final double sinφ = sin(φ);
+        final double tanφ = sinφ / cosφ;
+        final double A    = λ * cosφ;
+        final double A2   = A*A;
+        final double A3   = A*A2;
+        final double T    = tanφ * tanφ;
+        final double C    = eccentricitySquared * (cosφ*cosφ) / (1 - eccentricitySquared);
+        final double ν    = 1 / sqrt(1 - eccentricitySquared*(sinφ*sinφ));
+        final double y    = M(φ) - M0 + ν*tanφ*(A2/2 + (5 - T + 6*C)*(A2*A2) / 24);
+        final double x    = ν*(A - T*A3/6 - (8 - T + 8*C)*T*(A3*A2) / 120);
+        dstPts[dstOff  ]  = x;
+        dstPts[dstOff+1]  = y;
+        return null;
+    }
+
+    private double M(final double φ) {
+        final double e2 = eccentricitySquared;
+        final double e4 = e2*e2;
+        final double e6 = e4*e2;
+        final double sin2φ = sin(2*φ);
+        final double sin4φ = sin(4*φ);
+        final double sin6φ = sin(6*φ);
+        return (1  - e2/4 - 3*e4/64  -  5*e6/256 ) *     φ
+                - (3*e2/8 + 3*e4/32  + 45*e6/1024) * sin2φ
+                +         (15*e4/256 + 45*e6/1024) * sin4φ
+                -                     (35*e6/3072) * sin6φ;
+    }
+
+    /**
+     * 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 e2     = eccentricitySquared;
+        final double e4     = e2*e2;
+        final double e6     = e4*e2;
+        final double e1     = (1 - sqrt(1 - eccentricitySquared)) / (1 + sqrt(1 - eccentricitySquared));
+        final double e12    = e1*e1;
+        final double e13    = e12*e1;
+        final double e14    = e12*e12;
+        final double M1     = M0 + y;
+        final double μ1     = M1 / (1 - eccentricitySquared/4 - 3*e4/64 - 5*e6/256);
+        final double sinμ1  = sin(μ1);
+        final double sin2μ1 = sin(2*μ1);
+        final double sin4μ1 = sin(4*μ1);
+        final double sin6μ1 = sin(6*μ1);
+        final double sin8μ1 = sin(8*μ1);
+        final double φ1     = μ1 + (3*e1/2 - 27*e13/32) * sin2μ1 + (21*e12/16 - 55*e14/32)*sin4μ1
+                                 + (151*e13/96)*sin6μ1 + (1097*e14/512)*sin8μ1;
+        final double sinφ1  = sin(φ1);
+        final double cosφ1  = cos(φ1);
+        final double tanφ1  = sinφ1 / cosφ1;
+        final double ν1     = 1/sqrt(1 - eccentricitySquared*(sinφ1*sinφ1));
+        final double ρ1     = (1 - eccentricitySquared) / pow(1 - eccentricitySquared*(sinφ1*sinφ1),
1.5);
+        final double D      = x/ν1;
+        final double D2     = D*D;
+        final double D4     = D2*D2;
+        final double T1     = tanφ1 * tanφ1;
+        final double φ      = φ1 - (ν1*tanφ1 / ρ1) * (D2/2 - (1 + 3*T1)*D4 / 24);
+        final double λ      = (D - T1*(D2*D)/3 + (1 + 3*T1)*T1*(D4*D)/15) / cosφ1;
+        dstPts[dstOff  ]    = λ;
+        dstPts[dstOff+1]    = φ;
+    }
+}
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 0e68378..f612bbd 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
@@ -93,6 +93,7 @@ public final strictfp class ProvidersTest extends TestCase {
             AlbersEqualArea.class,
             TransverseMercator.class,
             TransverseMercatorSouth.class,
+            CassiniSoldner.class,
             PolarStereographicA.class,
             PolarStereographicB.class,
             PolarStereographicC.class,


Mime
View raw message