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,
|