This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit afdf695c15841e6f6d2ad23c24833111b4d0331e Author: Martin Desruisseaux <[email protected]> AuthorDate: Mon Aug 31 14:50:29 2020 +0200 When there is no direct Bursa-Wolf parameters between two CRS, accepts an indirect transformation through a common hub (typically WGS 84). This is a partial revert of commit 504e3ce3bfd44c6aea1ff9f4d5ba298d4b0a22f2 (November 14, 2013) together with new code for updating positional accuracy. An arbitrary value of 100 metres is taken when the domains of validity intersect, and 3000 metres (presumed worst case scenario) otherwise. Mailing list: https://lists.apache.org/thread.html/r751d8fca4997faa59202a6cf4c85f1b1824bf06abc315995fe927719%40%3Cuser.sis.apache.org%3E --- .../sis/internal/referencing/AnnotatedMatrix.java | 140 +++++++++++++++++++++ .../sis/internal/referencing/ExtentSelector.java | 33 ++++- .../referencing/PositionalAccuracyConstant.java | 40 ++++-- .../referencing/datum/DefaultGeodeticDatum.java | 53 ++++++-- .../operation/CoordinateOperationFinder.java | 15 ++- .../datum/DefaultGeodeticDatumTest.java | 44 ++++++- .../operation/CoordinateOperationFinderTest.java | 74 ++++++++++- 7 files changed, 370 insertions(+), 29 deletions(-) diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/AnnotatedMatrix.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/AnnotatedMatrix.java new file mode 100644 index 0000000..8308ca4 --- /dev/null +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/AnnotatedMatrix.java @@ -0,0 +1,140 @@ +/* + * 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; + +import org.opengis.referencing.operation.Matrix; +import org.opengis.metadata.quality.PositionalAccuracy; + + +/** + * A matrix augmented with annotation about transformation accuracy. + * We use this class for passing additional information in methods that returns only a {@link Matrix}. + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.1 + * @since 1.1 + * @module + */ +public final class AnnotatedMatrix implements Matrix, Cloneable { + /** + * The matrix which contains the actual values. + */ + private final Matrix matrix; + + /** + * Accuracy associated with this matrix. + */ + public final PositionalAccuracy accuracy; + + /** + * Creates a new matrix wrapping the given matrix. + * + * @param matrix the matrix which contains the actual values. + * @param accuracy accuracy associated with this matrix. + */ + private AnnotatedMatrix(final Matrix matrix, final PositionalAccuracy accuracy) { + this.matrix = matrix; + this.accuracy = accuracy; + } + + /** + * Returns an {@link AnnotatedMatrix} associates with {@link PositionalAccuracyConstant#INDIRECT_SHIFT_APPLIED}. + * + * @param matrix the matrix to wrap. + * @param intersect whether an intersection has been found between the two datum. + * @return the annotated matrix. + */ + public static Matrix indirect(final Matrix matrix, final boolean intersect) { + return new AnnotatedMatrix(matrix, + intersect ? PositionalAccuracyConstant.INDIRECT_SHIFT_APPLIED + : PositionalAccuracyConstant.DATUM_SHIFT_OMITTED); + } + + /** + * Returns the number of rows in this matrix. + */ + @Override + public int getNumRow() { + return matrix.getNumRow(); + } + + /** + * Returns the number of columns in this matrix. + */ + @Override + public int getNumCol() { + return matrix.getNumCol(); + } + + /** + * Returns {@code true} if this matrix is an identity matrix. + */ + @Override + public boolean isIdentity() { + return matrix.isIdentity(); + } + + /** + * Retrieves the value at the specified row and column of this matrix. + */ + @Override + public double getElement(int row, int column) { + return matrix.getElement(row, column); + } + + /** + * Modifies the value at the specified row and column of this matrix. + */ + @Override + public void setElement(int row, int column, double value) { + matrix.setElement(row, column, value); + } + + /** + * Returns a clone of this matrix. + */ + @Override + @SuppressWarnings("CloneDoesntCallSuperClone") + public Matrix clone() { + return new AnnotatedMatrix(matrix.clone(), accuracy); + } + + /** + * Compares this matrix with the given object for equality. + */ + @Override + public boolean equals(final Object obj) { + if (obj instanceof AnnotatedMatrix) { + final AnnotatedMatrix other = (AnnotatedMatrix) obj; + return matrix.equals(other.matrix) && accuracy.equals(other.accuracy); + } + return false; + } + + /** + * Returns a hash code value for this matrix. + */ + @Override + public int hashCode() { + return matrix.hashCode(); + } + + @Override + public String toString() { + return matrix.toString(); + } +} diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java index 6201282..59b1a8f 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java @@ -27,7 +27,7 @@ import org.apache.sis.metadata.iso.extent.Extents; * This may be extended to other kind of extent in any future SIS version. * * @author Martin Desruisseaux (Geomatys) - * @version 0.4 + * @version 1.1 * * @param <T> the type of object to be selected. * @@ -37,8 +37,9 @@ import org.apache.sis.metadata.iso.extent.Extents; public final class ExtentSelector<T> { /** * The area of interest, or {@code null} if none. + * This is specified at construction time, but can be modified later. */ - private final GeographicBoundingBox areaOfInterest; + private GeographicBoundingBox areaOfInterest; /** * The best object found so far. @@ -62,6 +63,25 @@ public final class ExtentSelector<T> { } /** + * Returns the area of interest. + * + * @return area of interest, or {@code null} if none. + */ + public final GeographicBoundingBox getAreaOfInterest() { + return areaOfInterest; + } + + /** + * Sets the area of interest to the intersection of the two given arguments. + * + * @param a1 first area of interest as a bounding box, or {@code null}. + * @param a2 second area of interest as an extent, or {@code null}. + */ + public final void setAreaOfInterest(final GeographicBoundingBox a1, final Extent a2) { + areaOfInterest = Extents.intersection(a1, Extents.getGeographicBoundingBox(a2)); + } + + /** * Evaluates the given extent against the criteria represented by the {@code ExtentSelector}. * If the intersection between the given extent and the area of interest is greater than any * previous intersection, then the given extent and object are remembered as the best match @@ -97,4 +117,13 @@ public final class ExtentSelector<T> { public T best() { return best; } + + /** + * Returns {@code true} if an intersection has been found. + * + * @return whether an intersection has been found. + */ + public boolean hasIntersection() { + return largestArea > 0; + } } diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionalAccuracyConstant.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionalAccuracyConstant.java index a510fae..9622a3d 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionalAccuracyConstant.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/PositionalAccuracyConstant.java @@ -43,7 +43,7 @@ import org.apache.sis.util.resources.Vocabulary; * Pre-defined positional accuracy resulting from some coordinate operations. * * @author Martin Desruisseaux (Geomatys) - * @version 0.6 + * @version 1.1 * * @see org.opengis.referencing.operation.Transformation#getCoordinateOperationAccuracy() * @@ -80,6 +80,14 @@ public final class PositionalAccuracyConstant extends DefaultAbsoluteExternalPos private static final double DATUM_SHIFT_ACCURACY = 25; /** + * Default accuracy of datum shifts when using an intermediate datum (typically WGS 84). + * Since this is a concatenation of two datum shifts, we use twice {@link #DATUM_SHIFT_ACCURACY}. + * The result is multiplied by 2 again as a margin because we have no guarantees that the domain + * of validity of the two datum are close enough for making this concatenation valid. + */ + public static final double INDIRECT_SHIFT_ACCURACY = 100; + + /** * Indicates that a {@linkplain org.opengis.referencing.operation.Transformation transformation} * requires a datum shift and some method has been applied. Datum shift methods often use * {@linkplain org.apache.sis.referencing.datum.BursaWolfParameters Bursa Wolf parameters}, @@ -94,11 +102,19 @@ public final class PositionalAccuracyConstant extends DefaultAbsoluteExternalPos * been found. Such datum shifts are approximations and may have 1 kilometer error. */ public static final PositionalAccuracy DATUM_SHIFT_OMITTED; + + /** + * Indicates that a {@linkplain org.opengis.referencing.operation.Transformation transformation} + * requires a datum shift, but only an indirect method has been found. The indirect method uses + * an intermediate datum, typically WGS 84. + */ + public static final PositionalAccuracy INDIRECT_SHIFT_APPLIED; static { final InternationalString desc = Vocabulary.formatInternational(Vocabulary.Keys.TransformationAccuracy); final InternationalString eval = Resources .formatInternational(Resources.Keys.ConformanceMeansDatumShift); - DATUM_SHIFT_APPLIED = new PositionalAccuracyConstant(desc, eval, true); - DATUM_SHIFT_OMITTED = new PositionalAccuracyConstant(desc, eval, false); + DATUM_SHIFT_APPLIED = new PositionalAccuracyConstant(desc, eval, true); + DATUM_SHIFT_OMITTED = new PositionalAccuracyConstant(desc, eval, false); + INDIRECT_SHIFT_APPLIED = new PositionalAccuracyConstant(desc, eval, true); } /** @@ -122,8 +138,9 @@ public final class PositionalAccuracyConstant extends DefaultAbsoluteExternalPos * @throws ObjectStreamException if the serialized object defines an unknown data type. */ private Object readResolve() throws ObjectStreamException { - if (equals(DATUM_SHIFT_APPLIED)) return DATUM_SHIFT_APPLIED; - if (equals(DATUM_SHIFT_OMITTED)) return DATUM_SHIFT_OMITTED; + if (equals(DATUM_SHIFT_APPLIED)) return DATUM_SHIFT_APPLIED; + if (equals(DATUM_SHIFT_OMITTED)) return DATUM_SHIFT_OMITTED; + if (equals(INDIRECT_SHIFT_APPLIED)) return INDIRECT_SHIFT_APPLIED; return this; } @@ -194,11 +211,14 @@ public final class PositionalAccuracyConstant extends DefaultAbsoluteExternalPos * about the return values chosen. */ if (operation instanceof Transformation) { - if (accuracies.contains(DATUM_SHIFT_APPLIED)) { - return DATUM_SHIFT_ACCURACY; - } - if (accuracies.contains(DATUM_SHIFT_OMITTED)) { - return UNKNOWN_ACCURACY; + for (final PositionalAccuracy element : accuracies) { + /* + * Really need identity comparisons, not Object.equals(Object), because the later + * does not distinguish between DATUM_SHIFT_APPLIED and INDIRECT_SHIFT_APPLIED. + */ + if (element == DATUM_SHIFT_APPLIED) return DATUM_SHIFT_ACCURACY; + if (element == DATUM_SHIFT_OMITTED) return UNKNOWN_ACCURACY; + if (element == INDIRECT_SHIFT_APPLIED) return INDIRECT_SHIFT_ACCURACY; } } /* diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultGeodeticDatum.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultGeodeticDatum.java index 47cf530..4d72af4 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultGeodeticDatum.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultGeodeticDatum.java @@ -27,17 +27,20 @@ import org.opengis.util.GenericName; import org.opengis.util.InternationalString; import org.opengis.metadata.Identifier; import org.opengis.metadata.extent.Extent; +import org.opengis.metadata.extent.GeographicBoundingBox; import org.opengis.referencing.crs.GeodeticCRS; import org.opengis.referencing.datum.Ellipsoid; import org.opengis.referencing.datum.PrimeMeridian; import org.opengis.referencing.datum.GeodeticDatum; import org.opengis.referencing.operation.Matrix; import org.apache.sis.referencing.operation.matrix.Matrices; +import org.apache.sis.referencing.operation.matrix.MatrixSIS; import org.apache.sis.referencing.operation.matrix.NoninvertibleMatrixException; import org.apache.sis.metadata.iso.extent.Extents; import org.apache.sis.internal.referencing.WKTKeywords; import org.apache.sis.internal.metadata.NameToIdentifier; import org.apache.sis.internal.metadata.MetadataUtilities; +import org.apache.sis.internal.referencing.AnnotatedMatrix; import org.apache.sis.internal.referencing.CoordinateOperations; import org.apache.sis.internal.referencing.ExtentSelector; import org.apache.sis.internal.util.CollectionsExt; @@ -122,7 +125,7 @@ import static org.apache.sis.internal.referencing.WKTUtilities.toFormattable; * constants. * * @author Martin Desruisseaux (IRD, Geomatys) - * @version 0.7 + * @version 1.1 * * @see DefaultEllipsoid * @see DefaultPrimeMeridian @@ -375,7 +378,7 @@ public class DefaultGeodeticDatum extends AbstractDatum implements GeodeticDatum * 1053 – <cite>Time-dependent Position Vector transformation</cite>. * * <p>If this datum and the given {@code targetDatum} do not use the same {@linkplain #getPrimeMeridian() prime meridian}, - * then it is caller's responsibility to to apply longitude rotation before to use the matrix returned by this method. + * then it is caller's responsibility to apply longitude rotation before to use the matrix returned by this method. * The target prime meridian should be Greenwich (see {@linkplain #DefaultGeodeticDatum(Map, Ellipsoid, PrimeMeridian) * constructor javadoc}), in which case the datum shift should be applied in a geocentric coordinate system having * Greenwich as the prime meridian.</p> @@ -439,23 +442,47 @@ public class DefaultGeodeticDatum extends AbstractDatum implements GeodeticDatum Logging.unexpectedException(Logging.getLogger(Loggers.COORDINATE_OPERATION), DefaultGeodeticDatum.class, "getPositionVectorTransformation", e); } + /* + * No direct tranformation found. Search for a path through an intermediate datum. + * First, search if there is some BursaWolfParameters for the same target in both + * `source` and `target` datum. If such an intermediate is found, ask for path: + * + * source → [common datum] → target + * + * A consequence of such indirect path is that it may connect unrelated datums + * if [common datum] is a world datum such as WGS84. We do not have a solution + * for preventing that. + */ + if (bursaWolf != null) { + final GeographicBoundingBox bbox = selector.getAreaOfInterest(); + for (final BursaWolfParameters toPivot : bursaWolf) { + selector.setAreaOfInterest(bbox, toPivot.getDomainOfValidity()); + candidate = ((DefaultGeodeticDatum) targetDatum).select(toPivot.getTargetDatum(), selector); + if (candidate != null) { + final Matrix step1 = createTransformation(toPivot, areaOfInterest); + final Matrix step2 = createTransformation(candidate, areaOfInterest); + /* + * MatrixSIS.multiply(MatrixSIS) is equivalent to AffineTransform.concatenate(…): + * First transform by the supplied transform and then transform the result by the + * original transform. + */ + try { + Matrix m = MatrixSIS.castOrCopy(step2).inverse().multiply(step1); + return AnnotatedMatrix.indirect(m, selector.hasIntersection()); + } catch (NoninvertibleMatrixException e) { + Logging.unexpectedException(Logging.getLogger(Loggers.COORDINATE_OPERATION), + DefaultGeodeticDatum.class, "getPositionVectorTransformation", e); + } + } + } + } } - /* - * In a previous version, we were used to search for a transformation path through a common datum: - * - * source → [common datum] → target - * - * This has been removed, because it was dangerous (many paths may be possible - we are better to rely on - * the EPSG database, which do define some transformation paths explicitly). Especially since our javadoc - * now said that associating BursaWolfParameters to GeodeticDatum is not recommended except in a few special - * cases, this method does not have a picture complete enough for attempting anything else than a direct path. - */ return null; } /** * Invokes {@link BursaWolfParameters#getPositionVectorTransformation(Date)} for a date calculated from - * the temporal elements on the given extent. This method chooses an instant located midway between the + * the temporal elements on the given extent. This method chooses an instant located midway between the * start and end time. */ private static Matrix createTransformation(final BursaWolfParameters bursaWolf, final Extent areaOfInterest) { 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 0ca99e2..f31c746 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 @@ -32,10 +32,12 @@ import org.opengis.referencing.crs.*; import org.opengis.referencing.datum.*; import org.opengis.referencing.operation.*; import org.opengis.metadata.Identifier; +import org.opengis.metadata.quality.PositionalAccuracy; import org.opengis.metadata.extent.GeographicBoundingBox; import org.opengis.parameter.ParameterValueGroup; import org.opengis.parameter.ParameterDescriptorGroup; import org.apache.sis.internal.referencing.AxisDirections; +import org.apache.sis.internal.referencing.AnnotatedMatrix; import org.apache.sis.internal.referencing.CoordinateOperations; import org.apache.sis.internal.referencing.EllipsoidalHeightCombiner; import org.apache.sis.internal.referencing.ReferencingUtilities; @@ -685,7 +687,18 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry { transform = mtFactory.createConcatenatedTransform(transform, after); } } - return asList(createFromMathTransform(properties(identifier), sourceCRS, targetCRS, transform, method, parameters, null)); + /* + * Adjust the accuracy information if the datum shift has been computed by an indirect path. + * The indirect path uses a third datum (typically WGS84) as an intermediate between the two + * specified datum. + */ + final Map<String, Object> properties = properties(identifier); + if (datumShift instanceof AnnotatedMatrix) { + properties.put(CoordinateOperation.COORDINATE_OPERATION_ACCURACY_KEY, new PositionalAccuracy[] { + ((AnnotatedMatrix) datumShift).accuracy + }); + } + return asList(createFromMathTransform(properties, sourceCRS, targetCRS, transform, method, parameters, null)); } /** diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/datum/DefaultGeodeticDatumTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/datum/DefaultGeodeticDatumTest.java index b6ac4fc..f1dfa77 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/datum/DefaultGeodeticDatumTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/datum/DefaultGeodeticDatumTest.java @@ -27,6 +27,9 @@ import org.opengis.test.Validators; import org.apache.sis.measure.Units; import org.apache.sis.xml.Namespaces; import org.apache.sis.io.wkt.Convention; +import org.apache.sis.referencing.operation.matrix.Matrix4; +import org.apache.sis.internal.referencing.AnnotatedMatrix; +import org.apache.sis.internal.referencing.PositionalAccuracyConstant; import org.apache.sis.metadata.iso.extent.DefaultExtent; import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox; import org.apache.sis.test.xml.TestCase; @@ -43,7 +46,7 @@ import static org.apache.sis.referencing.GeodeticObjectVerifier.*; * Tests the {@link DefaultGeodeticDatum} class. * * @author Martin Desruisseaux (Geomatys) - * @version 0.8 + * @version 1.1 * @since 0.4 * @module */ @@ -185,6 +188,45 @@ public final strictfp class DefaultGeodeticDatumTest extends TestCase { } /** + * Tests {@link DefaultGeodeticDatum#getPositionVectorTransformation(GeodeticDatum, Extent)} + * going through an indirect transformation. The main purpose of this test is to verify that + * the matrix is associated with {@link PositionalAccuracyConstant#INDIRECT_SHIFT_APPLIED}. + */ + @Test + @DependsOnMethod("testGetPositionVectorTransformation") + public void testIndirectTransformation() { + final Map<String,Object> properties = new HashMap<>(); + assertNull(properties.put(DefaultGeodeticDatum.NAME_KEY, "Invalid dummy datum")); + assertNull(properties.put(DefaultGeodeticDatum.BURSA_WOLF_KEY, BursaWolfParametersTest.createWGS72_to_WGS84())); + final DefaultGeodeticDatum global = new DefaultGeodeticDatum(properties, + GeodeticDatumMock.WGS72.getEllipsoid(), + GeodeticDatumMock.WGS72.getPrimeMeridian()); + /* + * Create a datum valid only in a specific region of the world and with no direct transformation to WGS72. + * However an indirect transformation to WGS72 is available through WGS84. + */ + properties.put(DefaultGeodeticDatum.BURSA_WOLF_KEY, BursaWolfParametersTest.createED87_to_WGS84()); + final DefaultGeodeticDatum local = new DefaultGeodeticDatum(properties, + GeodeticDatumMock.ED50.getEllipsoid(), + GeodeticDatumMock.ED50.getPrimeMeridian()); + /* + * Main test: verify that the transformation found is associated with accuracy information. + */ + final Matrix m = local.getPositionVectorTransformation(global, null); + assertInstanceOf("Should have accuracy information.", AnnotatedMatrix.class, m); + assertSame(PositionalAccuracyConstant.INDIRECT_SHIFT_APPLIED, ((AnnotatedMatrix) m).accuracy); + /* + * Following is an anti-regression test only (no authoritative values). + * Verified only opportunistically. + */ + assertMatrixEquals("getPositionVectorTransformation", new Matrix4( + 1, 7.961E-7, 7.287E-7, -82.981, + -7.961E-7, 1, 2.461E-6, -99.719, + -7.287E-7, -2.461E-6, 1, -115.209, + 0, 0, 0, 1), m, 0.01); + } + + /** * Tests {@link DefaultGeodeticDatum#toWKT()}. */ @Test diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java index b381f62..3ab3dea 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java @@ -38,12 +38,14 @@ import org.opengis.referencing.operation.Transformation; import org.opengis.referencing.operation.TransformException; import org.opengis.referencing.operation.ConcatenatedOperation; import org.opengis.referencing.operation.OperationNotFoundException; +import org.apache.sis.internal.referencing.PositionalAccuracyConstant; import org.apache.sis.referencing.operation.transform.LinearTransform; import org.apache.sis.referencing.operation.transform.MathTransforms; import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.crs.DefaultCompoundCRS; import org.apache.sis.referencing.crs.DefaultDerivedCRS; import org.apache.sis.referencing.CommonCRS; +import org.apache.sis.referencing.CRS; import org.apache.sis.io.wkt.WKTFormat; import static org.apache.sis.internal.referencing.Formulas.LINEAR_TOLERANCE; @@ -70,7 +72,7 @@ import static org.apache.sis.test.Assert.*; * Contrarily to {@link CoordinateOperationRegistryTest}, tests in this class are run without EPSG geodetic dataset. * * @author Martin Desruisseaux (Geomatys) - * @version 1.0 + * @version 1.1 * @since 0.7 * @module */ @@ -529,7 +531,7 @@ public final strictfp class CoordinateOperationFinderTest extends MathTransformT * @see <a href="https://issues.apache.org/jira/browse/SIS-364";>SIS-364</a> */ static String AGD66() { - return "PROJCS[“AGD66 / AMG zone 49”," + return "PROJCS[“AGD66 / AMG zone 49”, " + "GEOGCS[“AGD66”, " + "DATUM[“Australian_Geodetic_Datum_1966”, " + "SPHEROID[“Australian National Spheroid”,6378160, 298.25, AUTHORITY[“EPSG”,“7003”]]," @@ -550,6 +552,74 @@ public final strictfp class CoordinateOperationFinderTest extends MathTransformT } /** + * Tests a transformation between two CRS for which no direct bursa-wolf parameters are defined. + * However a transformation should still be possible indirectly, through WGS 84. + * + * @throws ParseException if a CRS used in this test can not be parsed. + * @throws FactoryException if the operation can not be created. + * @throws TransformException if an error occurred while converting the test points. + */ + @Test + public void testIndirectDatumShift() throws ParseException, FactoryException, TransformException { + final CoordinateReferenceSystem sourceCRS = parse( + "PROJCS[“RGF93 / Lambert-93”, " + + "GEOGCS[“RGF93”, " + + "DATUM[“Reseau Geodesique Francais 1993”, " + + "SPHEROID[“GRS 1980”, 6378137, 298.257222101], " + + "TOWGS84[0,0,0,0,0,0,0]], " + + "PRIMEM[“Greenwich”,0], " + + "UNIT[“degree”, 0.0174532925199433]], " + + "PROJECTION[“Lambert_Conformal_Conic_2SP”], " + + "PARAMETER[“standard_parallel_1”, 49], " + + "PARAMETER[“standard_parallel_2”, 44], " + + "PARAMETER[“latitude_of_origin”, 46.5], " + + "PARAMETER[“central_meridian”, 3], " + + "PARAMETER[“false_easting”, 700000], " + + "PARAMETER[“false_northing”, 6600000], " + + "UNIT[“metre”,1], " + + "AUTHORITY[“EPSG”,“2154”]]"); + + final CoordinateReferenceSystem targetCRS = parse( + "PROJCS[“Amersfoort / RD New”, " + + "GEOGCS[“Amersfoort”, " + + "DATUM[“Amersfoort”, " + + "SPHEROID[“Bessel 1841”, 6377397.155, 299.1528128], " + + "TOWGS84[565.417, 50.3319, 465.552, -0.398957, 0.343988, -1.8774, 4.0725]], " + + "PRIMEM[“Greenwich”,0], " + + "UNIT[“degree”,0.0174532925199433]], " + + "PROJECTION[“Oblique_Stereographic”], " + + "PARAMETER[“latitude_of_origin”, 52.15616055555555], " + + "PARAMETER[“central_meridian”, 5.38763888888889], " + + "PARAMETER[“scale_factor”, 0.9999079], " + + "PARAMETER[“false_easting”, 155000], " + + "PARAMETER[“false_northing”, 463000], " + + "UNIT[“metre”,1], " + + "AUTHORITY[“EPSG”,“28992”]]"); + /* + * Transform a point as a way to verify that a datum shift is applied. + * If no datum shift is applied, the point will be at 191 metres from + * expected value. + */ + final CoordinateOperation operation = finder.createOperation(sourceCRS, targetCRS); + tolerance = LINEAR_TOLERANCE; + transform = operation.getMathTransform(); + verifyTransform(new double[] {926713.702, 7348947.026}, + new double[] {220798.684, 577583.801}); // With datum shift through WGS84. + // 220762.487, 577396.040 // Without datum shift. + validate(); + /* + * The accuracy should tell that the datum shift is indirect (through WGS 84). + * However the value may differ depending on whether EPSG database has been + * used or not, because it depends on whether the datum have been completed + * with domain of validity. + */ + final double accuracy = CRS.getLinearAccuracy(operation); + if (accuracy != PositionalAccuracyConstant.UNKNOWN_ACCURACY) { + assertEquals("accuracy", PositionalAccuracyConstant.INDIRECT_SHIFT_ACCURACY, accuracy, STRICT); + } + } + + /** * Tests that an exception is thrown on attempt to grab a transformation between incompatible vertical CRS. * * @throws FactoryException if an exception other than the expected one occurred.
