sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Add tests and a GridGeometry.subgrid(areaOfInterest, targetResolution) method.
Date Fri, 28 Dec 2018 19:05:25 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


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new a5ad10a  Add tests and a GridGeometry.subgrid(areaOfInterest, targetResolution) method.
a5ad10a is described below

commit a5ad10ae13a984689b4550d56de6f12f19d9233c
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri Dec 28 20:04:58 2018 +0100

    Add tests and a GridGeometry.subgrid(areaOfInterest, targetResolution) method.
---
 .../org/apache/sis/coverage/grid/GridChange.java   |  84 +++++++++---
 .../org/apache/sis/coverage/grid/GridExtent.java   |  61 ++++++---
 .../org/apache/sis/coverage/grid/GridGeometry.java | 148 +++++++++++++++++----
 .../apache/sis/coverage/grid/PixelTranslation.java |   4 +-
 .../apache/sis/coverage/grid/GridChangeTest.java   | 104 +++++++++++++++
 .../apache/sis/coverage/grid/GridExtentTest.java   |  32 +++--
 .../apache/sis/coverage/grid/GridGeometryTest.java |  58 ++++++++
 .../org/apache/sis/test/suite/RasterTestSuite.java |   1 +
 8 files changed, 422 insertions(+), 70 deletions(-)

diff --git a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
index 4b42eb4..7c5faf1 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
@@ -29,10 +29,15 @@ import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.geometry.Envelopes;
+import org.apache.sis.util.collection.DefaultTreeTable;
+import org.apache.sis.util.collection.TableColumn;
+import org.apache.sis.util.collection.TreeTable;
 import org.apache.sis.util.resources.Vocabulary;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.util.CharSequences;
 import org.apache.sis.util.Classes;
+import org.apache.sis.util.Debug;
 
 
 /**
@@ -71,7 +76,7 @@ public class GridChange implements Serializable {
     private static final long serialVersionUID = -7819047186271172885L;
 
     /**
-     * The conversions from source grid to target grid, mapping pixel corners or pixel centers.
+     * The conversions from source grid to target grid, mapping cell corners or cell centers.
      */
     private final MathTransform mapCorners, mapCenters;
 
@@ -177,8 +182,8 @@ public class GridChange implements Serializable {
                 throw new TransformException(Errors.format(Errors.Keys.CanNotTransformEnvelope), e);
             }
             final GridExtent domain = source.getExtent();                             // May throw IncompleteGridGeometryException.
-            mapCorners  = path(source, crsChange, target, PixelInCell.CELL_CORNER);
-            mapCenters  = path(source, crsChange, target, PixelInCell.CELL_CENTER);
+            mapCorners = path(source, crsChange, target, PixelInCell.CELL_CORNER);
+            mapCenters = path(source, crsChange, target, PixelInCell.CELL_CENTER);
             final GridExtent extent = new GridExtent(domain.toCRS(mapCorners, mapCenters), rounding, margin, target.extent, null);
             targetExtent = extent.equals(target.extent) ? target.extent : extent.equals(domain) ? domain : extent;
             /*
@@ -203,7 +208,7 @@ public class GridChange implements Serializable {
      * @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 pixel corner or pixel center.
+     * @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
@@ -321,8 +326,9 @@ public class GridChange implements Serializable {
      * {@link #getTargetExtent()} as below for each dimension <var>i</var>:
      *
      * <ul>
-     *   <li>The {@linkplain GridExtent#getLow(int) low} coordinate is unchanged.</li>
-     *   <li>The {@linkplain GridExtent#getSize(int) size} is divided by {@code strides[i]}.</li>
+     *   <li>The {@linkplain GridExtent#getLow(int)  low}  is divided by {@code strides[i]}, rounded toward zero.</li>
+     *   <li>The {@linkplain GridExtent#getSize(int) size} is divided by {@code strides[i]}, rounded toward zero.</li>
+     *   <li>The {@linkplain GridExtent#getHigh(int) high} is recomputed from above low and size.</li>
      * </ul>
      *
      * The {@linkplain GridGeometry#getGridToCRS(PixelInCell) grid to CRS} transform is scaled accordingly
@@ -338,6 +344,8 @@ public class GridChange implements Serializable {
         double[] factors = null;
         Matrix toGiven = null;
         if (strides != null) {
+            // Validity of the strides values will be verified by GridExtent constructor below.
+            final GridExtent unscaled = extent;
             final int dimension = extent.getDimension();
             for (int i = Math.min(dimension, strides.length); --i >= 0;) {
                 final int s = strides[i];
@@ -350,10 +358,12 @@ public class GridChange implements Serializable {
                             toGiven = Matrices.createIdentity(dimension + 1);
                         }
                     }
-                    factors[i] = s;
+                    final double sd = s;
+                    factors[i] = sd;
                     if (toGiven != null) {
-                        toGiven.setElement(i, i, s);
-                        toGiven.setElement(i, dimension, (1-s) * extent.getLow(i));
+                        final long low = unscaled.getLow(i);
+                        toGiven.setElement(i, i, sd);
+                        toGiven.setElement(i, dimension, low - (low / s) * sd);     // (low / s) must be consistent with GridExtent.
                     }
                 }
             }
@@ -406,6 +416,54 @@ public class GridChange implements Serializable {
     }
 
     /**
+     * Returns a tree representation of this grid change. The tree representation
+     * is for debugging purpose only and may change in any future SIS version.
+     *
+     * @param  locale  the locale to use for textual labels.
+     * @return a tree representation of this grid change.
+     */
+    @Debug
+    private TreeTable toTree(final Locale locale) {
+        final TableColumn<CharSequence> column = TableColumn.VALUE_AS_TEXT;
+        final TreeTable tree = new DefaultTreeTable(column);
+        final TreeTable.Node root = tree.getRoot();
+        root.setValue(column, Classes.getShortClassName(this));
+        /*
+         * GridChange (example)
+         *   └─Target range
+         *       ├─Dimension 0: [ 2000 … 5475] (3476 cells)
+         *       └─Dimension 1: [-1000 … 7999] (9000 cells)
+         */
+        TreeTable.Node section = root.newChild();
+        section.setValue(column, "Target range");
+        final StringBuilder buffer = new StringBuilder(256);
+        getTargetExtent().appendTo(buffer, Vocabulary.getResources(locale));
+        for (final CharSequence line : CharSequences.splitOnEOL(buffer)) {
+            String text = line.toString().trim();
+            if (!text.isEmpty()) {
+                section.newChild().setValue(column, text);
+            }
+        }
+        /*
+         * GridChange (example)
+         *   └─Target strides
+         *       ├─{50, 300}
+         *       └─Global ≈ 175.0
+         */
+        buffer.setLength(0);
+        buffer.append('{');
+        for (int s : getTargetStrides()) {
+            if (buffer.length() > 1) buffer.append(", ");
+            buffer.append(s);
+        }
+        section = root.newChild();
+        section.setValue(column, "Target strides");
+        section.newChild().setValue(column, buffer.append('}').toString()); buffer.setLength(0);
+        section.newChild().setValue(column, buffer.append("Global ≈ ").append((float) getGlobalScale()).toString());
+        return tree;
+    }
+
+    /**
      * Returns a string representation of this grid change for debugging purpose.
      * The returned string is implementation dependent and may change in any future version.
      *
@@ -413,12 +471,6 @@ public class GridChange implements Serializable {
      */
     @Override
     public String toString() {
-        final String lineSeparator = System.lineSeparator();
-        final StringBuilder buffer = new StringBuilder(256)
-                .append(Classes.getShortClassName(this)).append(lineSeparator)
-                .append("└ Scale factor ≈ ").append((float) getGlobalScale()).append(lineSeparator)
-                .append("Target range").append(lineSeparator);
-        targetExtent.appendTo(buffer, Vocabulary.getResources((Locale) null));
-        return buffer.toString();
+        return toTree(null).toString();
     }
 }
diff --git a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
index f9f823c..0e17166 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
@@ -211,10 +211,38 @@ public class GridExtent implements Serializable {
     }
 
     /**
-     * Creates a new grid extent with the same lower values than the given extent, and high values adjusted by dividing
-     * the {@linkplain #getSize(int) size} by the given strides.  The policy of keeping the lower coordinates unchanged
-     * is consistent with {@link GridChange#getTargetGeometry(int...)} policy of keeping lower coordinates invariant
-     * under the scales.
+     * Creates a new grid extent with low and high values adjusted by dividing the {@linkplain #getSize(int) size}
+     * by the given strides. The policy of dividing the lower coordinates by the stride shall be kept consistent
+     * with {@link GridGeometry#subgrid(Envelope, double...)} and {@link GridChange#getTargetGeometry(int...)}
+     * computation of grid to CRS.
+     *
+     * <div class="note"><b>Note:</b>
+     * if a division does not produce an integer, then the size is rounded toward 0 (or toward negative infinity since
+     * sizes are always positive numbers). But the envelope computed by {@link GridChange#getTargetGeometry(int...)}
+     * will nevertheless be larger. This counter-intuitive effect is because the "grid to CRS" transform is scaled by
+     * the same factors, which results in larger cells. Assuming:
+     *
+     * <ul>
+     *   <li>All {@linkplain #getLow(int) low coordinates} = 0.</li>
+     *   <li>{@linkplain #getHigh(int) High coordinates} before scaling denoted <var>h₁</var>.</li>
+     *   <li>High coordinates after scaling denoted <var>h₂</var>.</li>
+     *   <li>Scale factor (or sub-sampling) <var>s</var> is an integer ≧ 1.</li>
+     * </ul>
+     *
+     * Then the envelope upper bounds <var>x</var> is:
+     *
+     * <ul>
+     *   <li>x = (h₁ + 1) × c</li>
+     *   <li>x = (h₂ + f) × c⋅s   which implies   h₂ = h₁/s   and   f = 1/s</li>
+     * </ul>
+     *
+     * If we modify the later equation for integer division instead than real numbers, we have:
+     *
+     * <blockquote>x = (h₂ + f) × c⋅s   where   h₂ = floor(h₁/s)   and   f = ((h₁ mod s) + 1)/s</blockquote>
+     *
+     * Because <var>s</var> ≧ 1, then <var>f</var> ≦ 1. But the <var>f</var> actually used by {@link #toCRS(MathTransform,
+     * MathTransform)} is hard-coded to 1 since it assumes that all cells are whole, i.e. it does not take in account that
+     * the last cell may actually be fraction of a cell. Since 1 ≧ <var>f</var>, the computed envelope may be larger.</div>
      *
      * @param  extent   the extent from which to compute a new extent.
      * @param  strides  the strides. Length shall be equal to the number of dimension and all values shall be greater than zero.
@@ -222,19 +250,21 @@ public class GridExtent implements Serializable {
      */
     GridExtent(final GridExtent extent, final int[] strides) {
         types = extent.types;
-        coordinates = extent.coordinates.clone();
-        final int m = getDimension();
+        final int m = extent.getDimension();
+        coordinates = new long[m << 1];
         for (int i=0; i<m; i++) {
             final int s = strides[i];
             if (s <= 0) {
                 throw new IllegalArgumentException(Errors.format(Errors.Keys.ValueNotGreaterThanZero_2, "strides[" + i + ']', s));
             }
             final int j = i + m;
-            final long size = coordinates[j] - coordinates[i] + 1;                             // Result is an unsigned number.
+            long low  = extent.coordinates[i];
+            long size = extent.coordinates[j] - low + 1;                                       // Result is an unsigned number.
             if (size == 0) throw new ArithmeticException("long overflow");
             long r = Long.divideUnsigned(size, s);
             if (r*s == size) r--;                           // Make inclusive if the division did not already rounded toward 0.
-            coordinates[j] = coordinates[i] + r;
+            coordinates[i] = low /= s;
+            coordinates[j] = low + r;
         }
     }
 
@@ -627,11 +657,11 @@ public class GridExtent implements Serializable {
 
     /**
      * Transforms this grid extent to a "real world" envelope using the given transform.
-     * The transform shall map <em>pixel corner</em> to real world coordinates.
+     * The transform shall map <em>cell corner</em> to real world coordinates.
      *
-     * @param  cornerToCRS  a transform from <em>pixel corners</em> to real world coordinates.
+     * @param  cornerToCRS  a transform from <em>cell corners</em> to real world coordinates.
      * @param  gridToCRS    the transform specified by the user. May be the same as {@code cornerToCRS}.
-     *                      If different, then this is assumed to map pixel centers instead than pixel corners.
+     *                      If different, then this is assumed to map cell centers instead than cell corners.
      * @return this grid extent in real world coordinates.
      *
      * @see #GridExtent(AbstractEnvelope, GridRoundingMode, int[], GridExtent, int[])
@@ -640,6 +670,7 @@ public class GridExtent implements Serializable {
         final int dimension = getDimension();
         GeneralEnvelope envelope = new GeneralEnvelope(dimension);
         for (int i=0; i<dimension; i++) {
+            // The +1.0 is for making high coordinate exclusive. See GridExtent(GridExtent, int...) for a discussion.
             envelope.setRange(i, coordinates[i], coordinates[i + dimension] + 1.0);
         }
         envelope = Envelopes.transform(cornerToCRS, envelope);
@@ -647,8 +678,8 @@ public class GridExtent implements Serializable {
             /*
              * If the envelope contains some NaN values, try to replace them by constant values inferred from the math transform.
              * We must use the MathTransform specified by the user (gridToCRS), not necessarily the cornerToCRS transform, because
-             * because inferring a 'cornerToCRS' by translation a 'centertoCRS' by 0.5 pixels increase the amount of NaN values in
-             * the matrix. For giving a chance to TransformSeparator to perform its work, we need the minimal amount of NaN values.
+             * inferring a 'cornerToCRS' by translating a 'centerToCRS' by 0.5 cell increase the amount of NaN values in the matrix.
+             * For giving a chance to TransformSeparator to perform its work, we need the minimal amount of NaN values.
              */
             final boolean isCenter = (gridToCRS != cornerToCRS);
             TransformSeparator separator = null;
@@ -674,8 +705,8 @@ public class GridExtent implements Serializable {
                             final double value = component.getElement(j, component.getNumCol() - 1);
                             /*
                              * Replace only the envelope NaN values by the translation term (non-NaN values are left unchanged).
-                             * If the gridToCRS map pixel corners, then we update only the lower bound since the transform maps
-                             * lower-left corner; the upper bound is unknown. If the gridToCRS maps pixel center, then we update
+                             * If the gridToCRS map cell corners, then we update only the lower bound since the transform maps
+                             * lower-left corner; the upper bound is unknown. If the gridToCRS maps cell center, then we update
                              * both lower and upper bounds to a value assumed to be in the center; the span is set to zero.
                              */
                             if (isCenter) {
diff --git a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
index a6adb7f..2447b22 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
@@ -55,6 +55,7 @@ import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.logging.Logging;
 import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.CharSequences;
+import org.apache.sis.util.ArraysExt;
 import org.apache.sis.util.Classes;
 import org.apache.sis.util.Debug;
 import org.apache.sis.io.TableAppender;
@@ -159,7 +160,7 @@ public class GridGeometry implements Serializable {
 
     /**
      * The conversion from grid indices to "real world" coordinates, or {@code null} if unknown.
-     * If non-null, the conversion shall map {@linkplain PixelInCell#CELL_CENTER pixel center}.
+     * If non-null, the conversion shall map {@linkplain PixelInCell#CELL_CENTER cell center}.
      * This conversion is usually, but not necessarily, affine.
      *
      * @see #CRS
@@ -169,7 +170,7 @@ public class GridGeometry implements Serializable {
     protected final MathTransform gridToCRS;
 
     /**
-     * Same conversion than {@link #gridToCRS} but from {@linkplain PixelInCell#CELL_CORNER pixel corner}
+     * Same conversion than {@link #gridToCRS} but from {@linkplain PixelInCell#CELL_CORNER cell corner}
      * instead than center. This transform is preferable to {@code gridToCRS} for transforming envelopes.
      *
      * @serial This field is serialized because it may be a value specified explicitly at construction time,
@@ -214,15 +215,23 @@ public class GridGeometry implements Serializable {
 
     /**
      * Creates a new grid geometry derived from the given grid geometry with a new extent and a modified transform.
-     * This constructor can be used for creating a grid geometry over a subregion, for example with the grid extent
-     * computed by {@link #getExtent(Envelope)}.
+     * This constructor is used for creating a grid geometry over a subregion (for example with the grid extent
+     * computed by {@link #getExtent(Envelope)}) or grid geometry for a sub-sampled raster.
      *
      * <p>If {@code toOther} is non-null, it should be a transform from the given {@code extent} coordinates to the
-     * {@code other} grid coordinates. The {@link #gridToCRS} transform of the new grid geometry will be set to the
-     * following concatenation:</p>
+     * {@code other} grid coordinates. That transform should be merely a {@linkplain MathTransforms#scale(double...)
+     * scale} and {@linkplain MathTransforms#translation(double...) translation} even if more complex transforms are
+     * accepted. The {@link #gridToCRS} transform of the new grid geometry will be set to the following concatenation:</p>
      *
      * <blockquote>{@code this.gridToCRS} = {@code toOther} → {@code other.gridToCRS}</blockquote>
      *
+     * The new {@linkplain #getEnvelope() grid geometry envelope} will be {@linkplain GeneralEnvelope#intersect(Envelope)
+     * clipped} to the envelope of the other grid geometry. This is for preventing the envelope to become larger under the
+     * effect of sub-sampling (because each cell become larger).
+     *
+     * <p>This constructor is for subclasses that wish to specialize their {@link #subgrid(Envelope, double...)} or related
+     * methods. Users should invoke {@code GridGeometry} member methods or use {@link GridChange} instead.</p>
+     *
      * @param  other    the other grid geometry to copy.
      * @param  extent   the new extent for the grid geometry to construct, or {@code null} if none.
      * @param  toOther  transform from this grid coordinates to {@code other} grid coordinates, or {@code null} if none.
@@ -230,8 +239,9 @@ public class GridGeometry implements Serializable {
      * @throws TransformException if the math transform can not compute the geospatial envelope from the grid extent.
      *
      * @see #getExtent(Envelope)
+     * @see #subgrid(Envelope, double...)
      */
-    public GridGeometry(final GridGeometry other, final GridExtent extent, final MathTransform toOther) throws TransformException {
+    protected GridGeometry(final GridGeometry other, final GridExtent extent, final MathTransform toOther) throws TransformException {
         ArgumentChecks.ensureNonNull("other", other);
         final int dimension = other.getDimension();
         this.extent = extent;
@@ -249,7 +259,12 @@ public class GridGeometry implements Serializable {
             resolution  = resolution(gridToCRS, extent);
             nonLinears  = findNonLinearTargets(gridToCRS);
         }
-        envelope = computeEnvelope(gridToCRS, other.envelope != null ? other.envelope.getCoordinateReferenceSystem() : null);
+        ImmutableEnvelope envelope = other.envelope;            // We will share the same instance if possible.
+        ImmutableEnvelope computed = computeEnvelope(gridToCRS, (envelope != null) ? envelope.getCoordinateReferenceSystem() : null, envelope);
+        if (computed == null || !computed.equals(envelope)) {
+            envelope = computed;
+        }
+        this.envelope = envelope;
         if (envelope == null && gridToCRS == null) {
             ArgumentChecks.ensureNonNull("extent", extent);
         }
@@ -263,7 +278,7 @@ public class GridGeometry implements Serializable {
      */
 
     /**
-     * Creates a new grid geometry from a grid extent and a mapping from pixel coordinates to "real world" coordinates.
+     * Creates a new grid geometry from a grid extent and a mapping from cell coordinates to "real world" coordinates.
      * At least one of {@code extent}, {@code gridToCRS} or {@code crs} arguments shall be non-null.
      * If {@code gridToCRS} is non-null, then {@code anchor} shall be non-null too with one of the following values:
      *
@@ -312,8 +327,8 @@ public class GridGeometry implements Serializable {
         this.extent      = extent;
         this.gridToCRS   = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CENTER);
         this.cornerToCRS = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CORNER);
-        this.envelope    = computeEnvelope(gridToCRS, crs);     // 'gridToCRS' specified by the user, not 'this.gridToCRS'.
-        this.resolution  = resolution(gridToCRS, extent);       // 'gridToCRS' or 'cornerToCRS' does not matter here.
+        this.envelope    = computeEnvelope(gridToCRS, crs, null);   // 'gridToCRS' specified by the user, not 'this.gridToCRS'.
+        this.resolution  = resolution(gridToCRS, extent);           // 'gridToCRS' or 'cornerToCRS' does not matter here.
         this.nonLinears  = findNonLinearTargets(gridToCRS);
     }
 
@@ -323,28 +338,36 @@ public class GridGeometry implements Serializable {
      *
      * @param  specified  the transform specified by the user. This is not necessarily {@link #gridToCRS}.
      * @param  crs        the coordinate reference system to declare in the envelope.
+     * @param  limits     if non-null, intersect with that envelope. The CRS must be the same than {@code crs}.
      */
-    private ImmutableEnvelope computeEnvelope(final MathTransform specified, final CoordinateReferenceSystem crs) throws TransformException {
-        GeneralEnvelope env = null;
+    private ImmutableEnvelope computeEnvelope(final MathTransform specified, final CoordinateReferenceSystem crs,
+            final Envelope limits) throws TransformException
+    {
+        final GeneralEnvelope env;
         if (extent != null && cornerToCRS != null) {
             env = extent.toCRS(cornerToCRS, specified);
             env.setCoordinateReferenceSystem(crs);
+            if (limits != null) {
+                env.intersect(limits);
+            }
         } else if (crs != null) {
             env = new GeneralEnvelope(crs);
             env.setToNaN();
+        } else {
+            return null;
         }
-        return (env != null) ? new ImmutableEnvelope(env) : null;
+        return new ImmutableEnvelope(env);
     }
 
     /**
-     * Creates a new grid geometry from a geospatial envelope and a mapping from pixel coordinates to "real world" coordinates.
+     * Creates a new grid geometry from a geospatial envelope and a mapping from cell coordinates to "real world" coordinates.
      * At least one of {@code gridToCRS} or {@code envelope} arguments shall be non-null.
      * If {@code gridToCRS} is non-null, then {@code anchor} shall be non-null too with one of the values documented in the
      * {@link #GridGeometry(GridExtent, PixelInCell, MathTransform, CoordinateReferenceSystem) constructor expecting a grid
      * extent}.
      *
-     * <p>The given envelope shall encompass all cell surfaces, from the left border of leftmost pixel to the right border
-     * of the rightmost pixel and similarly along other axes. This constructor tries to store a geospatial envelope close
+     * <p>The given envelope shall encompass all cell surfaces, from the left border of leftmost cell to the right border
+     * of the rightmost cell and similarly along other axes. This constructor tries to store a geospatial envelope close
      * to the specified envelope, but there is no guarantees that the envelope returned by {@link #getEnvelope()} will be
      * equal to the given envelope. The envelope stored in the new {@code GridGeometry} may be slightly smaller, larger or
      * shifted because the floating point values used in geospatial envelope can not always be mapped to the integer
@@ -518,8 +541,10 @@ public class GridGeometry implements Serializable {
      * Returns the bounding box of "real world" coordinates for this grid geometry.
      * This envelope is computed from the {@linkplain #getExtent() grid extent}, which is
      * {@linkplain #getGridToCRS(PixelInCell) transformed} to the "real world" coordinate system.
-     * The returned envelope encompasses all cell surfaces, from the left border of leftmost pixel
-     * to the right border of the rightmost pixel and similarly along other axes.
+     * The initial envelope encompasses all cell surfaces, from the left border of leftmost cell
+     * to the right border of the rightmost cell and similarly along other axes.
+     * If this grid geometry is a {@linkplain #subgrid(Envelope, double...) subgrid}, then the envelope is also
+     * {@linkplain GeneralEnvelope#intersect(Envelope) clipped} to the envelope of the original (non sub-sampled) grid geometry.
      *
      * @return the bounding box in "real world" coordinates (never {@code null}).
      * @throws IncompleteGridGeometryException if this grid geometry has no envelope —
@@ -568,6 +593,7 @@ public class GridGeometry implements Serializable {
      * @throws IncompleteGridGeometryException if this grid geometry has no extent or no "grid to CRS" transform.
      * @throws TransformException if an error occurred while converting the envelope coordinates to grid coordinates.
      *
+     * @see #subgrid(Envelope, double...)
      * @see #GridGeometry(GridGeometry, GridExtent, MathTransform)
      */
     public GridExtent getExtent(final Envelope areaOfInterest) throws IncompleteGridGeometryException, TransformException {
@@ -613,16 +639,10 @@ public class GridGeometry implements Serializable {
         return sub.equals(extent) ? extent : sub;
     }
 
-    /*
-     * Do not provide a convenience 'getGridToCRS()' method without PixelInCell or PixelOrientation argument.
-     * Experience shows that 0.5 pixel offset in image localization is a recurrent problem. We really want to
-     * force developers to think about whether they want a gridToCRS transform locating pixel corner or center.
-     */
-
     /**
      * Returns the conversion from grid coordinates to "real world" coordinates.
      * The conversion is often an affine transform, but not necessarily.
-     * Conversions from pixel indices to geospatial coordinates can be performed for example as below:
+     * Conversions from cell indices to geospatial coordinates can be performed for example as below:
      *
      * {@preformat java
      *     MathTransform  gridToCRS     = gridGeometry.getGridToCRS(PixelInCell.CELL_CENTER);
@@ -630,7 +650,7 @@ public class GridGeometry implements Serializable {
      *     DirectPosition aPixelCenter  = gridToCRS.transform(indicesOfCell, null);
      * }
      *
-     * Callers must specify whether they want the "real world" coordinates of pixel center or pixel corner.
+     * Callers must specify whether they want the "real world" coordinates of cell center or cell corner.
      * The cell corner is the one for which all grid indices have the smallest values (closest to negative infinity).
      * As a rule of thumb:
      *
@@ -666,6 +686,12 @@ public class GridGeometry implements Serializable {
         throw incomplete(GRID_TO_CRS, Resources.Keys.UnspecifiedTransform);
     }
 
+    /*
+     * Do not provide a convenience 'getGridToCRS()' method without PixelInCell or PixelOrientation argument.
+     * Experience shows that 0.5 pixel offset in image localization is a recurrent problem. We really want to
+     * force developers to think about whether they want a gridToCRS transform locating pixel corner or center.
+     */
+
     /**
      * Returns an <em>estimation</em> of the grid resolution, in units of the coordinate reference system axes.
      * The length of the returned array is the number of CRS dimensions, with {@code resolution[0]}
@@ -912,6 +938,74 @@ public class GridGeometry implements Serializable {
     }
 
     /**
+     * Returns a grid geometry over a sub-region of this grid geometry and optionally with sub-sampling.
+     * The given envelope can be expressed in any coordinate reference system (CRS) accepted by {@link #getExtent(Envelope)}.
+     * The target resolution, if provided, shall be in the units of CRS axes, in same order.
+     *
+     * @param  areaOfInterest    the desired spatiotemporal region in any CRS (transformations will be applied as needed),
+     *                           or {@code null} for not restricting the sub-grid to a sub-area.
+     * @param  targetResolution  the desired resolution in units of the CRS axes, or {@code null} or an empty array
+     *                           if no sub-sampling is desired.
+     * @return a grid geometry over the specified sub-region of this grid geometry.
+     * @throws IncompleteGridGeometryException if this grid geometry has no extent or no "grid to CRS" transform.
+     * @throws TransformException if an error occurred while converting the envelope coordinates to grid coordinates.
+     *
+     * @see #getExtent(Envelope)
+     */
+    public GridGeometry subgrid(final Envelope areaOfInterest, double... targetResolution) throws TransformException {
+        GridExtent domain = extent;
+        if (areaOfInterest != null) {
+            domain = getExtent(areaOfInterest);
+        }
+        MathTransform mt = null;
+        if (targetResolution != null && targetResolution.length != 0) {
+            /*
+             * Before to compute the strides, make sure we have all required information.
+             * One intent of following calls to getExtent() and getGridToCRS(…) is to get
+             * IncompleteGridGeometryException thrown if an information is missing.
+             */
+            if (domain == null) domain = getExtent();
+            final MathTransform gridToCRS = getGridToCRS(PixelInCell.CELL_CENTER);
+            /*
+             * Convert the target resolutions to sub-samplings as number of grid cells.
+             * The sub-sampling will be scale factors for the new "grid to CRS" transforms.
+             */
+            final int dimension = getDimension();
+            targetResolution = ArraysExt.resize(targetResolution, dimension);
+            Matrix m = gridToCRS.derivative(new DirectPositionView.Double(domain.getPointOfInterest()));
+            final double[] scales = Matrices.inverse(m).multiply(targetResolution);
+            /*
+             * Creates the matrix to pre-concatenate to "grid to CRS" transform. The (low /si) term in translation
+             * below must be computed in the same way than the "low" coordinates in GridExtent(GridExtent, int...)
+             * constructor.
+             */
+            m = Matrices.createIdentity(dimension + 1);
+            final int[] strides = new int[dimension];
+            Arrays.fill(strides, 1);
+            boolean isIdentity = true;
+            for (int i=0; i<dimension; i++) {
+                final double s = Math.floor(Math.abs(scales[i]));
+                if (s > 1) {                                                // Also for skipping NaN values.
+                    final int  si  = (int) s;
+                    final long low = domain.getLow(i);
+                    m.setElement(i, i, s);
+                    m.setElement(i, dimension, low - (low / si) * s);       // (low / si) must be consistent with GridExtent.
+                    strides[i] = si;
+                    isIdentity = false;
+                }
+            }
+            if (!isIdentity) {
+                mt = MathTransforms.linear(m);
+                domain = new GridExtent(domain, strides);
+            }
+        }
+        if (mt == null && Objects.equals(domain, extent)) {
+            return this;
+        }
+        return new GridGeometry(this, domain, mt);
+    }
+
+    /**
      * Returns a hash value for this grid geometry. This value need not remain
      * consistent between different implementations of the same class.
      */
diff --git a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/PixelTranslation.java b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/PixelTranslation.java
index 4dfd58e..a6f69f5 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/PixelTranslation.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/PixelTranslation.java
@@ -43,7 +43,7 @@ import org.apache.sis.referencing.operation.transform.MathTransforms;
  * </ul>
  *
  * This class provides also a few {@code translate(…)} convenience methods,
- * which apply the pixel translation on a given {@link MathTransform} instance.
+ * which apply the translation on a given {@link MathTransform} instance.
  *
  * <div class="note"><b>Example:</b>
  * if the following code snippet, {@code gridToCRS} is an {@link java.awt.geom.AffineTransform} from
@@ -165,7 +165,7 @@ public final class PixelTranslation extends Static implements Serializable {
     }
 
     /**
-     * Returns the position relative to the pixel center.
+     * Returns the position relative to the cell center.
      * This method is typically used for <var>n</var>-dimensional grids, where the number of dimension is unknown.
      * The translation is determined from the following table, with the same value applied to all dimensions:
      *
diff --git a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridChangeTest.java b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridChangeTest.java
new file mode 100644
index 0000000..00cc475
--- /dev/null
+++ b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridChangeTest.java
@@ -0,0 +1,104 @@
+/*
+ * 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.coverage.grid;
+
+import org.opengis.referencing.datum.PixelInCell;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.matrix.Matrix3;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.test.DependsOn;
+import org.apache.sis.test.TestCase;
+import org.junit.Test;
+
+import static org.apache.sis.test.ReferencingAssert.*;
+import static org.apache.sis.coverage.grid.GridExtentTest.assertExtentEquals;
+
+
+/**
+ * Tests {@link GridChange}.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since   1.0
+ * @module
+ */
+@DependsOn({GridExtentTest.class, GridGeometryTest.class})
+public final strictfp class GridChangeTest extends TestCase {
+    /**
+     * Creates a grid geometry with the given extent and scale for testing purpose.
+     * An arbitrary translation of (2,3) is added to the "grid to CRS" conversion.
+     */
+    private static GridGeometry grid(int xmin, int ymin, int xmax, int ymax, int xScale, int yScale) throws TransformException {
+        GridExtent extent = new GridExtent(null, new long[] {xmin, ymin}, new long[] {xmax, ymax}, true);
+        Matrix3 gridToCRS = new Matrix3();
+        gridToCRS.m00 = xScale;
+        gridToCRS.m11 = yScale;
+        gridToCRS.m02 = 200;            // Arbitrary translation.
+        gridToCRS.m12 = 500;
+        return new GridGeometry(extent, PixelInCell.CELL_CORNER, MathTransforms.linear(gridToCRS), null);
+    }
+
+    /**
+     * Tests the construction from grid geometries having a linear "grid to CRS" conversion.
+     *
+     * @throws TransformException if an error occurred while computing the grid geometry.
+     */
+    @Test
+    public void testLinear2D() throws TransformException {
+        GridGeometry source = grid(  10,   -20,  110,  180, 100, -300);     // Envelope x: [1200 … 11300]   y: [-53800 … 6500]
+        GridGeometry target = grid(2000, -1000, 9000, 8000,   2,   -1);     // Envelope x: [4200 … 18202]   y: [ -7501 … 1500]
+        GridChange   change = new GridChange(source, target);               // Envelope x: [4200 … 11300]   y: [ -7501 … 1500]
+        GridExtent   extent = change.getTargetExtent();
+        assertExtentEquals(extent, 0,  2000, 5549);                         // Subrange of target extent.
+        assertExtentEquals(extent, 1, -1000, 8000);
+        assertArrayEquals("strides", new int[] {50, 300}, change.getTargetStrides());       // s = scaleSource / scaleTarget
+        /*
+         * Scale factors in following matrix shall be the same than above strides.
+         * Translation appears only with PixelInCell different than the one used at construction.
+         */
+        Matrix3 c = new Matrix3();
+        c.m00 =  50;
+        c.m11 = 300;
+        assertMatrixEquals("CELL_CORNER", c, MathTransforms.getMatrix(change.getConversion(PixelInCell.CELL_CORNER)), STRICT);
+        c.m02 =  24.5;
+        c.m12 = 149.5;
+        assertMatrixEquals("CELL_CENTER", c, MathTransforms.getMatrix(change.getConversion(PixelInCell.CELL_CENTER)), STRICT);
+        /*
+         * If we do not ask for strides, the 'gridToCRS' transforms shall be the same than the 'target' geometry.
+         * The envelope is the intersection of the envelopes of 'source' and 'target' geometries, documented above.
+         */
+        GridGeometry tg = change.getTargetGeometry();
+        assertSame("extent",      extent, tg.getExtent());
+        assertSame("CELL_CORNER", target.getGridToCRS(PixelInCell.CELL_CORNER), tg.getGridToCRS(PixelInCell.CELL_CORNER));
+        assertSame("CELL_CENTER", target.getGridToCRS(PixelInCell.CELL_CENTER), tg.getGridToCRS(PixelInCell.CELL_CENTER));
+        GeneralEnvelope expected = new GeneralEnvelope(2);
+        expected.setRange(0,  4200, 11300);
+        expected.setRange(1, -7501,  1500);
+        assertEnvelopeEquals(expected, tg.getEnvelope(), STRICT);
+        /*
+         * If we ask for strides, then the envelope should be approximately the same or smaller. Note that without
+         * the clipping documented in GridExtent(GridExtent, int...) constructor, the envelope could be larger.
+         */
+        tg = change.getTargetGeometry(50, 300);
+        assertEnvelopeEquals(expected, tg.getEnvelope(), STRICT);
+        assertMatrixEquals("gridToCRS", new Matrix3(
+                100,    0, 200,
+                  0, -300, 600,
+                  0,    0,   1), MathTransforms.getMatrix(tg.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
+    }
+}
diff --git a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java
index 3519a8d..a6a4d29 100644
--- a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java
+++ b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java
@@ -38,14 +38,35 @@ import static org.apache.sis.test.Assert.*;
  */
 public final strictfp class GridExtentTest extends TestCase {
     /**
+     * Creates a three-dimensional grid extent to be shared by different tests.
+     */
+    private static GridExtent create3D() {
+        return new GridExtent(
+                new DimensionNameType[] {DimensionNameType.COLUMN, DimensionNameType.ROW, DimensionNameType.TIME},
+                new long[] {100, 200, 40}, new long[] {500, 800, 50}, false);
+    }
+
+    /**
      * Verifies the low and high values in the specified dimension of the given extent
      */
-    private static void assertExtentEquals(final GridExtent extent, final int dimension, final int low, final int high) {
+    static void assertExtentEquals(final GridExtent extent, final int dimension, final int low, final int high) {
         assertEquals("low",  low,  extent.getLow (dimension));
         assertEquals("high", high, extent.getHigh(dimension));
     }
 
     /**
+     * Tests the {@link GridExtent#GridExtent(GridExtent, int[])} constructor.
+     */
+    @Test
+    public void testCreateFromStrides() {
+        GridExtent extent = create3D();
+        extent = new GridExtent(extent, new int[] {4, 3, 9});
+        assertExtentEquals(extent, 0, 25, 124);                 // 100 cells
+        assertExtentEquals(extent, 1, 66, 265);                 // 200 cells
+        assertExtentEquals(extent, 2,  4,   5);                 //   2 cells
+    }
+
+    /**
      * Tests the {@link GridExtent#GridExtent(AbstractEnvelope, GridRoundingMode, int[], GridExtent, int[])} constructor.
      */
     @Test
@@ -98,15 +119,6 @@ public final strictfp class GridExtentTest extends TestCase {
     }
 
     /**
-     * Creates a three-dimensional grid extent to be shared by different tests.
-     */
-    private static GridExtent create3D() {
-        return new GridExtent(
-                new DimensionNameType[] {DimensionNameType.COLUMN, DimensionNameType.ROW, DimensionNameType.TIME},
-                new long[] {100, 200, 40}, new long[] {500, 800, 50}, false);
-    }
-
-    /**
      * Tests {@link GridExtent#subExtent(int, int)}.
      */
     @Test
diff --git a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
index c6ebd8f..bb651e7 100644
--- a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
+++ b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
@@ -143,6 +143,31 @@ public final strictfp class GridGeometryTest extends TestCase {
     }
 
     /**
+     * Tests the {@link GridGeometry#GridGeometry(GridGeometry, GridExtent, MathTransform)} constructor.
+     *
+     * @throws TransformException if an error occurred while using the "grid to CRS" transform.
+     */
+    @Test
+    public void testFromOther() throws TransformException {
+        long[]        low       = new long[] {  1,   3, 2};
+        long[]        high      = new long[] {101, 203, 4};
+        GridExtent    extent    = new GridExtent(null, low, high, false);
+        MathTransform gridToCRS = MathTransforms.translation(5, 7, 8);
+        GridGeometry  grid      = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, null);
+
+        low    = new long[] { 11,  35, 20};
+        high   = new long[] {120, 250, 39};
+        extent = new GridExtent(null, low, high, false);
+        grid   = new GridGeometry(grid, extent, MathTransforms.scale(2, 1, 3));
+        assertSame(extent, grid.getExtent());
+        assertMatrixEquals("gridToCRS", new Matrix4(
+                2, 0, 0, 5,
+                0, 1, 0, 7,     // Combination of above scales (diagonal) and translation (last column).
+                0, 0, 3, 8,
+                0, 0, 0, 1), MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CENTER)), STRICT);
+    }
+
+    /**
      * Tests construction from a <cite>grid to CRS</cite> having a 0.5 pixel translation.
      * This translation happens in transform mapping <cite>pixel center</cite> when the
      * corresponding <cite>pixel corner</cite> transformation is identity.
@@ -290,4 +315,37 @@ public final strictfp class GridGeometryTest extends TestCase {
         assertExtentEquals(new long[] { 56, 69, 2},
                            new long[] {130, 73, 4}, actual);
     }
+
+    /**
+     * Tests {@link GridGeometry#subgrid(Envelope, double...)}.
+     *
+     * @throws TransformException if an error occurred during computation.
+     */
+    @Test
+    @DependsOnMethod({"testFromGeospatialEnvelope", "testGetExtent"})
+    public void testSubgrid() throws TransformException {
+        final GeneralEnvelope envelope = new GeneralEnvelope(HardCodedCRS.WGS84_φλ);
+        envelope.setRange(0, -70, +80);
+        envelope.setRange(1,   5,  15);
+        final MathTransform gridToCRS = MathTransforms.linear(new Matrix3(
+            0,   0.5, -90,
+            0.5, 0,  -180,
+            0,   0,     1));
+        GridGeometry grid = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, envelope, GridRoundingMode.NEAREST);
+        assertExtentEquals(new long[] {370, 40}, new long[] {389, 339}, grid.getExtent());
+        assertEnvelopeEquals(envelope, grid.getEnvelope(), STRICT);
+        /*
+         * Set a sub-region. The grid extent and "grid to CRS" transform shall be adjusted
+         * in such a way that envelope computed from the new grid geometry is the same.
+         */
+        envelope.setRange(0, -50, +30);
+        envelope.setRange(1,   8,  12);
+        grid = grid.subgrid(envelope, 1, 2);
+        assertExtentEquals(new long[] {94, 40}, new long[] {95, 119}, grid.getExtent());
+        assertEnvelopeEquals(envelope, grid.getEnvelope(), STRICT);
+        assertMatrixEquals("gridToCRS", new Matrix3(
+                  0, 1,  -90,
+                  2, 0, -180,
+                  0, 0,    1), MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
+    }
 }
diff --git a/core/sis-raster/src/test/java/org/apache/sis/test/suite/RasterTestSuite.java b/core/sis-raster/src/test/java/org/apache/sis/test/suite/RasterTestSuite.java
index 2da0f9f..1aac8e2 100644
--- a/core/sis-raster/src/test/java/org/apache/sis/test/suite/RasterTestSuite.java
+++ b/core/sis-raster/src/test/java/org/apache/sis/test/suite/RasterTestSuite.java
@@ -34,6 +34,7 @@ import org.junit.BeforeClass;
     org.apache.sis.coverage.grid.PixelTranslationTest.class,
     org.apache.sis.coverage.grid.GridExtentTest.class,
     org.apache.sis.coverage.grid.GridGeometryTest.class,
+    org.apache.sis.coverage.grid.GridChangeTest.class,
     org.apache.sis.coverage.CategoryTest.class,
     org.apache.sis.coverage.CategoryListTest.class,
     org.apache.sis.coverage.SampleDimensionTest.class,


Mime
View raw message