You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2023/04/14 18:28:09 UTC

[sis] 02/02: Relax the restriction that all inputs to `BandAggregateGridCoverage` have the same grid geometry. With this commit, they are allowed to have different translations but not yet more complex changes.

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 9f8a4020152877c9a3dc82a162080a6c4d28f894
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Fri Apr 14 20:27:09 2023 +0200

    Relax the restriction that all inputs to `BandAggregateGridCoverage` have the same grid geometry.
    With this commit, they are allowed to have different translations but not yet more complex changes.
---
 .../coverage/grid/BandAggregateGridCoverage.java   |   9 +-
 .../coverage/grid/CoordinateOperationFinder.java   |   5 +-
 .../apache/sis/coverage/grid/DefaultEvaluator.java |   2 +-
 .../sis/coverage/grid/GridCoverageProcessor.java   |   2 +
 .../org/apache/sis/coverage/grid/GridExtent.java   |  22 +-
 .../org/apache/sis/image/MultiSourceLayout.java    |  17 +-
 .../sis/internal/coverage/CommonDomainFinder.java  | 372 +++++++++++++++++++++
 .../sis/internal/coverage/MultiSourceArgument.java |  39 ++-
 .../org/apache/sis/internal/feature/Resources.java |   5 +
 .../sis/internal/feature/Resources.properties      |   1 +
 .../sis/internal/feature/Resources_fr.properties   |   1 +
 .../grid/BandAggregateGridCoverageTest.java        | 137 ++++++++
 .../apache/sis/test/suite/FeatureTestSuite.java    |   1 +
 .../apache/sis/storage/aggregate/GridSlice.java    |  64 +---
 .../sis/storage/aggregate/GridSliceLocator.java    |   3 +-
 .../sis/storage/aggregate/GroupByTransform.java    |   4 +-
 16 files changed, 592 insertions(+), 92 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java
index 9e671073a9..11a98fe7a4 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java
@@ -64,6 +64,12 @@ final class BandAggregateGridCoverage extends GridCoverage {
      */
     private final int numBands;
 
+    /**
+     * Translations in units of grid cells to apply on an extent of this grid coverage
+     * for getting a "grid to CRS" transform compatible with each source.
+     */
+    private final long[][] gridTranslations;
+
     /**
      * Index of a sources having the same "grid to CRS" than this grid coverage, or -1 if none.
      */
@@ -94,6 +100,7 @@ final class BandAggregateGridCoverage extends GridCoverage {
         this.sources           = aggregate.sources();
         this.bandsPerSource    = aggregate.bandsPerSource(true);
         this.numBands          = aggregate.numBands();
+        this.gridTranslations  = aggregate.gridTranslations();
         this.sourceOfGridToCRS = aggregate.sourceOfGridToCRS();
         this.processor         = processor;
         this.dataType          = sources[0].getBandType();
@@ -138,7 +145,7 @@ final class BandAggregateGridCoverage extends GridCoverage {
         }
         final RenderedImage[] images = new RenderedImage[sources.length];
         for (int i=0; i<images.length; i++) {
-            images[i] = sources[i].render(sliceExtent);
+            images[i] = sources[i].render(sliceExtent.translate(gridTranslations[i]));
         }
         return processor.aggregateBands(images, bandsPerSource);
     }
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 e5f653f6d5..ee1ab98896 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
@@ -233,10 +233,11 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
     }
 
     /**
-     * Verifies whether the presence of a CRS considered mandatory, unless the CRS of opposite grid
-     * is also missing.
+     * Verifies the presence of a CRS considered mandatory,
+     * unless the CRS of the opposite grid is also missing.
      *
      * @param  rs  {@code true} is source CRS is mandatory, {@code false} if target CRS is mandatory.
+     * @throws IncompleteGridGeometryException if a mandatory CRS is missing.
      */
     final void verifyPresenceOfCRS(final boolean rs) {
         if ((rs ? target : source).isDefined(GridGeometry.CRS)) {
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java
index 66a7572630..65467914e6 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java
@@ -639,7 +639,7 @@ class DefaultEvaluator implements GridCoverage.Evaluator {
             }
         }
         // Modify fields only after everything else succeeded.
-        position     = new FractionalGridCoordinates.Position(crsToGrid.getTargetDimensions());
+        position    = new FractionalGridCoordinates.Position(crsToGrid.getTargetDimensions());
         inputCRS    = crs;
         inputToGrid = crsToGrid;
     }
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java
index 9394a7c603..d04f92c113 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java
@@ -773,6 +773,8 @@ public class GridCoverageProcessor implements Cloneable {
      *   <li>All coverages shall use the same data type in their rendered image.</li>
      * </ul>
      *
+     * Some of those restrictions may be relaxed in future versions.
+     *
      * @param  sources  coverages whose bands shall be aggregated, in order. At least one coverage must be provided.
      * @param  bandsPerSource  bands to use for each source coverage, in order. May contain {@code null} elements.
      * @return the aggregated coverage, or one of the sources if it can be used directly.
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
index de45ba8a24..99da3df96e 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
@@ -1724,18 +1724,18 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable
         ArgumentChecks.ensureNonNull("translation", translation);
         final int m = getDimension();
         final int length = Math.min(m, translation.length);
-        if (!isZero(translation, length)) {
-            final GridExtent translated = new GridExtent(this);
-            final long[] c = translated.coordinates;
-            for (int i=0; i < length; i++) {
-                final int  j = i + m;
-                final long t = translation[i];
-                c[i] = Math.addExact(c[i], t);
-                c[j] = Math.addExact(c[j], t);
-            }
-            return translated;
+        if (isZero(translation, length)) {
+            return this;
+        }
+        final GridExtent translated = new GridExtent(this);
+        final long[] c = translated.coordinates;
+        for (int i=0; i < length; i++) {
+            final int  j = i + m;
+            final long t = translation[i];
+            c[i] = Math.addExact(c[i], t);
+            c[j] = Math.addExact(c[j], t);
         }
-        return this;
+        return translated;
     }
 
     /**
diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java b/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java
index a49f466d96..d4eae377bc 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java
@@ -34,6 +34,7 @@ import org.apache.sis.internal.coverage.j2d.ImageLayout;
 import org.apache.sis.internal.coverage.j2d.ImageUtilities;
 import org.apache.sis.internal.coverage.j2d.ColorModelFactory;
 import org.apache.sis.internal.coverage.MultiSourceArgument;
+import org.apache.sis.internal.coverage.CommonDomainFinder;
 import org.apache.sis.coverage.grid.DisjointExtentException;
 
 
@@ -174,12 +175,18 @@ final class MultiSourceLayout extends ImageLayout {
                 }
                 /*
                  * Get the domain of the combined image to create.
-                 * TODO: current implementation computes the intersection of all sources.
-                 * But a future version should allow users to specify if they want intersection,
-                 * union or strict mode instead. A "strict" mode would prevent the combination of
-                 * images using different domains (i.e. raise an error if domains are not the same).
+                 * Current implementation computes the intersection of all sources.
                  */
-                ImageUtilities.clipBounds(source, domain);
+                if (CommonDomainFinder.INTERSECTION) {
+                    ImageUtilities.clipBounds(source, domain);
+                } else if (!domain.equals(ImageUtilities.getBounds(source))) {
+                    /*
+                     * TODO: a future version should allow users to specify if they want intersection,
+                     * union or strict mode instead. A "strict" mode would prevent the combination of
+                     * images using different domains (i.e. raise an error if domains are not the same).
+                     */
+                    throw new IllegalArgumentException();
+                }
                 if (domain.isEmpty()) {
                     throw new DisjointExtentException(Resources.format(Resources.Keys.SourceImagesDoNotIntersect));
                 }
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/CommonDomainFinder.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/CommonDomainFinder.java
new file mode 100644
index 0000000000..f6d28d8698
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/CommonDomainFinder.java
@@ -0,0 +1,372 @@
+/*
+ * 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.coverage;
+
+import java.util.Map;
+import java.util.LinkedHashMap;
+import java.util.NoSuchElementException;
+import org.opengis.util.FactoryException;
+import org.opengis.geometry.MismatchedDimensionException;
+import org.opengis.referencing.datum.PixelInCell;
+import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.operation.NoninvertibleTransformException;
+import org.apache.sis.referencing.CRS;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.IllegalGridGeometryException;
+import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
+import org.apache.sis.internal.feature.Resources;
+import org.apache.sis.internal.util.Numerics;
+import org.apache.sis.util.Numbers;
+
+
+/**
+ * Helper class for building a combined domain from a list of grid geometries.
+ * After construction, one of the following methods shall be invoked exactly once.
+ *
+ * <ul>
+ *   <li>{@link #setFromGridAligned(GridGeometry...)}</li>
+ * </ul>
+ *
+ * Then, the result can be obtained by the given methods:
+ *
+ * <ul>
+ *   <li>{@link #result()}</li>
+ *   <li>{@link #gridTranslations()}</li>
+ *   <li>{@link #sourceOfGridToCRS()}</li>
+ * </ul>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.4
+ * @since   1.4
+ */
+public final class CommonDomainFinder {
+    /**
+     * Whether to compute intersection (true) or union (false) of all grid coverages.
+     * This is not yet configurable in current version.
+     *
+     * <p>We use this constant for tracking the code to update when we will want to provide an option for using
+     * union or strict equality instead (the latter would be a mode that fails if all extents are not identical).
+     * Note that in the case of unions, it would be possible to specify coverages with no intersection.
+     * Whether we should accept that or raise an exception is still an open question.</p>
+     */
+    public static final boolean INTERSECTION = true;
+
+    /**
+     * The grid geometry taken as the reference for computing the translations of all grid geometry items.
+     * Values of {@link #gridTranslations} and {@link #itemTranslations} are relative to that reference.
+     * At first this is an arbitrary grid geometry, for example the first encountered one.
+     * After {@code CommonDomainFinder} completed its task, this is updated to the final result.
+     *
+     * @see #result()
+     */
+    private GridGeometry reference;
+
+    /**
+     * Coordinate reference system of the reference, or {@code null} if not yet known.
+     * It will be set to the CRS of the first grid geometry where the CRS is defined.
+     */
+    private CoordinateReferenceSystem crs;
+
+    /**
+     * The inverse of the "grid to CRS" transform of the grid geometry taken as a reference.
+     */
+    private MathTransform crsToGrid;
+
+    /**
+     * The convention to use for fetching the "grid to CRS" transforms.
+     */
+    private final PixelInCell anchor;
+
+    /**
+     * The combined extent, as the union or intersection of all grid extents.
+     */
+    private GridExtent extent;
+
+    /**
+     * The translation in units of grid cells from the {@linkplain #reference} grid geometry
+     * to the grid geometry in the key.
+     */
+    private final Map<GridGeometry,long[]> gridTranslations;
+
+    /**
+     * Translations in units of grid cells from the {@linkplain #reference} grid geometry to each item.
+     * For each index <var>i</var>, {@code itemTranslations[i]} is a value from {@link #gridTranslations}
+     * map and may be reused at more than one index <var>i</var>.
+     *
+     * @see #gridTranslations()
+     */
+    private long[][] itemTranslations;
+
+    /**
+     * If one of the grid geometries has the same "grid to CRS" than the common grid geometry, the index.
+     * Otherwise -1.
+     */
+    private int sourceOfGridToCRS;
+
+    /**
+     * Creates a new common domain finder.
+     *
+     * @param  anchor  the convention to use for fetching the "grid to CRS" transforms.
+     */
+    CommonDomainFinder(final PixelInCell anchor) {
+        this.anchor = anchor;
+        gridTranslations = new LinkedHashMap<>();
+    }
+
+    /**
+     * Computes a common grid geometry from the given items.
+     * All items shall share be aligned on the same grid.
+     * Items may be translated relative to each other,
+     * but the translations shall be an integer number of grid cells.
+     *
+     * <h4>Coordinate reference system</h4>
+     * If the CRS of a grid geometry is undefined, it is assumed the same than other grid geometries.
+     *
+     * @param  items  the grid geometries for which to compute a common grid geometry.
+     * @throws IllegalGridGeometryException if the specified item is not compatible with the reference grid geometry.
+     */
+    final void setFromGridAligned(final GridGeometry... items) {
+        itemTranslations = new long[items.length][];
+        for (int i=0; i<items.length; i++) {
+            itemTranslations[i] = gridTranslations.computeIfAbsent(items[i], this::itemToCommon);
+        }
+        /*
+         * Change the reference grid geometry for matching more closely the desired grid extent.
+         * If one item has exactly the desired grid extent, use it. Otherwise search for an item
+         * having the same origin. This criterion is arbitrary and may change in future version.
+         */
+        GridGeometry fallback = null;
+        for (final Map.Entry<GridGeometry,long[]> entry : gridTranslations.entrySet()) {
+            final GridGeometry item   = entry.getKey();
+            final GridExtent actual   = item.getExtent();
+            final GridExtent expected = extent.translate(entry.getValue());
+            if (actual.equals(expected)) {
+                setGridToCRS(items, item);      // Must be before `reference` assignation.
+                reference = item;
+                return;
+            }
+            if (fallback == null && extent.getLow().equals(item.getExtent().getLow())) {
+                fallback = item;
+            }
+        }
+        if (fallback == null) {
+            fallback = reference;
+        }
+        setGridToCRS(items, fallback);
+        try {
+            reference = fallback.relocate(extent);
+        } catch (TransformException e) {
+            throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IncompatibleGridGeometries), e);
+        }
+    }
+
+    /**
+     * Given a grid geometry with the desired "grid to CRS", saves its index in {@link #sourceOfGridToCRS}.
+     * This method updates all previously computed translations for making them relative to the new reference.
+     *
+     * Note: updating values of the {@link #gridTranslations} map indirectly update all
+     * {@link #itemTranslations} array elements.
+     */
+    private void setGridToCRS(final GridGeometry[] items, final GridGeometry item) {
+        sourceOfGridToCRS = indexOf(items, item);
+        if (item == reference) {                    // Quick check for a common case.
+            return;
+        }
+        final long[] oldReference = itemTranslations[indexOf(items, reference)];
+        final long[] newReference = itemTranslations[sourceOfGridToCRS];
+        final long[] change = new long[newReference.length];
+        for (int i=0; i < newReference.length; i++) {
+            change[i] = Math.subtractExact(newReference[i], oldReference[i]);
+        }
+        for (final long[] offset : gridTranslations.values()) {
+            for (int i=0; i < offset.length; i++) {
+                offset[i] = Math.subtractExact(offset[i], change[i]);
+            }
+        }
+    }
+
+    /**
+     * Returns the index of the given grid geometry.
+     * This method is invoked when the instance should always exist in the array.
+     */
+    private static int indexOf(final GridGeometry[] items, final GridGeometry item) {
+        for (int i=0; i<items.length; i++) {
+            if (items[i] == item) {
+                return i;
+            }
+        }
+        throw new NoSuchElementException();         // Should never happen.
+    }
+
+    /**
+     * Returns the common grid geometry computed from all specified items.
+     *
+     * @return a grid geometry which is the union or intersection of all specified items.
+     */
+    final GridGeometry result() {
+        if (crs != null && !reference.isDefined(GridGeometry.CRS)) {
+            reference = new GridGeometry(extent, anchor, reference.getGridToCRS(anchor), crs);
+        }
+        return reference;
+    }
+
+    /**
+     * Returns the translations (in units of grid cells) from the common grid geometry to all items.
+     * The items are the arguments given to {@link #setFromGridAligned(GridGeometry...)}, in order.
+     * The common grid geometry is the value returned by {@link #result()}.
+     *
+     * <p>The returned array should not be modified because it is not cloned.</p>
+     *
+     * @return translation from the common grid geometry to all items. This array is not cloned.
+     */
+    @SuppressWarnings("ReturnOfCollectionOrArrayField")
+    final long[][] gridTranslations() {
+        return itemTranslations;
+    }
+
+    /**
+     * If one items has the same "grid to CRS" than the common grid geometry, returns its index.
+     *
+     * @return index of an items having the desired "grid to CRS", or -1 if none.
+     */
+    final int sourceOfGridToCRS() {
+        return sourceOfGridToCRS;
+    }
+
+    /**
+     * Computes the translation in units of grid cells from the common grid geometry to the given item.
+     * This method also opportunistically computes the union or intersection of all grid extents.
+     *
+     * @param  item  the grid geometry for which to get a translation from the common grid geometry.
+     * @return translation in unis of grid cells. Note that the caller may reuse this array for many grid geometries.
+     * @throws IllegalGridGeometryException if the specified item is not compatible with the reference grid geometry.
+     */
+    private long[] itemToCommon(final GridGeometry item) {
+        /*
+         * Compute the change ourselves instead of invoking `GridGeometry.createTransformTo(…)`
+         * because we do not want wraparound handling when we search for a simple translation.
+         */
+        MathTransform change = item.getGridToCRS(anchor);
+        try {
+            if (crsToGrid == null) {
+                crsToGrid = change.inverse();
+                reference = item;
+            }
+            if (item.isDefined(GridGeometry.CRS)) {
+                final CoordinateReferenceSystem src = item.getCoordinateReferenceSystem();
+                if (crs == null) {
+                    crs = src;
+                } else {
+                    /*
+                     * Ask for a change of CRS without specifying an area of interest (AOI) on the assumption
+                     * that if the transformation is only a translation, the AOI would not make a difference.
+                     * It save not only the AOI computation cost, but also make easier for `findOperation(…)`
+                     * to use its cache.
+                     */
+                    change = MathTransforms.concatenate(change, CRS.findOperation(src, crs, null).getMathTransform());
+                }
+            }
+            change = MathTransforms.concatenate(change, crsToGrid);
+        } catch (FactoryException | NoninvertibleTransformException | MismatchedDimensionException e) {
+            throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IncompatibleGridGeometries), e);
+        }
+        final long[] offset = integerTranslation(MathTransforms.getMatrix(change), null);
+        if (offset == null) {
+            throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IncompatibleGridGeometries));
+        }
+        /*
+         * The grid geometry has been accepted as valid. Now compute the combined extent,
+         * taking offset in account. At this point this is an offset TO the common grid geometry.
+         * It will be converted to an offset FROM the common grid geometry after the extent update.
+         */
+        final GridExtent e = item.getExtent().translate(offset);
+        if (extent == null) {
+            extent = e;
+        } else if (INTERSECTION) {
+            extent = extent.intersect(e);
+        } else if (!extent.equals(e)) {
+            throw new IllegalGridGeometryException();
+        }
+        for (int i=0; i<offset.length; i++) {
+            offset[i] = Math.negateExact(offset[i]);
+        }
+        return offset;
+    }
+
+    /**
+     * If the given matrix is the identity matrix except for translation terms, returns the translation.
+     * Translation terms must be integer values and will be stored in the given {@code offset} array.
+     * If the matrix is not an integer translation, this method return {@code null} without modifying
+     * the given {@code offset} array.
+     *
+     * @param  change  conversion between two grid geometries, or {@code null}.
+     * @param  offset  where to store translation terms if the change is an integer translation, or {@code null}.
+     * @return the translation terms, or {@code null} if the given matrix does not met the conditions.
+     *
+     * @see org.apache.sis.referencing.operation.matrix.Matrices#isTranslation(Matrix)
+     */
+    public static long[] integerTranslation(final Matrix change, long[] offset) {
+        if (change == null) {
+            return null;
+        }
+        final int numRows = change.getNumRow();
+        final int numCols = change.getNumCol();
+        for (int j=0; j<numRows; j++) {
+            for (int i=0; i<numCols; i++) {
+                double tolerance = Numerics.COMPARISON_THRESHOLD;
+                double e = change.getElement(j, i);
+                if (i == j) {
+                    e--;
+                } else if (i == numCols - 1) {
+                    final double a = Math.abs(e);
+                    if (a > 1) {
+                        tolerance = Math.min(tolerance*a, 0.125);
+                    }
+                    e -= Math.rint(e);
+                }
+                if (!(Math.abs(e) <= tolerance)) {      // Use `!` for catching NaN.
+                    return null;
+                }
+            }
+        }
+        /*
+         * Store the translation terms after we have determined that they are integers.
+         * It must be an "all or nothing" operation (unless an exception is thrown).
+         */
+        if (offset == null) {
+            offset = new long[numRows - 1];
+        }
+        final int i = numCols - 1;
+        if (change instanceof ExtendedPrecisionMatrix) {
+            final var epm = (ExtendedPrecisionMatrix) change;
+            for (int j=0; j<offset.length; j++) {
+                final Number e = epm.getElementOrNull(j, i);
+                offset[j] = (e != null) ? Numbers.round(e) : 0;
+            }
+        } else {
+            for (int j=0; j<offset.length; j++) {
+                offset[j] = Math.round(change.getElement(j, i));
+            }
+        }
+        return offset;
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java
index 520ab3479a..3f3c1fe100 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java
@@ -25,6 +25,7 @@ import java.util.Objects;
 import java.util.function.Consumer;
 import java.util.function.Function;
 import java.util.function.ToIntFunction;
+import org.opengis.referencing.datum.PixelInCell;
 import org.apache.sis.coverage.SampleDimension;
 import org.apache.sis.coverage.grid.GridGeometry;
 import org.apache.sis.coverage.grid.IllegalGridGeometryException;
@@ -32,11 +33,10 @@ import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.ArraysExt;
-import org.apache.sis.util.ComparisonMode;
 
 
 /**
- * Helper class for building a list of sample dimensions from aggregated sources.
+ * Helper class for building a combined domain or range from aggregated sources.
  * This helper class is shared for aggregation operations on different sources:
  * rendered images, grid coverages and resources.
  *
@@ -112,6 +112,12 @@ public final class MultiSourceArgument<S> {
      */
     private List<SampleDimension> ranges;
 
+    /**
+     * Translations in units of grid cells to apply for obtaining a grid geometry
+     * compatible with the "grid to CRS" transform of a source.
+     */
+    private long[][] gridTranslations;
+
     /**
      * Index of a source having the same "grid to CRS" transform than the grid geometry
      * returned by {@link #domain(Function)}. If there is none, then this value is -1.
@@ -510,19 +516,28 @@ next:   for (int i=0; i<sources.length; i++) {          // `sources.length` may
      * @param  getter  the method to invoke for getting grid geometry from a source.
      * @return intersection of all grid geometries.
      * @throws IllegalGridGeometryException if a grid geometry is not compatible with the others.
-     *
-     * @todo Current implementation requires that all grid geometry are equal. We need to relax that.
      */
     public GridGeometry domain(final Function<S, GridGeometry> getter) {
         checkValidationState(true);
-        GridGeometry intersection = getter.apply(sources[0]);
-        for (int i=1; i < validatedSourceCount; i++) {
-            if (!intersection.equals(getter.apply(sources[i]), ComparisonMode.IGNORE_METADATA)) {
-                throw new IllegalGridGeometryException("Not yet supported on coverages with different grid geometries.");
-            }
-        }
-        sourceOfGridToCRS = 0;      // TODO: to be computed when different grid geometries will be allowed. Prefer widest extent.
-        return intersection;
+        final var finder = new CommonDomainFinder(PixelInCell.CELL_CORNER);
+        finder.setFromGridAligned(Arrays.stream(sources).map(getter).toArray(GridGeometry[]::new));
+        sourceOfGridToCRS = finder.sourceOfGridToCRS();
+        gridTranslations  = finder.gridTranslations();
+        return finder.result();
+    }
+
+    /**
+     * Returns the translations in units of grid cells to apply for obtaining a grid geometry
+     * compatible with the "grid to CRS" transform of a source.
+     *
+     * <p>The returned array should not be modified because it is not cloned.</p>
+     *
+     * @return translations from the common grid geometry to all items. This array is not cloned.
+     */
+    @SuppressWarnings("ReturnOfCollectionOrArrayField")
+    public long[][] gridTranslations() {
+        checkValidationState(true);
+        return gridTranslations;
     }
 
     /**
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java
index caf2a3ec8e..1c89d38a7d 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java
@@ -256,6 +256,11 @@ public final class Resources extends IndexedResourceBundle {
          */
         public static final short IncompatibleColorModel = 82;
 
+        /**
+         * At least two coverages have mutually incompatible grid geometries.
+         */
+        public static final short IncompatibleGridGeometries = 87;
+
         /**
          * The ({0}, {1}) tile has an unexpected size, number of bands or sample layout.
          */
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties
index 0a1c6c2018..853ab78e0a 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties
@@ -58,6 +58,7 @@ ImageAllowsTransparency           = Image allows transparency.
 ImageHasAlphaChannel              = Image has alpha channel.
 ImageIsOpaque                     = Image is opaque.
 IncompatibleColorModel            = Color model is incompatible with sample model.
+IncompatibleGridGeometries        = At least two coverages have mutually incompatible grid geometries.
 IncompatibleTile_2                = The ({0}, {1}) tile has an unexpected size, number of bands or sample layout.
 InsufficientBufferCapacity_3      = Data buffer capacity is insufficient for a grid of {0} cells \u00d7 {1} bands. Missing {2} elements.
 IterationIsFinished               = Iteration is finished.
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties
index 356dc1cecf..0ae9c2ce3b 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties
@@ -63,6 +63,7 @@ ImageAllowsTransparency           = L\u2019image permet la transparence.
 ImageHasAlphaChannel              = L\u2019image a un canal alpha.
 ImageIsOpaque                     = L\u2019image est opaque.
 IncompatibleColorModel            = Le mod\u00e8le de couleurs est incompatible.
+IncompatibleGridGeometries        = Au moins deux couvertures de donn\u00e9es ont des g\u00e9om\u00e9tries de grilles mutuellement incompatibles.
 IncompatibleTile_2                = La tuile ({0}, {1}) a une taille, un nombre de bandes ou une disposition des valeurs inattendu.
 InsufficientBufferCapacity_3      = La capacit\u00e9 du buffer est insuffisante pour une grille de {0} cellules \u00d7 {1} bandes. Il lui manque {2} \u00e9l\u00e9ments.
 IterationIsFinished               = L\u2019it\u00e9ration est termin\u00e9e.
diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/BandAggregateGridCoverageTest.java b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/BandAggregateGridCoverageTest.java
new file mode 100644
index 0000000000..ade6c2caa5
--- /dev/null
+++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/BandAggregateGridCoverageTest.java
@@ -0,0 +1,137 @@
+/*
+ * 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 java.util.Arrays;
+import java.awt.image.Raster;
+import java.awt.image.DataBufferInt;
+import org.opengis.referencing.datum.PixelInCell;
+import org.apache.sis.coverage.SampleDimension;
+import org.apache.sis.referencing.crs.HardCodedCRS;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.test.TestCase;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+
+/**
+ * Tests {@link BandAggregateGridCoverage}.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.4
+ * @since   1.4
+ */
+public final class BandAggregateGridCoverageTest extends TestCase {
+    /**
+     * Width and height of images created for tests.
+     */
+    private static final int WIDTH = 3, HEIGHT = 2;
+
+    /**
+     * The processor to use for creating the aggregated grid coverage.
+     */
+    private final GridCoverageProcessor processor;
+
+    /**
+     * Creates a new test case.
+     */
+    public BandAggregateGridCoverageTest() {
+        processor = new GridCoverageProcessor();
+    }
+
+    /**
+     * Tests aggregation with two coverages having the same grid geometry.
+     */
+    @Test
+    public void testSameGridGeometry() {
+        final GridCoverage c1 = createCoverage(-2, 4, 3, -1, 100, 200);
+        final GridCoverage c2 = createCoverage(-2, 4, 3, -1, 300);
+        final GridCoverage cr = processor.aggregateRanges(c1, c2);
+        assertEquals(c1.getGridGeometry(), cr.getGridGeometry());
+        assertEquals(c2.getGridGeometry(), cr.getGridGeometry());
+        assertEquals(3, cr.getSampleDimensions().size());
+        assertPixelsEqual(cr, 100, 200, 300, 101, 201, 301, 102, 202, 302,
+                              103, 203, 303, 104, 204, 304, 105, 205, 305);
+    }
+
+    /**
+     * Tests aggregation with two coverages having a translation in their grid extents.
+     * Their "grid to CRS" transforms are the same, which implies that the "real world"
+     * coordinates are different. The intersection is not equal to any source extent.
+     */
+    @Test
+    public void testNoCommonGridGeometry() {
+        final GridCoverage c1 = createCoverage(-2, 4, 3, -1, 100, 200);
+        final GridCoverage c2 = createCoverage(-1, 3, 3, -1, 300);
+        final GridCoverage cr = processor.aggregateRanges(c1, c2);
+        final GridExtent extent = cr.getGridGeometry().getExtent();
+        assertNotEquals(c1.getGridGeometry().getExtent(), extent);
+        assertNotEquals(c2.getGridGeometry().getExtent(), extent);
+        assertEquals(extent(-1, 4, 1, 5), extent);
+        assertEquals(3, cr.getSampleDimensions().size());
+        assertPixelsEqual(cr, 101, 201, 303, 102, 202, 304);
+    }
+
+    /**
+     * Returns a two-dimensional grid extents with the given bounding box.
+     * The maximal coordinates are exclusive.
+     */
+    private static GridExtent extent(final int minX, final int minY, final int maxX, final int maxY) {
+        return new GridExtent(null, new long[] {minX, minY}, new long[] {maxX, maxY}, false);
+    }
+
+    /**
+     * Creates a new grid coverage with bands starting with the given sample values.
+     * The length of the {@code bandValues} array is the number of bands to create.
+     * In a given band <var>b</var>, all pixels have the {@code bandValues[b]}.
+     *
+     * @param  minX        minimal <var>x</var> coordinate value of the grid extent.
+     * @param  minY        minimal <var>y</var> coordinate value of the grid extent.
+     * @param  translateX  <var>x</var> component of the "grid to CRS" translation.
+     * @param  translateY  <var>y</var> component of the "grid to CRS" translation.
+     * @param  bandValues  sample values for the first pixel.
+     * @return a coverage with an image where all pixels have the specified sample values.
+     */
+    private static GridCoverage createCoverage(final int minX, final int minY,
+            final int translateX, final int translateY, final int... bandValues)
+    {
+        final int numBands = bandValues.length;
+        final var samples  = new SampleDimension[numBands];
+        final var builder  = new SampleDimension.Builder();
+        for (int i=0; i<numBands; i++) {
+            samples[i] = builder.setName("Band with value " + bandValues[i]).build();
+        }
+        final int[] data = new int[WIDTH * HEIGHT * numBands];
+        for (int i=0; i<data.length; i++) {
+            data[i] = bandValues[i % numBands] + i / numBands;
+        }
+        final var values = new DataBufferInt(data, data.length);
+        final var domain = new GridGeometry(extent(minX, minY, minX+WIDTH, minY+HEIGHT),
+                PixelInCell.CELL_CORNER, MathTransforms.translation(translateX, translateY), HardCodedCRS.WGS84);
+        return new BufferedGridCoverage(domain, Arrays.asList(samples), values);
+    }
+
+    /**
+     * Asserts that all pixels from the given coverage have the expected values.
+     */
+    private static void assertPixelsEqual(final GridCoverage cr, final int... expected) {
+        final Raster r = cr.render(null).getData();
+        final int[] data = r.getPixels(r.getMinX(), r.getMinY(), r.getWidth(), r.getHeight(), (int[]) null);
+        assertArrayEquals(expected, data);
+    }
+}
diff --git a/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java b/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java
index 3d09e913e5..4f233d0744 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java
@@ -124,6 +124,7 @@ import org.junit.runners.Suite;
     org.apache.sis.coverage.grid.TranslatedGridCoverageTest.class,
     org.apache.sis.coverage.grid.ResampledGridCoverageTest.class,
     org.apache.sis.coverage.grid.DimensionalityReductionTest.class,
+    org.apache.sis.coverage.grid.BandAggregateGridCoverageTest.class,
 
     // Index and processing
     org.apache.sis.index.tree.PointTreeNodeTest.class,
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java
index cc8bff034f..502cf6878a 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java
@@ -29,10 +29,9 @@ import org.apache.sis.storage.DataStoreException;
 import org.apache.sis.coverage.grid.GridCoverage;
 import org.apache.sis.coverage.grid.GridGeometry;
 import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.internal.coverage.CommonDomainFinder;
 import org.apache.sis.internal.storage.MemoryGridResource;
-import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.internal.util.Strings;
-import org.apache.sis.util.Numbers;
 
 
 /**
@@ -48,6 +47,11 @@ import org.apache.sis.util.Numbers;
  * @since   1.3
  */
 final class GridSlice {
+    /**
+     * The "pixel in cell" value used for all "grid to CRS" computations.
+     */
+    static final PixelInCell CELL_ANCHOR = PixelInCell.CELL_CORNER;
+
     /**
      * The resource associated to this slice.
      */
@@ -87,57 +91,6 @@ final class GridSlice {
         offset   = new long[geometry.getDimension()];
     }
 
-    /**
-     * Sets the {@link #offset} terms to the values of the translation columns of the given matrix.
-     * This method shall be invoked if and only if {@link #isIntegerTranslation(Matrix)} returned {@code true}.
-     *
-     * @param  groupToSlice  conversion from source coordinates of {@link GroupByTransform#gridToCRS}
-     *                       to grid coordinates of {@link #geometry}.
-     * @throws ArithmeticException if a translation term is NaN or overflows {@code long} integer capacity.
-     *
-     * @see #getOffset(Map)
-     */
-    private void setOffset(final MatrixSIS groupToSlice) {
-        final int i = groupToSlice.getNumCol() - 1;
-        for (int j=0; j<offset.length; j++) {
-            offset[j] = Numbers.round(groupToSlice.getNumber(j, i));
-        }
-    }
-
-    /**
-     * Returns {@code true} if the given matrix is the identity matrix except for translation terms.
-     * Translation terms must be integer values.
-     *
-     * @param  groupToSlice  conversion from {@link GroupByTransform#gridToCRS}
-     *         source coordinates to {@link #gridToCRS} source coordinates.
-     * @return whether the matrix is identity, ignoring integer translation.
-     *
-     * @see org.apache.sis.referencing.operation.matrix.Matrices#isTranslation(Matrix)
-     */
-    private static boolean isIntegerTranslation(final Matrix groupToSlice) {
-        final int numRows = groupToSlice.getNumRow();
-        final int numCols = groupToSlice.getNumCol();
-        for (int j=0; j<numRows; j++) {
-            for (int i=0; i<numCols; i++) {
-                double tolerance = Numerics.COMPARISON_THRESHOLD;
-                double e = groupToSlice.getElement(j, i);
-                if (i == j) {
-                    e--;
-                } else if (i == numCols - 1) {
-                    final double a = Math.abs(e);
-                    if (a > 1) {
-                        tolerance = Math.min(tolerance*a, 0.125);
-                    }
-                    e -= Math.rint(e);
-                }
-                if (!(Math.abs(e) <= tolerance)) {      // Use `!` for catching NaN.
-                    return false;
-                }
-            }
-        }
-        return true;
-    }
-
     /**
      * Returns the group of objects associated to the CRS and "grid to CRS" transform.
      * The CRS comparisons ignore metadata and transform comparisons ignore integer translations.
@@ -152,14 +105,13 @@ final class GridSlice {
     final GroupByTransform getList(final List<GroupByCRS<GroupByTransform>> groups, final MergeStrategy strategy)
             throws NoninvertibleTransformException
     {
-        final MathTransform gridToCRS = geometry.getGridToCRS(PixelInCell.CELL_CORNER);
+        final MathTransform gridToCRS = geometry.getGridToCRS(CELL_ANCHOR);
         final MathTransform crsToGrid = gridToCRS.inverse();
         final List<GroupByTransform> transforms = GroupByCRS.getOrAdd(groups, geometry).members;
         synchronized (transforms) {
             for (final GroupByTransform c : transforms) {
                 final Matrix groupToSlice = c.linearTransform(crsToGrid);
-                if (groupToSlice != null && isIntegerTranslation(groupToSlice)) {
-                    setOffset(MatrixSIS.castOrCopy(groupToSlice));
+                if (CommonDomainFinder.integerTranslation(groupToSlice, offset) != null) {
                     c.strategy = strategy;
                     return c;
                 }
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java
index 6cecc715f3..137809f741 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java
@@ -22,7 +22,6 @@ import java.util.Arrays;
 import java.util.HashMap;
 import java.util.function.Function;
 import org.opengis.util.InternationalString;
-import org.opengis.referencing.datum.PixelInCell;
 import org.opengis.metadata.spatial.DimensionNameType;
 import org.apache.sis.coverage.grid.GridExtent;
 import org.apache.sis.coverage.grid.GridGeometry;
@@ -137,7 +136,7 @@ final class GridSliceLocator {
             return base;
         }
         extent = new GridExtent(axes, low, high, true);
-        return new GridGeometry(extent, PixelInCell.CELL_CORNER, base.getGridToCRS(PixelInCell.CELL_CORNER),
+        return new GridGeometry(extent, GridSlice.CELL_ANCHOR, base.getGridToCRS(GridSlice.CELL_ANCHOR),
                                 base.isDefined(GridGeometry.CRS) ? base.getCoordinateReferenceSystem() : null);
     }
 
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java
index 540483160e..d8884efe15 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java
@@ -51,7 +51,7 @@ final class GroupByTransform extends Group<GridSlice> {
     private final GridGeometry geometry;
 
     /**
-     * Value or {@code geometry.getGridToCRS(PixelInCell.CELL_CORNER)}.
+     * Value of {@code geometry.getGridToCRS(GridSlice.CELL_ANCHOR)}.
      * All {@linkplain #members} of this group use this transform,
      * possibly with integer differences in translation terms.
      */
@@ -68,7 +68,7 @@ final class GroupByTransform extends Group<GridSlice> {
      * Creates a new group of objects associated to the given transform.
      *
      * @param  geometry   geometry of the grid coverage or resource.
-     * @param  gridToCRS  value or {@code geometry.getGridToCRS(PixelInCell.CELL_CORNER)}.
+     * @param  gridToCRS  value of {@code geometry.getGridToCRS(GridSlice.CELL_ANCHOR)}.
      * @param  strategy   algorithm to apply when more than one grid coverage can be found at the same grid index.
      */
     GroupByTransform(final GridGeometry geometry, final MathTransform gridToCRS, final MergeStrategy strategy) {