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(
|