sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/03: GridDerivation.subgrid(GridGeometry) should support coverage spanning the anti-meridian.
Date Tue, 15 Sep 2020 12:51:56 GMT
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 c34da9e127731686de41ce3487deb64f2817fdfb
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Tue Sep 15 14:49:27 2020 +0200

    GridDerivation.subgrid(GridGeometry) should support coverage spanning the anti-meridian.
---
 .../coverage/grid/CoordinateOperationFinder.java   | 23 +++++++++-
 .../apache/sis/coverage/grid/GridDerivation.java   | 34 ++------------
 .../org/apache/sis/coverage/grid/GridGeometry.java |  4 ++
 .../sis/coverage/grid/GridDerivationTest.java      | 53 +++++++++++++++++++++-
 .../internal/referencing/WraparoundTransform.java  | 21 +++++++--
 .../sis/internal/storage/query/CoverageSubset.java |  3 --
 .../internal/storage/query/GridResourceMock.java   |  6 +--
 7 files changed, 100 insertions(+), 44 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
index fea754d..c92cc6e 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
@@ -67,7 +67,8 @@ import org.apache.sis.internal.system.Modules;
  *   <li>Use the area of interest and grid resolution for refining the coordinate operation
between two CRS.</li>
  * </ul>
  *
- * <b>Note:</b> this class does not provide the complete chain of operations
from grid to grid.
+ * <b>Note:</b> except for the {@link #gridToGrid(PixelInCell)} convenience method,
+ * this class does not provide the complete chain of operations from grid to grid.
  * It provides only the operation from <em>cell indices</em> in source grid to
<em>coordinates in the CRS</em>
  * of destination grid. Callers must add the last step (conversion from target CRS to cell
indices) themselves.
  *
@@ -197,6 +198,26 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     }
 
     /**
+     * Computes the transform from “grid coordinates of the source” to “grid coordinates
of the target”.
+     * This is a concatenation of {@link #gridToCRS(PixelInCell)} with target "CRS to grid"
transform.
+     *
+     * @param  anchor  whether the operation is between cell centers or cell corners.
+     * @return operation from source grid indices to target grid indices.
+     * @throws FactoryException if no operation can be found between the source and target
CRS.
+     * @throws TransformException if some coordinates can not be transformed to the specified
target.
+     * @throws IncompleteGridGeometryException if required CRS or a "grid to CRS" information
is missing.
+     */
+    final MathTransform gridToGrid(final PixelInCell anchor) throws FactoryException, TransformException
{
+        final MathTransform step1 = gridToCRS(anchor);
+        final MathTransform step2 = target.getGridToCRS(anchor);
+        if (step1.equals(step2)) {                                          // Optimization
for a common case.
+            return MathTransforms.identity(step1.getSourceDimensions());
+        } else {
+            return MathTransforms.concatenate(step1, step2.inverse());
+        }
+    }
+
+    /**
      * Computes the transform from “grid coordinates of the source” to “geospatial
coordinates of the target”.
      * It may be the identity operation. We try to take envelopes in account because the
operation choice may
      * depend on the geographic area.
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
index 30ccc29..52a5856 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
@@ -28,7 +28,6 @@ import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.operation.CoordinateOperation;
-import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.referencing.operation.transform.TransformSeparator;
@@ -249,7 +248,7 @@ public class GridDerivation {
      * @throws IllegalStateException if {@link #subgrid(Envelope, double...)} or {@link #slice(DirectPosition)}
      *         has already been invoked.
      *
-     * @see GridExtent#grow(boolean, boolean, long...)
+     * @see GridExtent#expand(long...)
      */
     public GridDerivation margin(final int... cellCounts) {
         ArgumentChecks.ensureNonNull("cellCounts", cellCounts);
@@ -426,10 +425,9 @@ public class GridDerivation {
             final MathTransform mapCorners, mapCenters;
             final GridExtent domain = gridOfInterest.getExtent();                  // May
throw IncompleteGridGeometryException.
             try {
-                final CoordinateOperation crsChange;
-                crsChange  = Envelopes.findOperation(gridOfInterest.envelope, base.envelope);
      // Any envelope may be null.
-                mapCorners = path(gridOfInterest, crsChange, base, PixelInCell.CELL_CORNER);
-                mapCenters = path(gridOfInterest, crsChange, base, PixelInCell.CELL_CENTER);
+                final CoordinateOperationFinder finder = new CoordinateOperationFinder(gridOfInterest,
base);
+                mapCorners = finder.gridToGrid(PixelInCell.CELL_CORNER);
+                mapCenters = finder.gridToGrid(PixelInCell.CELL_CENTER);
                 clipExtent(domain.toCRS(mapCorners, mapCenters, null));
             } catch (FactoryException | TransformException e) {
                 throw new IllegalGridGeometryException(e, "gridOfInterest");
@@ -452,30 +450,6 @@ public class GridDerivation {
     }
 
     /**
-     * Returns the concatenation of all transformation steps from the given source to the
given target.
-     * The transform maps grid coordinates (not envelopes).
-     *
-     * @param  source     the source grid geometry.
-     * @param  crsChange  the change of coordinate reference system, or {@code null} if none.
-     * @param  target     the target grid geometry.
-     * @param  anchor     whether we want the transform for cell corner or cell center.
-     */
-    private static MathTransform path(final GridGeometry source, final CoordinateOperation
crsChange,
-            final GridGeometry target, final PixelInCell anchor) throws NoninvertibleTransformException
-    {
-        MathTransform step1 = source.getGridToCRS(anchor);
-        MathTransform step2 = target.getGridToCRS(anchor);
-        if (crsChange != null) {
-            step1 = MathTransforms.concatenate(step1, crsChange.getMathTransform());
-        }
-        if (step1.equals(step2)) {                                          // Optimization
for a common case.
-            return MathTransforms.identity(step1.getSourceDimensions());
-        } else {
-            return MathTransforms.concatenate(step1, step2.inverse());
-        }
-    }
-
-    /**
      * Requests a grid geometry over a sub-region of the base grid geometry and optionally
with subsampling.
      * The given envelope does not need to be expressed in the same coordinate reference
system (CRS)
      * than {@linkplain GridGeometry#getCoordinateReferenceSystem() the CRS of the base grid
geometry};
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
index 0823ead..bdecc9a 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
@@ -1216,6 +1216,10 @@ public class GridGeometry implements LenientComparable, Serializable
{
     public MathTransform createTransformTo(final GridGeometry target, final PixelInCell anchor)
throws TransformException {
         ArgumentChecks.ensureNonNull("target", target);
         ArgumentChecks.ensureNonNull("anchor", anchor);
+        /*
+         * Inverse `source` and `target` because `CoordinateOperationFinder.inverse(…)`
does an
+         * effort for using `WraparoundTransform` only if needed (contrarily to `gridToCRS(…)`).
+         */
         final CoordinateOperationFinder finder = new CoordinateOperationFinder(target, this);
         MathTransform tr;
         try {
diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
index d84e827..bc58662 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
@@ -26,12 +26,15 @@ import org.opengis.referencing.datum.PixelInCell;
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.geometry.Envelope2D;
+import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+import org.apache.sis.referencing.cs.AxesConvention;
 import org.apache.sis.referencing.crs.DefaultCompoundCRS;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.matrix.Matrix3;
 import org.apache.sis.referencing.operation.matrix.Matrix4;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.HardCodedConversions;
 import org.apache.sis.referencing.crs.HardCodedCRS;
 import org.apache.sis.geometry.DirectPosition2D;
 import org.apache.sis.geometry.GeneralDirectPosition;
@@ -392,7 +395,9 @@ public final strictfp class GridDerivationTest extends TestCase {
     }
 
     /**
-     * Tests deriving a grid geometry with an envelope crossing the antimeridian.
+     * Tests deriving a grid geometry with a request crossing the antimeridian.
+     * The {@link GridGeometry} crossing the anti-meridian is the one given in
+     * argument to {@link GridDerivation#subgrid(GridGeometry)}.
      */
     @Test
     public void testSubgridCrossingAntiMeridian() {
@@ -414,6 +419,52 @@ public final strictfp class GridDerivationTest extends TestCase {
     }
 
     /**
+     * Tests deriving a grid geometry from a grid crossing the antimeridian.
+     * The {@link GridGeometry} crossing the anti-meridian is the one given
+     * to {@link GridDerivation} constructor.
+     */
+    @Test
+    public void testBasegridCrossingAntiMeridian() {
+        /*
+         * Longitudes from 100°E to 240°E (in WGS84 geographic CRS), which is equivalent
to 100°E to 120°W.
+         * That [100 … 240]°E range is compatible with the [0 … 360]° longitude range
declared in the CRS.
+         * Latitude range is from 21°S to 60°N, but this is secondary for this test.
+         */
+        final GridGeometry grid = new GridGeometry(
+                new GridExtent(null, null, new long[] {8400, 4860}, true), PixelInCell.CELL_CENTER,
+                MathTransforms.linear(new Matrix3(
+                        0.016664682775860015,  0,  100.00833234138793,
+                    0,  0.016663238016868958,      -20.991668380991566,
+                    0, 0, 1)),
+                HardCodedCRS.WGS84.forConvention(AxesConvention.POSITIVE_RANGE));
+        /*
+         * 180°W to 180″E (the world) and 80°S to 80°N in Mercator projection.
+         * Only the longitudes are of interest for this test.
+         */
+        final GridGeometry areaOfInterest = new GridGeometry(
+                new GridExtent(null, null, new long[] {256, 256}, false), PixelInCell.CELL_CORNER,
+                MathTransforms.linear(new Matrix3(
+                        156543.03392804097,  0,  -2.0037508342789244E7,
+                    0, -121066.95890409155,       1.5496570739723722E7,
+                    0, 0, 1)),
+                HardCodedConversions.mercator());
+        /*
+         * Since the area of interest covers the world (in longitude), the intersection should
+         * have the same range of longitudes than the source grid, which is [100 … 240]°E.
+         * Latitude range is also unchanged: 21°S to 60°N.
+         */
+        final GridGeometry result = grid.derive().subgrid(areaOfInterest).build();
+        assertEnvelopeEquals(new Envelope2D(null, 100, -21, 240 - 100, 60 - -21),
+                             result.getEnvelope(), Formulas.ANGULAR_TOLERANCE);
+        /*
+         * Following is empirical values provided as an anti-regression test.
+         * The longitude span is about (240° − 100°) ÷ 360° × 256 px  ≈  99.56 px.
+         * The latitude span is more complicated because of adjustement for resolution.
+         */
+        assertExtentEquals(new long[2], new long[] {100, 73}, result.getExtent());
+    }
+
+    /**
      * Tests deriving a grid geometry from an area of interest shifted by 360° before or
after the grid valid area.
      */
     @Test
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
index 95b6e6d..19411d6 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
@@ -319,14 +319,27 @@ public final class WraparoundTransform extends AbstractMathTransform
{
     }
 
     /**
+     * Whether to allow discontinuities in transform derivatives. Strictly speaking the derivative
+     * at 1, 2, 3… should be infinite because the coordinate value jumps "instantaneously"
from an
+     * integer value to 0.  However in practice we use derivatives as linear approximations
around
+     * small regions, not for calculations requiring strict mathematical values. An infinite
value
+     * goes against the approximation goal.  Furthermore whether a source coordinate is an
integer
+     * value or not is subject to rounding errors, which may cause unpredictable behavior
if
+     * discontinuities are allowed.
+     */
+    private static final boolean ALLOW_DISCONTINUITIES = false;
+
+    /**
      * Gets the derivative of this transform at a point.
      */
     @Override
     public Matrix derivative(final DirectPosition point) {
         final MatrixSIS derivative = Matrices.createIdentity(dimension);
-        final double v = point.getOrdinate(wraparoundDimension);
-        if (v == Math.floor(v)) {
-            derivative.setElement(wraparoundDimension, wraparoundDimension, Double.NEGATIVE_INFINITY);
+        if (ALLOW_DISCONTINUITIES) {
+            final double v = point.getOrdinate(wraparoundDimension);
+            if (v == Math.floor(v)) {
+                derivative.setElement(wraparoundDimension, wraparoundDimension, Double.NEGATIVE_INFINITY);
+            }
         }
         return derivative;
     }
@@ -349,7 +362,7 @@ public final class WraparoundTransform extends AbstractMathTransform {
             return null;
         }
         final MatrixSIS derivative = Matrices.createIdentity(dimension);
-        if (v == 0) {
+        if (ALLOW_DISCONTINUITIES && v == 0) {
             derivative.setElement(wraparoundDimension, wraparoundDimension, Double.NEGATIVE_INFINITY);
         }
         return derivative;
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/query/CoverageSubset.java
b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/query/CoverageSubset.java
index 3926441..4795275 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/query/CoverageSubset.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/query/CoverageSubset.java
@@ -39,9 +39,6 @@ import org.apache.sis.internal.util.UnmodifiableArrayList;
  * This implementation merges the domain and range specified by the query with
  * arguments of {@link GridCoverageResource#read(GridGeometry, int...)} method.
  *
- * <p>Current version does not yet support non-zero value for
- * {@link CoverageQuery#getSourceDomainExpansion()}.</p>
- *
  * @author  Martin Desruisseaux (Geomatys)
  * @version 1.1
  * @since   1.1
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/internal/storage/query/GridResourceMock.java
b/storage/sis-storage/src/test/java/org/apache/sis/internal/storage/query/GridResourceMock.java
index bb98f4a..b984874 100644
--- a/storage/sis-storage/src/test/java/org/apache/sis/internal/storage/query/GridResourceMock.java
+++ b/storage/sis-storage/src/test/java/org/apache/sis/internal/storage/query/GridResourceMock.java
@@ -24,7 +24,6 @@ import org.apache.sis.coverage.grid.GridCoverage;
 import org.apache.sis.coverage.grid.GridCoverage2D;
 import org.apache.sis.coverage.grid.GridExtent;
 import org.apache.sis.coverage.grid.GridGeometry;
-import org.apache.sis.coverage.grid.GridRoundingMode;
 import org.apache.sis.internal.storage.AbstractGridResource;
 
 import static org.junit.Assert.*;
@@ -94,10 +93,7 @@ final strictfp class GridResourceMock extends AbstractGridResource {
         if (domain == null) {
             domain = gridGeometry;
         } else {
-            domain = gridGeometry.derive()
-                    .rounding(GridRoundingMode.ENCLOSING)
-                    .subgrid(domain)
-                    .build();
+            domain = gridGeometry.derive().subgrid(domain).build();
         }
         final GridExtent extent = domain.getExtent();
         final BufferedImage img = new BufferedImage(


Mime
View raw message