sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch master updated: 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 cherry-pick from afdf695c15841e6f6d2ad23c24833111b4d0331e.
Date Mon, 31 Aug 2020 18:06:42 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/master by this push:
     new 902c9a0  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 cherry-pick from
afdf695c15841e6f6d2ad23c24833111b4d0331e.
902c9a0 is described below

commit 902c9a0df26a93457c9c588db296f021643fccd9
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Aug 31 20:05:26 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 cherry-pick from afdf695c15841e6f6d2ad23c24833111b4d0331e.
    
    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 7c6f601..eed95ef 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.extent.Extent;
 import org.opengis.referencing.ReferenceIdentifier;
+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 dd64de3..d0c7d77 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;
@@ -680,7 +682,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 c3fb178..add570b 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
  */
@@ -530,7 +532,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”]],"
@@ -551,6 +553,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.


Mime
View raw message