sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/03: Allow ResampledGridCoverage to work with grids having more than 2 dimensions.
Date Wed, 27 May 2020 17:07:09 GMT
This is an automated email from the ASF dual-hosted git repository.

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

commit 4bfbfb24ef754fc0f373b73e0cbeef03ed61c73b
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed May 27 14:58:08 2020 +0200

    Allow ResampledGridCoverage to work with grids having more than 2 dimensions.
    
    This fix requires improvement in sis-referencing module for allowing the creation of a CoordinateOperations
    even when the target CRS has more dimensions than the source CRS, provided that constant values (determined
    from the GridGeometry) can be assigned to supplemental dimensions.
    
    This commit contains work from Alexis Manin.
---
 .../coverage/grid/CoordinateOperationFinder.java   | 202 ++++++++++
 .../apache/sis/coverage/grid/GridCoverage2D.java   |  28 +-
 .../org/apache/sis/coverage/grid/GridExtent.java   |   2 +-
 .../sis/coverage/grid/ResampledGridCoverage.java   | 127 +++----
 .../coverage/grid/ResampledGridCoverageTest.java   | 420 ++++++++++++++++++++-
 .../java/org/apache/sis/image/TiledImageMock.java  |  15 +-
 .../java/org/apache/sis/test/FeatureAssert.java    |  41 ++
 .../internal/referencing/CoordinateOperations.java |  10 +-
 .../operation/CoordinateOperationContext.java      |  34 +-
 .../operation/CoordinateOperationFinder.java       | 100 +++--
 .../operation/CoordinateOperationRegistry.java     |   5 +
 .../operation/DefaultPassThroughOperation.java     |   2 +-
 .../referencing/operation/SubOperationInfo.java    | 348 +++++++++++++----
 .../apache/sis/referencing/crs/HardCodedCRS.java   |  24 ++
 .../org/apache/sis/util/resources/Vocabulary.java  |   5 +
 .../sis/util/resources/Vocabulary.properties       |   1 +
 .../sis/util/resources/Vocabulary_fr.properties    |   1 +
 17 files changed, 1107 insertions(+), 258 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
new file mode 100644
index 0000000..ab1ec74
--- /dev/null
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
@@ -0,0 +1,202 @@
+/*
+ * 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.util.function.Supplier;
+import org.opengis.geometry.Envelope;
+import org.opengis.util.FactoryException;
+import org.opengis.referencing.datum.PixelInCell;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.operation.CoordinateOperation;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.util.collection.BackingStoreException;
+import org.apache.sis.internal.referencing.CoordinateOperations;
+import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
+import org.apache.sis.geometry.Envelopes;
+import org.apache.sis.referencing.CRS;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+
+
+/**
+ * Finds a transform from points expressed in the CRS of a source coverage to points in the CRS of a target coverage.
+ * This class differs from {@link CRS#findOperation CRS#findOperation(…)} because of the gridded aspect of inputs.
+ * {@linkplain GridGeometry grid geometries} give more information about how referencing is applied on datasets.
+ * With them, we can detect dimensions where target coordinates are constrained to constant values
+ * because the {@linkplain GridExtent#getSize(int) grid size} is only one cell in those dimensions.
+ * This is an important difference because it allows us to find operations normally impossible,
+ * where we still can produce an operation to a target CRS even if some dimensions have no corresponding source CRS.
+ *
+ * <p><b>Note:</b> this class does not provide complete chain of transformation from grid to grid.
+ * It provides only the operation from the CRS of source to the CRS of destination.</p>
+ *
+ * @author  Alexis Manin (Geomatys)
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+final class CoordinateOperationFinder implements Supplier<double[]> {
+    /**
+     * Whether the operation is between cell centers or cell corners.
+     */
+    private PixelInCell anchor;
+
+    /**
+     * The grid geometry which is the source/target of the coordinate operation to find.
+     */
+    private final GridGeometry source, target;
+
+    /**
+     * The target coordinate values, computed only if needed.
+     */
+    private double[] coordinates;
+
+    /**
+     * The coordinate operation from source to target CRS, computed when first needed.
+     */
+    private CoordinateOperation operation;
+
+    /**
+     * Creates a new finder.
+     *
+     * @param  source  the grid geometry which is the source of the coordinate operation to find.
+     * @param  target  the grid geometry which is the target of the coordinate operation to find.
+     */
+    CoordinateOperationFinder(final GridGeometry source, final GridGeometry target) {
+        this.source = source;
+        this.target = target;
+    }
+
+    /**
+     * Returns the CRS of the source grid geometry. If neither the source and target grid geometry
+     * define a CRS, then this method returns {@code null}.
+     *
+     * @throws IncompleteGridGeometryException if the target grid geometry has a CRS but the source
+     *         grid geometry has none. Note that the converse is allowed, in which case the target
+     *         CRS is assumed the same than the source.
+     */
+    private CoordinateReferenceSystem getSourceCRS() {
+        return source.isDefined(GridGeometry.CRS) ||
+               target.isDefined(GridGeometry.CRS) ? source.getCoordinateReferenceSystem() : null;
+    }
+
+    /**
+     * Returns the target of the "corner to CRS" transform.
+     * May be {@code null} if the neither the source and target grid geometry define a CRS.
+     *
+     * @throws IncompleteGridGeometryException if the target grid geometry has a CRS but the source
+     *         grid geometry has none. Note that the converse is allowed, in which case the target
+     *         CRS is assumed the same than the source.
+     */
+    final CoordinateReferenceSystem getTargetCRS() {
+        return (operation != null) ? operation.getTargetCRS() : getSourceCRS();
+    }
+
+    /**
+     * Computes the transform from the grid coordinates of the source to geospatial coordinates of the target.
+     * It may be the identity operation. We try to take envelopes in account because the operation choice may
+     * depend on the geographic area.
+     *
+     * @param  anchor  whether the operation is between cell centers or cell corners.
+     * @return operation from source CRS to target CRS, or {@code null} if CRS are unspecified or equivalent.
+     * @throws FactoryException if no operation can be found between the source and target CRS.
+     * @throws TransformException if some coordinates can not be transformed to the specified target.
+     * @throws IncompleteGridGeometryException if required CRS or a "grid to CRS" information is missing.
+     */
+    final MathTransform gridToCRS(final PixelInCell anchor) throws FactoryException, TransformException {
+        /*
+         * If `coordinates` is non-null, it means that the `get()` method has been invoked in previous call
+         * to this `cornerToCRS(…)` method, which implies that `operation` depends on the `anchor` value.
+         * In such case we need to discard the previous `operation` value and recompute it.
+         */
+        this.anchor = anchor;
+        if (operation == null || coordinates != null) {
+            operation   = null;
+            coordinates = null;
+            final CoordinateReferenceSystem sourceCRS = getSourceCRS();
+            final CoordinateReferenceSystem targetCRS = target.isDefined(GridGeometry.CRS) ?
+                                                        target.getCoordinateReferenceSystem() : sourceCRS;
+            final Envelope sourceEnvelope = source.envelope;
+            final Envelope targetEnvelope = target.envelope;
+            try {
+                CoordinateOperations.CONSTANT_COORDINATES.set(this);
+                if (sourceEnvelope != null && targetEnvelope != null) {
+                    operation = Envelopes.findOperation(sourceEnvelope, targetEnvelope);
+                }
+                if (operation == null && sourceCRS != null) {
+                    DefaultGeographicBoundingBox areaOfInterest = null;
+                    if (sourceEnvelope != null || targetEnvelope != null) {
+                        areaOfInterest = new DefaultGeographicBoundingBox();
+                        areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope : sourceEnvelope);
+                    }
+                    operation = CRS.findOperation(sourceCRS, targetCRS, areaOfInterest);
+                }
+            } finally {
+                CoordinateOperations.CONSTANT_COORDINATES.remove();
+            }
+        }
+        /*
+         * The following line may throw IncompleteGridGeometryException, which is desired because
+         * if that transform is missing, we can not continue (we have no way to guess it).
+         */
+        MathTransform tr = source.getGridToCRS(anchor);
+        if (operation != null) {
+            tr = MathTransforms.concatenate(tr, operation.getMathTransform());
+        }
+        return tr;
+    }
+
+    /**
+     * Invoked when the target CRS have some dimensions that the source CRS does not have.
+     * For example this is invoked during the conversion from (<var>x</var>, <var>y</var>)
+     * coordinates to (<var>x</var>, <var>y</var>, <var>t</var>). If constant values can
+     * be given to the missing dimensions, than those values are returned. Otherwise this
+     * method returns {@code null}.
+     */
+    @Override
+    @SuppressWarnings("ReturnOfCollectionOrArrayField")
+    public double[] get() {
+        if (coordinates == null && target.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS)) {
+            final MathTransform gridToCRS = target.getGridToCRS(anchor);
+            coordinates = new double[gridToCRS.getTargetDimensions()];
+            double[] gc = new double[gridToCRS.getSourceDimensions()];
+            Arrays.fill(gc, Double.NaN);
+            final GridExtent extent = target.getExtent();
+            for (int i=0; i<gc.length; i++) {
+                final long low = extent.getLow(i);
+                if (low == extent.getHigh(i)) {
+                    gc[i] = low;
+                }
+            }
+            /*
+             * At this point, the only grid coordinates with finite values are the ones where the
+             * grid size is one cell (i.e. conversion to target CRS can produce only one value).
+             * After conversion with `gridToCRS`, the corresponding target dimensions will have
+             * non-NaN coordinate values only if they really do not depend on any dimension other
+             * than the one having a grid size of 1.
+             */
+            try {
+                gridToCRS.transform(gc, 0, coordinates, 0, 1);
+            } catch (TransformException e) {
+                throw new BackingStoreException(e);
+            }
+        }
+        return coordinates;
+    }
+}
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
index ef953b5..8462b97 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
@@ -160,24 +160,20 @@ public class GridCoverage2D extends GridCoverage {
     /**
      * Creates a new grid coverage for the resampling of specified source coverage.
      *
-     * @param  source       the coverage containing source values.
-     * @param  domain       the grid extent, CRS and conversion from cell indices to CRS.
-     * @param  extent       the {@code domain.getExtent()} value.
-     * @param  data         the sample values as a {@link RenderedImage}, with one band for each sample dimension.
-     * @param  xDimension   index of extent dimension for image <var>x</var> coordinates.
-     * @param  yDimension   index of extent dimension for image <var>y</var> coordinates.
+     * @param  source  the coverage containing source values.
+     * @param  domain  the grid extent, CRS and conversion from cell indices to CRS.
+     * @param  extent  the {@code domain.getExtent()} value.
+     * @param  data    the sample values as a {@link RenderedImage}, with one band for each sample dimension.
      */
-    GridCoverage2D(final GridCoverage source, final GridGeometry domain, final GridExtent extent,
-                   RenderedImage data, final int xDimension, final int yDimension)
-    {
+    GridCoverage2D(final GridCoverage source, final GridGeometry domain, final GridExtent extent, RenderedImage data) {
         super(source, domain);
-        data = unwrapIfSameSize(data);
-        this.data       = data;
-        this.xDimension = xDimension;
-        this.yDimension = yDimension;
-        gridToImageX    = subtractExact(data.getMinX(), extent.getLow(xDimension));
-        gridToImageY    = subtractExact(data.getMinY(), extent.getLow(yDimension));
-        gridGeometry2D  = new AtomicReference<>();
+        final int[] imageAxes = extent.getSubspaceDimensions(BIDIMENSIONAL);
+        xDimension     = imageAxes[0];
+        yDimension     = imageAxes[1];
+        this.data      = data = unwrapIfSameSize(data);
+        gridToImageX   = subtractExact(data.getMinX(), extent.getLow(xDimension));
+        gridToImageY   = subtractExact(data.getMinY(), extent.getLow(yDimension));
+        gridGeometry2D = new AtomicReference<>();
     }
 
     /**
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 5b258f6..fa98774 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
@@ -885,7 +885,7 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable
      *
      * @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 cell centers instead than cell corners.
+     *                      If different, then this is assumed to map cell centers instead of cell corners.
      * @param  fallback     bounds to use if some values still NaN, or {@code null} if none.
      * @return this grid extent in real world coordinates.
      * @throws TransformException if the envelope can not be computed with the given transform.
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
index 8ca6ebb..6d67c06 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
@@ -25,7 +25,6 @@ import org.opengis.util.FactoryException;
 import org.opengis.coverage.CannotEvaluateException;
 import org.opengis.referencing.datum.PixelInCell;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
-import org.opengis.referencing.operation.CoordinateOperation;
 import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.operation.MathTransform;
 import org.apache.sis.geometry.Envelopes;
@@ -34,12 +33,11 @@ import org.apache.sis.coverage.SampleDimension;
 import org.apache.sis.geometry.GeneralEnvelope;
 import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
-import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
 import org.apache.sis.referencing.operation.transform.LinearTransform;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.transform.TransformSeparator;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 import org.apache.sis.referencing.operation.matrix.Matrices;
-import org.apache.sis.referencing.CRS;
 import org.apache.sis.util.ComparisonMode;
 import org.apache.sis.util.Utilities;
 
@@ -80,13 +78,6 @@ final class ResampledGridCoverage extends GridCoverage {
     private final ImageProcessor imageProcessor;
 
     /**
-     * Indices of extent dimensions corresponding to image <var>x</var> and <var>y</var> coordinates.
-     * Typical values are 0 for {@code xDimension} and 1 for {@code yDimension}, but different values
-     * are allowed. This class select the dimensions having largest size.
-     */
-    private final int xDimension, yDimension;
-
-    /**
      * Creates a new grid coverage which will be the resampling of the given source.
      *
      * @param  source          the coverage to resample.
@@ -104,26 +95,6 @@ final class ResampledGridCoverage extends GridCoverage {
         this.source         = source;
         this.toSourceCorner = toSourceCorner;
         this.toSourceCenter = toSourceCenter;
-        final GridExtent extent = domain.getExtent();
-        long size1 = 0; int idx1 = 0;
-        long size2 = 0; int idx2 = 1;
-        final int dimension = extent.getDimension();
-        for (int i=0; i<dimension; i++) {
-            final long size = extent.getSize(i);
-            if (size > size1) {
-                size2 = size1; idx2 = idx1;
-                size1 = size;  idx1 = i;
-            } else if (size > size2) {
-                size2 = size;  idx2 = i;
-            }
-        }
-        if (idx1 < idx2) {          // Keep (x,y) dimensions in the order they appear.
-            xDimension = idx1;
-            yDimension = idx2;
-        } else {
-            xDimension = idx2;
-            yDimension = idx1;
-        }
         /*
          * Get fill values from background values declared for each band, if any.
          * If no background value is declared, default is 0 for integer data or
@@ -173,10 +144,10 @@ final class ResampledGridCoverage extends GridCoverage {
                 if (sourceGG.equals(targetGG, ComparisonMode.APPROXIMATE)) {
                     return source;
                 }
-                return new GridCoverage2D(source, targetGG, extent, source.render(null), xDimension, yDimension);
+                return new GridCoverage2D(source, targetGG, extent, source.render(null));
             }
         }
-        return new GridCoverage2D(source, gridGeometry, extent, render(null), xDimension, yDimension);
+        return new GridCoverage2D(source, gridGeometry, extent, render(null));
     }
 
     /**
@@ -248,42 +219,15 @@ final class ResampledGridCoverage extends GridCoverage {
     static GridCoverage create(final GridCoverage source, final GridGeometry target, final ImageProcessor processor)
             throws FactoryException, TransformException
     {
-        final CoordinateReferenceSystem sourceCRS = source.getCoordinateReferenceSystem();
-        final CoordinateReferenceSystem targetCRS = target.isDefined(GridGeometry.CRS) ?
-                                                    target.getCoordinateReferenceSystem() : sourceCRS;
-        /*
-         * Get the coordinate operation from source CRS to target CRS. It may be the identity operation,
-         * or null only if there is not enough information for determining the operation. We try to take
-         * envelopes in account because the operation choice may depend on the geographic area.
-         */
-        CoordinateOperation changeOfCRS = null;
         final GridGeometry sourceGG = source.getGridGeometry();
-        if (sourceGG.isDefined(GridGeometry.ENVELOPE) && target.isDefined(GridGeometry.ENVELOPE)) {
-            changeOfCRS = Envelopes.findOperation(sourceGG.getEnvelope(), target.getEnvelope());
-        }
-        if (changeOfCRS == null && sourceCRS != null && targetCRS != null) try {
-            DefaultGeographicBoundingBox areaOfInterest = null;
-            if (sourceGG.isDefined(GridGeometry.ENVELOPE)) {
-                areaOfInterest = new DefaultGeographicBoundingBox();
-                areaOfInterest.setBounds(sourceGG.getEnvelope());
-            }
-            changeOfCRS = CRS.findOperation(sourceCRS, targetCRS, areaOfInterest);
-        } catch (IncompleteGridGeometryException e) {
-            // Happen if the source GridCoverage does not define a CRS.
-            GridCoverageProcessor.recoverableException("resample", e);
-        }
+        final CoordinateOperationFinder changeOfCRS = new CoordinateOperationFinder(sourceGG, target);
         /*
          * Compute the transform from source pixels to target CRS (to be completed to target pixels later).
-         * The following line may throw IncompleteGridGeometryException, which is desired because if that
+         * The following lines may throw IncompleteGridGeometryException, which is desired because if that
          * transform is missing, we can not continue (we have no way to guess it).
          */
-        MathTransform sourceCornerToCRS = sourceGG.getGridToCRS(PixelInCell.CELL_CORNER);
-        MathTransform sourceCenterToCRS = sourceGG.getGridToCRS(PixelInCell.CELL_CENTER);
-        if (changeOfCRS != null) {
-            final MathTransform tr = changeOfCRS.getMathTransform();
-            sourceCornerToCRS = MathTransforms.concatenate(sourceCornerToCRS, tr);
-            sourceCenterToCRS = MathTransforms.concatenate(sourceCenterToCRS, tr);
-        }
+        final MathTransform sourceCornerToCRS = changeOfCRS.gridToCRS(PixelInCell.CELL_CORNER);
+        final MathTransform sourceCenterToCRS = changeOfCRS.gridToCRS(PixelInCell.CELL_CENTER);
         /*
          * Compute the transform from target grid to target CRS. This transform may be unspecified,
          * in which case we need to compute a default transform trying to preserve resolution at the
@@ -424,6 +368,7 @@ final class ResampledGridCoverage extends GridCoverage {
         GridGeometry resampled = target;
         ComparisonMode mode = ComparisonMode.IGNORE_METADATA;
         if (!target.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS | GridGeometry.CRS)) {
+            final CoordinateReferenceSystem targetCRS = changeOfCRS.getTargetCRS();
             resampled = new GridGeometry(targetExtent, PixelInCell.CELL_CENTER, targetCenterToCRS, targetCRS);
             mode = ComparisonMode.APPROXIMATE;
             if (target.isDefined(GridGeometry.ENVELOPE)) {
@@ -478,37 +423,53 @@ final class ResampledGridCoverage extends GridCoverage {
      */
     @Override
     public RenderedImage render(GridExtent sliceExtent) {
-        final RenderedImage values;         // Source image providing the values to resample.
-        final MathTransform toSource;       // From resampled image pixels to source image pixels.
+        if (sliceExtent == null) {
+            sliceExtent = gridGeometry.getExtent();
+        }
         final Rectangle     bounds;         // Bounds (in pixel coordinates) of resampled image.
+        final MathTransform toSource;       // From resampled image pixels to source image pixels.
+        final GridExtent    sourceExtent;
         try {
-            GridExtent sourceExtent = null;
-            if (sliceExtent != null) {
-                final GeneralEnvelope envelope = sliceExtent.toCRS(toSourceCorner, toSourceCenter, null);
-                sourceExtent = new GridExtent(envelope, GridRoundingMode.ENCLOSING, null, null, null);
-            }
-            values = source.render(sourceExtent);
-            if (sliceExtent  == null) sliceExtent  = gridGeometry.getExtent();
-            if (sourceExtent == null) sourceExtent = source.getGridGeometry().getExtent();
-            bounds = new Rectangle(Math.toIntExact(sliceExtent.getSize(xDimension)),
-                                   Math.toIntExact(sliceExtent.getSize(yDimension)));
+            final GeneralEnvelope sourceBounds = sliceExtent.toCRS(toSourceCorner, toSourceCenter, null);
+            sourceExtent = new GridExtent(sourceBounds, GridRoundingMode.ENCLOSING, null, null, null);
+            final int[] resampledDimensions = sliceExtent.getSubspaceDimensions(BIDIMENSIONAL);
+            final int[] sourceDimensions  =  sourceExtent.getSubspaceDimensions(BIDIMENSIONAL);
+            bounds = new Rectangle(Math.toIntExact(sliceExtent.getSize(resampledDimensions[0])),
+                                   Math.toIntExact(sliceExtent.getSize(resampledDimensions[1])));
+            /*
+             * The transform needs to be two-dimensional. The `toSourceCenter` transform does not met that condition,
+             * otherwise `specialize(…)` would have replaced this ResampledGridCoverage by a GridCoverage2D instance.
+             * Try to extract a two-dimensional part operating only on the slice dimensions having an extent larger
+             * than one cell. The choice of dimensions may vary between different calls to this `render(…)` method,
+             * depending on `sliceExtent` value.
+             */
+            final TransformSeparator sep = new TransformSeparator(toSourceCenter);
+            sep.addSourceDimensions(resampledDimensions);
+            sep.addTargetDimensions(sourceDimensions);
+            final MathTransform toSourceSlice = sep.separate();
             /*
              * `this.toSource` is a transform from source cell coordinates to target cell coordinates.
              * We need a transform from source pixel coordinates to target pixel coordinates (in images).
              * An offset may exist between cell coordinates and pixel coordinates.
              */
-            final MathTransform pixelsToTransform = MathTransforms.translation(
-                    sliceExtent.getLow(xDimension),
-                    sliceExtent.getLow(yDimension));
+            final MathTransform resampledToGrid = MathTransforms.translation(
+                    sliceExtent.getLow(resampledDimensions[0]),
+                    sliceExtent.getLow(resampledDimensions[1]));
 
-            final MathTransform transformToPixels = MathTransforms.translation(
-                    Math.negateExact(sourceExtent.getLow(xDimension)),
-                    Math.negateExact(sourceExtent.getLow(yDimension)));
+            final MathTransform gridToSource = MathTransforms.translation(
+                    Math.negateExact(sourceExtent.getLow(sourceDimensions[0])),
+                    Math.negateExact(sourceExtent.getLow(sourceDimensions[1])));
 
-            toSource = MathTransforms.concatenate(pixelsToTransform, toSourceCenter, transformToPixels);
-        } catch (TransformException | ArithmeticException e) {
+            toSource = MathTransforms.concatenate(resampledToGrid, toSourceSlice, gridToSource);
+        } catch (FactoryException | TransformException | ArithmeticException e) {
             throw new CannotEvaluateException(e.getLocalizedMessage(), e);
         }
+        /*
+         * Following call is potentially costly, depending on `source` implementation.
+         * For example it may cause loading of tiles from a file. For this reason we
+         * call this method only here, when remaining operations are unlikely to fail.
+         */
+        final RenderedImage values = source.render(sourceExtent);
         return imageProcessor.resample(bounds, toSource, values);
     }
 }
diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
index b08d22d..8e47a12 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
@@ -16,33 +16,51 @@
  */
 package org.apache.sis.coverage.grid;
 
+import java.util.Arrays;
 import java.util.Random;
+import java.awt.Color;
+import java.awt.Rectangle;
+import java.awt.image.BufferedImage;
 import java.awt.image.DataBuffer;
+import java.awt.image.Raster;
 import java.awt.image.RenderedImage;
+import java.awt.image.WritableRaster;
+import org.opengis.referencing.cs.AxisDirection;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
+import org.opengis.util.FactoryException;
 import org.apache.sis.geometry.Envelope2D;
 import org.apache.sis.geometry.ImmutableEnvelope;
 import org.apache.sis.image.Interpolation;
 import org.apache.sis.image.TiledImageMock;
+import org.apache.sis.internal.coverage.j2d.TiledImage;
 import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.referencing.CommonCRS;
 import org.apache.sis.referencing.crs.HardCodedCRS;
 import org.apache.sis.referencing.operation.HardCodedConversions;
+import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.transform.TransformSeparator;
 import org.apache.sis.test.TestUtilities;
 import org.apache.sis.test.DependsOn;
+import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.TestCase;
 import org.junit.Test;
 
 import static org.opengis.referencing.datum.PixelInCell.CELL_CENTER;
-import static org.opengis.test.Assert.*;
+import static org.apache.sis.test.FeatureAssert.*;
 
 
 /**
  * Tests the {@link ResampledGridCoverage} implementation.
- * The tests in this class does not verify interpolation values
+ * The tests in this class does not verify interpolation results
  * (this is {@link org.apache.sis.image.ResampledImageTest} job).
  * Instead it focus on the grid geometry inferred by the operation.
  *
  * @author  Martin Desruisseaux (Geomatys)
+ * @author  Alexis Manin (Geomatys)
  * @version 1.1
  * @since   1.1
  * @module
@@ -50,19 +68,37 @@ import static org.opengis.test.Assert.*;
 @DependsOn(org.apache.sis.image.ResampledImageTest.class)
 public final strictfp class ResampledGridCoverageTest extends TestCase {
     /**
+     * The random number generator used for generating some grid coverage values.
+     * Created only if needed.
+     */
+    private Random random;
+
+    /**
+     * Arbitrary non-zero pixel coordinates for image origin.
+     */
+    private int minX = -2, minY = -2;
+
+    /**
+     * Arbitrary non-zero grid coordinate for the <var>z</var> and <var>t</var> dimensions.
+     */
+    private int gridZ, gridT;
+
+    /**
      * Creates a small grid coverage with arbitrary data. The rendered image will
      * have only one tile since testing tiling is not the purpose of this class.
+     * This simple coverage is two-dimensional.
      */
-    private static GridCoverage2D createGridCoverage() {
-        final Random random  = TestUtilities.createRandomNumberGenerator();
+    private GridCoverage2D createCoverage2D() {
+        random = TestUtilities.createRandomNumberGenerator();
         final int width  = random.nextInt(8) + 3;
         final int height = random.nextInt(8) + 3;
+        minX = random.nextInt(32) - 10;
+        minY = random.nextInt(32) - 10;
         final TiledImageMock image = new TiledImageMock(
                 DataBuffer.TYPE_USHORT, 2,      // dataType and numBands
-                random.nextInt(32) - 10,        // minX
-                random.nextInt(32) - 10,        // minY
-                width, height,
-                width, height,
+                minX,  minY,
+                width, height,                  // Image size
+                width, height,                  // Tile size
                 random.nextInt(32) - 10,        // minTileX
                 random.nextInt(32) - 10);       // minTileY
         image.validate();
@@ -76,6 +112,125 @@ public final strictfp class ResampledGridCoverageTest extends TestCase {
     }
 
     /**
+     * Size of a quadrant in the coverage created by {@link #createCoverageND(boolean)}.
+     * The total image width and height are {@code 2*Q}.
+     */
+    private static final int QS = 3;
+
+    /**
+     * Creates a coverage in {@linkplain HardCodedCRS#WGS84_3D OGC:CRS:84 + elevation} reference system.
+     * If the {@code withTime} argument is {@code true}, then the coverage will also include a temporal
+     * dimension. The grid coverage characteristics are:
+     * <ul>
+     *   <li>Dimension is 6×6.</li>
+     *   <li>Grid extent is zero-based.</li>
+     *   <li>Envelope is arbitrary but stable (no random values).</li>
+     *   <li>Display oriented (origin is in upper-left corner).</li>
+     *   <li>3 byte bands for RGB coloration.</li>
+     *   <li>Each quarter of the overall image is filled with a plain color:
+     *     <table style="color:white;border-collapse:collapse;">
+     *       <tbody style="border:none">
+     *         <tr>
+     *           <td style="width:50%; background-color:black">Black</td>
+     *           <td style="width:50%; background-color:red">Red</td>
+     *         </tr>
+     *         <tr>
+     *           <td style="width:50%; background-color:green">Green</td>
+     *           <td style="width:50%; background-color:blue">Blue</td>
+     *         </tr>
+     *       </tbody>
+     *     </table>
+     *   </li>
+     * </ul>
+     *
+     * @param  withTime  {@code false} for a three-dimensional coverage, or {@code true} for adding a temporal dimension.
+     * @return a new three- or four-dimensional RGB Grid Coverage.
+     */
+    private GridCoverage createCoverageND(final boolean withTime) {
+        random = TestUtilities.createRandomNumberGenerator();
+        final BufferedImage image = new BufferedImage(2*QS, 2*QS, BufferedImage.TYPE_3BYTE_BGR);
+        final int[] color = new int[QS*QS];
+        /* Upper-left  quarter */ // Keep default value, which is black.
+        /* Upper-right quarter */ Arrays.fill(color, Color.RED  .getRGB()); image.setRGB(QS,  0, QS, QS, color, 0, QS);
+        /* Lower-left  quarter */ Arrays.fill(color, Color.GREEN.getRGB()); image.setRGB( 0, QS, QS, QS, color, 0, QS);
+        /* Lower-right quarter */ Arrays.fill(color, Color.BLUE .getRGB()); image.setRGB(QS, QS, QS, QS, color, 0, QS);
+        /*
+         * Create an image with origin different than 0. The image has a single tile
+         * with pixel coordinates from (-2, -2) inclusive to (4, 4) exclusive.
+         */
+        minX = -2;
+        minY = -2;
+        GridGeometry gg = createGridGeometryND(withTime ? HardCodedCRS.WGS84_4D : HardCodedCRS.WGS84_3D, 0, 1, 2, 3, false);
+        final TiledImage shiftedImage = new TiledImage(
+                image.getColorModel(),
+                image.getWidth(), image.getHeight(),        // Image size
+                random.nextInt(32) - 10,                    // minTileX
+                random.nextInt(32) - 10,                    // minTileY
+                image.getRaster().createTranslatedChild(minX, minY));
+        return new GridCoverage2D(gg, null, shiftedImage);
+    }
+
+    /**
+     * Creates the grid geometry associated with {@link #createCoverageND(boolean)}, optionally with swapped
+     * horizontal axes and flipped Y axis. The given CRS shall have 3 or 4 dimensions.
+     *
+     * @param  crs    the coordinate reference system to assign to the grid geometry.
+     * @param  x      dimension of <var>x</var> coordinates (typically 0).
+     * @param  y      dimension of <var>y</var> coordinates (typically 1).
+     * @param  z      dimension of <var>z</var> coordinates (typically 2).
+     * @param  t      dimension of <var>t</var> coordinates (typically 3). Ignored if the CRS is not four-dimensional.
+     * @param  flipY  whether to flip the <var>y</var> axis.
+     */
+    private GridGeometry createGridGeometryND(final CoordinateReferenceSystem crs,
+            final int x, final int y, final int z, final int t, final boolean flipY)
+    {
+        final int dim = crs.getCoordinateSystem().getDimension();
+        final long[] lower = new long[dim];
+        final long[] upper = new long[dim];
+        Arrays.fill(upper, StrictMath.min(x,y), StrictMath.max(x,y)+1, 2*QS - 1);
+        final MatrixSIS gridToCRS = Matrices.createIdentity(dim + 1);
+        gridToCRS.setElement(x, x,    44./(2*QS));  // X scale
+        gridToCRS.setElement(x, dim, -50./(2*QS));  // X translation
+        gridToCRS.setElement(y, y,   -3.5);         // Y scale
+        gridToCRS.setElement(y, dim, -0.75);        // Y translation
+        gridToCRS.setElement(z, dim, -100);
+        lower[z] = upper[z] = gridZ = 7;            // Arbitrary non-zero position in the grid.
+        if (t < dim) {
+            gridToCRS.setElement(t, dim, 48055);
+            lower[t] = upper[t] = gridT = 12;
+        }
+        if (flipY) {
+            gridToCRS.setElement(y, y,     3.5);    // Inverse sign.
+            gridToCRS.setElement(y, dim, -18.25);
+        }
+        return new GridGeometry(new GridExtent(null, lower, upper, true),
+                    CELL_CENTER, MathTransforms.linear(gridToCRS), crs);
+    }
+
+    /**
+     * Verifies that the given target coverage has the same pixel values than the source coverage.
+     * This method opportunistically verifies that the target {@link GridCoverage} instance has a
+     * {@link GridCoverage#render(GridExtent)} implementation conforms to the specification, i.e.
+     * that requesting only a sub-area results in an image where pixel coordinate (0,0) corresponds
+     * to cell coordinates in the lower corner of specified {@code sliceExtent}.
+     */
+    private void assertContentEquals(final GridCoverage source, final GridCoverage target) {
+        final int tx = random.nextInt(3);
+        final int ty = random.nextInt(3);
+        final GridExtent sourceExtent = source.gridGeometry.getExtent();
+        final int newWidth   = (int) sourceExtent.getSize(0) - tx;
+        final int newHeight  = (int) sourceExtent.getSize(1) - ty;
+        GridExtent subExtent = new GridExtent(
+                (int) sourceExtent.getLow(0) + tx,
+                (int) sourceExtent.getLow(1) + ty,
+                newWidth,
+                newHeight
+        );
+        assertPixelsEqual(source.render(null), new Rectangle(tx, ty, newWidth, newHeight),
+                          target.render(subExtent), new Rectangle(newWidth, newHeight));
+    }
+
+    /**
      * Returns a resampled coverage using processor with default configuration.
      * We use processor instead than instantiating {@link ResampledGridCoverage} directly in order
      * to test {@link GridCoverageProcessor#resample(GridCoverage, GridGeometry)} method as well.
@@ -87,42 +242,48 @@ public final strictfp class ResampledGridCoverageTest extends TestCase {
     }
 
     /**
-     * Tests application of an identity transform from an explicitly specified grid geometry.
+     * Tests application of an identity transform computed from an explicitly given "grid to CRS" transform.
      * We expect the source coverage to be returned unchanged.
      *
      * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
      */
     @Test
     public void testExplicitIdentity() throws TransformException {
-        final GridCoverage2D source = createGridCoverage();
+        final GridCoverage2D source = createCoverage2D();
         GridGeometry gg = source.getGridGeometry();
         gg = new GridGeometry(null, CELL_CENTER, gg.getGridToCRS(CELL_CENTER), gg.getCoordinateReferenceSystem());
         final GridCoverage target = resample(source, gg);
         assertSame("Identity transform should result in same coverage.", source, target);
+        assertContentEquals(source, target);
     }
 
     /**
      * Tests application of an identity transform without specifying explicitly the desired grid geometry.
+     * This test is identical to {@link #testExplicitIdentity()} except that the "grid to CRS" transform
+     * specified to the {@code resample(…)} operation is null.
      *
      * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
      */
     @Test
+    @DependsOnMethod("testExplicitIdentity")
     public void testImplicitIdentity() throws TransformException {
-        final GridCoverage2D source = createGridCoverage();
+        final GridCoverage2D source = createCoverage2D();
         GridGeometry gg = source.getGridGeometry();
         gg = new GridGeometry(null, CELL_CENTER, null, gg.getCoordinateReferenceSystem());
         final GridCoverage target = resample(source, gg);
         assertSame("Identity transform should result in same coverage.", source, target);
+        assertContentEquals(source, target);
     }
 
     /**
-     * Tests application of an axis swapping.
+     * Tests application of axis swapping in a two-dimensional coverage.
+     * This test verifies the envelope of resampled coverage.
      *
      * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
      */
     @Test
     public void testAxisSwap() throws TransformException {
-        final GridCoverage2D source = createGridCoverage();
+        final GridCoverage2D source = createCoverage2D();
         GridGeometry gg = new GridGeometry(null, CELL_CENTER, null, HardCodedCRS.WGS84_φλ);
         final GridCoverage target = resample(source, gg);
         /*
@@ -153,13 +314,119 @@ public final strictfp class ResampledGridCoverageTest extends TestCase {
     }
 
     /**
-     * Tests resampling of a sub-region.
+     * Tests application of axis swapping in a three-dimensional coverage, together with an axis flip.
+     * This test verifies that the pixel values of resampled coverage are found in expected quadrant.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    @DependsOnMethod("testAxisSwap")
+    public void testAxisSwapAndFlip() throws TransformException {
+        final GridCoverage  source      = createCoverageND(false);
+        final GridGeometry  target      = createGridGeometryND(CommonCRS.WGS84.geographic3D(), 1, 0, 2, 3, true);
+        final GridCoverage  result      = resample(source, target);
+        final RenderedImage sourceImage = source.render(null);
+        final RenderedImage targetImage = result.render(null);
+        assertEquals(target, result.getGridGeometry());
+        assertEquals("minX", 0, sourceImage.getMinX());                     // As per GridCoverage.render(…) contract.
+        assertEquals("minY", 0, sourceImage.getMinY());
+        assertEquals("minX", 0, targetImage.getMinX());
+        assertEquals("minY", 0, targetImage.getMinY());
+        assertPixelsEqual(sourceImage, new Rectangle( 0, QS, QS, QS),
+                          targetImage, new Rectangle( 0,  0, QS, QS));      // Green should be top-left.
+        assertPixelsEqual(sourceImage, new Rectangle( 0,  0, QS, QS),
+                          targetImage, new Rectangle(QS,  0, QS, QS));      // Black should be upper-right.
+        assertPixelsEqual(sourceImage, new Rectangle(QS, QS, QS, QS),
+                          targetImage, new Rectangle( 0, QS, QS, QS));      // Blue should be lower-left.
+        assertPixelsEqual(sourceImage, new Rectangle(QS,  0, QS, QS),
+                          targetImage, new Rectangle(QS, QS, QS, QS));      // Red should be lower-right.
+    }
+
+    /**
+     * Tests an operation moving the dimension of temporal axis in a four-dimensional coverage.
      *
      * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
      */
     @Test
-    public void testSubArea() throws TransformException {
-        final GridCoverage2D source = createGridCoverage();           // Envelope2D(20, 15, 60, 62)
+    public void testTemporalAxisMoved() throws TransformException {
+        final GridCoverage source = createCoverageND(true);
+        final GridGeometry target = createGridGeometryND(HardCodedCRS.TIME_WGS84, 1, 2, 3, 0, false);
+        final GridCoverage result = resample(source, target);
+        assertAxisDirectionsEqual("Expected (t,λ,φ,H) axes.",
+                result.getGridGeometry().getCoordinateReferenceSystem().getCoordinateSystem(),
+                AxisDirection.FUTURE, AxisDirection.EAST, AxisDirection.NORTH, AxisDirection.UP);
+
+        assertPixelsEqual(source.render(null), null, result.render(null), null);
+    }
+
+    /**
+     * Tests resampling in a sub-region specified by a grid extent. This method uses a three-dimensional coverage,
+     * which implies that this method also tests the capability to identify which slice needs to be resampled.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    public void testSubGridExtent() throws TransformException {
+        final GridCoverage source     = createCoverageND(false);
+        final GridGeometry sourceGeom = source.getGridGeometry();
+        final GridGeometry targetGeom = new GridGeometry(
+                new GridExtent(null, new long[] {2, 2, gridZ}, new long[] {5, 5, gridZ}, true),
+                CELL_CENTER, sourceGeom.gridToCRS,
+                sourceGeom.getCoordinateReferenceSystem());
+
+        final GridCoverage result = resample(source, targetGeom);
+        assertEquals(targetGeom, result.getGridGeometry());
+        /*
+         * Verify that the target coverage contains all pixel values of the source coverage.
+         * Iteration over source pixels needs to be restricted to the `targetGeom` extent.
+         */
+        final RenderedImage sourceImage = source.render(null);
+        RenderedImage targetImage = result.render(null);
+        assertPixelsEqual(sourceImage, new Rectangle(2, 2, 4, 4),
+                          targetImage, null);
+        /*
+         * Verify GridCoverage.render(GridExtent) contract: the origin of the returned image
+         * shall be the lower-left corner of `sliceExtent`, which is (3,3) in this test.
+         */
+        targetImage = result.render(new GridExtent(3, 3, 2, 2));
+        assertPixelsEqual(sourceImage, new Rectangle (3, 3, 2, 2),
+                          targetImage, new Rectangle (0, 0, 2, 2));
+    }
+
+    /**
+     * Tests resampling in a sub-region specified by a grid extent spanning a single column.
+     * When trying to optimize resampling by dropping dimensions, it can happen that transform dimensions
+     * are reduced to 1D. However, it is a problem for image case which requires 2D coordinates.
+     * So we must ensure that resample conversion keeps at least two dimensions.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    public void testSubGridExtentColumnar() throws TransformException {
+        final GridCoverage2D source   = createCoverage2D();
+        final GridGeometry sourceGeom = source.getGridGeometry();
+        final GridExtent sourceExtent = sourceGeom.getExtent();
+        final GridExtent targetExtent = new GridExtent(null,
+                new long[] {sourceExtent.getLow(0), sourceExtent.getLow (1)},
+                new long[] {sourceExtent.getLow(0), sourceExtent.getHigh(1)}, true);
+        final GridGeometry targetGeom = new GridGeometry(
+                targetExtent, CELL_CENTER,
+                sourceGeom.getGridToCRS(CELL_CENTER),
+                sourceGeom.getCoordinateReferenceSystem());
+
+        final GridCoverage result = resample(source, targetGeom);
+        final int height = (int) targetExtent.getSize(1);
+        assertPixelsEqual(source.render(null), new Rectangle(0, 0, 1, height), result.render(null), null);
+    }
+
+    /**
+     * Tests resampling in a sub-region specified by an envelope.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    public void testSubGeographicArea() throws TransformException {
+        final GridCoverage2D source = createCoverage2D();             // Envelope2D(20, 15, 60, 62)
         GridGeometry gg = new GridGeometry(null, new Envelope2D(HardCodedCRS.WGS84, 18, 20, 17, 31));
         final GridCoverage target = resample(source, gg);
         final GridExtent sourceExtent = source.getGridGeometry().getExtent();
@@ -169,13 +436,13 @@ public final strictfp class ResampledGridCoverageTest extends TestCase {
     }
 
     /**
-     * Tests application of a reprojection.
+     * Tests application of a non-linear transform.
      *
      * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
      */
     @Test
     public void testReprojection() throws TransformException {
-        final GridCoverage2D source = createGridCoverage();
+        final GridCoverage2D source = createCoverage2D();
         GridGeometry gg = new GridGeometry(null, CELL_CENTER, null, HardCodedConversions.mercator());
         final GridCoverage target = resample(source, gg);
         assertTrue("GridExtent.startsAtZero", target.getGridGeometry().getExtent().startsAtZero());
@@ -187,4 +454,121 @@ public final strictfp class ResampledGridCoverageTest extends TestCase {
         assertEquals(sourceExtent.getSize(0),   targetExtent.getSize(0));
         assertTrue  (sourceExtent.getSize(1) <= targetExtent.getSize(1));
     }
+
+    /**
+     * Tests application of a three-dimensional transform which can not be reduced to a two-dimensional transform.
+     * It happens for example when transformation of <var>x</var> or <var>y</var> coordinate depends on <var>z</var>
+     * coordinate value. In such case we can not separate the 3D transform into (2D + 1D) transforms. This method
+     * verifies that {@link ResampledGridCoverage} nevertheless manages to do its work even in that situation.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    @org.junit.Ignore("Needs more development")
+    public void testNonSeparableGridToCRS() throws TransformException {
+        final GridCoverage source = createCoverageND(false);
+        final MatrixSIS nonSeparableMatrix = Matrices.createDiagonal(4, 4);
+        nonSeparableMatrix.setElement(0, 2, 1);     // Make X dependent of Z.
+        nonSeparableMatrix.setElement(1, 2, 1);     // Make Y dependent of Z.
+        final MathTransform nonSeparableG2C = MathTransforms.concatenate(
+                MathTransforms.linear(nonSeparableMatrix),
+                source.getGridGeometry().getGridToCRS(CELL_CENTER));
+        {
+            /*
+             * The test in this block is not a `ResampleGridCoverage` test, but rather a
+             * check for a condition that we need for the test performed in this method.
+             */
+            final TransformSeparator separator = new TransformSeparator(nonSeparableG2C);
+            separator.addSourceDimensions(0, 1);
+            separator.addTargetDimensions(0, 1);
+            try {
+                final MathTransform separated = separator.separate();
+                fail("Test requires a non-separable transform, but separation succeed: " + separated);
+            } catch (FactoryException e) {
+                // Successful check.
+            }
+        }
+        final GridGeometry targetGeom = new GridGeometry(
+                source.getGridGeometry().getExtent(),
+                CELL_CENTER, nonSeparableG2C,
+                source.getCoordinateReferenceSystem());
+        /*
+         * Real test is below (above code was only initialization).
+         */
+        final GridCoverage result = resample(source, targetGeom);
+        assertPixelsEqual(source.render(null), null, result.render(null), null);
+    }
+
+    /**
+     * Tests the addition of a temporal axis. The value to insert in the temporal coordinate can be computed
+     * from the four-dimensional "grid to CRS" transform given in argument to the {@code resample(…)} method,
+     * combined with the source grid extent.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    public void crs3D_to_crs4D() throws TransformException {
+        final GridCoverage source3D = createCoverageND(false);
+        final GridGeometry target4D = createGridGeometryND(HardCodedCRS.WGS84_4D, 0, 1, 2, 3, false);
+        final GridCoverage result   = resample(source3D, target4D);
+        assertEquals(target4D, result.getGridGeometry());
+        assertPixelsEqual(source3D.render(null), null, result.render(null), null);
+    }
+
+    /**
+     * Tests the removal of temporal axis.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
+     */
+    @Test
+    @org.junit.Ignore("To debug")
+    public void crs4D_to_crs3D() throws TransformException {
+        final GridGeometry target3D = createGridGeometryND(HardCodedCRS.WGS84_3D, 0, 1, 2, 3, false);
+        final GridCoverage source4D = createCoverageND(true);
+        final GridCoverage result   = resample(source4D, target3D);
+        assertEquals(target3D, result.getGridGeometry());
+        assertPixelsEqual(source4D.render(null), null, result.render(null), null);
+    }
+
+    /**
+     * Returns an image with only the queries part of the given image.
+     * This is an helper tools which can be invoked during debugging
+     * session in IDE capable to display images.
+     *
+     * <p><b>Usage:</b> Add a new watch calling this method on wanted image.</p>
+     *
+     * <p><b>Limitations:</b>
+     * <ul>
+     *   <li>If given image color-model is null, this method assumes 3 byte/RGB image.</li>
+     *   <li>Works only with single-tile images.</li>
+     * </ul>
+     *
+     * @param source  the image to display.
+     * @param extent  if non-null, crop rendering to the rectangle defined by given extent,
+     *                assuming extent low coordinate matches source image (0,0) coordinate.
+     * @return the image directly displayable through debugger.
+     */
+    private static BufferedImage debug(final RenderedImage source, final GridExtent extent) {
+        Raster tile = source.getTile(source.getMinTileX(), source.getMinTileY());
+        final int width, height;
+        if (extent == null) {
+            tile   = tile.createTranslatedChild(0, 0);
+            width  = tile.getWidth();
+            height = tile.getHeight();
+        } else {
+            width  = StrictMath.toIntExact(extent.getSize(0));
+            height = StrictMath.toIntExact(extent.getSize(1));
+            tile   = tile.createChild(0, 0, width, height, 0, 0, null);
+        }
+        final BufferedImage view;
+        if (source.getColorModel() == null) {
+            view = new BufferedImage(width, height, BufferedImage.TYPE_3BYTE_BGR);
+            view.getRaster().setRect(tile);
+        } else {
+            final WritableRaster wr = tile.createCompatibleWritableRaster(0, 0, width, height);
+            wr.setRect(tile);
+            view = new BufferedImage(source.getColorModel(), wr, false, null);
+        }
+        return view;
+    }
 }
diff --git a/core/sis-feature/src/test/java/org/apache/sis/image/TiledImageMock.java b/core/sis-feature/src/test/java/org/apache/sis/image/TiledImageMock.java
index 690d2fe..6c13ba8 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/image/TiledImageMock.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/image/TiledImageMock.java
@@ -28,6 +28,7 @@ import java.awt.image.WritableRaster;
 import java.awt.image.WritableRenderedImage;
 import java.util.Random;
 import java.util.concurrent.atomic.AtomicInteger;
+import org.apache.sis.internal.coverage.j2d.ImageUtilities;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.ArraysExt;
 
@@ -366,10 +367,18 @@ public final strictfp class TiledImageMock extends PlanarImage implements Writab
     }
 
     /**
-     * Not needed for the tests.
+     * Sets a rectangle of this image to the contents of given raster.
+     * The raster is assumed to be in the same coordinate space as this image.
+     * Current implementation can set raster covering only only one tile.
      */
     @Override
-    public void setData(Raster r) {
-        throw new UnsupportedOperationException();
+    public void setData(final Raster r) {
+        final int minX = r.getMinX();
+        final int minY = r.getMinY();
+        final int tx = ImageUtilities.pixelToTileX(this, minX);
+        final int ty = ImageUtilities.pixelToTileY(this, minY);
+        assertEquals("Unsupported operation.", tx, ImageUtilities.pixelToTileX(this, minX + r.getWidth()  - 1));
+        assertEquals("Unsupported operation.", ty, ImageUtilities.pixelToTileX(this, minY + r.getHeight() - 1));
+        tile(tx, ty, true).setRect(r);
     }
 }
diff --git a/core/sis-feature/src/test/java/org/apache/sis/test/FeatureAssert.java b/core/sis-feature/src/test/java/org/apache/sis/test/FeatureAssert.java
index c08b455..1f9ba4f 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/test/FeatureAssert.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/test/FeatureAssert.java
@@ -16,7 +16,12 @@
  */
 package org.apache.sis.test;
 
+import java.awt.Point;
+import java.awt.Rectangle;
 import java.awt.image.Raster;
+import java.awt.image.RenderedImage;
+import org.opengis.coverage.grid.SequenceType;
+import org.apache.sis.image.PixelIterator;
 
 import static org.junit.Assert.*;
 
@@ -38,6 +43,42 @@ public strictfp class FeatureAssert extends ReferencingAssert {
     }
 
     /**
+     * Verifies that sample values in the given image are equal to the expected integer values.
+     *
+     * @param  expected       the expected sample values.
+     * @param  boundsExpected bounds of the expected region, or {@code null} for the whole image.
+     * @param  actual         the image to verify.
+     * @param  boundsActual   bounds of the actual region, or {@code null} for the whole image.
+     */
+    public static void assertPixelsEqual(final RenderedImage expected, final Rectangle boundsExpected,
+                                         final RenderedImage actual,   final Rectangle boundsActual)
+    {
+        if (boundsExpected != null && boundsActual != null) {
+            assertEquals("width",  boundsExpected.width,  boundsActual.width);
+            assertEquals("height", boundsExpected.height, boundsActual.height);
+        }
+        final PixelIterator ie = new PixelIterator.Builder().setIteratorOrder(SequenceType.LINEAR).setRegionOfInterest(boundsExpected).create(expected);
+        final PixelIterator ia = new PixelIterator.Builder().setIteratorOrder(SequenceType.LINEAR).setRegionOfInterest(boundsActual).create(actual);
+        double[] ev = null;
+        double[] av = null;
+        while (ie.next()) {
+            assertTrue(ia.next());
+            ev = ie.getPixel(ev);
+            av = ia.getPixel(av);
+            assertEquals(ev.length, av.length);
+            for (int band=0; band<ev.length; band++) {
+                final double e = ev[band];
+                final double a = av[band];
+                if (Double.doubleToRawLongBits(a) != Double.doubleToRawLongBits(e)) {
+                    final Point p = ia.getPosition();
+                    fail(mismatchedSampleValue(p.x, p.y, p.x - actual.getMinX(), p.y - actual.getMinY(), band, e, a));
+                }
+            }
+        }
+        assertFalse(ia.next());
+    }
+
+    /**
      * Verifies that sample values in the given raster are equal to the expected integer values.
      *
      * @param  raster    the raster to verify.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/CoordinateOperations.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/CoordinateOperations.java
index 21397ee..4bfff39 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/CoordinateOperations.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/CoordinateOperations.java
@@ -20,6 +20,7 @@ import java.util.Set;
 import java.util.Map;
 import java.util.HashMap;
 import java.util.Collections;
+import java.util.function.Supplier;
 import javax.measure.UnitConverter;
 import javax.measure.IncommensurableException;
 import org.opengis.referencing.cs.CSFactory;
@@ -50,7 +51,7 @@ import org.apache.sis.util.collection.Containers;
  * until first needed. Contains also utility methods related to coordinate operations.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.7
  * @module
  */
@@ -81,6 +82,13 @@ public final class CoordinateOperations extends SystemListener {
     public static final String OPERATION_TYPE_KEY = "operationType";
 
     /**
+     * Value of {@link org.apache.sis.referencing.operation.CoordinateOperationContext#getConstantCoordinates()}.
+     * This thread-local is used as a workaround for the fact that we do not yet provide a public API for this
+     * functionality. This workaround should be deleted after a public API is defined.
+     */
+    public static final ThreadLocal<Supplier<double[]>> CONSTANT_COORDINATES = new ThreadLocal<>();
+
+    /**
      * Cached values or {@link #wrapAroundChanges wrapAroundChanges(…)}, created when first needed.
      * Indices are bit masks computed by {@link #changes changes(…)}. Since the most common "wrap around" axes
      * are longitude at dimension 0 or 1, and some measurement of time (in climatology) at dimension 2 or 3,
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationContext.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationContext.java
index 49b3855..1d5e76c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationContext.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationContext.java
@@ -17,16 +17,20 @@
 package org.apache.sis.referencing.operation;
 
 import java.io.Serializable;
+import java.util.function.Supplier;
 import java.util.function.Predicate;
 import org.opengis.metadata.extent.Extent;
 import org.opengis.metadata.extent.GeographicBoundingBox;
 import org.opengis.referencing.operation.CoordinateOperation;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.internal.referencing.CoordinateOperations;
 import org.apache.sis.metadata.iso.extent.DefaultExtent;
 import org.apache.sis.metadata.iso.extent.Extents;
 import org.apache.sis.internal.util.CollectionsExt;
 import org.apache.sis.measure.Latitude;
 import org.apache.sis.measure.Longitude;
 import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.util.collection.BackingStoreException;
 
 
 /**
@@ -54,7 +58,7 @@ import org.apache.sis.util.ArgumentChecks;
  * late binding implementations.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.7
  * @module
  *
@@ -219,4 +223,32 @@ public class CoordinateOperationContext implements Serializable {
     Predicate<CoordinateOperation> getOperationFilter() {
         return null;
     }
+
+    /**
+     * Invoked when some coordinates in the target CRS can not be computed from coordinates in the source CRS.
+     * For example if the source CRS has (<var>x</var>, <var>y</var>) axes and the target CRS has (<var>x</var>,
+     * <var>y</var>, <var>t</var>) axes, then this method is invoked for determining which value to assign to the
+     * <var>t</var> coordinate. In some cases the user can tell that the coordinate should be set to a constant value.
+     *
+     * <p>If this method returns {@code null} (which is the default), then the {@link CoordinateOperationFinder} caller
+     * will throw an {@link org.opengis.referencing.operation.OperationNotFoundException}. Otherwise the returned array
+     * should have a length equals to the number of dimensions in the full (usually compound) target CRS.
+     * Only coordinate values in dimensions without source (the <var>t</var> dimension in above example) will be used.
+     * All other coordinate values will be ignored.
+     *
+     * @return coordinate values to take as constants for the specified target component, or {@code null} if none.
+     * @throws TransformException if the coordinates can not be computed. This exception may occur when the constant
+     *         coordinate values are the results of performing a coordinate operation in advance.
+     *
+     * @todo Non-public API for now, pending more feedback from experience. A public method would be non-static.
+     */
+    static double[] getConstantCoordinates() throws TransformException {
+        final Supplier<double[]> f = CoordinateOperations.CONSTANT_COORDINATES.get();
+        if (f != null) try {
+            return f.get();
+        } catch (BackingStoreException e) {
+            throw e.unwrapOrRethrow(TransformException.class);
+        }
+        return null;
+    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
index a0bfc03..30016a5 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
@@ -112,7 +112,7 @@ import static org.apache.sis.util.Utilities.equalsIgnoreMetadata;
  * </ul>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see DefaultCoordinateOperationFactory#createOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, CoordinateOperationContext)
  *
@@ -733,7 +733,7 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
         /*
          * Transform from ellipsoidal height to the height requested by the caller.
          * This operation requires the horizontal components (φ,λ) of source CRS,
-         * unless the user asked for ellipsooidal height (which strictly speaking
+         * unless the user asked for ellipsoidal height (which strictly speaking
          * is not allowed by ISO 19111). Those horizontal components are given by
          * the interpolation CRS.
          *
@@ -897,36 +897,32 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
             final CoordinateReferenceSystem targetCRS, final List<? extends SingleCRS> targetComponents)
             throws FactoryException
     {
-        final SubOperationInfo[] infos = new SubOperationInfo[targetComponents.size()];
-        final boolean[]   sourceIsUsed = new boolean[sourceComponents.size()];
-        final CoordinateReferenceSystem[] stepComponents = new CoordinateReferenceSystem[infos.length];
         /*
-         * Operations found are stored in 'infos', but are not yet wrapped in PassThroughOperations.
-         * We need to know first if some coordinate values need reordering for matching the target CRS
-         * order. We also need to know if any source coordinates should be dropped.
+         * Operations found are stored in the `infos` array, but are not yet wrapped in PassThroughOperations.
+         * We need to know first if some coordinate values need reordering for matching the target CRS order.
+         * We also need to know if any source coordinates should be dropped.
          */
-        for (int i=0; i<infos.length; i++) {
-            if ((infos[i] = SubOperationInfo.create(this, sourceIsUsed, sourceComponents, targetComponents.get(i))) == null) {
-                throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS));
-            }
-            stepComponents[i] = infos[i].operation.getSourceCRS();
+        final SubOperationInfo[] infos;
+        try {
+            infos = SubOperationInfo.createSteps(this, sourceComponents, targetComponents);
+        } catch (TransformException e) {
+            throw new FactoryException(notFoundMessage(sourceCRS, targetCRS), e);
+        }
+        if (infos == null) {
+            throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS));
         }
         /*
          * At this point, a coordinate operation has been found for all components of the target CRS.
          * However the CoordinateOperation.getSourceCRS() values are not necessarily in the same order
-         * than the components of the source CRS given to this method, and some dimensions may be dropped.
+         * than in the `sourceComponents` list given to this method, and some dimensions may be dropped.
          * The matrix computed by sourceToSelected(…) gives us the rearrangement needed for the coordinate
          * operations that we just found.
          */
-        int remainingSourceDimensions = 0;
-        for (final SubOperationInfo component : infos) {
-            remainingSourceDimensions += component.endAtDimension - component.startAtDimension;
-        }
-        final Matrix select = SubOperationInfo.sourceToSelected(
-                sourceCRS.getCoordinateSystem().getDimension(), remainingSourceDimensions, infos);
+        final CoordinateReferenceSystem[] stepComponents = SubOperationInfo.getSourceCRS(infos);
+        final Matrix select = SubOperationInfo.sourceToSelected(sourceCRS.getCoordinateSystem().getDimension(), infos);
         /*
-         * First, we need a CRS matching the above-cited rearrangement. That CRS will be named 'stepSourceCRS'
-         * and its components will be named 'stepComponents'. Then we will execute a loop in which each component
+         * First, we need a CRS matching the above-cited rearrangement. That CRS will be named `stepSourceCRS`
+         * and its components will be named `stepComponents`. Then we will execute a loop in which each component
          * is progressively (one by one) updated from a source component to a target component. A new step CRS is
          * recreated each time, since it will be needed for each PassThroughOperation.
          */
@@ -951,22 +947,22 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
          * SubOperationInfo, because those indices are not relative to the same CompoundCRS.
          */
         int endAtDimension = 0;
-        final int startOfIdentity = SubOperationInfo.startOfIdentity(infos);
+        int remainingSourceDimensions = select.getNumRow() - 1;
+        final int indexOfFinal = SubOperationInfo.indexOfFinal(infos);
         for (int i=0; i<stepComponents.length; i++) {
+            final SubOperationInfo          info   = infos[i];
             final CoordinateReferenceSystem source = stepComponents[i];
-            final CoordinateReferenceSystem target = targetComponents.get(i);
-            CoordinateOperation subOperation = infos[i].operation;
-            final MathTransform subTransform = subOperation.getMathTransform();
+            final CoordinateReferenceSystem target = targetComponents.get(info.targetComponentIndex);
             /*
-             * In order to compute 'stepTargetCRS', replace in-place a single element in 'stepComponents'.
-             * For each step except the last one, 'stepTargetCRS' is a mix of target and source CRS. Only
-             * after the loop finished, 'stepTargetCRS' will become the complete targetCRS definition.
+             * In order to compute `stepTargetCRS`, replace in-place a single element in `stepComponents`.
+             * For each step except the last one, `stepTargetCRS` is a mix of target CRS and source CRS.
+             * Only after the loop finished, `stepTargetCRS` will become the complete target definition.
              */
             final CoordinateReferenceSystem stepTargetCRS;
-            stepComponents[i] = target;
-            if (i >= startOfIdentity) {
+            stepComponents[info.targetComponentIndex] = target;
+            if (i >= indexOfFinal) {
                 stepTargetCRS = targetCRS;              // If all remaining transforms are identity, we reached the final CRS.
-            } else if (subTransform.isIdentity()) {
+            } else if (info.isIdentity()) {
                 stepTargetCRS = stepSourceCRS;          // In any identity transform, the source and target CRS are equal.
             } else if (stepComponents.length == 1) {
                 stepTargetCRS = target;                 // Slight optimization of the next block.
@@ -977,31 +973,14 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
             int delta = source.getCoordinateSystem().getDimension();
             final int startAtDimension = endAtDimension;
             endAtDimension += delta;
-            /*
-             * Constructs the pass through transform only if there is at least one coordinate to pass.
-             * Actually the code below would work inconditionally, but we perform this check anyway
-             * for avoiding the creation of intermediate objects.
-             */
-            if (!(startAtDimension == 0 && endAtDimension == remainingSourceDimensions)) {
-                final Map<String,?> properties = IdentifiedObjects.getProperties(subOperation);
-                /*
-                 * The DefaultPassThroughOperation constuctor expect a SingleOperation.
-                 * In most case, the 'subOperation' is already of this kind. However if
-                 * it is not, try to copy it in such object.
-                 */
-                final SingleOperation op;
-                if (subOperation instanceof SingleOperation) {
-                    op = (SingleOperation) subOperation;
-                } else {
-                    op = factorySIS.createSingleOperation(properties,
-                            subOperation.getSourceCRS(), subOperation.getTargetCRS(), null,
-                            new DefaultOperationMethod(subTransform), subTransform);
-                }
-                subOperation = new DefaultPassThroughOperation(properties, stepSourceCRS, stepTargetCRS,
-                        op, startAtDimension, remainingSourceDimensions - endAtDimension);
+            final int numTrailingCoordinates = remainingSourceDimensions - endAtDimension;
+            CoordinateOperation subOperation = info.operation;
+            if ((startAtDimension | numTrailingCoordinates) != 0) {
+                subOperation = new DefaultPassThroughOperation(IdentifiedObjects.getProperties(subOperation),
+                        stepSourceCRS, stepTargetCRS, subOperation, startAtDimension, numTrailingCoordinates);
             }
             /*
-             * Concatenate the operation with the ones we have found so far, and use the current 'stepTargetCRS'
+             * Concatenate the operation with the ones we have found so far, and use the current `stepTargetCRS`
              * as the source CRS for the next operation step. We also need to adjust the dimension indices,
              * since the previous operations may have removed some dimensions. Note that the delta may also
              * be negative in a few occasions.
@@ -1012,6 +991,17 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
             endAtDimension -= delta;
             remainingSourceDimensions -= delta;
         }
+        /*
+         * If there is some target dimensions that are set to constant values instead than computed from the
+         * source dimensions, add those constants as the last step. It never occur except is some particular
+         * contexts like when computing a transform between two `GridGeometry` instances.
+         */
+        if (stepComponents.length < infos.length) {
+            final Matrix m = SubOperationInfo.createConstantOperation(infos, stepComponents.length,
+                    stepSourceCRS.getCoordinateSystem().getDimension(),
+                        targetCRS.getCoordinateSystem().getDimension());
+            operation = concatenate(operation, createFromAffineTransform(CONSTANTS, stepSourceCRS, targetCRS, m));
+        }
         return asList(operation);
     }
 
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java
index 291af97..7827710 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationRegistry.java
@@ -113,6 +113,11 @@ class CoordinateOperationRegistry {
     static final Identifier IDENTITY = createIdentifier(Vocabulary.Keys.Identity);
 
     /**
+     * The identifier for an operation setting some coordinates to constant values.
+     */
+    static final Identifier CONSTANTS = createIdentifier(Vocabulary.Keys.Constants);
+
+    /**
      * The identifier for conversion using an affine transform for axis swapping and/or unit conversions.
      */
     static final Identifier AXIS_CHANGES = createIdentifier(Vocabulary.Keys.AxisChanges);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/DefaultPassThroughOperation.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/DefaultPassThroughOperation.java
index 15215c1..717129a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/DefaultPassThroughOperation.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/DefaultPassThroughOperation.java
@@ -71,7 +71,7 @@ public class DefaultPassThroughOperation extends AbstractCoordinateOperation imp
     private CoordinateOperation operation;
 
     /**
-     * Constructs a single operation from a set of properties.
+     * Constructs a pass-through operation from a set of properties.
      * The properties given in argument follow the same rules than for the
      * {@linkplain AbstractCoordinateOperation#AbstractCoordinateOperation(Map, CoordinateReferenceSystem,
      * CoordinateReferenceSystem, CoordinateReferenceSystem, MathTransform) super-class constructor}.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/SubOperationInfo.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/SubOperationInfo.java
index f797bae..0b80fd8 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/SubOperationInfo.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/SubOperationInfo.java
@@ -19,21 +19,27 @@ package org.apache.sis.referencing.operation;
 import java.util.List;
 import org.opengis.referencing.crs.*;
 import org.opengis.referencing.operation.CoordinateOperation;
-import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.OperationNotFoundException;
+import org.opengis.referencing.operation.TransformException;
 import org.opengis.util.FactoryException;
 import org.apache.sis.util.logging.Logging;
 import org.apache.sis.internal.system.Loggers;
 import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 
 
 /**
- * Information about the relationship between a source component and a target component
- * in {@code CompoundCRS} instances.
+ * Information about the operation from a source component to a target component in {@code CompoundCRS} instances.
+ * An instance of {@code SubOperationInfo} is created for each target CRS component. This class allows to collect
+ * information about all operation steps before to start the creation of pass-through operations. This separation
+ * is useful for applying reordering.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.7
- * @since   0.7
+ * @version 1.1
+ *
+ * @see CoordinateOperationFinder#createOperationStep(CoordinateReferenceSystem, List, CoordinateReferenceSystem, List)
+ *
+ * @since 0.7
  * @module
  */
 final class SubOperationInfo {
@@ -71,105 +77,219 @@ final class SubOperationInfo {
     }
 
     /**
-     * The coordinate operation between a source component and a target component.
+     * The coordinate operation between a source CRS component and a target CRS component.
+     * Exactly one of {@link #operation} or {@link #constants} shall be non-null.
      */
     final CoordinateOperation operation;
 
     /**
-     * Returns the first dimension (inclusive) where the source component CRS begins in the source compound CRS.
+     * The constant values to store in target coordinates, or {@code null} if none. This array is usually null.
+     * It may be non-null if no source CRS component has been found for a target CRS component. For example if
+     * the source CRS has (<var>x</var>, <var>y</var>) axes and the target CRS has (<var>x</var>, <var>y</var>,
+     * <var>t</var>) axes, then this array may be set to a non-null value for specifying the <var>t</var> value.
+     * The array length is the number of dimensions in the full (usually compound) target CRS, but only the
+     * coordinate values for sourceless dimensions are used. Other coordinates are ignored and can be NaN.
+     * Exactly one of {@link #operation} or {@link #constants} shall be non-null.
+     *
+     * @see CoordinateOperationContext#getConstantCoordinates()
+     */
+    private final double[] constants;
+
+    /**
+     * The first dimension (inclusive) and the last dimension (exclusive) where the {@link SingleCRS} starts/ends
+     * in the full (usually compound) CRS. It may be an index in source dimensions or in target dimensions,
+     * depending on the following rules:
+     *
+     * <ul>
+     *   <li>If {@link #operation} is non-null, then this is an index in the <em>source</em> dimensions.
+     *       This is the normal case.</li>
+     *   <li>Otherwise the operation is sourceless with coordinate results described by the {@link #constants} array.
+     *       Since there is no source, the index become an index in target dimensions instead of source dimensions.</li>
+     * </ul>
+     *
+     * @see #sourceToSelected(int, SubOperationInfo[])
      */
-    final int startAtDimension;
+    private final int startAtDimension, endAtDimension;
 
     /**
-     * Returns the last dimension (exclusive) where the source component CRS ends in the source compound CRS.
+     * Index of this instance in the array of {@code SubOperationInfo} instances,
+     * before the reordering applied by {@link #getSourceCRS(SubOperationInfo[])}.
      */
-    final int endAtDimension;
+    final int targetComponentIndex;
 
     /**
-     * Creates a new instance containing the given information.
+     * Creates a new instance wrapping the given coordinate operation or coordinate constants.
+     * Exactly one of {@code operation} or {@code constants} shall be non-null.
      */
-    private SubOperationInfo(final CoordinateOperation operation, final int startAtDimension, final int endAtDimension) {
-        this.operation        = operation;
-        this.startAtDimension = startAtDimension;
-        this.endAtDimension   = endAtDimension;
+    private SubOperationInfo(final int targetComponentIndex, final CoordinateOperation operation,
+                final double[] constants, final int startAtDimension, final int endAtDimension)
+    {
+        this.operation            = operation;
+        this.constants            = constants;
+        this.startAtDimension     = startAtDimension;
+        this.endAtDimension       = endAtDimension;
+        this.targetComponentIndex = targetComponentIndex;
+        assert (operation == null) != (constants == null);
     }
 
     /**
-     * Searches in given list of source components for an operation capable to convert or transform coordinates
-     * to the given target CRS. If no such operation can be found, then this method returns {@code null}.
+     * Searches in given list of source components for operations capable to transform coordinates to each target CRS.
+     * There is one {@code SubOperationInfo} per target CRS because we need to satisfy all target dimensions, while it
+     * is okay to ignore some source dimensions. If an operation can not be found, then this method returns {@code null}.
      *
-     * @param  caller        the object which is inferring a coordinate operation.
-     * @param  sourceIsUsed  flags for keeping trace of which source has been used.
-     * @param  sources       all components of the source CRS.
-     * @param  target        one component of the target CRS.
-     * @return information about a coordinate operation from a source CRS to the given target CRS, or {@code null}.
+     * @param  caller   the object which is inferring a coordinate operation.
+     * @param  sources  all components of the source CRS.
+     * @param  targets  all components of the target CRS.
+     * @return information about each coordinate operation from a source CRS to a target CRS, or {@code null}.
      * @throws FactoryException if an error occurred while grabbing a coordinate operation.
+     * @throws TransformException if an error occurred while computing the sourceless coordinate constants.
      */
-    static SubOperationInfo create(final CoordinateOperationFinder caller, final boolean[] sourceIsUsed,
-            final List<? extends SingleCRS> sources, final SingleCRS target) throws FactoryException
+    static SubOperationInfo[] createSteps(final CoordinateOperationFinder caller,
+                                          final List<? extends SingleCRS> sources,
+                                          final List<? extends SingleCRS> targets)
+            throws FactoryException, TransformException
     {
-        OperationNotFoundException failure = null;
-        final Class<?> targetType = type(target);
-        for (final Class<?>[] sourceTypes : COMPATIBLE_TYPES) {
-            if (sourceTypes[0].isAssignableFrom(targetType)) {
-                for (final Class<?> sourceType : sourceTypes) {
-                    int startAtDimension;
-                    int endAtDimension = 0;
-                    for (int i=0; i<sourceIsUsed.length; i++) {
-                        final SingleCRS source = sources.get(i);
-                        startAtDimension = endAtDimension;
-                        endAtDimension += source.getCoordinateSystem().getDimension();
-                        if (!sourceIsUsed[i] && sourceType.isAssignableFrom(type(source))) {
-                            final CoordinateOperation operation;
-                            try {
-                                operation = caller.createOperation(source, target);
-                            } catch (OperationNotFoundException exception) {
-                                if (failure == null) {
-                                    failure = exception;
-                                } else {
-                                    failure.addSuppressed(exception);
+        final SubOperationInfo[] infos = new SubOperationInfo[targets.size()];
+        final boolean[] sourceComponentIsUsed  =  new boolean[sources.size()];
+        /*
+         * Iterate over target CRS because all of them must have an operation.
+         * One the other hand, source CRS can be used zero or one (not two) time.
+         *
+         * Note: the loops on source CRS and on target CRS follow the same pattern.
+         *       The first 6 lines of code are the same with only `target` replaced
+         *       by `source`.
+         */
+        int targetStartAtDimension;
+        int targetEndAtDimension = 0;
+next:   for (int targetComponentIndex = 0; targetComponentIndex < infos.length; targetComponentIndex++) {
+            final SingleCRS target = targets.get(targetComponentIndex);
+            targetStartAtDimension = targetEndAtDimension;
+            targetEndAtDimension  += target.getCoordinateSystem().getDimension();
+
+            final Class<?> targetType = type(target);
+            OperationNotFoundException failure = null;
+            /*
+             * For each target CRS, search for a source CRS which has not yet been used.
+             * Some sources may be left unused after this method completion. Check only
+             * sources that may be compatible according the `COMPATIBLE_TYPES` array.
+             */
+            int sourceStartAtDimension;
+            int sourceEndAtDimension = 0;
+            for (int sourceComponentIndex = 0; sourceComponentIndex < sourceComponentIsUsed.length; sourceComponentIndex++) {
+                final SingleCRS source = sources.get(sourceComponentIndex);
+                sourceStartAtDimension = sourceEndAtDimension;
+                sourceEndAtDimension  += source.getCoordinateSystem().getDimension();
+                if (!sourceComponentIsUsed[sourceComponentIndex]) {
+                    for (final Class<?>[] sourceTypes : COMPATIBLE_TYPES) {
+                        if (sourceTypes[0].isAssignableFrom(targetType)) {
+                            for (final Class<?> sourceType : sourceTypes) {
+                                if (sourceType.isAssignableFrom(type(source))) {
+                                    final CoordinateOperation operation;
+                                    try {
+                                        operation = caller.createOperation(source, target);
+                                    } catch (OperationNotFoundException exception) {
+                                        if (failure == null) {
+                                            failure = exception;
+                                        } else {
+                                            failure.addSuppressed(exception);
+                                        }
+                                        continue;
+                                    }
+                                    /*
+                                     * Found an operation.  Exclude the source component from the list because each source
+                                     * should be used at most once by SourceComponent. Note that the same source may still
+                                     * be used again in another context if that source is also an interpolation CRS.
+                                     *
+                                     * EXAMPLE: consider a coordinate operation from (GeodeticCRS₁, VerticalCRS₁) source
+                                     * to (GeodeticCRS₂, VerticalCRS₂) target.  The source GeodeticCRS₁ should be mapped
+                                     * to exactly one target component (which is GeodeticCRS₂)  and  VerticalCRS₁ mapped
+                                     * to VerticalCRS₂.  But the operation on vertical coordinates may need GeodeticCRS₁
+                                     * for doing its work, so GeodeticCRS₁ is needed twice.  However when needed for the
+                                     * vertical coordinate operation, the GeodeticCRS₁ is used as an "interpolation CRS".
+                                     * Interpolation CRS are handled in other code paths; it is not the business of this
+                                     * SourceComponent class to care about them. From the point of view of this class,
+                                     * GeodeticCRS₁ is used only once.
+                                     */
+                                    sourceComponentIsUsed[sourceComponentIndex] = true;
+                                    infos[targetComponentIndex] = new SubOperationInfo(targetComponentIndex,
+                                            operation, null, sourceStartAtDimension, sourceEndAtDimension);
+
+                                    if (failure != null) {
+                                        Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION),
+                                                CoordinateOperationFinder.class, "decompose", failure);
+                                    }
+                                    continue next;
                                 }
-                                continue;
-                            }
-                            /*
-                             * Found an operation.  Exclude the source component from the list because each source
-                             * should be used at most once by SourceComponent. Note that the same source may still
-                             * be used again in another context if that source is also an interpolation CRS.
-                             *
-                             * EXAMPLE: consider a coordinate operation from (GeodeticCRS₁, VerticalCRS₁) source
-                             * to (GeodeticCRS₂, VerticalCRS₂) target.  The source GeodeticCRS₁ should be mapped
-                             * to exactly one target component (which is GeodeticCRS₂)  and  VerticalCRS₁ mapped
-                             * to VerticalCRS₂.  But the operation on vertical coordinates may need GeodeticCRS₁
-                             * for doing its work, so GeodeticCRS₁ is needed twice.  However when needed for the
-                             * vertical coordinate operation, the GeodeticCRS₁ is used as an "interpolation CRS".
-                             * Interpolation CRS are handled in other code paths; it is not the business of this
-                             * SourceComponent class to care about them. From the point of view of this class,
-                             * GeodeticCRS₁ is used only once.
-                             */
-                            sourceIsUsed[i] = true;
-                            if (failure != null) {
-                                Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION),
-                                        CoordinateOperationFinder.class, "decompose", failure);
                             }
-                            return new SubOperationInfo(operation, startAtDimension, endAtDimension);
                         }
                     }
                 }
             }
+            if (failure != null) {
+                throw failure;
+            }
+            /*
+             * If we reach this point, we have not been able to find a source CRS that we can map to the target CRS.
+             * Usually this is fatal; returning null will instruct the caller to throw `OperationNotFoundException`.
+             * However in some contexts (e.g. when searching for an operation between two `GridGeometry` instances)
+             * it is possible to assign a constant value to the target coordinates. Those values can not be guessed
+             * by `sis-referencing`; they must be provided by caller. If such constants are specified, then we will
+             * try to apply them.
+             */
+            final double[] constants = CoordinateOperationContext.getConstantCoordinates();
+            if (constants == null || constants.length < targetEndAtDimension) {
+                return null;
+            }
+            for (int i = targetStartAtDimension; i < targetEndAtDimension; i++) {
+                if (Double.isNaN(constants[i])) return null;
+            }
+            infos[targetComponentIndex] = new SubOperationInfo(targetComponentIndex,
+                    null, constants, targetStartAtDimension, targetEndAtDimension);
+        }
+        return infos;
+    }
+
+    /**
+     * Returns the source CRS of given operations. This method modifies the given array in-place by moving all
+     * sourceless operations last. Then an array is returned with the source CRS of only ordinary operations.
+     * Each CRS at index <var>i</var> in the returned array is the component from {@link #startAtDimension}
+     * inclusive to {@link #endAtDimension} exclusive in the complete (usually compound) source CRS analyzed
+     * by {@link CoordinateOperationFinder}.
+     *
+     * @param  selected  all operations from source to target {@link CompoundCRS}.
+     * @return source CRS of all ordinary operations (excluding operations producing constant values).
+     */
+    static CoordinateReferenceSystem[] getSourceCRS(final SubOperationInfo[] selected) {
+        int n = selected.length;
+        final int last = n - 1;
+        for (int i=0; i<n; i++) {
+            final SubOperationInfo component = selected[i];
+            if (component.operation == null) {
+                System.arraycopy(selected, i+1, selected, i, last - i);
+                selected[last] = component;
+                n--;
+            }
         }
-        if (failure != null) {
-            throw failure;
+        final CoordinateReferenceSystem[] stepComponents = new CoordinateReferenceSystem[n];
+        for (int i=0; i<n; i++) {
+            stepComponents[i] = selected[i].operation.getSourceCRS();
         }
-        return null;
+        return stepComponents;
     }
 
     /**
-     * Returns the dimension from which all remaining operations are identity.
+     * Returns the index of the last non-identity operation. This is used as a slight optimization for deciding when
+     * {@link CoordinateOperationFinder} can stop to create intermediate target {@link CompoundCRS} instances because
+     * all remaining operations leave target coordinates unchanged. It may help to skip a few operations for example
+     * when converting (<var>x</var>, <var>y</var>, <var>t</var>) coordinates where <var>t</var> value is unchanged.
+     *
+     * @param  selected  all operations from source to target {@link CompoundCRS}.
+     * @return index of the last non-identity operation, inclusive.
      */
-    static int startOfIdentity(final SubOperationInfo[] selected) {
+    static int indexOfFinal(final SubOperationInfo[] selected) {
         int n = selected.length;
         while (n != 0) {
-            if (!selected[--n].operation.getMathTransform().isIdentity()) {
+            if (!selected[--n].isIdentity()) {
                 break;
             }
         }
@@ -177,15 +297,49 @@ final class SubOperationInfo {
     }
 
     /**
-     * Returns a matrix for an affine transform from all source coordinates to the coordinates of the
-     * source components selected for participating in the coordinate operation.
+     * Returns a matrix for an affine transform moving coordinate values from their position in the source CRS to a
+     * position in the order {@link #operation}s are applied. This matrix is needed because {@link #operation} may
+     * select any source CRS in the list of {@link SingleCRS} given to the {@link #createSteps createSteps(…)} method;
+     * the source CRS are not necessarily picked in the same order as they appear in the list.
+     *
+     * <div class="note"><b>Example:</b>
+     * if the source CRS has (<var>x</var>, <var>y</var>, <var>t</var>) coordinates and the target CRS has
+     * (<var>t</var>, <var>x</var>, <var>y</var>) coordinates with some operation applied on <var>x</var>
+     * and <var>y</var>, then the operations will be applied in that order:
+     *
+     * <ol>
+     *   <li>An operation for <var>t</var>, because it is the first coordinate to appear in target CRS.</li>
+     *   <li>An operation for (<var>x</var>, <var>y</var>), because those coordinates are next in target CRS.</li>
+     * </ol>
+     *
+     * Since {@link DefaultPassThroughOperation} can not take coordinates before the "first affected coordinate"
+     * dimension and move them into the "trailing coordinates" dimension, we have to reorder coordinates before
+     * to create the pass-through operations. This is done by the following matrix:
+     *
+     * {@preformat math
+     *   ┌   ┐   ┌         ┐┌   ┐
+     *   │ t │   │ 0 0 1 0 ││ x │
+     *   │ x │ = │ 1 0 0 0 ││ y │
+     *   │ y │   │ 0 1 0 0 ││ t │
+     *   │ 1 │   │ 0 0 0 1 ││ 1 │
+     *   └   ┘   └         ┘└   ┘
+     * }
+     * </div>
      *
-     * @param sourceDimensions    number of dimension of the source {@code CompoundCRS}.
-     * @param selectedDimensions  number of source dimensions needed by the coordinate operations.
-     * @param selected            all {@code SourceComponent} instances needed for the target {@code CompoundCRS}.
+     * Furthermore some dimensions may be dropped,
+     * e.g. from (<var>x</var>, <var>y</var>, <var>t</var>) to (<var>x</var>, <var>y</var>).
+     *
+     * @param  sourceDimensions  number of dimensions in the source {@link CompoundCRS}.
+     * @param  selected          all operations from source to target {@link CompoundCRS}.
+     * @return mapping from source {@link CompoundCRS} to each {@link CoordinateOperation#getSourceCRS()}.
      */
-    static Matrix sourceToSelected(final int sourceDimensions, final int selectedDimensions, final SubOperationInfo[] selected) {
-        final Matrix select = Matrices.createZero(selectedDimensions + 1, sourceDimensions + 1);
+    static MatrixSIS sourceToSelected(final int sourceDimensions, final SubOperationInfo[] selected) {
+        int selectedDimensions = 0;
+        for (final SubOperationInfo component : selected) {
+            if (component.operation == null) break;
+            selectedDimensions += component.endAtDimension - component.startAtDimension;
+        }
+        final MatrixSIS select = Matrices.createZero(selectedDimensions + 1, sourceDimensions + 1);
         select.setElement(selectedDimensions, sourceDimensions, 1);
         int j = 0;
         for (final SubOperationInfo component : selected) {
@@ -195,4 +349,40 @@ final class SubOperationInfo {
         }
         return select;
     }
+
+    /**
+     * Returns {@code true} if the coordinate operation wrapped by this object is an identity transform.
+     */
+    final boolean isIdentity() {
+        return (operation != null) && operation.getMathTransform().isIdentity();
+    }
+
+    /**
+     * Returns the matrix of an operation setting some coordinates to a constant values.
+     *
+     * @param  selected  all operations from source to target {@link CompoundCRS}.
+     * @param  n         index of the first selected operation which describe a constant value.
+     * @param  srcDim    number of dimensions in the target CRS of previous operation step.
+     * @param  tgtDim    number of dimensions in the full (usually compound) target CRS.
+     */
+    static MatrixSIS createConstantOperation(final SubOperationInfo[] selected, int n,
+                                             final int srcDim, final int tgtDim)
+    {
+        final boolean[] targetDimensionIsUsed = new boolean[tgtDim];
+        final MatrixSIS m = Matrices.createZero(tgtDim + 1, srcDim + 1);
+        m.setElement(tgtDim, srcDim, 1);
+        do {
+            final SubOperationInfo component = selected[n];
+            for (int j = component.startAtDimension; j < component.endAtDimension; j++) {
+                m.setElement(j, srcDim, component.constants[j]);
+                targetDimensionIsUsed[j] = true;
+            }
+        } while (++n < selected.length);
+        for (int i=0,j=0; j<tgtDim; j++) {
+            if (!targetDimensionIsUsed[j]) {
+                m.setElement(j, i++, 1);
+            }
+        }
+        return m;
+    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/HardCodedCRS.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/HardCodedCRS.java
index 930bc9d..65bc107 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/HardCodedCRS.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/crs/HardCodedCRS.java
@@ -78,6 +78,16 @@ public final strictfp class HardCodedCRS {
             properties("WGS 84 (3D)", null), HardCodedDatum.WGS84, HardCodedCS.GEODETIC_3D);
 
     /**
+     * A four-dimensional geographic coordinate reference system using the WGS84 datum.
+     * This CRS uses (<var>longitude</var>, <var>latitude</var>, <var>height</var>, <var>time</var>)
+     * with the 3 first dimensions specified by {@link #WGS84_3D} and the fourth dimension specified
+     * by {@link #TIME}.
+     *
+     * @see #TIME_WGS84
+     */
+    public static final DefaultCompoundCRS WGS84_4D;
+
+    /**
      * A two-dimensional geographic coordinate reference system using the Paris prime meridian.
      * This CRS uses (<var>longitude</var>, <var>latitude</var>) coordinates with longitude values
      * increasing towards the East and latitude values increasing towards the North.
@@ -234,6 +244,19 @@ public final strictfp class HardCodedCRS {
             getProperties(HardCodedCS.DAYS), HardCodedDatum.MODIFIED_JULIAN, HardCodedCS.DAYS);
 
     /**
+     * A four-dimensional geographic coordinate reference system with time as the first axis.
+     * This CRS uses (<var>time</var>, <var>longitude</var>, <var>latitude</var>, <var>height</var>)
+     * with the first dimension specified by {@link #TIME} and the 3 last dimensions specified by {@link #WGS84_3D}.
+     * Such axis order is unusual but we use it as a way to verify that SIS is robust to arbitrary axis order.
+     */
+    public static final DefaultCompoundCRS TIME_WGS84 = new DefaultCompoundCRS(
+            properties("time + WGS 84 (3D)", null), TIME, WGS84_3D);;
+
+    static {
+        WGS84_4D = new DefaultCompoundCRS(properties("WGS 84 (3D) + time", null), WGS84_3D, TIME);
+    }
+
+    /**
      * A (λ,φ,H) CRS where <var>H</var> is the {@link #GRAVITY_RELATED_HEIGHT}.
      * This constant uses the "geoid" term as an approximation for the gravity related height.
      */
@@ -250,6 +273,7 @@ public final strictfp class HardCodedCRS {
     /**
      * A (H,λ,φ) CRS where <var>H</var> is the {@link #GRAVITY_RELATED_HEIGHT}.
      * This is the same CRS than {@link #GEOID_3D} but with height first.
+     * Such axis order is unusual but we use it as a way to verify that SIS is robust to arbitrary axis order.
      */
     public static final DefaultCompoundCRS GEOID_ZXY = new DefaultCompoundCRS(
             properties("height + WGS 84", null), GRAVITY_RELATED_HEIGHT, WGS84);
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
index 3e97eb0..e7e080e 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.java
@@ -225,6 +225,11 @@ public final class Vocabulary extends IndexedResourceBundle {
         public static final short ConstantPressureSurface = 32;
 
         /**
+         * Constants
+         */
+        public static final short Constants = 233;
+
+        /**
          * Container
          */
         public static final short Container = 33;
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
index 6b157c2..f0224e2 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary.properties
@@ -46,6 +46,7 @@ Code_1                  = {0} code
 Colors                  = Colors
 ColorIndex              = Color index
 Commands                = Commands
+Constants               = Constants
 ConstantPressureSurface = Constant pressure surface
 Container               = Container
 Conversion              = Conversion
diff --git a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
index 0c89b49..84dd01e 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
+++ b/core/sis-utility/src/main/java/org/apache/sis/util/resources/Vocabulary_fr.properties
@@ -53,6 +53,7 @@ Code_1                  = Code {0}
 Colors                  = Couleurs
 ColorIndex              = Indice de couleur
 Commands                = Commandes
+Constants               = Constantes
 ConstantPressureSurface = Surface \u00e0 pression constante
 Container               = Conteneur
 Conversion              = Conversion


Mime
View raw message