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 c8ed9467ae089b2d3a9327cb23c1a1be6a4cef65 Author: Martin Desruisseaux <[email protected]> AuthorDate: Fri Apr 29 18:55:53 2022 +0200 Port `LambertAzimuthalEqualArea` from Geotk. We rewrote the equations from EPSG guidance notes. https://issues.apache.org/jira/browse/SIS-228 --- .../provider/LambertAzimuthalEqualArea.java | 186 +++++++++++ .../LambertAzimuthalEqualAreaSpherical.java | 77 +++++ .../java/org/apache/sis/referencing/Builder.java | 14 +- .../operation/projection/EqualAreaProjection.java | 29 +- .../projection/LambertAzimuthalEqualArea.java | 283 ++++++++++++++++ ...g.opengis.referencing.operation.OperationMethod | 2 + .../referencing/provider/ProvidersTest.java | 2 + .../projection/LambertAzimuthalEqualAreaTest.java | 357 +++++++++++++++++++++ .../sis/test/suite/ReferencingTestSuite.java | 1 + 9 files changed, 938 insertions(+), 13 deletions(-) diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/LambertAzimuthalEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/LambertAzimuthalEqualArea.java new file mode 100644 index 0000000000..18f34bacd1 --- /dev/null +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/LambertAzimuthalEqualArea.java @@ -0,0 +1,186 @@ +/* + * 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.opengis.referencing.operation.PlanarProjection; +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.NormalizedProjection; + + +/** + * The provider for "<cite>Lambert Azimuthal Equal Area</cite>" projection (EPSG:9820). + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.2 + * + * @see <a href="http://geotiff.maptools.org/proj_list/lambert_azimuthal_equal_area.html";>GeoTIFF parameters for Lambert Azimuthal Equal Area</a> + * + * @since 1.2 + * @module + */ +@XmlTransient +public class LambertAzimuthalEqualArea extends MapProjection { + /** + * For cross-version compatibility. + */ + private static final long serialVersionUID = 3360095493388421774L; + + /** + * The EPSG identifier, to be preferred to the name when available. + * The {@code IDENTIFIER_OF_BASE} name is for avoiding confusion with + * the {@code IDENTIFIER} name used in subclasses. + */ + public static final String IDENTIFIER_OF_BASE = "9820"; + + /** + * The operation parameter descriptor for the <cite>Latitude of natural origin</cite> (φ₀) parameter value. + * Valid values range is [-90 … 90]°. The 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_center </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> CenterLat </td></tr> + * <tr><td> Proj4: </td><td> lat_0 </td></tr> + * </table> + */ + public static final ParameterDescriptor<Double> 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> longitude_of_center </td></tr> + * <tr><td> ESRI: </td><td> Central_Meridian </td></tr> + * <tr><td> NetCDF: </td><td> longitude_of_projection_origin </td></tr> + * <tr><td> GeoTIFF: </td><td> CenterLong </td></tr> + * <tr><td> Proj4: </td><td> lon_0 </td></tr> + * </table> + */ + public static final ParameterDescriptor<Double> LONGITUDE_OF_ORIGIN; + + /** + * 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 = Equirectangular.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 = Equirectangular.FALSE_NORTHING; + + /** + * The group of all parameters expected by this coordinate operation. + */ + private static final ParameterDescriptorGroup PARAMETERS; + static { + final ParameterBuilder builder = builder(); + LATITUDE_OF_ORIGIN = createLatitude(renameAlias(builder, Equirectangular.LATITUDE_OF_ORIGIN, + Citations.OGC, AlbersEqualArea.LATITUDE_OF_FALSE_ORIGIN), true); + + LONGITUDE_OF_ORIGIN = createLongitude(renameAlias(builder, Equirectangular.LONGITUDE_OF_ORIGIN, + Citations.OGC, AlbersEqualArea.LONGITUDE_OF_FALSE_ORIGIN)); + + PARAMETERS = builder + .addIdentifier(IDENTIFIER_OF_BASE) + .addName( "Lambert Azimuthal Equal Area") + .addName(Citations.OGC, "Lambert_Azimuthal_Equal_Area") + .addName(Citations.ESRI, "Lambert_Azimuthal_Equal_Area") + .addName(Citations.NETCDF, "LambertAzimuthalEqualArea") + .addName(Citations.GEOTIFF, "CT_LambertAzimEqualArea") + .addName(Citations.PROJ4, "laea") + .addIdentifier(Citations.GEOTIFF, "10") + .addIdentifier(Citations.MAP_INFO, "4") + .createGroupForMapProjection( + LATITUDE_OF_ORIGIN, + LONGITUDE_OF_ORIGIN, + FALSE_EASTING, + FALSE_NORTHING); + } + + /** + * Constructs a new provider. + */ + public LambertAzimuthalEqualArea() { + super(PARAMETERS); + } + + /** + * Constructs a math transform provider from a set of parameters. + * + * @param parameters the set of parameters (never {@code null}). + */ + LambertAzimuthalEqualArea(final ParameterDescriptorGroup parameters) { + super(parameters); + } + + /** + * Returns the operation type for this map projection. + * + * @return {@code PlanarProjection.class} + */ + @Override + public Class<PlanarProjection> getOperationType() { + return PlanarProjection.class; + } + + /** + * {@inheritDoc} + * + * @return the map projection created from the given parameter values. + */ + @Override + protected final NormalizedProjection createProjection(final Parameters parameters) { + return new org.apache.sis.referencing.operation.projection.LambertAzimuthalEqualArea(this, parameters); + } +} diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/LambertAzimuthalEqualAreaSpherical.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/LambertAzimuthalEqualAreaSpherical.java new file mode 100644 index 0000000000..6be0502cd2 --- /dev/null +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/LambertAzimuthalEqualAreaSpherical.java @@ -0,0 +1,77 @@ +/* + * 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.apache.sis.parameter.ParameterBuilder; +import org.opengis.parameter.ParameterDescriptor; +import org.opengis.parameter.ParameterDescriptorGroup; + + +/** + * The provider for "<cite>Lambert Azimuthal Equal Area (Spherical)</cite>" projection (EPSG:1027). + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.2 + * @since 1.2 + * @module + */ +@XmlTransient +public final class LambertAzimuthalEqualAreaSpherical extends LambertAzimuthalEqualArea { + /** + * For cross-version compatibility. + */ + private static final long serialVersionUID = -8356440804799174833L; + + /** + * The EPSG identifier, to be preferred to the name when available. + */ + public static final String IDENTIFIER = "1027"; + + /** + * The group of all parameters expected by this coordinate operation. + */ + static final ParameterDescriptorGroup PARAMETERS; + static { + final ParameterBuilder builder = builder(); + + final ParameterDescriptor<Double> latitudeOfOrigin = createLatitude(builder + .addNamesAndIdentifiers(LATITUDE_OF_ORIGIN).setDeprecated(true) + .addName("Spherical latitude of origin"), true); + + final ParameterDescriptor<Double> longitudeOfOrigin = createLongitude(builder + .addNamesAndIdentifiers(LONGITUDE_OF_ORIGIN).setDeprecated(true) + .addName("Spherical longitude of origin")); + + PARAMETERS = builder + .addIdentifier(IDENTIFIER) + .addName("Lambert Azimuthal Equal Area (Spherical)") + .setDeprecated(true).addIdentifier("9821").setDeprecated(false) + .createGroupForMapProjection( + latitudeOfOrigin, + longitudeOfOrigin, + FALSE_EASTING, + FALSE_NORTHING); + } + + /** + * Constructs a new provider. + */ + public LambertAzimuthalEqualAreaSpherical() { + super(PARAMETERS); + } +} diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/Builder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/Builder.java index ae12007d43..81e054b5c2 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/Builder.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/Builder.java @@ -985,18 +985,16 @@ public abstract class Builder<B extends Builder<B>> { } /** - * Sets whether the next {@code Identifier} or {@code IdentifiedObject}s to create shall be considered deprecated. - * Deprecated objects exist in some {@linkplain org.opengis.referencing.AuthorityFactory authority factories} like - * the EPSG database. - * - * <p>Note that this method does not apply to name and identifiers, which have their own - * {@code addDeprecatedFoo(…)} methods.</p> + * Sets whether the next {@code GenericName}s, {@code Identifier}s or {@code IdentifiedObject}s + * to create shall be considered deprecated. + * Deprecated objects exist in some {@linkplain org.opengis.referencing.AuthorityFactory authority factories} + * like the EPSG database. * * <p><b>Lifetime:</b> * Deprecation status is cleared after a {@code createXXX(…)} method has been invoked.</p> * - * @param deprecated {@code true} if the next names, identifiers and identified objects should be - * considered deprecated, or {@code false} otherwise. + * @param deprecated {@code true} if the next names, identifiers and identified + * objects should be considered deprecated, or {@code false} otherwise. * @return {@code this}, for method call chaining. * * @see AbstractIdentifiedObject#isDeprecated() diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java index 47a47558c4..3abebb6a42 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java @@ -32,7 +32,7 @@ import static org.apache.sis.math.MathFunctions.atanh; * hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct classes.</p> * * @author Martin Desruisseaux (Geomatys) - * @version 1.0 + * @version 1.2 * @since 0.8 * @module */ @@ -40,7 +40,7 @@ abstract class EqualAreaProjection extends NormalizedProjection { /** * For cross-version compatibility. */ - private static final long serialVersionUID = 2221537810038553082L; + private static final long serialVersionUID = 5880625564193782957L; /** * The threshold value of {@link #eccentricity} at which we consider the accuracy of the @@ -65,7 +65,7 @@ abstract class EqualAreaProjection extends NormalizedProjection { * * <div class="note"><b>Serialization note:</b> * we do not strictly need to serialize those fields since they could be computed after deserialization. - * Bu we serialize them anyway in order to simplify a little bit this class (it allows us to keep those + * But we serialize them anyway in order to simplify a little bit this class (it allows us to keep those * fields final) and because values computed after deserialization could be slightly different than the * ones computed after construction if a future version of the constructor uses the double-double values * provided by {@link Initializer}.</div> @@ -74,8 +74,21 @@ abstract class EqualAreaProjection extends NormalizedProjection { /** * Value of {@link #qm(double)} function (part of Snyder equation (3-12)) at pole (sinφ = 1). + * The <var>q</var><sub>p</sub> constant is defined by EPSG guidance note as: + * + * <blockquote>(1 – ℯ²)⋅([1 / (1 – ℯ²)] – {[1/(2ℯ)] ln [(1 – ℯ) / (1 + ℯ)]})</blockquote> + * + * But in this class, we omit the (1 – ℯ²) factor. + * + * <h4>Spherical case</h4> + * In the spherical case (ℯ=0), the value is exactly 2. */ - private final double qmPolar; + final double qmPolar; + + /** + * {@code true} if {@link #eccentricity} is zero. + */ + final boolean isSpherical; /** * {@code true} if the {@link #eccentricity} value is greater than or equals to the eccentricity threshold, @@ -90,6 +103,7 @@ abstract class EqualAreaProjection extends NormalizedProjection { */ EqualAreaProjection(final Initializer initializer) { super(initializer); + isSpherical = (eccentricitySquared == 0); final double e2 = eccentricitySquared; final double e4 = e2 * e2; final double e6 = e2 * e4; @@ -135,6 +149,7 @@ abstract class EqualAreaProjection extends NormalizedProjection { c4β = other.c4β; c6β = other.c6β; qmPolar = other.qmPolar; + isSpherical = other.isSpherical; useIterations = other.useIterations; } @@ -165,7 +180,7 @@ abstract class EqualAreaProjection extends NormalizedProjection { * Check for zero eccentricity is required because qm_ellipsoid(sinφ) would * simplify to sinφ + atanh(0) / 0 == sinφ + 0/0, thus producing NaN. */ - return (eccentricity == 0) ? 2*sinφ : qm_ellipsoid(sinφ); + return isSpherical ? 2*sinφ : qm_ellipsoid(sinφ); } /** @@ -202,6 +217,10 @@ abstract class EqualAreaProjection extends NormalizedProjection { * have a sufficient amount of terms for achieving the centimetric precision, so we "finish" it by the * iterative method. The series expansion is nevertheless useful for reducing the number of iterations. * + * <h4>Spherical case</h4> + * In the spherical case, this method returns {@code asin(y/2)}. This method does not use that + * simplification for the spherical case. This optimization is left to the caller if desired. + * * @param y in the cylindrical case, this is northing on the normalized ellipsoid. * @return the latitude in radians. */ diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java new file mode 100644 index 0000000000..8a4cac76ee --- /dev/null +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java @@ -0,0 +1,283 @@ +/* + * 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.referencing.operation.matrix.Matrix2; +import org.apache.sis.referencing.operation.matrix.MatrixSIS; +import org.apache.sis.referencing.operation.transform.ContextualParameters; +import org.apache.sis.internal.util.DoubleDouble; +import org.apache.sis.parameter.Parameters; +import org.apache.sis.util.Workaround; + +import static java.lang.Math.*; +import static org.apache.sis.internal.referencing.Formulas.fastHypot; +import static org.apache.sis.internal.referencing.provider.LambertAzimuthalEqualArea.*; + + +/** + * <cite>Lambert Azimuthal Equal Area</cite> projection (EPSG code 9820). + * See the following references for an overview: + * <ul> + * <li><a href="https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection";>Lambert Azimuthal projection on Wikipedia</a></li> + * <li><a href="https://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html";>Lambert Azimuthal Equal-Area projection on MathWorld</a></li> + * </ul> + * + * @author Martin Desruisseaux (Geomatys) + * @author Rémi Maréchal (Geomatys) + * @version 1.2 + * @since 1.2 + * @module + */ +public class LambertAzimuthalEqualArea extends EqualAreaProjection { + /** + * For cross-version compatibility. + */ + private static final long serialVersionUID = 7419101460426922558L; + + /** + * Sine and cosine of authalic latitude of origin. In the spherical case, + * the authalic latitude ß and the geodetic latitude φ are the same. + * + * <h4>Possible simplification</h4> + * In equatorial case, sin(ß₀)=0 and cos(ß₀)=1. We could do special cases with simplifications in the + * {@code transform(…)} formulas, but the result does not seem simpler enough to worth code branching. + * + * <p>In the polar case, sin(ß₀)=1 and cos(ß₀)=0. But the equations become indeterminate (we get 0/0) + * and a different set of equations must be used.</p> + */ + private final double sinß0, cosß0; + + /** + * {@code true} if the projection is at a pole. + * This implementation does not need to distinguish between North and South pole. + * Formulas in this class are for the South case only, and this class handles the + * North case by reverting the sign of φ before conversion and y after conversion. + */ + private final boolean polar; + + /** + * Creates a Lambert Azimuthal Equal Area projection from the given parameters. + * + * @param method description of the projection parameters. + * @param parameters the parameter values of the projection to create. + */ + public LambertAzimuthalEqualArea(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"). + */ + @SuppressWarnings("fallthrough") + @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.FALSE_EASTING, FALSE_EASTING); + roles.put(ParameterRole.FALSE_NORTHING, FALSE_NORTHING); + return new Initializer(method, parameters, roles, null); + } + + /** + * 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 LambertAzimuthalEqualArea(final Initializer initializer) { + super(initializer); + final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN)); + final double sinφ0 = sin(φ0); + final double cosφ0 = cos(φ0); + polar = abs(sinφ0) == 1; // sin(φ) == 1 implies a tolerance ∆φ≈1E-6°. + sinß0 = qm(sinφ0) / qmPolar; + cosß0 = sqrt(1 - sinß0*sinß0); + /* + * In the polar case we have cos(φ₀) ≈ 0 and cos(ß₀) ≈ 0, which cause D = 0/0. + * Trying to evaluate the indeterminate with L'Hôpital's rule produce infinity. + * Consequently a different set of formulas for the polar form must be used not + * only here but also in the `transform(…)` and `inverseTransform(…)` methods. + */ + final MatrixSIS denormalize = getContextualParameters().getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION); + if (polar) { + final DoubleDouble c = initializer.axisLengthRatio(); + denormalize.convertBefore(0, c, null); + denormalize.convertBefore(1, c, null); + if (φ0 > 0) { // North pole case: use South pole formulas with sign of φ and y reverted. + final MatrixSIS normalize = getContextualParameters().getMatrix(ContextualParameters.MatrixRole.NORMALIZATION); + normalize .convertBefore(1, -1, null); + denormalize.convertBefore(1, -1, null); + } + } else { + final double D = cosφ0 / (sqrt(1 - eccentricitySquared*(sinφ0*sinφ0)) * cosß0); + final double Rq2 = (1 - eccentricitySquared) * qmPolar/2; + denormalize.convertBefore(0, D, null); + denormalize.convertBefore(1, Rq2/D, null); + } + } + + /** + * Creates a new projection initialized to the same parameters than the given one. + */ + LambertAzimuthalEqualArea(final LambertAzimuthalEqualArea other) { + super(other); + sinß0 = other.sinß0; + cosß0 = other.cosß0; + polar = other.polar; + } + + /** + * Returns the names of additional internal parameters which need to be taken in account when + * comparing two {@code LambertAzimuthalEqualArea} projections or formatting them in debug mode. + */ + @Override + final String[] getInternalParameterNames() { + return new String[] {"ß₀"}; + } + + /** + * Returns the values of additional internal parameters which need to be taken in account when + * comparing two {@code LambertAzimuthalEqualArea} projections or formatting them in debug mode. + */ + @Override + final double[] getInternalParameterValues() { + return new double[] {asin(sinß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 sinλ = sin(λ); + final double cosλ = cos(λ); + final double sinφ = sin(φ); + final double qm = qm(sinφ); + if (!polar) { + final double sinß = qm / qmPolar; + final double cosß = sqrt(1 - sinß*sinß); + final double c = sinß0*sinß + cosß0*cosß*cosλ + 1; + /* + * The `c` factor is 0 when projecting the antipodal point of origin. + * The antipodal point is actually all points on a circle of radius 2. + * We can not return a single point, so NaN seems a safer value. + */ + double B = sqrt(2 / c); + if (B == Double.POSITIVE_INFINITY) { + B = Double.NaN; + } + final double x = cosß*sinλ; // (x,y)×B = (E,N) + final double y = cosß0*sinß - sinß0*cosß*cosλ; + if (dstPts != null) { + dstPts[dstOff ] = B * x; + dstPts[dstOff+1] = B * y; + } + /* + * End of map projection. Compute the derivative if it was requested. + */ + if (!derivate) { + return null; + } + final double dsinß_dφ = dqm_dφ(sinφ, cos(φ)) / qmPolar; + final double dcosß_dφ = -dsinß_dφ * (sinß/cosß); + final double cosλdcosß = dcosß_dφ * cosλ; + final double dy_dλ = sinß0*x; + final double db_dλ = cosß0*x / (2*c); + final double dy_dφ = (cosß0*dsinß_dφ - sinß0*cosλdcosß); + final double db_dφ = -(sinß0*dsinß_dφ + cosß0*cosλdcosß) / (2*c); + return new Matrix2( + B * (cosλ + db_dλ*sinλ)*cosß, + B * (dcosß_dφ + db_dφ*cosß)*sinλ, + B * (dy_dλ + db_dλ*y), + B * (dy_dφ + db_dφ*y)); + } + /* + * Polar case needs special code because formula for the oblique case become indeterminate. + * North pole case is handled with the same formula, but with sign of φ reversed by the normalize + * affine transform (which cause the sign of q to be also reversed). After the transform, sign of + * y will be reversed by the denormalize affine transform, so it should not be reversed here. + */ + double ρ = sqrt(qmPolar + qm); + if (sinφ == 1) { + ρ = Double.NaN; // Antipodal point. + } + final double x = ρ * sinλ; + final double y = ρ * cosλ; + if (dstPts != null) { + dstPts[dstOff ] = x; + dstPts[dstOff+1] = y; + } + if (!derivate) { + return null; + } + final double db_dφ = dqm_dφ(sinφ, cos(φ)) / (2*ρ); + return new Matrix2(y, db_dφ * sinλ, + -x, db_dφ * cosλ); + } + + /** + * Converts the specified (<var>x</var>,<var>y</var>) coordinates and stores the (λ,φ) result in {@code dstPts}. + * + * @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 + { + double x = srcPts[srcOff ]; + double y = srcPts[srcOff+1]; + double sinß_qp; // sin(ß) × qmPolar + if (polar) { + sinß_qp = (x*x + y*y) - qmPolar; + } else { + final double ρ = fastHypot(x, y); + final double C = 2*asin(ρ / 2); + final double cosC = cos(C); + final double sinC = sin(C); + final double ysinC = y*sinC; + sinß_qp = cosC*sinß0; + if (cosC != 1) { // cos(C) == 1 implies y/ρ = 0/0. + sinß_qp += ysinC*cosß0/ρ; + } + sinß_qp *= qmPolar; + y = ρ*cosC*cosß0 - ysinC*sinß0; + x *= sinC; + } + dstPts[dstOff ] = atan2(x, y); + dstPts[dstOff+1] = isSpherical ? asin(sinß_qp/2) : φ(sinß_qp); + } + + /* + * We do not provide a specialized sub-class for the spherical case + * because the simplifications are too small. + */ +} 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 0e11fa01f6..52f61edf69 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 @@ -42,6 +42,8 @@ org.apache.sis.internal.referencing.provider.LambertConformalBelgium org.apache.sis.internal.referencing.provider.LambertConformalMichigan org.apache.sis.internal.referencing.provider.LambertCylindricalEqualArea org.apache.sis.internal.referencing.provider.LambertCylindricalEqualAreaSpherical +org.apache.sis.internal.referencing.provider.LambertAzimuthalEqualArea +org.apache.sis.internal.referencing.provider.LambertAzimuthalEqualAreaSpherical org.apache.sis.internal.referencing.provider.AlbersEqualArea org.apache.sis.internal.referencing.provider.TransverseMercator org.apache.sis.internal.referencing.provider.TransverseMercatorSouth 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 4810aa159a..fb24e65e70 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 @@ -91,6 +91,8 @@ public final strictfp class ProvidersTest extends TestCase { LambertConformalMichigan.class, LambertCylindricalEqualArea.class, LambertCylindricalEqualAreaSpherical.class, + LambertAzimuthalEqualArea.class, + LambertAzimuthalEqualAreaSpherical.class, AlbersEqualArea.class, TransverseMercator.class, TransverseMercatorSouth.class, diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualAreaTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualAreaTest.java new file mode 100644 index 0000000000..72427c2dac --- /dev/null +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualAreaTest.java @@ -0,0 +1,357 @@ +/* + * 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.referencing.operation.transform.CoordinateDomain; +import org.apache.sis.referencing.operation.transform.MathTransformFactoryMock; +import org.apache.sis.internal.referencing.provider.MapProjection; +import org.apache.sis.internal.metadata.ReferencingServices; +import org.apache.sis.parameter.Parameters; +import org.junit.Test; + +import static java.lang.Double.NaN; +import static java.lang.StrictMath.*; +import static org.apache.sis.math.MathFunctions.SQRT_2; +import static org.apache.sis.internal.referencing.Formulas.LINEAR_TOLERANCE; +import static org.apache.sis.internal.referencing.Formulas.ANGULAR_TOLERANCE; +import static org.junit.Assert.*; + + +/** + * Tests the {@link LambertAzimuthalEqualArea} class. We test using various values + * of the latitude of origin, which is the only parameter impacting the internal + * coefficients of that class (except for the eccentricity). + * + * @author Martin Desruisseaux (Geomatys) + * @author Rémi Maréchal (Geomatys) + * @version 1.2 + * @since 1.2 + * @module + */ +public final strictfp class LambertAzimuthalEqualAreaTest extends MapProjectionTestCase { + /** + * The radius of the sphere used in sphere test cases. + */ + private static final double SPHERE_RADIUS = ReferencingServices.AUTHALIC_RADIUS; + + /** + * Returns the provider for the map projection to tesT. + */ + private static MapProjection provider(final boolean elliptical) { + return elliptical ? new org.apache.sis.internal.referencing.provider.LambertAzimuthalEqualArea() + : new org.apache.sis.internal.referencing.provider.LambertAzimuthalEqualAreaSpherical(); + } + + /** + * Creates and validates a new instance of {@link LambertAzimuthalEqualArea}. + * + * @param elliptical {@code false} for a sphere, or {@code true} for WGS84 ellipsoid. + * @param latitudeOfOrigin the latitude of origin. + * @param complete whether to create the full projection, working on degrees instead of radians. + * @throws FactoryException if an error occurred while creating the map projection. + */ + private void createProjection(final boolean elliptical, final double latitudeOfOrigin, final boolean complete) + throws FactoryException + { + final MapProjection provider = provider(elliptical); + final Parameters parameters = parameters(provider, elliptical); + parameters.parameter("latitude_of_origin").setValue(latitudeOfOrigin); + LambertAzimuthalEqualArea projection = new LambertAzimuthalEqualArea(provider, parameters); + if (complete) { + transform = projection.createMapProjection(new MathTransformFactoryMock(provider)); + } else { + transform = projection; + } + validate(); + } + + /** + * Tests oblique case with a point computed with PROJ. + * Command line was: + * + * <pre>cs2cs "EPSG:4326" +to +type=crs +proj=laea +no_defs +lat_0=90 +lon_0=0 +datum=WGS84 <<< "72 50"</pre> + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testPolar() throws FactoryException, TransformException { + tolerance = LINEAR_TOLERANCE; + createProjection(true, 90, true); + final double[] point = new double[] {50, 72}; + final double[] expected = new double[] {1533302.80, -1286593.82}; + verifyTransform(point, expected); + } + + /** + * Tests oblique case with a point computed with PROJ. + * Command line was: + * + * <pre>cs2cs "EPSG:4326" +to +type=crs +proj=laea +no_defs +lat_0=37 +lon_0=0 +datum=WGS84 <<< "25 -9"</pre> + * + * Note that {@link #runGeoapiTest()} performs a similar test. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testOblique() throws FactoryException, TransformException { + tolerance = LINEAR_TOLERANCE; + createProjection(true, 37, true); + final double[] point = new double[] {-9, 25}; + final double[] expected = new double[] {-911656.53, -1288191.92}; + verifyTransform(point, expected); + } + + /** + * Tests self-consistency of the unitary projection on a sphere. + * The projection works with radians on a sphere of radius 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testConsistencyOfUnitaryOnSphere() throws FactoryException, TransformException { + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + derivativeDeltas = new double[] {delta, delta}; + tolerance = ANGULAR_TOLERANCE; + for (int φ = -90; φ <= 90; φ += 15) { + createProjection(false, φ, false); + verifyInDomain(CoordinateDomain.GEOGRAPHIC_RADIANS, 862247543); + } + } + + /** + * Tests self-consistency of the unitary projection on an ellipse. + * The projection works with radians on an ellipsoid of semi-major axis length 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testConsistencyOfUnitaryOnEllipse() throws FactoryException, TransformException { + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + derivativeDeltas = new double[] {delta, delta}; + tolerance = ANGULAR_TOLERANCE; + for (int φ = -90; φ <= 90; φ += 15) { + createProjection(true, φ, false); + verifyInDomain(CoordinateDomain.GEOGRAPHIC_RADIANS, 484117986); + } + } + + /** + * Tests the projection at a few particular points in the oblique case. + * The particular points are the origin and the (almost) antipodal points. + * The projection works with radians on a sphere of radius 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testParticularPointsWithObliqueOnSphere() throws FactoryException, TransformException { + testParticularPointsWithOblique(false); + } + + /** + * Tests the projection at a few particular points in the oblique case. + * The particular points are the origin and the (almost) antipodal points. + * The projection works with degrees on an ellipsoid of semi-major axis length 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testParticularPointsWithObliqueOnEllipse() throws FactoryException, TransformException { + testParticularPointsWithOblique(true); + } + + /** + * Implementation of {@link #testParticularPointsWithObliqueOnSphere()} + * and {@link #testParticularPointsWithObliqueOnEllipse()}. + */ + private void testParticularPointsWithOblique(final boolean elliptical) throws FactoryException, TransformException { + tolerance = LINEAR_TOLERANCE; + createProjection(elliptical, 45, true); + + // Projects the origin. + final double[] point = new double[] {0, 45}; + final double[] expected = new double[] {0, 0}; + verifyTransform(point, expected); + + // Project the antipode. + point[0] = 180; + point[1] = -45; + transform.transform(point, 0, point, 0, 1); + assertEquals("E", NaN, point[0], tolerance); + assertEquals("N", NaN, point[1], tolerance); + + // Project the almost-antipode. + point[0] = 180; + point[1] = -44.999; + transform.transform(point, 0, point, 0, 1); + assertEquals("E", 0, point[0], tolerance); + assertEquals("N", 2*SPHERE_RADIUS, point[1], 0.002*SPHERE_RADIUS); + } + + /** + * Tests the projection at a few extreme points in the polar case. + * The projection works with radians on a sphere of radius 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testParticularPointsWithPolarOnSphere() throws FactoryException, TransformException { + testParticularPointsWithPolar(false); + } + + /** + * Tests the projection at a few extreme points in the polar case. + * The projection works with radians on an ellipsoid of semi-major axis length 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testParticularPointsWithPolarOnEllipse() throws FactoryException, TransformException { + testParticularPointsWithPolar(true); + } + + /** + * Implementation of {@link #testParticularPointsWithPolarOnSphere()} + * and {@link #testParticularPointsWithPolarOnEllipse()}. + */ + private void testParticularPointsWithPolar(final boolean elliptical) throws FactoryException, TransformException { + tolerance = LINEAR_TOLERANCE; + createProjection(elliptical, 90, true); + /* + * Project the origin. Result should be (0,0). Do not test the inverse projection + * because the longitude could be anything and still be the North pole. We test that + * by projecting again with an other longitude, set to 45°, and expect the same result. + */ + final double[] point = new double[] {0, 90}; + final double[] expected = new double[] {0, 0}; + isInverseTransformSupported = false; + verifyTransform(point, expected); + transform.inverse().transform(expected, 0, point, 0, 1); + assertEquals("λ", 0, point[0], 180); + assertEquals("φ", 90, point[1], ANGULAR_TOLERANCE); + /* + * Same point as above (i.e. projection origin), but using a different longitude value. + * Same expected result because longitude should have no effect at a pole. + */ + point[0] = 45; + point[1] = 90; + verifyTransform(point, expected); + transform.inverse().transform(expected, 0, point, 0, 1); + assertEquals("λ", 45, point[0], 180); + assertEquals("φ", 90, point[1], ANGULAR_TOLERANCE); + /* + * Project a point on the equator, at 0° and at 180° longitude. + * Result should be (0, √2) positive or negative depending on the longitude. + */ + point[0] = 0; expected[0] = 0; + point[1] = 0; expected[1] = -SQRT_2*SPHERE_RADIUS; + isInverseTransformSupported = true; + if (elliptical) tolerance = 0.5; // Because the use of `SPHERE_RADIUS` is approximate. + verifyTransform(point, expected); + + point[0] = 180; expected[0] = 0; + point[1] = 0; expected[1] = SQRT_2*SPHERE_RADIUS; + verifyTransform(point, expected); + /* + * Project the antipode. Result would be (0, -2) if this operation was allowed. + * Actually the formulas would work, but every points on a circle or radius 2 + * would be the pole so returning a single value may not be appropriate. + */ + point[0] = 0; + point[1] = -90; + transform.transform(point, 0, point, 0, 1); + assertEquals("E", NaN, point[0], tolerance); + assertEquals("N", NaN, point[1], tolerance); + } + + /** + * Tests consistency when converting a point forward then backward. + * The projection works with radians on a sphere of radius 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testRoundtrip() throws FactoryException, TransformException { + tolerance = ANGULAR_TOLERANCE; + createProjection(false, -75, false); + final double[] point = new double[] {toRadians(-90), toRadians(-8)}; + transform.transform(point, 0, point, 0, 1); + transform.inverse().transform(point, 0, point, 0, 1); + assertEquals(-90, toDegrees(point[0]), tolerance); + assertEquals( -8, toDegrees(point[1]), tolerance); + } + + /** + * Creates a projection and tests the derivatives at a few points. + * The projection works with radians on an ellipsoid of semi-major axis length 1. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + */ + @Test + public void testDerivative() throws FactoryException, TransformException { + tolerance = 1E-9; + final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres. + derivativeDeltas = new double[] {delta, delta}; + + // Polar projection. + createProjection(true, 90, false); + verifyDerivative(toRadians(-6), toRadians(80)); + + // Intentionally above the pole. + verifyDerivative(toRadians(-6), toRadians(100)); + + // Polar projection, spherical formulas. + createProjection(false, 90, false); + verifyDerivative(toRadians(-6), toRadians(85)); + + // Equatorial projection, spherical formulas. + createProjection(false, 0, false); + verifyDerivative(toRadians(3), toRadians(4)); + + // Oblique projection, ellipsoidal formulas. + createProjection(true, 8, false); + verifyDerivative(toRadians(-6), toRadians(2)); + + // Oblique projection, spherical formulas. + createProjection(false, 8, false); + verifyDerivative(toRadians(-6), toRadians(2)); + } + + /** + * Runs the test defined in the GeoAPI-conformance module. + * + * @throws FactoryException if an error occurred while creating the map projection. + * @throws TransformException if an error occurred while projecting a coordinate. + * + * @see org.opengis.test.referencing.ParameterizedTransformTest#testLambertAzimuthalEqualArea() + */ + @Test + public void runGeoapiTest() throws FactoryException, TransformException { + createGeoApiTest(provider(true)).testLambertAzimuthalEqualArea(); + } +} 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 22a252697f..9417060f15 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 @@ -192,6 +192,7 @@ import org.junit.BeforeClass; org.apache.sis.referencing.operation.projection.EqualAreaProjectionTest.class, org.apache.sis.referencing.operation.projection.CylindricalEqualAreaTest.class, org.apache.sis.referencing.operation.projection.AlbersEqualAreaTest.class, + org.apache.sis.referencing.operation.projection.LambertAzimuthalEqualAreaTest.class, org.apache.sis.referencing.operation.projection.MeridianArcTest.class, org.apache.sis.referencing.operation.projection.SinusoidalTest.class, org.apache.sis.referencing.operation.projection.PolyconicTest.class,
