Merge the addition of AlbersEqualArea projection from JDK7 branch.
git-svn-id: https://svn.apache.org/repos/asf/sis/branches/JDK6@1754380 13f79535-47bb-0310-9956-ffa450edef68
diff --git a/core/sis-feature/src/main/java/org/apache/sis/feature/FeatureFormat.java b/core/sis-feature/src/main/java/org/apache/sis/feature/FeatureFormat.java
index 972b077..7e7f3ab 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/feature/FeatureFormat.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/feature/FeatureFormat.java
@@ -32,6 +32,7 @@
import org.opengis.util.GenericName;
import org.apache.sis.io.TableAppender;
import org.apache.sis.io.TabularFormat;
+import org.apache.sis.util.CharSequences;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.resources.Vocabulary;
@@ -47,7 +48,6 @@
import org.opengis.feature.FeatureType;
import org.opengis.feature.FeatureAssociationRole;
import org.opengis.feature.Operation;
-import org.apache.sis.util.CharSequences;
/**
diff --git a/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AssociationRoleBuilder.java b/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AssociationRoleBuilder.java
index 15aa566..6a9ba54 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AssociationRoleBuilder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AssociationRoleBuilder.java
@@ -18,6 +18,8 @@
import org.opengis.util.GenericName;
import org.apache.sis.feature.DefaultAssociationRole;
+
+// Branch-dependent imports
import org.opengis.feature.FeatureType;
import org.opengis.feature.PropertyType;
import org.opengis.feature.FeatureAssociationRole;
diff --git a/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AttributeTypeBuilder.java b/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AttributeTypeBuilder.java
index 891dcf0..e5d1abf 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AttributeTypeBuilder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/feature/builder/AttributeTypeBuilder.java
@@ -32,7 +32,6 @@
// Branch-dependent imports
import org.apache.sis.internal.jdk7.Objects;
import org.opengis.feature.AttributeType;
-import org.opengis.feature.Feature;
import org.opengis.feature.PropertyType;
@@ -254,9 +253,6 @@
* @see #characteristics()
*/
public <C> CharacteristicTypeBuilder<C> addCharacteristic(final Class<C> type) {
- if (valueClass == Feature.class) {
- throw new UnsupportedOperationException(errors().getString(Errors.Keys.IllegalOperationForValueClass_1, valueClass));
- }
ensureNonNull("type", type);
final CharacteristicTypeBuilder<C> characteristic = new CharacteristicTypeBuilder<C>(this, type);
characteristics.add(characteristic);
diff --git a/core/sis-feature/src/main/java/org/apache/sis/feature/builder/PropertyTypeBuilder.java b/core/sis-feature/src/main/java/org/apache/sis/feature/builder/PropertyTypeBuilder.java
index da3eed1..acd4eb4 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/feature/builder/PropertyTypeBuilder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/feature/builder/PropertyTypeBuilder.java
@@ -18,6 +18,8 @@
import org.opengis.util.GenericName;
import org.apache.sis.util.resources.Errors;
+
+// Branch-dependent imports
import org.opengis.feature.AttributeType;
import org.opengis.feature.FeatureType;
import org.opengis.feature.PropertyType;
diff --git a/core/sis-feature/src/test/java/org/apache/sis/feature/builder/FeatureTypeBuilderTest.java b/core/sis-feature/src/test/java/org/apache/sis/feature/builder/FeatureTypeBuilderTest.java
index 3dd944b..16e1d5a 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/feature/builder/FeatureTypeBuilderTest.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/feature/builder/FeatureTypeBuilderTest.java
@@ -19,8 +19,11 @@
import java.util.Iterator;
import com.esri.core.geometry.Geometry;
import com.esri.core.geometry.Point;
+import org.apache.sis.internal.feature.AttributeConvention;
+import org.apache.sis.feature.DefaultFeatureTypeTest;
import org.apache.sis.referencing.crs.HardCodedCRS;
import org.apache.sis.test.DependsOnMethod;
+import org.apache.sis.test.TestUtilities;
import org.apache.sis.test.TestCase;
import org.junit.Test;
@@ -30,9 +33,6 @@
import org.opengis.feature.AttributeType;
import org.opengis.feature.FeatureType;
import org.opengis.feature.PropertyType;
-import org.apache.sis.feature.DefaultFeatureTypeTest;
-import org.apache.sis.internal.feature.AttributeConvention;
-import org.apache.sis.test.TestUtilities;
/**
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/CRS.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/CRS.java
index 7258591..f7323e1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/CRS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/CRS.java
@@ -39,6 +39,7 @@
import org.opengis.referencing.crs.TemporalCRS;
import org.opengis.referencing.crs.VerticalCRS;
import org.opengis.referencing.crs.EngineeringCRS;
+import org.opengis.referencing.operation.OperationNotFoundException;
import org.opengis.metadata.citation.Citation;
import org.opengis.metadata.extent.Extent;
import org.opengis.metadata.extent.GeographicBoundingBox;
@@ -282,11 +283,17 @@
* for checking if the operation has sufficient accuracy for caller's purpose.</li>
* </ul>
*
+ * If the source and target CRS are equivalent, then this method returns an operation backed by an
+ * {@linkplain org.apache.sis.referencing.operation.transform.AbstractMathTransform#isIdentity() identity}
+ * transform. If there is no known operation between the given pair of CRS, then this method throws an
+ * {@link OperationNotFoundException}.
+ *
* @param sourceCRS the CRS of source coordinates.
* @param targetCRS the CRS of target coordinates.
* @param areaOfInterest the area of interest, or {@code null} if none.
* @return the mathematical operation from {@code sourceCRS} to {@code targetCRS}.
- * @throws FactoryException if the operation can not be created.
+ * @throws OperationNotFoundException if no operation was found between the given pair of CRS.
+ * @throws FactoryException if the operation can not be created for another reason.
*
* @see DefaultCoordinateOperationFactory#createOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, CoordinateOperationContext)
*
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
index 2fe6697..60d88e5 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
@@ -356,8 +356,12 @@
final CoordinateOperation step1;
try {
step1 = inverse(sourceCRS.getConversionFromBase());
+ } catch (OperationNotFoundException exception) {
+ throw exception;
+ } catch (FactoryException exception) {
+ throw new OperationNotFoundException(canNotInvert(sourceCRS), exception);
} catch (NoninvertibleTransformException exception) {
- throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS), exception);
+ throw new OperationNotFoundException(canNotInvert(sourceCRS), exception);
}
return concatenate(step1, step2);
}
@@ -387,8 +391,12 @@
final CoordinateOperation step1;
try {
step1 = inverse(sourceCRS.getConversionFromBase());
+ } catch (OperationNotFoundException exception) {
+ throw exception;
+ } catch (FactoryException exception) {
+ throw new OperationNotFoundException(canNotInvert(sourceCRS), exception);
} catch (NoninvertibleTransformException exception) {
- throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS), exception);
+ throw new OperationNotFoundException(canNotInvert(sourceCRS), exception);
}
return concatenate(step1, step2, step3);
}
@@ -1064,4 +1072,15 @@
private static String notFoundMessage(final IdentifiedObject source, final IdentifiedObject target) {
return Errors.format(Errors.Keys.CoordinateOperationNotFound_2, CRSPair.label(source), CRSPair.label(target));
}
+
+ /**
+ * Returns an error message for "Can not invert operation XYZ.".
+ * This is used for the construction of {@link OperationNotFoundException}.
+ *
+ * @param crs the CRS having a conversion that can not be inverted.
+ * @return A default error message.
+ */
+ private static String canNotInvert(final GeneralDerivedCRS crs) {
+ return Errors.format(Errors.Keys.NonInvertibleOperation_1, crs.getConversionFromBase().getName().getCode());
+ }
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java
index 965b25a..0276e37 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java
@@ -500,24 +500,52 @@
/**
* Creates the inverse of the given single operation.
+ * If this operation succeed, then the returned coordinate operations has the following properties:
+ *
+ * <ul>
+ * <li>Its {@code sourceCRS} is the {@code targetCRS} of the given operation.</li>
+ * <li>Its {@code targetCRS} is the {@code sourceCRS} of the given operation.</li>
+ * <li>Its {@code interpolationCRS} is {@code null}.</li>
+ * <li>Its {@code MathTransform} is the
+ * {@linkplain org.apache.sis.referencing.operation.transform.AbstractMathTransform#inverse() inverse}
+ * of the {@code MathTransform} of this operation.</li>
+ * <li>Its domain of validity and accuracy is the same.</li>
+ * </ul>
+ *
+ * <div class="note"><b>Note:</b>
+ * in many cases, the inverse operation is numerically less accurate than the direct operation because it
+ * uses approximations like series expansions or iterative methods. However the numerical errors caused by
+ * those approximations are not of interest here, because they are usually much smaller than the inaccuracy
+ * due to the stochastic nature of coordinate transformations (not to be confused with coordinate conversions;
+ * see ISO 19111 for more information).</div>
*/
final CoordinateOperation inverse(final SingleOperation op) throws NoninvertibleTransformException, FactoryException {
final CoordinateReferenceSystem sourceCRS = op.getSourceCRS();
final CoordinateReferenceSystem targetCRS = op.getTargetCRS();
final MathTransform transform = op.getMathTransform().inverse();
+ final OperationMethod method = InverseOperationMethod.create(op.getMethod());
+ final Map<String,Object> properties = properties(INVERSE_OPERATION);
+ InverseOperationMethod.properties(op, properties);
+ /*
+ * Find a hint about whether the coordinate operation is a transformation or a conversion,
+ * but do not set any conversion subtype. In particular, do not specify a Projection type,
+ * because the inverse of a Projection does not implement the Projection interface.
+ */
Class<? extends CoordinateOperation> type = null;
if (op instanceof Transformation) type = Transformation.class;
else if (op instanceof Conversion) type = Conversion.class;
- final Map<String,Object> properties = properties(INVERSE_OPERATION);
- InverseOperationMethod.putMetadata(op, properties);
- InverseOperationMethod.putParameters(op, properties);
- return createFromMathTransform(properties, targetCRS, sourceCRS,
- transform, InverseOperationMethod.create(op.getMethod()), null, type);
+ return createFromMathTransform(properties, targetCRS, sourceCRS, transform, method, null, type);
}
/**
* Creates the inverse of the given operation, which may be single or compound.
*
+ * <p><b>Design note:</b>
+ * we do not provide a {@code AbstractCoordinateOperation.inverse()} method. If the user wants an inverse method,
+ * he should invoke {@code CRS.findOperation(targetCRS, sourceCRS, null)} or something equivalent. This is because
+ * a new query of EPSG database may be necessary, and if no explicit definition is found there is too many arbitrary
+ * values to set in a default inverse operation for making that API public.</p>
+ *
* @param operation The operation to invert, or {@code null}.
* @return The inverse of {@code operation}, or {@code null} if none.
* @throws NoninvertibleTransformException if the operation is not invertible.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/InverseOperationMethod.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/InverseOperationMethod.java
index 21af1bd..09b865e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/InverseOperationMethod.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/InverseOperationMethod.java
@@ -45,7 +45,7 @@
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.7
- * @version 0.7
+ * @version 0.8
* @module
*/
@XmlTransient
@@ -107,7 +107,15 @@
}
/**
- * Copies accuracy and domain of validity metadata from the given operation into the given properties map.
+ * Infers the properties to give to an inverse coordinate operation.
+ * The returned map will contains three kind of information:
+ *
+ * <ul>
+ * <li>Metadata (domain of validity, accuracy)</li>
+ * <li>Parameter values, if possible</li>
+ * </ul>
+ *
+ * This method copies accuracy and domain of validity metadata from the given operation.
* We presume that the inverse operation has the same accuracy than the direct operation.
*
* <div class="note"><b>Note:</b>
@@ -116,25 +124,25 @@
* those approximations are not of interest here, because they are usually much smaller than the inaccuracy
* due to the stochastic nature of coordinate transformations (not to be confused with coordinate conversions;
* see ISO 19111 for more information).</div>
+ *
+ * If the inverse of the given operation can be represented by inverting the sign of all numerical
+ * parameter values, then this method copies also those parameters in a {@code "parameters"} entry.
+ *
+ * @param source the operation for which to get the inverse parameters.
+ * @param target where to store the inverse parameters.
*/
- static void putMetadata(final SingleOperation source, final Map<String,Object> target) {
+ static void properties(final SingleOperation source, final Map<String,Object> target) {
target.put(SingleOperation.DOMAIN_OF_VALIDITY_KEY, source.getDomainOfValidity());
final Collection<PositionalAccuracy> accuracy = source.getCoordinateOperationAccuracy();
if (!Containers.isNullOrEmpty(accuracy)) {
target.put(SingleOperation.COORDINATE_OPERATION_ACCURACY_KEY,
accuracy.toArray(new PositionalAccuracy[accuracy.size()]));
}
- }
-
- /**
- * If the inverse of the given operation can be represented by inverting the sign of all numerical
- * parameter values, copies those parameters in a {@code "parameters"} entry in the given map.
- * Otherwise does nothing.
- *
- * @param source the operation for which to get the inverse parameters.
- * @param target where to store the inverse parameters.
- */
- static void putParameters(final SingleOperation source, final Map<String,Object> target) {
+ /*
+ * If the inverse of the given operation can be represented by inverting the sign of all numerical
+ * parameter values, copies those parameters in a "parameters" entry in the properties map.
+ * Otherwise does nothing.
+ */
final boolean isInvertible = isInvertible(source.getMethod());
final ParameterValueGroup parameters = source.getParameterValues();
final ParameterValueGroup copy = parameters.getDescriptor().createValue();
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..16fb4e3
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -0,0 +1,349 @@
+/*
+ * 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, ParameterDescriptor<Double>>(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];
+ /*
+ * Note: Synder suggests to reverse the sign of x, y and ρ₀ if n is negative. It should not done in Apache SIS
+ * implementation because (x,y) are premultiplied by n (by the normalization affine transform) before to enter
+ * in this method, so if n was negative those values have already their sign reverted.
+ */
+ 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 ef5c624..a582e93 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.
@@ -177,16 +170,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);
}
/**
@@ -195,7 +179,6 @@
CylindricalEqualArea(final CylindricalEqualArea other) {
super(other);
variant = other.variant;
- qmPolar = other.qmPolar;
}
/**
@@ -236,8 +219,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.
@@ -266,7 +249,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)
}
}
}
@@ -284,7 +267,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..71ffb1c 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;
@@ -50,6 +51,14 @@
*
* Note that since this boolean is static final, the compiler should exclude the code in the branch that is never
* executed (no need to comment-out that code).
+ *
+ *
+ * <p><b>BENCHMARK AND ANALYSIS:</b>
+ * as of July 2016, benchmarking shows no benefit in using trigonometric identities for {@code EqualAreaProjection}
+ * (contrarily to {@link ConformalProjection} where we did measured a benefit). This may be because in this class,
+ * the series expansion is unconditionally followed by iterative method in order to reach the centimetric precision.
+ * We observe that the original series expansion allows convergence in only one iteration, while the formulas using
+ * trigonometric identifies often requires two iterations. Consequently we disallow those modifications for now.</p>
*/
private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = false;
@@ -70,6 +79,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 +100,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 +116,7 @@
ci4 *= 8;
ci8 *= 64;
}
+ qmPolar = qm(1);
}
/**
@@ -113,6 +130,7 @@
ci2 = other.ci2;
ci4 = other.ci4;
ci8 = other.ci8;
+ qmPolar = other.qmPolar;
}
/**
@@ -120,6 +138,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 +150,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 +186,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
@@ -179,24 +220,62 @@
final double t4β = 0.5 - sin2_β; // = sin(4β) / ( 4⋅sin(2β))
final double t8β = (cos2_β - sin2_β)*(cos2_β*cos2_β - cos2_β + 1./8); // = sin(8β) / (32⋅sin(2β))
- assert identityEquals(t2β, sin(2*β) / ( 2 ));
- assert identityEquals(t4β, sin(4*β) / ( 8 * t2β));
- assert identityEquals(t8β, sin(8*β) / (64 * t2β));
+ assert ConformalProjection.identityEquals(t2β, sin(2*β) / ( 2 ));
+ assert ConformalProjection.identityEquals(t4β, sin(4*β) / ( 8 * t2β));
+ assert ConformalProjection.identityEquals(t8β, sin(8*β) / (64 * t2β));
- return (ci8*t8β + ci4*t4β + ci2) * t2β + β;
+ φ = (ci8*t8β + ci4*t4β + ci2) * t2β + β;
}
- }
-
- /**
- * Verifies if a trigonometric identity produced the expected value. This method is used in assertions only.
- * The tolerance threshold is approximatively 1.5E-12 (note that it still about 7000 time greater than
- * {@code Math.ulp(1.0)}).
- *
- * @see #ALLOW_TRIGONOMETRIC_IDENTITIES
- */
- private static boolean identityEquals(final double actual, final double expected) {
- // Use !(a > b) instead of (a <= b) in order to tolerate NaN.
- return !(abs(actual - expected) > (ANGULAR_TOLERANCE / 1000));
+ /*
+ * 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 φ;
+ }
+ }
+ /*
+ * In the Albers Equal Area discussion, Synder said that above algorithm does not converge if
+ *
+ * q = ±(1 - (1-ℯ²)/(2ℯ) ⋅ ln((1-ℯ)/(1+ℯ)))
+ *
+ * which we rewrite as
+ *
+ * q = ±(1 + (1-ℯ²)⋅atanh(ℯ)/ℯ)
+ *
+ * Given that y = q/(1-ℯ²) (see above comment), we rewrite as
+ *
+ * y = ±(1/(1-ℯ²) + atanh(ℯ)/ℯ) = ±qmPolar
+ *
+ * which implies sinβ = ±1. This is consistent with Synder discussion of Cylndrical Equal Area
+ * projection, where he said exactly that about the same formula (that it does not converge for
+ * β = ±90°). In both case, Synder said that the result is φ = β, with the same sign.
+ */
+ final double as = abs(sinβ);
+ if (abs(as - 1) < ANGULAR_TOLERANCE) {
+ return copySign(PI/2, y); // Value is at a pole.
+ }
+ if (as >= 1 || Double.isNaN(y)) {
+ return Double.NaN; // Value "after" the pole.
+ }
+ // Value should have converged but did not.
+ 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 9b2d962..a4447ae 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));
@@ -275,11 +275,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;
}
@@ -376,8 +376,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
@@ -561,7 +561,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 762e2ac..bc4e963 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..54535d0
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AlbersEqualAreaTest.java
@@ -0,0 +1,258 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import org.opengis.util.FactoryException;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.referencing.Formulas;
+
+import static java.lang.StrictMath.*;
+import static java.lang.Double.NaN;
+
+// 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 {
+ /**
+ * 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 {
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.AlbersEqualArea(),
+ 6370997, // Semi-major axis from Synder table 15
+ 6370997, // Semi-minor axis
+ 0, // Central meridian
+ 0, // Latitude of origin
+ 29.5, // Standard parallel 1 (from Synder table 15)
+ 45.5, // Standard parallel 2 (from Synder table 15)
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
+
+ final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres.
+ derivativeDeltas = new double[] {delta, delta};
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ 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}, // Minimal input ordinate values
+ new double[] {+20, 50}, // Maximal input ordinate values
+ new int[] { 5, 5}, // Number of points to test
+ 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 {
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.AlbersEqualArea(),
+ 6378206.4, // Semi-major axis from Synder table 15
+ 6356583.8, // Semi-minor axis
+ 0, // Central meridian
+ 0, // Latitude of origin
+ 29.5, // Standard parallel 1 (from Synder table 15)
+ 45.5, // Standard parallel 2 (from Synder table 15)
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
+
+ final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres.
+ derivativeDeltas = new double[] {delta, delta};
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ 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}, // Minimal input ordinate values
+ new double[] {+20, 50}, // Maximal input ordinate values
+ new int[] { 5, 5}, // Number of points to test
+ 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
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.AlbersEqualArea(),
+ 6400000, // Semi-major axis
+ 6400000, // Semi-minor axis
+ 0, // Central meridian
+ 0, // Latitude of origin
+ 0, // Standard parallel 1
+ 2, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
+
+ verifyTransform(new double[] {2, 1}, new double[] {223334.085, 111780.432});
+
+ // Ellipsoidal case
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.AlbersEqualArea(),
+ 6378137, // Semi-major axis
+ 6356752.314140347, // Semi-minor axis
+ 0, // Central meridian
+ 0, // Latitude of origin
+ 0, // Standard parallel 1
+ 2, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
+ verifyTransform(new double[] {2, 1}, new double[] {222571.609, 110653.327});
+ }
+
+ /**
+ * Tests a few "special" points which need special care in inverse projection algorithm.
+ *
+ * @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 testSingularity() throws FactoryException, TransformException {
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.AlbersEqualArea(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0, // Central meridian
+ 0, // Latitude of origin
+ 0, // Standard parallel 1
+ 2, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
+
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ verifyTransform(new double[] {0, 0,
+ 0, +90,
+ 0, -90},
+ new double[] {0, 0,
+ 0, +6420271.594575703, // Computed empirically with SIS (not from an external source).
+ 0, -6309429.217});
+ }
+
+ /**
+ * Tests conversion of random points with non-zero central meridian, standard parallel
+ * and false easting/northing.
+ *
+ * @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 testRandomPoints() throws FactoryException, TransformException {
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.AlbersEqualArea(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 12, // Central meridian
+ NaN, // Latitude of origin (none)
+ 24, // Standard parallel 1
+ 40, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 300, // False easting
+ 200); // False northing
+
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ derivativeDeltas = new double[] {100, 100};
+ final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres.
+ derivativeDeltas = new double[] {delta, delta};
+ verifyInDomain(new double[] {-40, 10}, // Minimal input ordinate values
+ new double[] {+40, 60}, // Maximal input ordinate values
+ new int[] { 5, 5}, // Number of points to test
+ TestUtilities.createRandomNumberGenerator());
+ }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/Benchmark.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/Benchmark.java
new file mode 100644
index 0000000..4dc1dbd
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/Benchmark.java
@@ -0,0 +1,281 @@
+/*
+ * 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.List;
+import java.util.Random;
+import java.io.IOException;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.NoninvertibleTransformException;
+import org.opengis.referencing.operation.TransformException;
+import org.opengis.util.FactoryException;
+import org.apache.sis.math.Statistics;
+import org.apache.sis.math.StatisticsFormat;
+import org.apache.sis.measure.Latitude;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.internal.util.Constants;
+import org.apache.sis.internal.referencing.provider.AbstractProvider;
+import org.apache.sis.referencing.operation.transform.LinearTransform;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.transform.MathTransformFactoryMock;
+
+
+/**
+ * Measures the performance of a given map projection implementation.
+ * This class can be used for comparing different implementation alternatives,
+ * for example with {@code ALLOW_TRIGONOMETRIC_IDENTITIES} flag on or off.
+ *
+ * <p><b>Usage:</b> modify the provider created in the {@code main} method if needed, and run. Change map projection
+ * implementation (for example by changing a {@code ALLOW_TRIGONOMETRIC_IDENTITIES} flag value) and run again.</p>
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @since 0.8
+ * @version 0.8
+ * @module
+ */
+@SuppressWarnings("UseOfSystemOutOrSystemErr")
+public final class Benchmark {
+ /**
+ * Runs the benchmark and prints the time result to the standard output.
+ * Edit this method for measuring the performance of a different map projection implementation.
+ *
+ * @param args ignored.
+ * @throws Exception if an error occurred while creating the map projection, projecting a point, <i>etc.</i>
+ */
+ public static void main(String[] args) throws Exception {
+ final Benchmark benchmark = new Benchmark( // Put on next line the provider of the projection to benchmark.
+ new org.apache.sis.internal.referencing.provider.LambertConformal2SP(),
+ 8, // Central meridian
+ 25, // Standard parallel 1
+ 40); // Standard parallel 2
+
+ for (int i=0; i<10; i++) {
+ benchmark.run(true);
+ System.gc();
+ Thread.sleep(1000);
+ }
+ benchmark.printStatistics();
+ }
+
+ /**
+ * Number of points to project. Can be modified freely.
+ */
+ private static final int NUM_POINTS = 1000000;
+
+ /**
+ * Number of dimension (not modifiable).
+ */
+ private static final int DIMENSION = 2;
+
+ /**
+ * The math transforms to use for the forward projection.
+ */
+ private final Transforms forward;
+
+ /**
+ * The math transforms to use for the inverse projection.
+ */
+ private final Transforms inverse;
+
+ /**
+ * The coordinates to project, filled with random points.
+ */
+ private final double[] coordinates;
+
+ /**
+ * The result of map projections.
+ */
+ private final double[] result;
+
+ /**
+ * Difference between inverse projections and the original coordinates.
+ */
+ private final Statistics errors;
+
+ /**
+ * Random number generator for the points to project.
+ */
+ private final Random random;
+
+ /**
+ * Prepares benchmarking for the map projection created by the given provider.
+ */
+ private Benchmark(final AbstractProvider provider,
+ final double centralMeridian,
+ final double standardParallel1,
+ final double standardParallel2) throws FactoryException, NoninvertibleTransformException
+ {
+ random = new Random();
+ final Parameters values = Parameters.castOrWrap(provider.getParameters().createValue());
+ values.parameter(Constants.SEMI_MAJOR) .setValue(MapProjectionTestCase.WGS84_A);
+ values.parameter(Constants.SEMI_MINOR) .setValue(MapProjectionTestCase.WGS84_B);
+ values.parameter(Constants.CENTRAL_MERIDIAN) .setValue(centralMeridian);
+ values.parameter(Constants.STANDARD_PARALLEL_1).setValue(standardParallel1);
+ values.parameter(Constants.STANDARD_PARALLEL_2).setValue(standardParallel2);
+ forward = new Transforms("Forward", new MathTransformFactoryMock(provider).createParameterizedTransform(values));
+ inverse = new Transforms("Inverse", forward.projection.inverse());
+ coordinates = new double[NUM_POINTS * DIMENSION];
+ final double λmin = centralMeridian - 40;
+ final double Δλ = 80;
+ final double φmin = standardParallel1 * 0.75;
+ final double Δφ = Math.min(standardParallel2 * 1.25, Latitude.MAX_VALUE) - φmin;
+ for (int i=0; i<coordinates.length;) {
+ coordinates[i++] = random.nextDouble() * Δλ + λmin; // Longitude
+ coordinates[i++] = random.nextDouble() * Δφ + φmin; // Latitude
+ }
+ result = new double[coordinates.length];
+ errors = new Statistics("Errors (cm)");
+ }
+
+ /**
+ * Decomposition of the {@code MathTransform} that perform map projections.
+ */
+ private static final class Transforms {
+ /**
+ * The original (full) map projection.
+ */
+ final MathTransform projection;
+
+ /**
+ * The transform that performs the non-linear part of the map projection.
+ * This is the transform that we want to benchmark.
+ */
+ private final MathTransform kernel;
+
+ /**
+ * The normalization performed before the map projection.
+ */
+ private final LinearTransform normalize;
+
+ /**
+ * The denormalization performed after the map projection.
+ */
+ private final LinearTransform denormalize;
+
+ /**
+ * Statistics about the time needed for performing the map projection.
+ */
+ private final Statistics performance;
+
+ /**
+ * Creates a decomposition of the given map projection.
+ *
+ * @param label a label for display purpose.
+ * @param projection the map projection to benchmark.
+ */
+ private Transforms(final String label, final MathTransform projection) {
+ this.projection = projection;
+ performance = new Statistics(label);
+ final List<MathTransform> steps = MathTransforms.getSteps(projection);
+ int kernelIndex = -1;
+ for (int i=steps.size(); --i >= 0;) {
+ if (!(steps.get(i) instanceof LinearTransform)) {
+ if (kernelIndex >= 0) {
+ throw new IllegalArgumentException("Found more than one non-linear kernel.");
+ }
+ kernelIndex = i;
+ }
+ }
+ if (kernelIndex < 0) {
+ throw new IllegalArgumentException("Non-linear kernel not found.");
+ }
+ kernel = steps.get(kernelIndex);
+ normalize = singleton(steps.subList(0, kernelIndex));
+ denormalize = singleton(steps.subList(kernelIndex + 1, steps.size()));
+ }
+
+ /**
+ * Returns the single elements of the given list, or an identity transform if none.
+ */
+ private static LinearTransform singleton(final List<MathTransform> components) {
+ switch (components.size()) {
+ case 0: return MathTransforms.identity(DIMENSION);
+ case 1: return (LinearTransform) components.get(0);
+ default: throw new IllegalArgumentException("Unexpected number of components.");
+ }
+ }
+
+ /**
+ * Runs the benchmark on the complete map projection, including the linear parts.
+ */
+ final void runComplete(final double[] sources, final double[] targets) throws TransformException {
+ long time = System.nanoTime();
+ projection.transform(sources, 0, targets, 0, NUM_POINTS);
+ time = System.nanoTime() - time;
+ final double seconds = time * NANOS_TO_SECONDS;
+ System.out.printf("%s time: %1.4f%n", performance.name(), seconds);
+ performance.accept(seconds);
+ }
+
+ /**
+ * Runs the benchmark only on the non-linear part of the map projection,
+ * ignoring the linear parts.
+ */
+ final void runKernel(final double[] sources, final double[] targets) throws TransformException {
+ normalize.transform(sources, 0, targets, 0, NUM_POINTS);
+ long time = System.nanoTime();
+ kernel.transform(targets, 0, targets, 0, NUM_POINTS);
+ time = System.nanoTime() - time;
+ denormalize.transform(targets, 0, targets, 0, NUM_POINTS);
+ final double seconds = time * NANOS_TO_SECONDS;
+ System.out.printf("%s time: %1.4f%n", performance.name(), seconds);
+ performance.accept(seconds);
+ }
+
+ /**
+ * For reporting the time measured by {@link System#nanoTime()}.
+ */
+ private static final double NANOS_TO_SECONDS = 1E-9;
+ }
+
+ /**
+ * Runs the benchmark.
+ *
+ * @param kernelOnly {@code true} for measuring the performance of only the non-linear part,
+ * or {@code false} for measuring the performance of the whole projection.
+ */
+ private void run(final boolean kernelOnly) throws TransformException {
+ if (kernelOnly) {
+ forward.runKernel(coordinates, result);
+ inverse.runKernel(result, result);
+ } else {
+ forward.runComplete(coordinates, result);
+ inverse.runComplete(result, result);
+ }
+ for (int i=0; i<NUM_POINTS;) {
+ final double dx = result[i] - coordinates[i]; coordinates[i++] += random.nextDouble() - 0.5;
+ final double dy = result[i] - coordinates[i]; coordinates[i++] += random.nextDouble() - 0.5;
+ errors.accept(Math.hypot(dx*DEGREES_TO_CENTIMETRES, dy*DEGREES_TO_CENTIMETRES));
+ }
+ }
+
+ /**
+ * For reporting the errors in more convenient units.
+ * This is an approximative conversion based on the nautical mile length.
+ */
+ private static final double DEGREES_TO_CENTIMETRES = 60*1852*100;
+
+ /**
+ * Prints statistics about measured time.
+ */
+ private void printStatistics() throws IOException {
+ System.out.println();
+ StatisticsFormat.getInstance().format(new Statistics[] {forward.performance, inverse.performance}, System.out);
+ System.out.printf("%nAverage error is %1.2E cm (standard deviation %1.1E).%n", errors.mean(), errors.standardDeviation(false));
+ System.out.flush();
+ }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
index e9498ae..98246fa 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
@@ -178,7 +178,7 @@
/**
* Tests the {@link NormalizedProjection#dy_dφ(double, double)} method.
*
- * @throws TransformException Should never happen.
+ * @throws TransformException if an error occurred while projecting a point.
*/
@Test
@DependsOnMethod("testExpOfNorthing")
@@ -235,7 +235,7 @@
* In theory only the [-90° … +90°] range needs to be tested. However the function is still
* consistent in the [-90° … +270°] range so we test that range for tracking this fact.
*
- * @throws ProjectionException Should never happen.
+ * @throws ProjectionException if an error occurred while projecting a point.
*/
@Test
@DependsOnMethod("testExpOfNorthing")
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalEqualAreaTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalEqualAreaTest.java
index f71db67..3e9267b 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalEqualAreaTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/CylindricalEqualAreaTest.java
@@ -22,10 +22,12 @@
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.referencing.provider.LambertCylindricalEqualArea;
import org.apache.sis.internal.referencing.provider.LambertCylindricalEqualAreaSpherical;
+import org.apache.sis.referencing.operation.transform.CoordinateDomain;
import org.apache.sis.test.DependsOnMethod;
import org.junit.Test;
import static java.lang.StrictMath.*;
+import static java.lang.Double.NaN;
/**
@@ -38,16 +40,6 @@
*/
public final strictfp class CylindricalEqualAreaTest extends MapProjectionTestCase {
/**
- * Creates a map projection.
- */
- private void createCompleteProjection(final boolean ellipse,
- final double centralMeridian, final double standardParallel) throws FactoryException
- {
- createCompleteProjection(new LambertCylindricalEqualArea(),
- ellipse, centralMeridian, 0, standardParallel, 1, 0, 0);
- }
-
- /**
* Tests the derivatives at a few points. This method compares the derivatives computed by
* the projection with an estimation of derivatives computed by the finite differences method.
*
@@ -69,7 +61,16 @@
*/
@Test
public void testEllipsoidal() throws FactoryException, TransformException {
- createCompleteProjection(true, 0, 0);
+ createCompleteProjection(new LambertCylindricalEqualArea(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0, // Central meridian
+ NaN, // Latitude of origin
+ 0, // Standard parallel 1
+ NaN, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
toleranceModifier = ToleranceModifier.PROJECTION;
final double λ = 2;
@@ -90,7 +91,16 @@
@Test
@DependsOnMethod("testEllipsoidal")
public void testSpherical() throws FactoryException, TransformException {
- createCompleteProjection(false, 0, 0);
+ createCompleteProjection(new LambertCylindricalEqualArea(),
+ 6371007, // Semi-major axis length
+ 6371007, // Semi-minor axis length
+ 0, // Central meridian
+ NaN, // Latitude of origin
+ 0, // Standard parallel 1
+ NaN, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
toleranceModifier = ToleranceModifier.PROJECTION;
final double λ = 2;
@@ -114,7 +124,16 @@
@Test
@DependsOnMethod("testSpherical")
public void testSphericalWithConformalSphereRadius() throws FactoryException, TransformException {
- createCompleteProjection(new LambertCylindricalEqualAreaSpherical(), true, 0, 0, 0, 1, 0, 0);
+ createCompleteProjection(new LambertCylindricalEqualAreaSpherical(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0, // Central meridian
+ NaN, // Latitude of origin
+ 0, // Standard parallel 1
+ NaN, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 0, // False easting
+ 0); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
toleranceModifier = ToleranceModifier.PROJECTION;
final double λ = 2;
@@ -125,4 +144,32 @@
new double[] {x, y, -x, y, x, -y, -x, -y});
testDerivative();
}
+
+ /**
+ * Tests conversion of random points with non-zero central meridian, standard parallel
+ * and false easting/northing.
+ *
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a point.
+ */
+ @Test
+ public void testRandomPoints() throws FactoryException, TransformException {
+ createCompleteProjection(new LambertCylindricalEqualArea(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 12, // Central meridian
+ NaN, // Latitude of origin (none)
+ 24, // Standard parallel 1
+ NaN, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 300, // False easting
+ 200); // False northing
+
+ tolerance = Formulas.LINEAR_TOLERANCE;
+ toleranceModifier = ToleranceModifier.PROJECTION;
+ derivativeDeltas = new double[] {100, 100};
+ final double delta = toRadians(100.0 / 60) / 1852; // Approximatively 100 metres.
+ derivativeDeltas = new double[] {delta, delta};
+ verifyInDomain(CoordinateDomain.GEOGRAPHIC_SAFE, 0);
+ }
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EquirectangularTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EquirectangularTest.java
index 3b03f5e..3a0da83 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EquirectangularTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EquirectangularTest.java
@@ -27,6 +27,7 @@
import org.apache.sis.test.ReferencingAssert;
import org.junit.Test;
+import static java.lang.Double.NaN;
import static java.lang.StrictMath.toRadians;
import static org.apache.sis.internal.metadata.ReferencingServices.AUTHALIC_RADIUS;
@@ -38,7 +39,7 @@
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.6
- * @version 0.6
+ * @version 0.8
* @module
*/
public final strictfp class EquirectangularTest extends MapProjectionTestCase {
@@ -61,7 +62,7 @@
* an affine transform, the WKT formatter should handle this projection in a special way and shows the
* projection parameters instead than the affine transform parameters (except in "show internal" mode).
*
- * @throws FactoryException should never happen.
+ * @throws FactoryException if an error occurred while creating the map projection.
*/
@Test
public void testWKT() throws FactoryException {
@@ -90,19 +91,19 @@
/**
* Tests a simple transform on a sphere.
*
- * @throws FactoryException should never happen.
- * @throws TransformException should never happen.
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a point.
*/
@Test
public void testSimpleTransform() throws FactoryException, TransformException {
createCompleteProjection();
verifyTransform(
- new double[] { // (λ,φ) coordinates in degrees to project.
+ new double[] { // (λ,φ) coordinates in degrees to project.
0, 0,
2, 0,
0, 3
},
- new double[] { // Expected (x,y) results in metres.
+ new double[] { // Expected (x,y) results in metres.
0, 0,
AUTHALIC_RADIUS*toRadians(2), 0,
0, AUTHALIC_RADIUS*toRadians(3)
@@ -113,18 +114,21 @@
* Tests conversion of random points. This test is actually of limited interest since the Equirectangular
* projection is implemented by an affine transform, which has been tested elsewhere.
*
- * @throws FactoryException should never happen.
- * @throws TransformException should never happen.
+ * @throws FactoryException if an error occurred while creating the map projection.
+ * @throws TransformException if an error occurred while projecting a point.
*/
@Test
public void testRandomPoints() throws FactoryException, TransformException {
- createCompleteProjection(new Equirectangular(), true,
- 0.5, // Central meridian
- 0, // Latitude of origin (none)
- 20, // Standard parallel
- 1, // Scale factor (none)
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new Equirectangular(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0.5, // Central meridian
+ 0, // Latitude of origin (none)
+ 20, // Standard parallel 1
+ NaN, // Standard parallel 2
+ NaN, // Scale factor (none)
+ 200, // False easting
+ 100); // False northing
tolerance = Formulas.LINEAR_TOLERANCE; // Not NORMALIZED_TOLERANCE since this is not a NormalizedProjection.
derivativeDeltas = new double[] {100, 100};
verifyInDomain(CoordinateDomain.GEOGRAPHIC, 0);
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
index 51918c4..67e69a8 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
@@ -49,7 +49,7 @@
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @since 0.6
- * @version 0.7
+ * @version 0.8
* @module
*/
@DependsOn(ConformalProjectionTest.class)
@@ -112,12 +112,12 @@
/**
* Tests the projection at some special latitudes (0, ±π/2, NaN and others).
*
- * @throws ProjectionException Should never happen.
+ * @throws ProjectionException if an error occurred while projecting a point.
*/
@Test
public void testSpecialLatitudes() throws ProjectionException {
- if (transform == null) { // May have been initialized by 'testSphericalCase'.
- createNormalizedProjection(true, 40); // Elliptical case
+ if (transform == null) { // May have been initialized by 'testSphericalCase'.
+ createNormalizedProjection(true, 40); // Elliptical case
}
final double INF = POSITIVE_INFINITY;
assertEquals ("Not a number", NaN, transform(NaN), NORMALIZED_TOLERANCE);
@@ -161,7 +161,7 @@
* Tests the derivatives at a few points. This method compares the derivatives computed by
* the projection with an estimation of derivatives computed by the finite differences method.
*
- * @throws TransformException Should never happen.
+ * @throws TransformException if an error occurred while projecting a point.
*/
@Test
@DependsOnMethod("testSpecialLatitudes")
@@ -246,28 +246,34 @@
@Test
@DependsOnMethod("testLambertConicConformal1SP")
public void testLambertConicConformalWestOrientated() throws FactoryException, TransformException {
- createCompleteProjection(new LambertConformal1SP(), false,
- 0.5, // Central meridian
- 40, // Latitude of origin
- 0, // Standard parallel (none)
- 0.997, // Scale factor
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new LambertConformal1SP(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0.5, // Central meridian
+ 40, // Latitude of origin
+ NaN, // Standard parallel 1
+ NaN, // Standard parallel 2
+ 0.997, // Scale factor
+ 200, // False easting
+ 100); // False northing
final MathTransform reference = transform;
- createCompleteProjection(new LambertConformalWest(), false,
- 0.5, // Central meridian
- 40, // Latitude of origin
- 0, // Standard parallel (none)
- 0.997, // Scale factor
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new LambertConformalWest(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0.5, // Central meridian
+ 40, // Latitude of origin
+ NaN, // Standard parallel 1
+ NaN, // Standard parallel 2
+ 0.997, // Scale factor
+ 200, // False easting
+ 100); // False northing
final Random random = TestUtilities.createRandomNumberGenerator();
final double[] sources = new double[20];
for (int i=0; i<sources.length;) {
- sources[i++] = 20 * random.nextDouble(); // Longitude
- sources[i++] = 10 * random.nextDouble() + 35; // Latitude
+ sources[i++] = 20 * random.nextDouble(); // Longitude
+ sources[i++] = 10 * random.nextDouble() + 35; // Latitude
}
final double[] expected = new double[sources.length];
reference.transform(sources, 0, expected, 0, sources.length/2);
@@ -326,13 +332,16 @@
@Test
@DependsOnMethod("testSphericalCase")
public void compareEllipticalWithSpherical() throws FactoryException, TransformException {
- createCompleteProjection(new LambertConformal1SP(), false,
- 0.5, // Central meridian
- 40, // Latitude of origin
- 0, // Standard parallel (none)
- 0.997, // Scale factor
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new LambertConformal1SP(),
+ 6371007, // Semi-major axis length
+ 6371007, // Semi-minor axis length
+ 0.5, // Central meridian
+ 40, // Latitude of origin
+ NaN, // Standard parallel 1
+ NaN, // Standard parallel 2
+ 0.997, // Scale factor
+ 200, // False easting
+ 100); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
compareEllipticalWithSpherical(CoordinateDomain.GEOGRAPHIC_SAFE, 0);
}
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..619a8e8 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
@@ -16,9 +16,9 @@
*/
package org.apache.sis.referencing.operation.projection;
-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,8 +28,10 @@
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.Double.isNaN;
import static java.lang.StrictMath.*;
import static org.junit.Assert.*;
@@ -39,10 +41,20 @@
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.6
- * @version 0.6
+ * @version 0.8
* @module
*/
-strictfp class MapProjectionTestCase extends MathTransformTestCase {
+abstract strictfp class MapProjectionTestCase extends MathTransformTestCase {
+ /**
+ * Semi-major axis length of WGS84 ellipsoid.
+ */
+ static final double WGS84_A = 6378137;
+
+ /**
+ * Semi-minor axis length of WGS84 ellipsoid.
+ */
+ static final double WGS84_B = 6356752.314245179;
+
/**
* Tolerance level for comparing formulas on the unitary sphere or ellipsoid.
*/
@@ -55,11 +67,22 @@
}
/**
- * Returns the parameters to use for instantiating the projection to test.
+ * Instantiates the object to use for running GeoAPI 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.
+ * @return the GeoAPI test class using the given provider.
+ */
+ static ParameterizedTransformTest createGeoApiTest(final MapProjection provider) {
+ return new ParameterizedTransformTest(new MathTransformFactoryMock(provider));
+ }
+
+ /**
+ * Returns the parameters to use for instantiating the projection to test.
+ * The parameters are initialized with the ellipse semi-axis lengths.
+ *
+ * @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());
@@ -73,45 +96,57 @@
}
/**
- * 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.
- */
- static ParameterizedTransformTest createGeoApiTest(final MapProjection provider) {
- return new ParameterizedTransformTest(new MathTransformFactoryMock(provider));
- }
-
- /**
* Initializes a complete projection (including conversion from degrees to radians) for the given provider.
- * This method uses arbitrary central meridian, scale factor, false easting and false northing for increasing
- * the chances to detect a mismatch. The result is stored in the {@link #transform} field.
*/
- final void createCompleteProjection(final DefaultOperationMethod provider, final boolean ellipse,
+ final void createCompleteProjection(final DefaultOperationMethod provider,
+ final double semiMajor,
+ final double semiMinor,
final double centralMeridian,
final double latitudeOfOrigin,
- final double standardParallel,
+ final double standardParallel1,
+ final double standardParallel2,
final double scaleFactor,
final double falseEasting,
final double falseNorthing) throws FactoryException
{
- final Parameters parameters = parameters(provider, ellipse);
- if (centralMeridian != 0) parameters.parameter(Constants.CENTRAL_MERIDIAN) .setValue(centralMeridian, NonSI.DEGREE_ANGLE);
- if (latitudeOfOrigin != 0) parameters.parameter(Constants.LATITUDE_OF_ORIGIN) .setValue(latitudeOfOrigin);
- if (standardParallel != 0) parameters.parameter(Constants.STANDARD_PARALLEL_1).setValue(standardParallel);
- if (scaleFactor != 1) parameters.parameter(Constants.SCALE_FACTOR) .setValue(scaleFactor);
- if (falseEasting != 0) parameters.parameter(Constants.FALSE_EASTING) .setValue(falseEasting);
- if (falseNorthing != 0) parameters.parameter(Constants.FALSE_NORTHING) .setValue(falseNorthing);
- transform = new MathTransformFactoryMock(provider).createParameterizedTransform(parameters);
+ final Parameters values = Parameters.castOrWrap(provider.getParameters().createValue());
+ values.parameter(Constants.SEMI_MAJOR).setValue(semiMajor);
+ values.parameter(Constants.SEMI_MINOR).setValue(semiMinor);
+ if (semiMajor == WGS84_A && semiMinor == WGS84_B) {
+ values.parameter(Constants.INVERSE_FLATTENING).setValue(298.257223563);
+ }
+ if (!isNaN(centralMeridian)) values.parameter(Constants.CENTRAL_MERIDIAN) .setValue(centralMeridian);
+ if (!isNaN(latitudeOfOrigin)) values.parameter(Constants.LATITUDE_OF_ORIGIN) .setValue(latitudeOfOrigin);
+ if (!isNaN(standardParallel1)) values.parameter(Constants.STANDARD_PARALLEL_1).setValue(standardParallel1);
+ if (!isNaN(standardParallel2)) values.parameter(Constants.STANDARD_PARALLEL_2).setValue(standardParallel2);
+ if (!isNaN(scaleFactor)) values.parameter(Constants.SCALE_FACTOR) .setValue(scaleFactor);
+ if (!isNaN(falseEasting)) values.parameter(Constants.FALSE_EASTING) .setValue(falseEasting);
+ if (!isNaN(falseNorthing)) values.parameter(Constants.FALSE_NORTHING) .setValue(falseNorthing);
+ transform = new MathTransformFactoryMock(provider).createParameterizedTransform(values);
validate();
}
/**
+ * 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 {
@@ -119,7 +154,7 @@
coordinate[1] = φ;
((NormalizedProjection) transform).transform(coordinate, 0, coordinate, 0, false);
final double y = coordinate[1];
- if (!Double.isNaN(y) && !Double.isInfinite(y)) {
+ if (!isNaN(y) && !Double.isInfinite(y)) {
assertEquals(0, coordinate[0], tolerance);
}
return y;
@@ -129,8 +164,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 {
@@ -138,7 +173,7 @@
coordinate[1] = y;
((NormalizedProjection) transform).inverseTransform(coordinate, 0, coordinate, 0);
final double φ = coordinate[1];
- if (!Double.isNaN(φ)) {
+ if (!isNaN(φ)) {
/*
* Opportunistically verify that the longitude is still zero. However the longitude value is meaningless
* at poles. We can not always use coordinate[0] for testing if we are at a pole because its calculation
@@ -159,9 +194,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/referencing/operation/projection/MercatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java
index 002a8e7..0f15fe5 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorTest.java
@@ -44,7 +44,7 @@
* @author Simon Reynard (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @since 0.6
- * @version 0.7
+ * @version 0.8
* @module
*/
@DependsOn(ConformalProjectionTest.class)
@@ -70,7 +70,7 @@
* the ellipsoid eccentricity. We expect nothing else because all other parameters are used
* by the (de)normalization affine transforms instead than the {@link Mercator} class itself.
*
- * @throws NoninvertibleTransformException should never happen.
+ * @throws NoninvertibleTransformException if the transform can not be inverted.
*
* @see LambertConicConformalTest#testNormalizedWKT()
*/
@@ -90,18 +90,21 @@
* Tests WKT of a complete map projection.
*
* @throws FactoryException if an error occurred while creating the map projection.
- * @throws NoninvertibleTransformException should never happen.
+ * @throws NoninvertibleTransformException if the transform can not be inverted.
*/
@Test
@DependsOnMethod("testNormalizedWKT")
public void testCompleteWKT() throws FactoryException, NoninvertibleTransformException {
- createCompleteProjection(new Mercator1SP(), true,
- 0.5, // Central meridian
- 0, // Latitude of origin (none)
- 0, // Standard parallel (none)
- 0.997, // Scale factor
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new Mercator1SP(),
+ WGS84_A, // Semi-major axis length
+ WGS84_B, // Semi-minor axis length
+ 0.5, // Central meridian
+ NaN, // Latitude of origin (none)
+ NaN, // Standard parallel 1 (none)
+ NaN, // Standard parallel 2 (none)
+ 0.997, // Scale factor
+ 200, // False easting
+ 100); // False northing
assertWktEquals("PARAM_MT[“Mercator_1SP”,\n" +
" PARAMETER[“semi_major”, 6378137.0],\n" +
@@ -289,13 +292,16 @@
@Test
@DependsOnMethod("testSphericalCase")
public void compareEllipticalWithSpherical() throws FactoryException, TransformException {
- createCompleteProjection(new Mercator1SP(), false,
- 0.5, // Central meridian
- 0, // Latitude of origin (none)
- 0, // Standard parallel (none)
- 0.997, // Scale factor
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new Mercator1SP(),
+ 6371007, // Semi-major axis length
+ 6371007, // Semi-minor axis length
+ 0.5, // Central meridian
+ NaN, // Latitude of origin (none)
+ NaN, // Standard parallel 1 (none)
+ NaN, // Standard parallel 2 (none)
+ 0.997, // Scale factor
+ 200, // False easting
+ 100); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
compareEllipticalWithSpherical(CoordinateDomain.GEOGRAPHIC_SAFE, 0);
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java
index 66ea07a..d415658 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/PolarStereographicTest.java
@@ -32,6 +32,7 @@
import org.apache.sis.test.DependsOn;
import org.junit.Test;
+import static java.lang.Double.NaN;
import static java.lang.StrictMath.*;
@@ -40,7 +41,7 @@
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.6
- * @version 0.6
+ * @version 0.8
* @module
*/
@DependsOn(NormalizedProjectionTest.class)
@@ -169,11 +170,14 @@
private void compareEllipticalWithSpherical(final CoordinateDomain domain, final double latitudeOfOrigin,
final long randomSeed) throws FactoryException, TransformException
{
- createCompleteProjection(new PolarStereographicA(), false,
- 0.5, // Central meridian
- latitudeOfOrigin, // Latitude of origin
- 0, // Standard parallel (none)
- 0.994, // Scale factor
+ createCompleteProjection(new PolarStereographicA(),
+ 6371007, // Semi-major axis length
+ 6371007, // Semi-minor axis length
+ 0.5, // Central meridian
+ latitudeOfOrigin, // Latitude of origin
+ NaN, // Standard parallel 1 (none)
+ NaN, // Standard parallel 2 (none)
+ 0.994, // Scale factor
200, // False easting
100); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
index bf09139..5952a8a 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/TransverseMercatorTest.java
@@ -27,6 +27,7 @@
import org.apache.sis.test.DependsOn;
import org.junit.Test;
+import static java.lang.Double.NaN;
import static java.lang.StrictMath.toRadians;
import static org.apache.sis.test.Assert.*;
@@ -96,13 +97,16 @@
@Test
@DependsOnMethod("testTransverseMercator")
public void compareEllipticalWithSpherical() throws FactoryException, TransformException {
- createCompleteProjection(new org.apache.sis.internal.referencing.provider.TransverseMercator(), false,
- 0.5, // Central meridian
- 2.5, // Latitude of origin
- 0, // Standard parallel (none)
- 0.997, // Scale factor
- 200, // False easting
- 100); // False northing
+ createCompleteProjection(new org.apache.sis.internal.referencing.provider.TransverseMercator(),
+ 6371007, // Semi-major axis length
+ 6371007, // Semi-minor axis length
+ 0.5, // Central meridian
+ 2.5, // Latitude of origin
+ NaN, // Standard parallel 1 (none)
+ NaN, // Standard parallel 1 (none)
+ 0.997, // Scale factor
+ 200, // False easting
+ 100); // False northing
tolerance = Formulas.LINEAR_TOLERANCE;
compareEllipticalWithSpherical(CoordinateDomain.RANGE_10, 0);
}
@@ -110,7 +114,7 @@
/**
* Creates a projection and derivates a few points.
*
- * @throws TransformException Should never happen.
+ * @throws TransformException if an error occurred while projecting a point.
*/
@Test
public void testSphericalDerivative() throws TransformException {
@@ -127,7 +131,7 @@
/**
* Creates a projection and derivates a few points.
*
- * @throws TransformException Should never happen.
+ * @throws TransformException if an error occurred while projecting a point.
*/
@Test
public void testEllipsoidalDerivative() 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,
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
index 7dc9a57..efd82d0 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
@@ -807,6 +807,11 @@
public static final short NonInvertibleMatrix_2 = 81;
/**
+ * Can not invert the “{0}” operation.
+ */
+ public static final short NonInvertibleOperation_1 = 232;
+
+ /**
* Transform is not invertible.
*/
public static final short NonInvertibleTransform = 82;
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
index 222180f..714ccb5 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
@@ -166,6 +166,7 @@
NonEquilibratedParenthesis_2 = Missing a \u2018{1}\u2019 parenthesis in \u201c{0}\u201d.
NonInvertibleConversion = Conversion is not invertible.
NonInvertibleMatrix_2 = Non invertible {0}\u00d7{1} matrix.
+NonInvertibleOperation_1 = Can not invert the \u201c{0}\u201d operation.
NonInvertibleTransform = Transform is not invertible.
NonAngularUnit_1 = \u201c{0}\u201d is not an angular unit.
NonLinearUnit_1 = \u201c{0}\u201d is not a linear unit.
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
index 4c43ae0..a3145ad 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
@@ -163,6 +163,7 @@
NonEquilibratedParenthesis_2 = Il manque une parenth\u00e8se \u2018{1}\u2019 dans \u00ab\u202f{0}\u202f\u00bb.
NonInvertibleConversion = La conversion n\u2019est pas inversible.
NonInvertibleMatrix_2 = Matrice {0}\u00d7{1} non inversible.
+NonInvertibleOperation_1 = Ne peut pas inverser l\u2019op\u00e9ration \u00ab\u202f{0}\u202f\u00bb.
NonInvertibleTransform = La transformation n\u2019est pas inversible.
NonAngularUnit_1 = \u00ab\u202f{0}\u202f\u00bb n\u2019est pas une unit\u00e9 d\u2019angles.
NonLinearUnit_1 = \u00ab\u202f{0}\u202f\u00bb n\u2019est pas une unit\u00e9 de longueurs.
