Initial implementation of AlbersEqualArea. Share some more code with CylindricalEqualArea when applicable.
git-svn-id: https://svn.apache.org/repos/asf/sis/branches/JDK8@1753540 13f79535-47bb-0310-9956-ffa450edef68
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AlbersEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AlbersEqualArea.java
new file mode 100644
index 0000000..d09c984
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AlbersEqualArea.java
@@ -0,0 +1,197 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing.provider;
+
+import javax.xml.bind.annotation.XmlTransient;
+import org.opengis.parameter.ParameterDescriptor;
+import org.opengis.parameter.ParameterDescriptorGroup;
+import org.opengis.referencing.operation.ConicProjection;
+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>Albers Equal Area</cite>" projection (EPSG:9822).
+ *
+ * @author Rueben Schulz (UBC)
+ * @author Martin Desruisseaux (Geomatys)
+ * @since 0.8
+ * @version 0.8
+ * @module
+ *
+ * @see <a href="http://www.remotesensing.org/geotiff/proj_list/albers_equal_area_conic.html">Albers Equal-Area Conic on RemoteSensing.org</a>
+ */
+@XmlTransient
+public final class AlbersEqualArea extends MapProjection {
+ /**
+ * For cross-version compatibility.
+ */
+ private static final long serialVersionUID = -7489679528438418778L;
+
+ /**
+ * The operation parameter descriptor for the <cite>Latitude of false origin</cite> (φf) parameter value.
+ * Valid values range is [-90 … 90]° and default value is 0°.
+ */
+ public static final ParameterDescriptor<Double> LATITUDE_OF_FALSE_ORIGIN;
+
+ /**
+ * The operation parameter descriptor for the <cite>Longitude of false origin</cite> (λf) parameter value.
+ * Valid values range is [-180 … 180]° and default value is 0°.
+ */
+ public static final ParameterDescriptor<Double> LONGITUDE_OF_FALSE_ORIGIN;
+
+ /**
+ * The operation parameter descriptor for the <cite>Latitude of 1st standard parallel</cite> parameter value.
+ * Valid values range is [-90 … 90]° and default value is the value given to the {@link #LATITUDE_OF_FALSE_ORIGIN}
+ * parameter.
+ */
+ public static final ParameterDescriptor<Double> STANDARD_PARALLEL_1;
+
+ /**
+ * The operation parameter descriptor for the <cite>Latitude of 2nd standard parallel</cite> parameter value.
+ * Valid values range is [-90 … 90]° and default value is the value given to the {@link #STANDARD_PARALLEL_2}
+ * parameter.
+ */
+ public static final ParameterDescriptor<Double> STANDARD_PARALLEL_2;
+
+ /**
+ * The operation parameter descriptor for the <cite>Easting at false origin</cite> (Ef) parameter value.
+ * Valid values range is unrestricted and default value is 0 metre.
+ */
+ public static final ParameterDescriptor<Double> EASTING_AT_FALSE_ORIGIN;
+
+ /**
+ * The operation parameter descriptor for the <cite>Northing at false origin</cite> (Nf) parameter value.
+ * Valid values range is unrestricted and default value is 0 metre.
+ */
+ public static final ParameterDescriptor<Double> NORTHING_AT_FALSE_ORIGIN;
+
+ /**
+ * The group of all parameters expected by this coordinate operation.
+ */
+ static final ParameterDescriptorGroup PARAMETERS;
+ static {
+ final ParameterBuilder builder = builder();
+ /*
+ * EPSG: Latitude of false origin
+ * OGC: latitude_of_center
+ * ESRI: Latitude_Of_Origin
+ * NetCDF: latitude_of_projection_origin
+ * GeoTIFF: NatOriginLat
+ */
+ LATITUDE_OF_FALSE_ORIGIN = createLatitude(
+ renameAlias(LambertConformal2SP.LATITUDE_OF_FALSE_ORIGIN, Citations.GEOTIFF, Mercator1SP.LATITUDE_OF_ORIGIN, builder)
+ .rename(Citations.OGC, "latitude_of_center"), true);
+ /*
+ * EPSG: Longitude of false origin
+ * OGC: longitude_of_center
+ * ESRI: Central_Meridian
+ * NetCDF: longitude_of_central_meridian
+ * GeoTIFF: NatOriginLong
+ */
+ LONGITUDE_OF_FALSE_ORIGIN = createLongitude(
+ renameAlias(LambertConformal2SP.LONGITUDE_OF_FALSE_ORIGIN, Citations.GEOTIFF, Mercator1SP.LONGITUDE_OF_ORIGIN, builder)
+ .rename(Citations.OGC, "longitude_of_center"));
+ /*
+ * EPSG: Latitude of 1st standard parallel
+ * OGC: standard_parallel_1
+ * ESRI: Standard_Parallel_1
+ * NetCDF: standard_parallel
+ * GeoTIFF: StdParallel1
+ *
+ * Special case: default value shall be the value of LATITUDE_OF_FALSE_ORIGIN.
+ */
+ STANDARD_PARALLEL_1 = LambertConformal2SP.STANDARD_PARALLEL_1;
+ /*
+ * EPSG: Latitude of 2nd standard parallel
+ * OGC: standard_parallel_2
+ * ESRI: Standard_Parallel_2
+ * NetCDF: standard_parallel
+ * GeoTIFF: StdParallel2
+ *
+ * Special case: default value shall be the value of STANDARD_PARALLEL_1.
+ */
+ STANDARD_PARALLEL_2 = LambertConformal2SP.STANDARD_PARALLEL_2;
+ /*
+ * EPSG: Easting at false origin
+ * OGC: false_easting
+ * ESRI: False_Easting
+ * NetCDF: false_easting
+ * GeoTIFF: FalseEasting
+ */
+ EASTING_AT_FALSE_ORIGIN = createShift(
+ renameAlias(LambertConformal2SP.EASTING_AT_FALSE_ORIGIN, Citations.GEOTIFF, LambertConformal1SP.FALSE_EASTING, builder));
+ /*
+ * EPSG: Northing at false origin
+ * OGC: false_northing
+ * ESRI: False_Northing
+ * NetCDF: false_northing
+ * GeoTIFF: FalseNorthing
+ */
+ NORTHING_AT_FALSE_ORIGIN = createShift(
+ renameAlias(LambertConformal2SP.NORTHING_AT_FALSE_ORIGIN, Citations.GEOTIFF, LambertConformal1SP.FALSE_NORTHING, builder));
+
+ PARAMETERS = builder
+ .addIdentifier( "9822")
+ .addName( "Albers Equal Area")
+ .addName(Citations.OGC, "Albers_Conic_Equal_Area")
+ .addName(Citations.ESRI, "Albers_Equal_Area_Conic")
+ .addName(Citations.ESRI, "Albers")
+ .addName(Citations.NETCDF, "AlbersEqualArea")
+ .addName(Citations.GEOTIFF, "CT_AlbersEqualArea")
+ .addName(Citations.PROJ4, "aea")
+ .addIdentifier(Citations.GEOTIFF, "11")
+ .addIdentifier(Citations.MAP_INFO, "9")
+ .addIdentifier(Citations.S57, "1")
+ .createGroupForMapProjection(
+ LATITUDE_OF_FALSE_ORIGIN,
+ LONGITUDE_OF_FALSE_ORIGIN,
+ STANDARD_PARALLEL_1,
+ STANDARD_PARALLEL_2,
+ EASTING_AT_FALSE_ORIGIN,
+ NORTHING_AT_FALSE_ORIGIN);
+ }
+
+ /**
+ * Constructs a new provider.
+ */
+ public AlbersEqualArea() {
+ super(PARAMETERS);
+ }
+
+ /**
+ * Returns the operation type for this map projection.
+ *
+ * @return {@code ConicProjection.class}
+ */
+ @Override
+ public final Class<ConicProjection> getOperationType() {
+ return ConicProjection.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.AlbersEqualArea(this, parameters);
+ }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
new file mode 100644
index 0000000..b2e601a
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -0,0 +1,343 @@
+/*
+ * 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.internal.referencing.Formulas;
+import org.apache.sis.internal.util.DoubleDouble;
+import org.apache.sis.measure.Latitude;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.util.Workaround;
+import org.apache.sis.util.resources.Errors;
+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 static java.lang.Math.*;
+import static org.apache.sis.internal.referencing.provider.AlbersEqualArea.*;
+
+
+/**
+ * <cite>Albers Equal Area</cite> projection (EPSG code 9822).
+ * See the <a href="http://mathworld.wolfram.com/AlbersEqual-AreaConicProjection.html">Albers Equal-Area
+ * Conic projection on MathWorld</a> for an overview.
+ *
+ * <p>The {@code "standard_parallel_2"} parameter is optional and will be given the same value as
+ * {@code "standard_parallel_1"} if not set.</p>
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @author Rémi Maréchal (Geomatys)
+ * @since 0.8
+ * @version 0.8
+ * @module
+ */
+public class AlbersEqualArea extends EqualAreaProjection {
+ /**
+ * For cross-version compatibility.
+ */
+ private static final long serialVersionUID = -3024658742514888646L;
+
+ /**
+ * Internal coefficients for computation, depending only on eccentricity and values of standards parallels.
+ * This is defined as {@literal n = (m₁² – m₂²) / (α₂ – α₁)} in §1.3.13 of IOGP Publication 373-7-2 (april 2015).
+ *
+ * <p>In Apache SIS implementation, we use modified formulas in which the (1 - ℯ²) factor is omitted in
+ * {@link #qm(double)} calculation. Consequently what we get is a modified value <var>nm</var> which is
+ * related to Synder's <var>n</var> value by {@literal n = nm / (1 - ℯ²)}. The omitted (1 - ℯ²) factor
+ * is either taken in account by the (de)normalization matrix, or cancels with other (1 - ℯ²) factors
+ * when we develop the formulas.</p>
+ *
+ * <p>Note that in the spherical case, <var>nm</var> = Synder's <var>n</var>.</p>
+ */
+ final double nm;
+
+ /**
+ * Internal coefficients for computation, depending only on values of standards parallels.
+ * This is defined as {@literal C = m₁² + (n⋅α₁)} in §1.3.13 of IOGP Publication 373-7-2 –
+ * Geomatics Guidance Note number 7, part 2 – April 2015.
+ */
+ final double C;
+
+ /**
+ * Creates an Albers 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 AlbersEqualArea(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.FALSE_EASTING, EASTING_AT_FALSE_ORIGIN);
+ roles.put(ParameterRole.FALSE_NORTHING, NORTHING_AT_FALSE_ORIGIN);
+ roles.put(ParameterRole.CENTRAL_MERIDIAN, LONGITUDE_OF_FALSE_ORIGIN);
+ return new Initializer(method, parameters, roles, (byte) 0);
+ }
+
+ /**
+ * Work around for RFE #4093999 in Sun's bug database
+ * ("Relax constraint on placement of this()/super() call in constructors").
+ */
+ @Workaround(library="JDK", version="1.7")
+ private AlbersEqualArea(final Initializer initializer) {
+ super(initializer);
+ double φ0 = initializer.getAndStore(LATITUDE_OF_FALSE_ORIGIN);
+ double φ1 = initializer.getAndStore(STANDARD_PARALLEL_1, φ0);
+ double φ2 = initializer.getAndStore(STANDARD_PARALLEL_2, φ1);
+ if (abs(φ1 + φ2) < Formulas.ANGULAR_TOLERANCE) {
+ throw new IllegalArgumentException(Errors.format(Errors.Keys.LatitudesAreOpposite_2,
+ new Latitude(φ1), new Latitude(φ2)));
+ }
+ final boolean secant = (abs(φ1 - φ2) >= Formulas.ANGULAR_TOLERANCE);
+ φ0 = toRadians(φ0);
+ φ1 = toRadians(φ1);
+ φ2 = toRadians(φ2);
+ final double sinφ0 = sin(φ0);
+ final double sinφ1 = sin(φ1);
+ final double cosφ1 = cos(φ1);
+ final double sinφ2 = sin(φ2);
+ final double cosφ2 = cos(φ2);
+ final double m1 = initializer.scaleAtφ(sinφ1, cosφ1); // = cos(φ₁) / √(1 – ℯ²sin²φ₁)
+ final double α1 = qm(sinφ1); // Omitted ×(1-ℯ²)
+ if (secant) {
+ final double m2 = initializer.scaleAtφ(sinφ2, cosφ2); // = cos(φ₂) / √(1 – ℯ²sin²φ₂)
+ final double α2 = qm(sinφ2); // Omitted ×(1-ℯ²)
+ nm = (m1*m1 - m2*m2) / (α2 - α1); // n = nm / (1-ℯ²)
+ } else {
+ nm = sinφ1;
+ }
+ C = m1*m1 + nm*α1; // Omitted (1-ℯ²) term in nm cancels with omitted (1-ℯ²) term in α₁.
+ /*
+ * Compute rn = (1-ℯ²)/nm, which is the reciprocal of the "real" n used in Synder and EPSG guidance note.
+ * We opportunistically use double-double arithmetic since the MatrixSIS operations use them anyway, but
+ * we do not really have that accuracy because of the limited precision of 'nm'. The intend is rather to
+ * increase the chances term cancellations happen during concatenation of coordinate operations.
+ */
+ final DoubleDouble rn = DoubleDouble.verbatim(1);
+ rn.subtract(initializer.eccentricitySquared);
+ rn.divide(nm, 0);
+ /*
+ * Compute ρ₀ = √(C - n⋅q(sinφ₀))/n with multiplication by a omitted because already taken in account
+ * by the denormalization matrix. Omitted (1-ℯ²) term in nm cancels with omitted (1-ℯ²) term in qm(…).
+ * See above note about double-double arithmetic usage.
+ */
+ final DoubleDouble ρ0 = DoubleDouble.verbatim(C - nm*qm(sinφ0));
+ ρ0.sqrt();
+ ρ0.multiply(rn);
+ /*
+ * At this point, all parameters have been processed. Now process to their
+ * validation and the initialization of (de)normalize affine transforms.
+ */
+ final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
+ final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
+ denormalize.convertBefore(0, rn, null); rn.negate();
+ denormalize.convertBefore(1, rn, ρ0); rn.inverseDivide(-1, 0);
+ normalize.convertAfter(0, rn, null);
+ super.computeCoefficients();
+ }
+
+ /**
+ * Creates a new projection initialized to the same parameters than the given one.
+ */
+ AlbersEqualArea(final AlbersEqualArea other) {
+ super(other);
+ nm = other.nm;
+ C = other.C;
+ }
+
+ /**
+ * Returns the names of additional internal parameters which need to be taken in account when
+ * comparing two {@code AlbersEqualArea} projections or formatting them in debug mode.
+ */
+ @Override
+ final String[] getInternalParameterNames() {
+ return new String[] {"n", "C"};
+ }
+
+ /**
+ * Returns the values of additional internal parameters which need to be taken in account when
+ * comparing two {@code AlbersEqualArea} projections or formatting them in debug mode.
+ */
+ @Override
+ final double[] getInternalParameterValues() {
+ return new double[] {nm / (1 - eccentricitySquared), C};
+ }
+
+ /**
+ * 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>.
+ *
+ * <p>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.</p>
+ *
+ * @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 {
+ AlbersEqualArea 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
+ {
+ final double θ = srcPts[srcOff ]; // θ = n⋅λ
+ final double φ = srcPts[srcOff+1];
+ final double cosθ = cos(θ);
+ final double sinθ = sin(θ);
+ final double sinφ = sin(φ);
+ final double ρ = sqrt(C - nm*qm_ellipsoid(sinφ));
+ if (dstPts != null) {
+ dstPts[dstOff ] = ρ * sinθ;
+ dstPts[dstOff+1] = ρ * cosθ;
+ }
+ if (!derivate) {
+ return null;
+ }
+ /*
+ * End of map projection. Now compute the derivative.
+ */
+ final double me = 1 - eccentricitySquared;
+ final double dρ_dφ = -0.5 * nm*dqm_dφ(sinφ, cos(φ)*me) / (me*ρ);
+ return new Matrix2(cosθ*ρ, dρ_dφ*sinθ, // ∂x/∂λ, ∂x/∂φ
+ -sinθ*ρ, dρ_dφ*cosθ); // ∂y/∂λ, ∂y/∂φ
+ }
+
+ /**
+ * 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
+ {
+ final double x = srcPts[srcOff ];
+ final double y = srcPts[srcOff+1];
+ dstPts[dstOff ] = atan2(x, y);
+ dstPts[dstOff+1] = φ((C - (x*x + y*y)) / nm);
+ /*
+ * Note: Synder 14-19 gives q = (C - ρ²n²/a²)/n where ρ = √(x² + (ρ₀ - y)²).
+ * But in Apache SIS implementation, ρ₀ has already been subtracted by the matrix before we reach this point.
+ * So we can simplify by ρ² = x² + y². Furthermore the matrix also divided x and y by a (the semi-major axis
+ * length) before this method, and multiplied by n. so what we have is actually (ρ⋅n/a)² = x² + y².
+ * So the formula become:
+ *
+ * q = (C - (x² + y²)) / n
+ *
+ * We divide by nm instead of n, so a (1-ℯ²) term is missing. But that missing term will be cancelled with
+ * the missing (1-ℯ²) term in qmPolar (the divisor applied by the φ(double) method that we invoke).
+ */
+ }
+
+
+ /**
+ * Provides the transform equations for the spherical case of the Albers Equal Area projection.
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @author Rémi Maréchal (Geomatys)
+ * @since 0.8
+ * @version 0.8
+ * @module
+ */
+ static final class Spherical extends AlbersEqualArea {
+ /**
+ * For cross-version compatibility.
+ */
+ private static final long serialVersionUID = 9090765015127854096L;
+
+ /**
+ * 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 AlbersEqualArea other) {
+ super(other);
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @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]; // θ = n⋅λ
+ final double φ = srcPts[srcOff+1];
+ final double cosθ = cos(θ);
+ final double sinθ = sin(θ);
+ final double sinφ = sin(φ);
+ final double ρ = sqrt(C - 2*nm*sinφ); // Synder 14-3 with radius and division by n omitted.
+ if (dstPts != null) {
+ dstPts[dstOff ] = ρ * sinθ; // Synder 14-1
+ dstPts[dstOff+1] = ρ * cosθ; // Synder 14-2
+ }
+ if (!derivate) {
+ return null;
+ }
+ final double dρ_dφ = -nm*cos(φ) / ρ;
+ return new Matrix2(cosθ*ρ, dρ_dφ*sinθ, // ∂x/∂λ, ∂x/∂φ
+ -sinθ*ρ, dρ_dφ*cosθ); // ∂y/∂λ, ∂y/∂φ
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @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];
+ dstPts[dstOff ] = atan2(x, y); // Part of Synder 14-11
+ dstPts[dstOff+1] = asin((C - (x*x + y*y)) / (nm*2)); // Synder 14-8 modified
+ }
+ }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
index 9c2d16b..ae74e73 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
@@ -97,13 +97,6 @@
private final byte variant;
/**
- * Value of {@link #qm(double)} function (part of Snyder equation (3-12)) at pole (sinφ = 1).
- *
- * @see #computeCoefficients()
- */
- private transient double qmPolar;
-
- /**
* Creates a Cylindrical Equal Area projection from the given parameters.
*
* @param method Description of the projection parameters.
@@ -176,16 +169,7 @@
ik.divide(k0);
denormalize.convertAfter(0, k0, null);
denormalize.convertAfter(1, ik, null);
- computeCoefficients();
- }
-
- /**
- * Invoked at construction time or on deserialization for computing the transient fields.
- */
- @Override
- final void computeCoefficients() {
super.computeCoefficients();
- qmPolar = qm(1);
}
/**
@@ -194,7 +178,6 @@
CylindricalEqualArea(final CylindricalEqualArea other) {
super(other);
variant = other.variant;
- qmPolar = other.qmPolar;
}
/**
@@ -235,8 +218,8 @@
final double φ = srcPts[srcOff+1];
final double sinφ = sin(φ);
if (dstPts != null) {
- dstPts[dstOff ] = srcPts[srcOff]; // Multiplication by k₀ will be applied by the denormalization matrix.
- dstPts[dstOff+1] = qm(sinφ); // Multiplication by (1-ℯ²)/(2k₀) will be applied by the denormalization matrix.
+ dstPts[dstOff ] = srcPts[srcOff]; // Multiplication by k₀ will be applied by the denormalization matrix.
+ dstPts[dstOff+1] = qm_ellipsoid(sinφ); // Multiplication by (1-ℯ²)/(2k₀) will be applied by the denormalization matrix.
}
/*
* End of map projection. Now compute the derivative, if requested.
@@ -265,7 +248,7 @@
dstOff--;
while (--numPts >= 0) {
final double φ = dstPts[dstOff += 2]; // Same as srcPts[srcOff + 1].
- dstPts[dstOff] = qm(sin(φ)); // Part of Synder equation (10-15)
+ dstPts[dstOff] = qm_ellipsoid(sin(φ)); // Part of Synder equation (10-15)
}
}
}
@@ -283,7 +266,7 @@
{
final double y = srcPts[srcOff+1]; // Must be before writing x.
dstPts[dstOff ] = srcPts[srcOff ]; // Must be before writing y.
- dstPts[dstOff+1] = φ(y / qmPolar);
+ dstPts[dstOff+1] = φ(y);
/*
* Equation 10-26 of Synder gives β = asin(2y⋅k₀/(a⋅qPolar)).
* In our case it simplifies to sinβ = (y/qmPolar) because:
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 d6a101c..0c17bf6 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
@@ -18,6 +18,7 @@
import java.io.IOException;
import java.io.ObjectInputStream;
+import org.apache.sis.util.resources.Errors;
import static java.lang.Math.*;
import static org.apache.sis.math.MathFunctions.atanh;
@@ -70,6 +71,13 @@
private transient double ci2, ci4, ci8;
/**
+ * Value of {@link #qm(double)} function (part of Snyder equation (3-12)) at pole (sinφ = 1).
+ *
+ * @see #computeCoefficients()
+ */
+ private transient double qmPolar;
+
+ /**
* Creates a new normalized projection from the parameters computed by the given initializer.
*
* @param initializer The initializer for computing map projection internal parameters.
@@ -84,8 +92,8 @@
*/
void computeCoefficients() {
final double e2 = eccentricitySquared;
- final double e4 = e2 * e2;
- final double e6 = e2 * e4;
+ final double e4 = e2 * e2;
+ final double e6 = e2 * e4;
ci2 = 517/5040. * e6 + 31/180. * e4 + 1/3. * e2;
ci4 = 251/3780. * e6 + 23/360. * e4;
ci8 = 761/45360. * e6;
@@ -100,6 +108,7 @@
ci4 *= 8;
ci8 *= 64;
}
+ qmPolar = qm(1);
}
/**
@@ -113,6 +122,7 @@
ci2 = other.ci2;
ci4 = other.ci4;
ci8 = other.ci8;
+ qmPolar = other.qmPolar;
}
/**
@@ -120,6 +130,8 @@
* In order to get the <var>q</var> function, this method output must be multiplied
* by <code>(1 - {@linkplain #eccentricitySquared})</code>.
*
+ * <p>The <var>q</var> variable is named <var>α</var> in EPSG guidance notes.</p>
+ *
* <p>This equation has the following properties:</p>
*
* <ul>
@@ -130,13 +142,28 @@
* <li>q(0) = 0</li>
* </ul>
*
- * In the spherical case, <var>q</var> = 2⋅sinφ. It is caller responsibility to ensure that this
- * method is not invoked in the spherical case, since this implementation does not work in such case.
+ * In the spherical case, <var>q</var> = 2⋅sinφ.
*
* @param sinφ sine of the latitude <var>q</var> is calculated for.
* @return <var>q</var> from Snyder equation (3-12).
*/
final double qm(final double sinφ) {
+ /*
+ * 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φ);
+ }
+
+ /**
+ * Same as {@link #qm(double)} but without check about whether the map projection is a spherical case.
+ * It is caller responsibility to ensure that this method is not invoked in the spherical case, since
+ * this implementation does not work in such case.
+ *
+ * @param sinφ sine of the latitude <var>q</var> is calculated for.
+ * @return <var>q</var> from Snyder equation (3-12).
+ */
+ final double qm_ellipsoid(final double sinφ) {
final double ℯsinφ = eccentricity * sinφ;
return sinφ / (1 - ℯsinφ*ℯsinφ) + atanh(ℯsinφ) / eccentricity;
}
@@ -151,23 +178,29 @@
* @return the {@code qm} derivative at the specified latitude.
*/
final double dqm_dφ(final double sinφ, final double cosφ) {
- final double ℯsinφ2 = eccentricitySquared * (sinφ*sinφ);
- return (cosφ / (1 - ℯsinφ2)) * (1 + ((1 + ℯsinφ2) / (1 - ℯsinφ2)));
+ final double t = 1 - eccentricitySquared*(sinφ*sinφ);
+ return 2*cosφ / (t*t);
}
/**
- * Computes the latitude using equation 3-18 from Synder.
+ * Computes the latitude using equation 3-18 from Synder, followed by iterative resolution of Synder 3-16.
+ * If theory, the series expansion given by equation 3-18 (φ ≈ c₂⋅sin(2β) + c₄⋅sin(4β) + c₈⋅sin(8β)) should
+ * be used in replacement of the iterative method. However in practice the series expansion seems to not
+ * 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.
*
- * @param sinβ see Synder equation 10-26.
+ * @param y in the cylindrical case, this is northing on the normalized ellipsoid.
* @return the latitude in radians.
*/
- final double φ(final double sinβ) {
+ final double φ(final double y) throws ProjectionException {
+ final double sinβ = y / qmPolar;
final double β = asin(sinβ);
+ double φ;
if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
- return ci8 * sin(8*β)
- + ci4 * sin(4*β)
- + ci2 * sin(2*β)
- + β; // Synder 3-18
+ φ = ci8 * sin(8*β)
+ + ci4 * sin(4*β)
+ + ci2 * sin(2*β)
+ + β; // Synder 3-18
} else {
/*
* Same formula than above, but rewriten using trigonometric identities in order to avoid
@@ -183,8 +216,33 @@
assert identityEquals(t4β, sin(4*β) / ( 8 * t2β));
assert identityEquals(t8β, sin(8*β) / (64 * t2β));
- return (ci8*t8β + ci4*t4β + ci2) * t2β + β;
+ φ = (ci8*t8β + ci4*t4β + ci2) * t2β + β;
}
+ /*
+ * At this point φ is close to the desired value, but may have an error of a few centimetres.
+ * Use the iterative method for reaching the last part of missing accuracy. Usually this loop
+ * will perform exactly one iteration, no more, because φ is already quite close to the result.
+ *
+ * Mathematical note: Synder 3-16 gives q/(1-ℯ²) instead of y in the calculation of Δφ below.
+ * For Cylindrical Equal Area projection, Synder 10-17 gives q = (qPolar⋅sinβ), which simplifies
+ * as y.
+ *
+ * For Albers Equal Area projection, Synder 14-19 gives q = (C - ρ²n²/a²)/n, which we rewrite
+ * as q = (C - ρ²)/n (see comment in AlbersEqualArea.inverseTransform(…) for the mathematic).
+ * The y value given to this method is y = (C - ρ²) / (n⋅(1-ℯ²)) = q/(1-ℯ²), the desired value.
+ */
+ for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
+ final double sinφ = sin(φ);
+ final double cosφ = cos(φ);
+ final double ℯsinφ = eccentricity * sinφ;
+ final double ome = 1 - ℯsinφ*ℯsinφ;
+ final double Δφ = ome*ome/(2*cosφ) * (y - sinφ/ome - atanh(ℯsinφ)/eccentricity);
+ φ += Δφ;
+ if (abs(Δφ) <= ITERATION_TOLERANCE) {
+ return φ;
+ }
+ }
+ throw new ProjectionException(Errors.format(Errors.Keys.NoConvergence));
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
index 237afba..ec356e1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
@@ -295,7 +295,7 @@
*
* <blockquote>ν = 1 / √(1 - ℯ²⋅sin²φ)</blockquote>
*
- * This method returns 1/ν².
+ * This method returns 1/ν², which is the (1 - ℯ²⋅sin²φ) part of above equation.
*
* <div class="section">Relationship with Snyder</div>
* This is related to functions (14-15) from Snyder (used for computation of scale factors
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
index a289459..4d6b8af 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
@@ -117,7 +117,7 @@
}
/**
- * Internal coefficients for computation, depending only on values of standards parallels.
+ * Internal coefficients for computation, depending only on eccentricity and values of standards parallels.
* This is defined as {@literal n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂)} in §1.3.1.1 of
* IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015.
*
@@ -150,8 +150,8 @@
* <li><cite>"Lambert Conic Conformal (2SP Michigan)"</cite>.</li>
* </ul>
*
- * @param method Description of the projection parameters.
- * @param parameters The parameter values of the projection to create.
+ * @param method description of the projection parameters.
+ * @param parameters the parameter values of the projection to create.
*/
public LambertConicConformal(final OperationMethod method, final Parameters parameters) {
this(initializer(method, parameters));
@@ -274,11 +274,11 @@
* for reducing the amount of calls to the logarithmic function. Note that this equation
* tends toward 0/0 if φ₁ ≈ φ₂, which force us to do a special check for the SP1 case.
*/
- if (abs(φ1 - φ2) >= ANGULAR_TOLERANCE) { // Should be 'true' for 2SP case.
+ if (abs(φ1 - φ2) >= ANGULAR_TOLERANCE) { // Should be 'true' for 2SP case.
final double sinφ2 = sin(φ2);
final double m2 = initializer.scaleAtφ(sinφ2, cos(φ2));
final double t2 = expOfNorthing(φ2, eccentricity*sinφ2);
- n = log(m1/m2) / log(t1/t2); // Tend toward 0/0 if φ1 ≈ φ2.
+ n = log(m1/m2) / log(t1/t2); // Tend toward 0/0 if φ1 ≈ φ2.
} else {
n = -sinφ1;
}
@@ -375,8 +375,8 @@
* <p>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.</p>
*
- * @param factory The factory to use for creating the transform.
- * @return The map projection from (λ,φ) to (<var>x</var>,<var>y</var>) coordinates.
+ * @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
@@ -560,7 +560,7 @@
double x = srcPts[srcOff ];
double y = srcPts[srcOff+1];
final double ρ = hypot(x, y);
- dstPts[dstOff ] = atan2(x, y); // Really (x,y), not (y,x);
+ dstPts[dstOff ] = atan2(x, y); // Really (x,y), not (y,x);
dstPts[dstOff+1] = PI/2 - 2*atan(pow(1/ρ, 1/n));
}
}
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 e93d3d7..3a0f820 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
@@ -27,13 +27,14 @@
org.apache.sis.internal.referencing.provider.PseudoMercator
org.apache.sis.internal.referencing.provider.RegionalMercator
org.apache.sis.internal.referencing.provider.MillerCylindrical
-org.apache.sis.internal.referencing.provider.LambertCylindricalEqualArea
-org.apache.sis.internal.referencing.provider.LambertCylindricalEqualAreaSpherical
org.apache.sis.internal.referencing.provider.LambertConformal1SP
org.apache.sis.internal.referencing.provider.LambertConformal2SP
org.apache.sis.internal.referencing.provider.LambertConformalWest
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.AlbersEqualArea
org.apache.sis.internal.referencing.provider.TransverseMercator
org.apache.sis.internal.referencing.provider.TransverseMercatorSouth
org.apache.sis.internal.referencing.provider.PolarStereographicA
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 cbe2c51..31e571b 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
@@ -79,13 +79,14 @@
PseudoMercator.class,
RegionalMercator.class,
MillerCylindrical.class,
- LambertCylindricalEqualArea.class,
- LambertCylindricalEqualAreaSpherical.class,
LambertConformal1SP.class,
LambertConformal2SP.class,
LambertConformalWest.class,
LambertConformalBelgium.class,
LambertConformalMichigan.class,
+ LambertCylindricalEqualArea.class,
+ LambertCylindricalEqualAreaSpherical.class,
+ AlbersEqualArea.class,
TransverseMercator.class,
TransverseMercatorSouth.class,
PolarStereographicA.class,
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
new file mode 100644
index 0000000..75265e2
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
@@ -0,0 +1,178 @@
+/*
+ * 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.MathTransform;
+import org.opengis.referencing.operation.MathTransformFactory;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.internal.referencing.provider.MapProjection;
+import org.apache.sis.internal.system.DefaultFactories;
+import org.apache.sis.parameter.Parameters;
+
+import static java.lang.StrictMath.*;
+
+// Test dependencies
+import org.opengis.test.ToleranceModifier;
+import org.apache.sis.test.DependsOnMethod;
+import org.apache.sis.test.DependsOn;
+import org.apache.sis.test.TestUtilities;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+
+/**
+ * Tests the {@link AlbersEqualArea} class. We test using various values of standard parallels.
+ * We do not test with various values of the latitude of origin, because its only effect is to
+ * modify the translation term on the <var>y</var> axis.
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @author Rémi Maréchal (Geomatys)
+ * @since 0.8
+ * @version 0.8
+ * @module
+ */
+@DependsOn(CylindricalEqualAreaTest.class)
+public final strictfp class AlbersEqualAreaTest extends MapProjectionTestCase {
+ /**
+ * Creates a new map projection. See the class javadoc for an explanation about
+ * why we ask only for the standard parallels and not the latitude of origin.
+ *
+ * @param a semi-major axis length.
+ * @param b semi-minor axis length.
+ * @param φ1 first standard parallel.
+ * @param φ2 second standard parallel.
+ * @return newly created projection.
+ * @throws FactoryException if an error occurred while creating the map projection.
+ */
+ private static MathTransform create(final double a, final double b, double φ1, double φ2) throws FactoryException {
+ final MapProjection method = new org.apache.sis.internal.referencing.provider.AlbersEqualArea();
+ final Parameters values = Parameters.castOrWrap(method.getParameters().createValue());
+ values.parameter("semi_major").setValue(a);
+ values.parameter("semi_minor").setValue(b);
+ values.parameter("standard_parallel_1").setValue(φ1);
+ values.parameter("standard_parallel_2").setValue(φ2);
+ return method.createMathTransform(DefaultFactories.forBuildin(MathTransformFactory.class), values);
+ }
+
+ /**
+ * Returns whether the given projection is the spherical implementation.
+ */
+ private static boolean isSpherical(final AlbersEqualArea transform) {
+ return transform instanceof AlbersEqualArea.Spherical;
+ }
+
+ /**
+ * Tests the unitary projection on a sphere.
+ *
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a point.
+ */
+ @Test
+ public void testSphere() throws FactoryException, TransformException {
+ final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres.
+ derivativeDeltas = new double[] {delta, delta};
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ transform = create(6370997, 6370997, 29.5, 45.5); // Standard parallels from Synder table 15.
+ final AlbersEqualArea kernel = (AlbersEqualArea) getKernel();
+ assertTrue("isSpherical", isSpherical(kernel));
+ assertEquals("n", 0.6028370, kernel.nm, 0.5E-7); // Expected 'n' value from Synder table 15.
+ /*
+ * When stepping into the AlbersEqualArea.Sphere.transform(…) method with a debugger, the
+ * expected value of 6370997*ρ/n is 6910941 (value taken from ρ column in Synder table 15).
+ */
+ verifyTransform(new double[] {0, 50}, new double[] {0, 5373933.180});
+ /*
+ * Expect 6370997*ρ/n ≈ 8022413 (can be verified only with the debugger)
+ */
+ verifyTransform(new double[] {0, 40}, new double[] {0, 4262461.266});
+ /*
+ * Expect 6370997*ρ/n ≈ 9695749 (can be verified only with the debugger)
+ */
+ verifyTransform(new double[] {0, 25}, new double[] {0, 2589125.654});
+ /*
+ * Verify consistency with random points.
+ */
+ verifyInDomain(new double[] {-20, 20}, new double[] {20, 50}, new int[] {5, 5},
+ TestUtilities.createRandomNumberGenerator());
+ }
+
+ /**
+ * Tests the unitary projection on an ellipse.
+ *
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a point.
+ */
+ @Test
+ @DependsOnMethod("testSphere")
+ public void testEllipse() throws FactoryException, TransformException {
+ final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres.
+ derivativeDeltas = new double[] {delta, delta};
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ transform = create(6378206.4, 6356583.8, 29.5, 45.5); // Standard parallels from Synder table 15.
+ final AlbersEqualArea kernel = (AlbersEqualArea) getKernel();
+ assertFalse("isSpherical", isSpherical(kernel));
+ /*
+ * Expected 'n' value from Synder table 15. The division by (1-ℯ²) is because Apache SIS omits this factor
+ * in its calculation of n (we rather take it in account in (de)normalization matrices and elsewhere).
+ */
+ assertEquals("n", 0.6029035, kernel.nm / (1 - kernel.eccentricitySquared), 0.5E-7);
+ /*
+ * When stepping into the AlbersEqualArea.Sphere.transform(…) method with a debugger, the expected
+ * value of 6378206.4*ρ/(nm/(1-ℯ²)) is 6931335 (value taken from ρ column in Synder table 15).
+ */
+ verifyTransform(new double[] {0, 50}, new double[] {0, 5356698.435});
+ /*
+ * Expect 6378206.4*ρ/(nm/(1-ℯ²)) ≈ 8042164 (can be verified only with the debugger)
+ */
+ verifyTransform(new double[] {0, 40}, new double[] {0, 4245869.390});
+ /*
+ * Expect 6378206.4*ρ/(nm/(1-ℯ²)) ≈ 9710969 (can be verified only with the debugger)
+ */
+ verifyTransform(new double[] {0, 25}, new double[] {0, 2577064.350});
+ /*
+ * Verify consistency with random points.
+ */
+ verifyInDomain(new double[] {-20, 20}, new double[] {20, 50}, new int[] {5, 5},
+ TestUtilities.createRandomNumberGenerator());
+ }
+
+ /**
+ * Uses Proj.4 test point has a reference.
+ *
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a point.
+ */
+ @Test
+ @DependsOnMethod("testEllipse")
+ public void compareWithProj4() throws FactoryException, TransformException {
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ tolerance = Formulas.LINEAR_TOLERANCE;
+
+ // Spherical case
+ transform = create(6400000, 6400000, 0, 2);
+ verifyTransform(new double[] {2, 1}, new double[] {223334.085, 111780.432});
+
+ // Ellipsoidal case
+ transform = create(6378137, 6356752.314140347, 0, 2);
+ verifyTransform(new double[] {2, 1}, new double[] {222571.609, 110653.327});
+ }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
index ed95188..8ee3fbb 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
@@ -19,6 +19,7 @@
import javax.measure.unit.NonSI;
import org.opengis.util.FactoryException;
import org.opengis.referencing.datum.Ellipsoid;
+import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import org.opengis.test.referencing.ParameterizedTransformTest;
import org.apache.sis.parameter.Parameters;
@@ -28,6 +29,7 @@
import org.apache.sis.referencing.operation.transform.CoordinateDomain;
import org.apache.sis.referencing.operation.transform.MathTransformTestCase;
import org.apache.sis.referencing.operation.transform.MathTransformFactoryMock;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.datum.GeodeticDatumMock;
import static java.lang.StrictMath.*;
@@ -39,7 +41,7 @@
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.6
- * @version 0.6
+ * @version 0.8
* @module
*/
strictfp class MapProjectionTestCase extends MathTransformTestCase {
@@ -57,9 +59,9 @@
/**
* Returns the parameters to use for instantiating the projection to test.
*
- * @param provider The provider of the projection to test.
- * @param ellipse {@code false} for a sphere, or {@code true} for WGS84 ellipsoid.
- * @return The parameters to use for instantiating the projection.
+ * @param provider the provider of the projection to test.
+ * @param ellipse {@code false} for a sphere, or {@code true} for WGS84 ellipsoid.
+ * @return the parameters to use for instantiating the projection.
*/
static Parameters parameters(final DefaultOperationMethod provider, final boolean ellipse) {
final Parameters parameters = Parameters.castOrWrap(provider.getParameters().createValue());
@@ -75,8 +77,8 @@
/**
* Instantiates the object to use for running GeoAPI test.
*
- * @param provider The provider of the projection to test.
- * @return The GeoAPI test class using the given provider.
+ * @param provider the provider of the projection to test.
+ * @return the GeoAPI test class using the given provider.
*/
static ParameterizedTransformTest createGeoApiTest(final MapProjection provider) {
return new ParameterizedTransformTest(new MathTransformFactoryMock(provider));
@@ -107,11 +109,26 @@
}
/**
+ * Returns the {@code NormalizedProjection} component of the current transform.
+ */
+ final NormalizedProjection getKernel() {
+ NormalizedProjection kernel = null;
+ for (final MathTransform component : MathTransforms.getSteps(transform)) {
+ if (component instanceof NormalizedProjection) {
+ assertNull("Found more than one kernel.", kernel);
+ kernel = (NormalizedProjection) component;
+ }
+ }
+ assertNotNull("Kernel not found.", kernel);
+ return kernel;
+ }
+
+ /**
* Projects the given latitude value. The longitude is fixed to zero.
* This method is useful for testing the behavior close to poles in a simple case.
*
- * @param φ The latitude.
- * @return The northing.
+ * @param φ the latitude.
+ * @return the northing.
* @throws ProjectionException if the projection failed.
*/
final double transform(final double φ) throws ProjectionException {
@@ -129,8 +146,8 @@
* Inverse projects the given northing value. The easting is fixed to zero.
* This method is useful for testing the behavior close to poles in a simple case.
*
- * @param y The northing.
- * @return The latitude.
+ * @param y the northing.
+ * @return the latitude.
* @throws ProjectionException if the projection failed.
*/
final double inverseTransform(final double y) throws ProjectionException {
@@ -159,9 +176,9 @@
* The spherical formulas are arbitrarily selected as the reference implementation because they are simpler,
* so less bug-prone, than the elliptical formulas.
*
- * @param domain The domain of the numbers to be generated.
- * @param randomSeed The seed for the random number generator, or 0 for choosing a random seed.
- * @throws TransformException If a conversion, transformation or derivative failed.
+ * @param domain the domain of the numbers to be generated.
+ * @param randomSeed the seed for the random number generator, or 0 for choosing a random seed.
+ * @throws TransformException if a conversion, transformation or derivative failed.
*/
final void compareEllipticalWithSpherical(final CoordinateDomain domain, final long randomSeed)
throws TransformException
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 3a5fdc9..081e003 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
@@ -165,6 +165,7 @@
org.apache.sis.referencing.operation.projection.PolarStereographicTest.class,
org.apache.sis.referencing.operation.projection.ObliqueStereographicTest.class,
org.apache.sis.referencing.operation.projection.CylindricalEqualAreaTest.class,
+ org.apache.sis.referencing.operation.projection.AlbersEqualAreaTest.class,
// Coordinate operation and derived Coordinate Reference Systems (cyclic dependency).
org.apache.sis.referencing.operation.DefaultTransformationTest.class,
