sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch master updated: Cherry-pick the https://issues.apache.org/jira/browse/SIS-409 fix from geoapi-4.0 branch.
Date Mon, 17 Feb 2020 10:32:33 GMT
This is an automated email from the ASF dual-hosted git repository.

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


The following commit(s) were added to refs/heads/master by this push:
     new 34fb39f  Cherry-pick the https://issues.apache.org/jira/browse/SIS-409 fix from geoapi-4.0 branch.
34fb39f is described below

commit 34fb39faacf137eb511c0aa0bd027bad21714797
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Feb 17 11:14:12 2020 +0100

    Cherry-pick the https://issues.apache.org/jira/browse/SIS-409 fix from geoapi-4.0 branch.
---
 .../apache/sis/internal/referencing/RTreeNode.java | 356 +++++++++++
 .../apache/sis/internal/referencing/Resources.java |  10 +
 .../sis/internal/referencing/Resources.properties  |   2 +
 .../internal/referencing/Resources_fr.properties   |   2 +
 .../referencing/j2d/IntervalRectangle.java         |  59 +-
 .../apache/sis/internal/referencing/j2d/Tile.java  | 678 +++++++++++++++++++++
 .../internal/referencing/j2d/TileOrganizer.java    | 515 ++++++++++++++++
 .../internal/referencing/j2d/TileTranslation.java  |  75 +++
 .../sis/internal/referencing/j2d/package-info.java |   2 +-
 .../referencing/provider/AbstractProvider.java     |  16 +-
 .../provider/DatumShiftGridCompressed.java         |  13 +-
 .../referencing/provider/DatumShiftGridFile.java   | 324 +++++++++-
 .../referencing/provider/DatumShiftGridGroup.java  | 306 ++++++++++
 .../referencing/provider/DatumShiftGridLoader.java |  14 +-
 .../provider/FranceGeocentricInterpolation.java    |   7 +-
 .../referencing/provider/GeocentricAffine.java     |  11 +-
 .../internal/referencing/provider/Molodensky.java  |   5 +-
 .../sis/internal/referencing/provider/NADCON.java  |   5 +-
 .../sis/internal/referencing/provider/NTv1.java    |  90 +++
 .../sis/internal/referencing/provider/NTv2.java    | 557 ++++++++++++-----
 .../sis/referencing/datum/DatumShiftGrid.java      |  82 ++-
 .../operation/builder/LinearTransformBuilder.java  |   9 +-
 .../operation/matrix/AffineTransforms2D.java       |  14 +-
 .../operation/transform/MathTransforms.java        |   7 +-
 .../transform/SpecializableTransform.java          | 303 ++++-----
 ...g.opengis.referencing.operation.OperationMethod |   1 +
 .../FranceGeocentricInterpolationTest.java         |   8 +-
 .../internal/referencing/provider/NADCONTest.java  |  20 +-
 .../internal/referencing/provider/NTv2Test.java    | 120 +++-
 .../referencing/provider/ProvidersTest.java        |   1 +
 .../java/org/apache/sis/internal/util/Strings.java |  32 +-
 .../org/apache/sis/util/resources/Vocabulary.java  |   5 +
 .../sis/util/resources/Vocabulary.properties       |   1 +
 .../sis/util/resources/Vocabulary_fr.properties    |   1 +
 .../src/test/java/org/apache/sis/test/Assume.java  |   2 +-
 .../internal/netcdf/SatelliteGroundTrackTest.java  |   2 +-
 36 files changed, 3179 insertions(+), 476 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/RTreeNode.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/RTreeNode.java
new file mode 100644
index 0000000..dcb6e5f
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/RTreeNode.java
@@ -0,0 +1,356 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.function.Consumer;
+import org.opengis.geometry.DirectPosition;
+import org.opengis.geometry.Envelope;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.io.wkt.Formatter;
+import org.apache.sis.util.Utilities;
+import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.collection.TreeTable;
+import org.apache.sis.util.collection.TableColumn;
+import org.apache.sis.util.collection.DefaultTreeTable;
+
+
+/**
+ * Placeholder for a RTree. This simple implementation is not designed for scalability or performance.
+ * This class is there for minimal service, to be replaced in some future Apache SIS version by a more
+ * sophisticated tree. Current version of this class is okay for small trees where big nodes are added
+ * before small nodes.
+ *
+ * <p>A {@code RTreeNode} instance is used as the root of the tree. Children nodes are stored as a linked list
+ * (the list is implemented by the {@link #sibling} field, which reference the next element in the list).</p>
+ *
+ * <div class="note"><b>Possible evolution:</b>
+ * a future version could avoid extending {@link GeneralEnvelope}. Instead we could provide abstract
+ * {@code contains(…)} methods and let subclasses define them, with possibly more efficient implementations.
+ * We would still need an implementation that delegate to {@link GeneralEnvelope} since that class has the
+ * advantage of handling envelopes crossing the anti-meridian.</div>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public class RTreeNode extends GeneralEnvelope {
+    /**
+     * For cross-version compatibility.
+     */
+    private static final long serialVersionUID = -6544217991652682694L;
+
+    /**
+     * The parent, or {@code null} if none.
+     *
+     * @see #getParent()
+     */
+    private RTreeNode parent;
+
+    /**
+     * The first node fully included in this code, or {@code null} if none. If non-null, that node and all
+     * its {@linkplain #sibling}s shall be fully included in this {@code RTreeNode}. The child may contain
+     * other children, thus forming a tree from this wider node to smaller nodes.
+     *
+     * @see #getChildren()
+     */
+    private RTreeNode firstChild;
+
+    /**
+     * The next node having the same parent than this node. This is used for creating a linked list of nodes
+     * that are the children of the {@linkplain #parent}.
+     *
+     * <div class="note"><b>Design note:</b>
+     * an {@code RTreeNode children} array instead of {@link #firstChild} + {@link #sibling} would be more intuitive.
+     * But the use of linked list avoid one level of indirection and is one less object to create for each node.
+     * The gain may be negligible with a few hundreds nodes, but a future version of this class may target much
+     * more numerous nodes.</div>
+     */
+    private RTreeNode sibling;
+
+    /**
+     * Creates a new node for the given envelope.
+     *
+     * @param  area  the region covered by this node.
+     */
+    public RTreeNode(final Envelope area) {
+        super(area);
+    }
+
+    /**
+     * Returns the parent of this node, or {@code null} if none.
+     *
+     * @return the parent of this node, or {@code null} if none.
+     */
+    public final RTreeNode getParent() {
+        return parent;
+    }
+
+    /**
+     * Returns the immediate children of this node, or an empty list if none.
+     *
+     * @return the immediate children of this node.
+     */
+    public final List<RTreeNode> getChildren() {
+        final List<RTreeNode> children = new ArrayList<>();
+        for (RTreeNode child = firstChild; child != null; child = child.sibling) {
+            assert child.parent == this : child;
+            children.add(child);
+        }
+        return children;
+    }
+
+    /**
+     * Executes the given action on the given node, all its children and all its siblings.
+     *
+     * @param node    the node on which to execute the action, or {@code null} if none.
+     * @param action  the action to execute.
+     */
+    public static void walk(RTreeNode node, final Consumer<? super RTreeNode> action) {
+        while (node != null) {
+            action.accept(node);
+            walk(node.firstChild, action);
+            node = node.sibling;
+        }
+    }
+
+    /**
+     * Adds the given node to the tree having this node as the root. This method assumes that this
+     * node or its siblings are likely to contain the nodes given to this method. This method will
+     * work even if this assumption does not hold, but the tree will be inefficient in such case.
+     *
+     * <p>A "professional" R-Tree would make a better effort for creating a balanced tree here.
+     * Current version of this class uses a simple algorithm, okay for small trees where big nodes
+     * are added before small nodes. This is not a good general purpose R-Tree.</p>
+     *
+     * @param  node  the node to add to this tree. If this node has sibling, they will be added too.
+     *
+     * @see #finish()
+     */
+    public final void addNode(RTreeNode node) {
+detach: for (RTreeNode next; node != null; node = next) {
+            next = node.sibling;
+            node.sibling = null;                    // Detach the siblings for adding each node individually.
+            RTreeNode tail, receiver = this;
+            do {
+                tail = receiver;
+                if (receiver.tryAddChild(node)) {
+                    continue detach;                // Node accepted by a child, process the next node if any.
+                }
+                receiver = receiver.sibling;
+            } while (receiver != null);
+            tail.sibling = node;                    // Not a child of this node, but assume common parent.
+            node.parent = parent;
+        }
+    }
+
+    /**
+     * Tries to add a child. This method checks if the given candidate is fully included in this {@code RTreeNode}.
+     * The given node may have children but shall not have any {@linkplain #sibling}, because this method does not
+     * check if the siblings would also be included in this node.
+     *
+     * @param  candidate  the single node (without sibling) to add as a child if possible.
+     * @return whether the given node has been added.
+     */
+    private boolean tryAddChild(final RTreeNode candidate) {
+        assert candidate.sibling == null : candidate;
+        if (contains(candidate)) {
+            RTreeNode child = firstChild;
+            if (child == null) {
+                firstChild = candidate;                     // This node had no children before this method call.
+            } else {
+                RTreeNode lastChild;
+                do {
+                    lastChild = child;
+                    if (child.tryAddChild(candidate)) {
+                        return true;                        // Given node added to a child instead than to this node.
+                    }
+                    child = child.sibling;
+                } while (child != null);
+                lastChild.sibling = candidate;              // Add last in the linked list of this node children.
+            }
+            candidate.parent = this;
+            return true;
+        }
+        return false;
+    }
+
+    /**
+     * Finishes the construction of the tree. This methods sets the CRS of all nodes to a common value.
+     * An exception is thrown if incompatible CRS are found. This method does not verify the number of
+     * dimensions; this check should have been done by the caller.
+     *
+     * <p>This method should be invoked only on the instance on which {@link #addNode(RTreeNode)} has been
+     * invoked.</p>
+     *
+     * @return the root of the tree (may be {@code this}).
+     */
+    public final RTreeNode finish() {
+        final Uniformizer action = new Uniformizer();
+        walk(this, action);
+        action.set = true;
+        walk(this, action);
+        RTreeNode next = sibling;
+        if (next == null) {
+            return this;
+        }
+        /*
+         * If there is more than one node, create a synthetic parent containing them all.
+         * We do not need to traverse all node for this task; only the immediate children
+         * of the new parent.
+         */
+        parent = new RTreeNode(this);
+        parent.firstChild = this;
+        do {
+            parent.add(next);           // Compute union of envelopes (not an addition of node).
+            next.parent = parent;
+            next = next.sibling;
+        } while (next != null);
+        return parent;
+    }
+
+    /**
+     * The action for getting a common CRS for all nodes, then for setting the CRS of all nodes.
+     * This action will be executed on the whole tree, with recursive traversal of children.
+     */
+    private static final class Uniformizer implements Consumer<RTreeNode> {
+        /** The CRS common to all nodes. */
+        private CoordinateReferenceSystem common;
+
+        /** {@code false} for collecting information, or {@code true} for setting all node CRS. */
+        boolean set;
+
+        /** Invoked for each node of the tree. */
+        @Override public void accept(final RTreeNode node) {
+            if (set) {
+                node.setCoordinateReferenceSystem(common);
+            } else {
+                final CoordinateReferenceSystem crs = node.getCoordinateReferenceSystem();
+                if (common == null) {
+                    common = crs;
+                } else if (crs != null && !Utilities.equalsIgnoreMetadata(common, crs)) {
+                    throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedCRS));
+                }
+            }
+        }
+    }
+
+    /**
+     * Returns the node that contains the given position, or {@code null} if none.
+     * The given node does not need to be the root of the tree; it can be any node.
+     * This method will check that node first. If it does not contain the position,
+     * then this method goes up in the parent nodes.
+     *
+     * <p>If consecutive positions are close to each other, then a good strategy is
+     * to give to this method the most recently used node.</p>
+     *
+     * @param  node  any node in the tree. Should be node most likely to contain the position.
+     * @param  pos   the position of the node to locate.
+     * @return the smallest node containing the given position.
+     */
+    public static RTreeNode locate(RTreeNode node, final DirectPosition pos) {
+        RTreeNode skip = null;      // For avoiding to invoke `contains(pos)` twice on the same mode.
+        do {
+            if (node.contains(pos)) {
+                /*
+                 * Before to return the node we just found, check if a child contains the point.
+                 * If we find a child containing the point, then `node` is updated to that child.
+                 */
+                RTreeNode candidate = node.firstChild;
+                while (candidate != null) {
+                    if (candidate != skip && candidate.contains(pos)) {
+                        node = candidate;
+                        candidate = node.firstChild;        // Continue checking children of the child.
+                    } else {
+                        candidate = candidate.sibling;      // Check other children of the parent node.
+                    }
+                }
+                return node;
+            }
+            skip = node;
+            node = node.parent;
+        } while (node != null);
+        return null;
+    }
+
+    /**
+     * Returns a hash code value based on the this node and its children,
+     * ignoring the parent and the siblings.
+     *
+     * @return a hash code value for this node and its children.
+     */
+    @Override
+    public int hashCode() {
+        return super.hashCode() + 37 * getChildren().hashCode();
+    }
+
+    /**
+     * Compares this node with the given object for equality, ignoring the parent and the siblings.
+     *
+     * @param  obj  the object to compare with this node.
+     * @return whether the two objects are of the same class, have equal envelope and equal children list.
+     */
+    @Override
+    public boolean equals(final Object obj) {
+        return super.equals(obj) && getChildren().equals(((RTreeNode) obj).getChildren());
+    }
+
+    /**
+     * Formats this envelope as a "{@code DOMAIN}" element (non-standard).
+     */
+    @Override
+    protected String formatTo(final Formatter formatter) {
+        super.formatTo(formatter);
+        return "Domain";
+    }
+
+    /**
+     * Returns a string representation of this node and its children as a tree.
+     * This method is for debugging purposes.
+     *
+     * @return a string representation of this node and all its children.
+     */
+    @Override
+    public String toString() {
+        return toTree().toString();
+    }
+
+    /**
+     * Returns a representation of this node and its children as a tree.
+     * This is mostly for debugging purposes.
+     *
+     * @return a tree representation of this node and all its children.
+     */
+    public final TreeTable toTree() {
+        final DefaultTreeTable tree = new DefaultTreeTable(TableColumn.VALUE);
+        toTree(tree.getRoot());
+        return tree;
+    }
+
+    /**
+     * Adds this node and its children to the given node. This method invokes itself recursively.
+     */
+    private void toTree(final TreeTable.Node addTo) {
+        addTo.setValue(TableColumn.VALUE, super.toString());
+        for (final RTreeNode child : getChildren()) {
+            child.toTree(addTo.newChild());
+        }
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
index c456c97..4d2916f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.java
@@ -290,6 +290,11 @@ public final class Resources extends IndexedResourceBundle {
         public static final short LoadingDatumShiftFile_1 = 32;
 
         /**
+         * Misaligned datum shift grid in “{0}”.
+         */
+        public static final short MisalignedDatumShiftGrid_1 = 94;
+
+        /**
          * The “{1}” parameter could have been omitted. But it has been given a value of {2} which does
          * not match the definition of the “{0}” ellipsoid.
          */
@@ -521,6 +526,11 @@ public final class Resources extends IndexedResourceBundle {
          * Parameter values have not been specified.
          */
         public static final short UnspecifiedParameterValues = 70;
+
+        /**
+         * Using datum shift grid from “{0}” to “{1}” created on {2} (updated on {3}).
+         */
+        public static final short UsingDatumShiftGrid_4 = 93;
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
index 854340c..78af8c7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources.properties
@@ -31,6 +31,8 @@ IgnoredServiceProvider_3          = More than one service provider of type \u201
 InverseOperationUsesSameSign      = Inverse operation uses the same parameter value.
 InverseOperationUsesOppositeSign  = Inverse operation uses this parameter value with opposite sign.
 LoadingDatumShiftFile_1           = Loading datum shift file \u201c{0}\u201d.
+UsingDatumShiftGrid_4             = Using datum shift grid from \u201c{0}\u201d to \u201c{1}\u201d created on {2} (updated on {3}).
+MisalignedDatumShiftGrid_1        = Misaligned datum shift grid in \u201c{0}\u201d.
 MismatchedEllipsoidAxisLength_3   = The \u201c{1}\u201d parameter could have been omitted. But it has been given a value of {2} which does not match the definition of the \u201c{0}\u201d ellipsoid.
 MismatchedOperationFactories_2    = No coordinate operation from \u201c{0}\u201d to \u201c{1}\u201d because of mismatched factories.
 MisnamedParameter_1               = Despite its name, this parameter is effectively \u201c{0}\u201d.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
index 7a32d59..09e9e73 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Resources_fr.properties
@@ -36,6 +36,8 @@ IgnoredServiceProvider_3          = Plusieurs fournisseurs de service de type \u
 InverseOperationUsesSameSign      = L\u2019op\u00e9ration inverse utilise la m\u00eame valeur pour ce param\u00e8tre.
 InverseOperationUsesOppositeSign  = L\u2019op\u00e9ration inverse utilise ce param\u00e8tre avec la valeur de signe oppos\u00e9.
 LoadingDatumShiftFile_1           = Chargement du fichier de changement de r\u00e9f\u00e9rentiel \u00ab\u202f{0}\u202f\u00bb.
+UsingDatumShiftGrid_4             = Utilise la grille de changement de r\u00e9f\u00e9rentiel de \u00ab\u202f{0}\u202f\u00bb vers \u00ab\u202f{1}\u202f\u00bb cr\u00e9\u00e9e le {2} (mise \u00e0 jour le {3}).
+MisalignedDatumShiftGrid_1        = Les grilles de changement de r\u00e9f\u00e9rentiel de \u00ab\u202f{0}\u202f\u00bb ne sont pas align\u00e9es.
 MismatchedEllipsoidAxisLength_3   = Le param\u00e8tre \u00ab\u202f{1}\u202f\u00bb aurait pu \u00eatre omis. Mais il lui a \u00e9t\u00e9 donn\u00e9 la valeur {2} qui ne correspond pas \u00e0 la d\u00e9finition de l\u2019ellipso\u00efde \u00ab\u202f{0}\u202f\u00bb.
 MismatchedOperationFactories_2    = Il n\u2019y a pas d\u2019op\u00e9rations allant de \u00ab\u202f{0}\u202f\u00bb vers \u00ab\u202f{1}\u202f\u00bb parce que ces derniers sont associ\u00e9s \u00e0 deux fabriques diff\u00e9rentes.
 MisnamedParameter_1               = Malgr\u00e9 son nom, ce param\u00e8tre produit en r\u00e9alit\u00e9 l\u2019effet d\u2019un \u00ab\u202f{0}\u202f\u00bb.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
index 346222a..caa06bf 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
@@ -41,7 +41,13 @@ import org.apache.sis.internal.util.Strings;
  *       {@link Envelope2D#contains(Rectangle2D)} or {@link Envelope2D#intersects(Rectangle2D)} methods.</li>
  * </ul>
  *
- * This class does <strong>not</strong> support by itself rectangles spanning the anti-meridian of a geographic CRS.
+ * This class provides the following additional methods which are not defined in {@link Rectangle2D}:
+ * <ul>
+ *   <li>{@link #containsInclusive(double, double)}</li>
+ *   <li>{@link #distanceSquared(double, double)}</li>
+ * </ul>
+ *
+ * This class does <strong>not</strong> support by itself rectangles crossing the anti-meridian of a geographic CRS.
  * However the {@link #getX()}, {@link #getY()}, {@link #getWidth()} and {@link #getHeight()} methods are defined in
  * the straightforward way expected by {@link Envelope2D#intersects(Rectangle2D)} and similar methods for computing
  * correct result if the given {@code Rectangle2D} crosses the anti-meridian.
@@ -53,7 +59,7 @@ import org.apache.sis.internal.util.Strings;
  * recommended approach, but for Apache SIS private classes this is a way to reduce pressure on garbage collector.</div>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.1
  * @since   0.8
  * @module
  */
@@ -370,6 +376,8 @@ public class IntervalRectangle extends Rectangle2D {
 
     /**
      * Tests if a specified coordinate is inside the boundary of this {@code Rectangle2D}.
+     * The maximal <var>x</var> and <var>y</var> values (i.e. the right and bottom edges
+     * of this rectangle) are exclusive.
      *
      * @param  x  the <var>x</var> coordinates to test.
      * @param  y  the <var>y</var> coordinates to test.
@@ -381,6 +389,53 @@ public class IntervalRectangle extends Rectangle2D {
     }
 
     /**
+     * Tests if a specified coordinate is inside the boundary of this {@code Rectangle2D}
+     * with all edges inclusive.
+     *
+     * @param  x  the <var>x</var> coordinates to test.
+     * @param  y  the <var>y</var> coordinates to test.
+     * @return {@code true} if the specified coordinates are inside this rectangle, all edges included.
+     */
+    public final boolean containsInclusive(final double x, final double y) {
+        return (x >= xmin && x <= xmax && y >= ymin && y <= ymax);
+    }
+
+    /**
+     * Returns the square of the minimal distance between a point and this rectangle.
+     * If the point is inside the rectangle or on the edge, then this method returns 0.
+     *
+     * @param  x  the <var>x</var> coordinates to test.
+     * @param  y  the <var>y</var> coordinates to test.
+     * @return square of minimal distance, or 0 if the point is inside this rectangle.
+     */
+    public final double distanceSquared(final double x, final double y) {
+        int outcode = 0;
+        double dx = java.lang.Double.POSITIVE_INFINITY;
+        double dy = java.lang.Double.POSITIVE_INFINITY;
+        double d;
+        if ((d = (x - xmax)) >= 0)           {dx = d; outcode =  1;}
+        if ((d = (y - ymax)) >= 0)           {dy = d; outcode |= 2;}
+        if ((d = (xmin - x)) >= 0 && d < dx) {dx = d; outcode |= 1;}
+        if ((d = (ymin - y)) >= 0 && d < dy) {dy = d; outcode |= 2;}
+        switch (outcode) {
+            case 1:  return dx*dx;                                  // Only x coordinate is outside.
+            case 2:  return dy*dy;                                  // Only y coordinate is outside.
+            case 3:  return dx*dx + dy*dy;                          // Both coordinates are outside.
+            case 0:  assert containsInclusive(x, y); return 0;      // No coordinate is outside.
+            default: throw new AssertionError(outcode);
+        }
+        /*
+         * Note: if we want non-squared distance in a future version, we can rely on the fact
+         * that `dx` and `dy` are guaranteed positives (no need to take the absolute value).
+         * So we would replace the 3 first cases as below:
+         *
+         *     case 1:  return dx;
+         *     case 2:  return dy;
+         *     case 3:  return Math.hypot(dx, dy);
+         */
+    }
+
+    /**
      * Determines where the specified coordinates lie with respect to this rectangle.
      * This method computes a binary OR of the appropriate mask values indicating,
      * for each side of this {@code Rectangle2D}, whether or not the specified coordinates
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Tile.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Tile.java
new file mode 100644
index 0000000..b2ec2ae
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/Tile.java
@@ -0,0 +1,678 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing.j2d;
+
+import java.util.Locale;
+import java.util.Collection;
+import java.util.Optional;
+import java.io.Writer;
+import java.io.StringWriter;
+import java.io.Serializable;
+import java.io.IOException;
+import java.awt.Point;
+import java.awt.Dimension;
+import java.awt.Rectangle;
+import java.awt.geom.AffineTransform;
+import javax.imageio.ImageReader;                           // For javadoc
+import org.opengis.metadata.spatial.PixelOrientation;
+import org.apache.sis.util.resources.Vocabulary;
+import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.util.Classes;
+import org.apache.sis.io.TableAppender;
+
+
+/**
+ * A tile identified by a location, a dimension and a subsampling.
+ * This class can be used for constructing a mosaic or a pyramid of images.
+ * While the Javadoc discusses image I/O operations, this {@code Tile} class is not restricted to imagery.
+ * This class is also used for managing tiles in a datum shift file encoded in NTv2 format.
+ *
+ * <p>Each tile contains the following:</p>
+ * <ul class="verbose">
+ *   <li><b>A format name or a provider of {@link ImageReader}</b> (optional).
+ *   The same format is typically used for every tiles, but this is not mandatory.
+ *   An {@linkplain ImageReader image reader} can be instantiated before a tile is read.</li>
+ *
+ *   <li><b>An image input</b> (optional), typically a {@link java.nio.file.Path} or {@link java.net.URL}.
+ *   The input is often different for every tile to be read, but this is not mandatory. For example tiles
+ *   could be stored at different {@linkplain #getImageIndex() image index} in the same file.</li>
+ *
+ *   <li><b>An image index</b> to be given to {@link ImageReader#read(int)} for reading the tile.
+ *   This index is often 0.</li>
+ *
+ *   <li><b>The upper-left corner</b> in the destination image as a {@link Point},
+ *   or the upper-left corner together with the image size as a {@link Rectangle}.
+ *   If the upper-left corner has been given as a point, then the
+ *   {@linkplain ImageReader#getWidth(int) width} and {@linkplain ImageReader#getHeight(int) height}
+ *   may be obtained from the image reader when first needed, which may have a slight performance cost.
+ *   If the upper-left corner has been given as a rectangle instead, then this performance cost is avoided
+ *   but the user is responsible for the accuracy of the information provided.
+ *
+ *     <div class="note"><b>NOTE:</b>
+ *     the upper-left corner is the {@linkplain #getLocation() location} of this tile in the
+ *     {@linkplain javax.imageio.ImageReadParam#setDestination destination image} when no
+ *     {@linkplain javax.imageio.ImageReadParam#setDestinationOffset destination offset} are specified.
+ *     If the user specified a destination offset, then the tile location will be translated accordingly
+ *     for the image being read.
+ *     </div></li>
+ *
+ *   <li><b>The subsampling relative to the tile having the best resolution.</b>
+ *   This is not the subsampling to apply when reading this tile, but rather the subsampling that we would
+ *   need to apply on the tile having the finest resolution in order to produce an image equivalent to this tile.
+ *   The subsampling is (1,1) for the tile having the finest resolution, (2,3) for an overview having
+ *   half the width and third of the height for the same geographic extent, <i>etc.</i>
+ *   (note that overviews are not required to have the same geographic extent - the above is just an example).
+ *
+ *     <div class="note"><b>NOTE 1:</b>
+ *     the semantic assumes that overviews are produced by subsampling, not by interpolation or pixel averaging.
+ *     The later are not prohibited, but doing so introduce some subsampling-dependent variations in images read,
+ *     which would not be what we would expect from a strictly compliant {@link ImageReader}.</div>
+ *
+ *     <div class="note"><b>NOTE 2:</b>
+ *     tile {@linkplain #getLocation() location} and {@linkplain #getRegion() region} coordinates should be
+ *     specified in the overview pixel units - they should <em>not</em> be pre-multiplied by subsampling.
+ *     This multiplication should be performed automatically by a {@code TileManager} when comparing regions
+ *     from tiles at different subsampling levels.
+ *     </div></li>
+ * </ul>
+ *
+ * The tiles are not required to be arranged on a regular grid, but performances may be better if they are.
+ * {@link TileOrganizer} is responsible for analyzing the layout of a collection of tiles.
+ *
+ * <h2>Multi-threading</h2>
+ * This class is thread-safe. In addition {@code Tile} instances can be considered as immutable after construction.
+ * However some properties may be available only after the tiles have been processed by a {@link TileOrganizer},
+ * or only after {@link #fetchSize()} has been invoked.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public class Tile implements Serializable {
+    /**
+     * For cross-version compatibility during serialization.
+     */
+    private static final long serialVersionUID = 1638238437701248681L;
+
+    /**
+     * The upper-left corner in the mosaic (destination image). Should be considered as final,
+     * since this class is supposed to be mostly immutable. However the value can be changed
+     * by {@link #translate(int, int)} before the {@code Tile} instance is made public.
+     *
+     * @see #getLocation()
+     * @see #getRegion()
+     */
+    private int x, y;
+
+    /**
+     * The size of the tile, or 0 if not yet computed.
+     *
+     * @see #getSize()
+     * @see #getRegion()
+     */
+    private int width, height;
+
+    /**
+     * The subsampling relative to the tile having the finest resolution. If this tile is the one with
+     * finest resolution, then the value shall be 1. Should never be 0 or negative, except if the value
+     * has not yet been computed.
+     *
+     * <p>This field should be considered as final. It is not final only because
+     * {@link TileOrganizer} may compute this value automatically.</p>
+     *
+     * @see #getSubsampling()
+     */
+    private int xSubsampling, ySubsampling;
+
+    /**
+     * The "grid to real world" transform, used by {@link TileOrganizer} in order to compute
+     * the {@linkplain #getRegion() region} for this tile. This field is set to {@code null} when
+     * {@link TileOrganizer}'s work is in progress, and set to a new value on completion.
+     *
+     * <p><b>Note:</b> {@link TileOrganizer} really needs a new instance for each tile.
+     * No caching allowed before {@link TileOrganizer} processing.
+     * Caching is allowed <em>after</em> {@link TileOrganizer} processing is completed.</p>
+     */
+    private AffineTransform gridToCRS;
+
+    /**
+     * Creates a tile for the given tile location. This constructor can be used when the size of
+     * the tile is unknown. This tile size will be fetched automatically by {@link #fetchSize()}
+     * when {@link #getSize()} or {@link #getRegion()} is invoked for the first time.
+     *
+     * @param location     the upper-left corner in the mosaic (destination image).
+     * @param subsampling  the subsampling relative to the tile having the finest resolution,
+     *                     or {@code null} if none. If non-null, width and height shall be strictly positive.
+     *                     This argument can be understood as pixel size relative to finest resolution.
+     */
+    public Tile(final Point location, final Dimension subsampling) {
+        ArgumentChecks.ensureNonNull("location", location);
+        x = location.x;
+        y = location.y;
+        setSubsampling(subsampling);
+    }
+
+    /**
+     * Creates a tile for the given region. This constructor should be used when the size of the tile is known.
+     * This information avoid the cost of fetching the size when {@link #getSize()} or {@link #getRegion()} is
+     * first invoked.
+     *
+     * @param region       the region (location and size) in the mosaic (destination image).
+     * @param subsampling  the subsampling relative to the tile having the finest resolution,
+     *                     or {@code null} if none. If non-null, width and height shall be strictly positive.
+     *                     This argument can be understood as pixel size relative to finest resolution.
+     * @throws IllegalArgumentException if the given region {@linkplain Rectangle#isEmpty() is empty}.
+     */
+    public Tile(final Rectangle region, final Dimension subsampling) {
+        ArgumentChecks.ensureNonNull("location", region);
+        x      = region.x;
+        y      = region.y;
+        width  = region.width;
+        height = region.height;
+        if (width <= 0 || height <= 0) {
+            throw new IllegalArgumentException(Errors.format(Errors.Keys.EmptyArgument_1, "region"));
+        }
+        setSubsampling(subsampling);
+    }
+
+    /**
+     * Creates a tile for the given region and <cite>"grid to real world"</cite> transform.
+     * This constructor can be used when the {@linkplain #getLocation() location} of the tile is unknown.
+     * The location and subsampling will be computed automatically when this tile will be processed by a
+     * {@link TileOrganizer}.
+     *
+     * <p>When using this constructor, the {@link #getLocation()}, {@link #getRegion()} and
+     * {@link #getSubsampling()} methods will throw an {@link IllegalStateException} until
+     * this tile has been processed by a {@link TileOrganizer}, which will compute those
+     * values automatically.</p>
+     *
+     * @param region     the tile region, or {@code null} if unknown.
+     *                   The (<var>x</var>,<var>y</var> location of this region is typically (0,0).
+     *                   The final location will be computed when this tile will be given to a {@link TileOrganizer}.
+     * @param gridToCRS  the <cite>"grid to real world"</cite> transform mapping pixel
+     *                   {@linkplain PixelOrientation#UPPER_LEFT upper left} corner.
+     */
+    public Tile(final Rectangle region, final AffineTransform gridToCRS) {
+        ArgumentChecks.ensureNonNull("gridToCRS", gridToCRS);
+        if (region != null) {
+            x      = region.x;
+            y      = region.y;
+            width  = Math.max(region.width,  0);                        // Empty region authorized.
+            height = Math.max(region.height, 0);
+        }
+        this.gridToCRS = new AffineTransform(gridToCRS);                // Really need a new instance - no cache
+    }
+
+    /**
+     * Creates a new tile for the given final transform.
+     * This is used for storing {@link TileOrganizer} results.
+     */
+    Tile(final AffineTransform gridToCRS, final Rectangle region) {
+        this.x         = region.x;
+        this.y         = region.y;
+        this.width     = region.width;
+        this.height    = region.height;
+        this.gridToCRS = gridToCRS;                                     // Should be an AffineTransform2D instance.
+        setSubsampling(null);
+    }
+
+    /**
+     * Checks if the location, region, and subsampling can be returned. Throws an exception if this
+     * tile has been created without location and not yet processed by {@link TileOrganizer}.
+     */
+    private void ensureDefined() throws IllegalStateException {
+        if (xSubsampling == 0 || ySubsampling == 0) {
+            throw new IllegalStateException();
+        }
+    }
+
+    /**
+     * Returns the tile upper-left corner coordinates in the mosaic.
+     *
+     * @return the tile upper-left corner.
+     * @throws IllegalStateException if this tile has been {@linkplain #Tile(Rectangle, AffineTransform)
+     *         created without location} and has not yet been processed by {@link TileOrganizer}.
+     *
+     * @see javax.imageio.ImageReadParam#setDestinationOffset(Point)
+     */
+    public synchronized Point getLocation() throws IllegalStateException {
+        ensureDefined();
+        return new Point(x, y);
+    }
+
+    /**
+     * Returns the image size. If this tile has been created with the {@linkplain #Tile(Rectangle, Dimension)
+     * constructor expecting a rectangle}, then the dimension of that rectangle is returned.
+     * Otherwise {@link #fetchSize()} is invoked and the result is cached for future usage.
+     *
+     * <p>At the difference of {@link #getLocation()} and {@link #getRegion()}, this method never
+     * throw {@link IllegalStateException} because the tile size does not depend on the processing
+     * performed by {@link TileOrganizer}.</p>
+     *
+     * @return the tile size.
+     * @throws IOException if an I/O operation was required for fetching the tile size and that operation failed.
+     * @throws IllegalStateException if this class does not have sufficient information for providing a tile size.
+     */
+    public synchronized Dimension getSize() throws IOException {
+        // No call to ensureDefined().
+        if ((width | height) == 0) {
+            final Dimension size = fetchSize();
+            width  = size.width;
+            height = size.height;
+        }
+        return new Dimension(width, height);
+    }
+
+    /**
+     * Invoked when the tile size need to be read or computed. The default implementation throws
+     * {@link IllegalStateException} since this base class has no information for computing a tile size.
+     * Subclasses can override and, for example, get the size with {@link ImageReader#getWidth(int)} and
+     * {@link ImageReader#getHeight(int)}.
+     *
+     * @return the tile size.
+     * @throws IOException if an I/O operation was required for fetching the tile size and that operation failed.
+     * @throws IllegalStateException if this class does not have sufficient information for providing a tile size.
+     */
+    protected Dimension fetchSize() throws IOException {
+        throw new IllegalStateException();
+    }
+
+    /**
+     * Returns the upper-left corner location in the mosaic together with the tile size.
+     * If this tile has been created with the {@linkplain #Tile(Rectangle, Dimension)
+     * constructor expecting a rectangle}, a copy of the specified rectangle is returned.
+     * Otherwise {@link #fetchSize()} is invoked and the result is cached for future usage.
+     *
+     * @return the region in the mosaic (destination image).
+     * @throws IOException if an I/O operation was required for fetching the tile size and that operation failed.
+     * @throws IllegalStateException if this tile has been {@linkplain #Tile(Rectangle, AffineTransform) created
+     *         without location} and has not yet been processed by {@link TileOrganizer}, of if this tile does
+     *         not have enough information for providing a tile size.
+     *
+     * @see javax.imageio.ImageReadParam#setSourceRegion(Rectangle)
+     */
+    public synchronized Rectangle getRegion() throws IllegalStateException, IOException {
+        ensureDefined();
+        if ((width | height) == 0) {
+            final Dimension size = fetchSize();
+            width  = size.width;
+            height = size.height;
+        }
+        return new Rectangle(x, y, width, height);
+    }
+
+    /**
+     * Returns the {@linkplain #getRegion() region} multiplied by the subsampling.
+     * This is this tile coordinates in the units of the tile having the finest resolution,
+     * as opposed to other methods which are always in units relative to this tile.
+     *
+     * @return the region in units relative to the tile having the finest resolution.
+     * @throws IOException if an I/O operation was required for fetching the tile size and that operation failed.
+     * @throws IllegalStateException if this tile has been {@linkplain #Tile(Rectangle, AffineTransform) created
+     *         without location} and has not yet been processed by {@link TileOrganizer}, of if this tile does
+     *         not have enough information for providing a tile size.
+     */
+    public synchronized Rectangle getAbsoluteRegion() throws IOException {
+        final Rectangle region = getRegion();
+        final int sx = xSubsampling;
+        final int sy = ySubsampling;
+        region.x      *= sx;
+        region.y      *= sy;
+        region.width  *= sx;
+        region.height *= sy;
+        return region;
+    }
+
+    /**
+     * Invoked by {@link TileOrganizer} only. No other caller allowed.
+     * {@link #setSubsampling(Dimension)} must be invoked prior this method.
+     *
+     * <p>Note that invoking this method usually invalidate {@link #gridToCRS}.
+     * Calls to this method should be followed by {@link #translate(int, int)}
+     * for fixing the "gridToCRS" value.</p>
+     *
+     * @param  region  the region to assign to this tile in units of tile having finest resolution.
+     * @throws ArithmeticException if {@link #setSubsampling(Dimension)} method has not be invoked.
+     */
+    final void setAbsoluteRegion(final Rectangle region) throws ArithmeticException {
+        assert Thread.holdsLock(this);
+        final int sx = xSubsampling;
+        final int sy = ySubsampling;
+        assert (region.width % sx) == 0 && (region.height % sy) == 0 : region;
+        x      = region.x      / sx;
+        y      = region.y      / sy;
+        width  = region.width  / sx;
+        height = region.height / sy;
+    }
+
+    /**
+     * Returns the subsampling relative to the tile having the finest resolution.
+     * The return value can be interpreted as "pixel size" relative to tiles having the finest resolution.
+     * This method never return {@code null}, and the width and height shall never be smaller than 1.
+     *
+     * @return the subsampling along <var>x</var> and <var>y</var> axes.
+     * @throws IllegalStateException if this tile has been {@linkplain #Tile(Rectangle, AffineTransform)
+     *         created without location} and has not yet been processed by {@link TileOrganizer}.
+     *
+     * @see javax.imageio.ImageReadParam#setSourceSubsampling(int, int, int, int)
+     */
+    public synchronized Dimension getSubsampling() throws IllegalStateException {
+        ensureDefined();
+        return new Dimension(xSubsampling, ySubsampling);
+    }
+
+    /**
+     * Sets the subsampling to the given dimension.
+     * Invoked by constructors and {@link TileOrganizer} only.
+     */
+    final void setSubsampling(final Dimension subsampling) throws IllegalStateException {
+        // No assert Thread.holdsLock(this) because invoked from constructors.
+        if ((xSubsampling | ySubsampling) != 0) {
+            throw new IllegalStateException();                          // Should never happen.
+        }
+        if (subsampling != null) {
+            ArgumentChecks.ensureBetween("width",  0, 0xFFFF, subsampling.width);
+            ArgumentChecks.ensureBetween("height", 0, 0xFFFF, subsampling.height);
+            xSubsampling = subsampling.width;
+            ySubsampling = subsampling.height;
+        } else {
+            xSubsampling = ySubsampling = 1;
+        }
+    }
+
+    /**
+     * If the user-supplied transform is waiting for processing by {@link TileOrganizer}, returns it.
+     * Otherwise returns {@code null}. This method is for internal usage by {@link TileOrganizer} only.
+     *
+     * <p>This method clears the {@link #gridToCRS} field before to return. This is a way to tell that
+     * processing is in progress, and also a safety against transform usage while it may become invalid.</p>
+     *
+     * @return the transform, or {@code null} if none. This method does not clone the returned value -
+     *         {@link TileOrganizer} will reference and modify directly that transform.
+     */
+    final synchronized AffineTransform getPendingGridToCRS() {
+        if ((xSubsampling | ySubsampling) != 0) {
+            // No transform waiting to be processed.
+            return null;
+        }
+        final AffineTransform at = gridToCRS;
+        gridToCRS = null;
+        return at;
+    }
+
+    /**
+     * Returns the <cite>"grid to real world"</cite> transform, or {@code null} if unknown.
+     * This transform is derived from the value given to the constructor, but may not be identical
+     * since it may have been {@linkplain AffineTransform#translate(double, double) translated}
+     * in order to get a uniform grid geometry for every tiles.
+     *
+     * <div class="note"><b>Tip:</b>
+     * the <a href="https://en.wikipedia.org/wiki/World_file">World File</a> coefficients of this tile
+     * (i.e. the <cite>grid to CRS</cite> transform that we would have if the pixel in the upper-left
+     * corner always had indices (0,0)) can be computed as below:
+     *
+     * {@preformat java
+     *     Point location = tile.getLocation();
+     *     AffineTransform gridToCRS = new AffineTransform(tile.getGridToCRS());
+     *     gridToCRS.translate(location.x, location.y);
+     * }
+     * </div>
+     *
+     * @return the <cite>"grid to real world"</cite> transform mapping pixel
+     *         {@linkplain PixelOrientation#UPPER_LEFT upper left} corner, or {@code null} if undefined.
+     * @throws IllegalStateException if this tile has been {@linkplain #Tile(Rectangle, AffineTransform)
+     *         created without location} and has not yet been processed by {@link TileOrganizer}.
+     */
+    public synchronized AffineTransform2D getGridToCRS() throws IllegalStateException {
+        ensureDefined();
+        /*
+         * The cast should not fail: if the `gridToCRS` is the one specified at construction time,
+         * then `ensureDefined()` should have thrown an IllegalStateException. Otherwise this tile
+         * have been processed by `TileOrganizer`, which has set an `AffineTransform2D` instance.
+         * If we get a ClassCastException below, then there is a bug in our pre/post conditions.
+         */
+        return (AffineTransform2D) gridToCRS;
+    }
+
+    /**
+     * Sets the new <cite>"grid to real world"</cite> transform to use after the translation performed by
+     * {@link #translate(int, int)}, if any. The given instance should be immutable; it will not be cloned.
+     *
+     * @param  at  the <cite>"grid to real world"</cite> transform mapping pixel
+     *             {@linkplain PixelOrientation#UPPER_LEFT upper left} corner.
+     * @throws IllegalStateException if another transform was already assigned to this tile.
+     */
+    final void setGridToCRS(final AffineTransform at) throws IllegalStateException {
+        assert Thread.holdsLock(this);
+        if (gridToCRS == null) {
+            gridToCRS = at;
+        } else if (!gridToCRS.equals(at)) {
+            throw new IllegalStateException();
+        }
+    }
+
+    /**
+     * Translates this tile. For internal usage by {@link TileOrganizer} only.
+     *
+     * <p>Reminder: {@link #setGridToCRS(AffineTransform)} should be invoked after this method.</p>
+     *
+     * @param  dx  the translation to apply on <var>x</var> values (often 0).
+     * @param  dy  the translation to apply on <var>y</var> values (often 0).
+     */
+    final void translate(final int dx, final int dy) {
+        assert Thread.holdsLock(this);
+        x += dx;
+        y += dy;
+        gridToCRS = null;
+    }
+
+    /**
+     * Returns a name for the tile format or tile input, or an empty value if none.
+     * The format name can be inferred for example from an {@link javax.imageio.spi.ImageReaderSpi}.
+     * The input name is typically (but not necessarily) a file name or URL.
+     *
+     * @param  input  {@code false} for the file format name, or {@code true} for the file input name.
+     * @return the format or input name.
+     */
+    public Optional<String> getName(final boolean input) {
+        return Optional.empty();
+    }
+
+    /**
+     * Returns the image index to be given to the image reader for reading this tile.
+     * The default implementation returns 0.
+     *
+     * @return the image index, numbered from 0.
+     *
+     * @see ImageReader#read(int)
+     */
+    public int getImageIndex() {
+        return 0;
+    }
+
+    /*
+     * Intentionally no implementation for `equals()` and `hashCode()`. Tile is an "almost immutable" class
+     * which can still be modified (only once) by MocaicCalculator, or by read operations during `getSize()`
+     * or `getRegion()` execution. This causes confusing behavior when used in an HashMap. We are better to
+     * rely on system identity. For example `DatumShiftGridGroup` rely on the capability to locate Tiles in
+     * HashMap before and after they have been processed by `TileOrganizer`.
+     */
+
+    /**
+     * Returns a string representation of this tile for debugging purposes.
+     *
+     * @return a string representation of this tile.
+     */
+    @Override
+    public String toString() {
+        final StringBuilder buffer = new StringBuilder(Classes.getShortClassName(this)).append('[');
+        if ((xSubsampling | ySubsampling) != 0) {
+            buffer.append("location=(");
+            if (width == 0 && height == 0) {
+                final Point location = getLocation();
+                buffer.append(location.x).append(',').append(location.y);
+            } else try {
+                final Rectangle region = getRegion();
+                buffer.append(region.x).append(',').append(region.y)
+                      .append("), size=(").append(region.width).append(',').append(region.height);
+            } catch (IOException e) {
+                /*
+                 * Should not happen since we checked that `getRegion()` should be easy.
+                 * If it happen anyway, put the exception message at the place where
+                 * coordinates were supposed to appear, so we can debug.
+                 */
+                buffer.append(e);
+            }
+            final Dimension subsampling = getSubsampling();
+            buffer.append("), subsampling=(").append(subsampling.width)
+                  .append(',').append(subsampling.height).append(')');
+        } else {
+            /*
+             * Location and subsampling not yet computed, so don't display it. We can not
+             * invoke `getRegion()` neither since it would throw an IllegalStateException.
+             * Since we have to read the fields directly, make sure that this instance is
+             * not a subclass, otherwise those values may be wrong.
+             */
+            if ((width != 0 || height != 0) && getClass() == Tile.class) {
+                buffer.append("size=(").append(width).append(',').append(height).append(')');
+            }
+        }
+        return buffer.append(']').toString();
+    }
+
+    /**
+     * Returns a string representation of a collection of tiles.
+     * The tiles are formatted in a table in iteration order.
+     *
+     * <p>This method is not public because it can consume a large amount of memory (the underlying
+     * {@link StringBuffer} can be quite large). Users are encouraged to use the method expecting a
+     * {@link Writer}, which may be expensive too but less than this method.</p>
+     *
+     * @param tiles    the tiles to format in a table.
+     * @param maximum  the maximum number of tiles to format. If there is more tiles, a message will be
+     *                 formatted below the table. A reasonable value like 5000 is recommended because
+     *                 attempt to format millions of tiles leads to {@link OutOfMemoryError}.
+     * @return a string representation of the given tiles as a table.
+     */
+    static String toString(final Collection<Tile> tiles, final int maximum) {
+        final StringWriter writer = new StringWriter();
+        try {
+            writeTable(tiles, writer, maximum);
+        } catch (IOException e) {
+            // Should never happen since we are writing to a StringWriter.
+            throw new AssertionError(e);
+        }
+        return writer.toString();
+    }
+
+    /**
+     * Formats a collection of tiles in a table.
+     * The tiles are appended in iteration order.
+     *
+     * @param tiles    the tiles to format in a table.
+     * @param out      where to write the table.
+     * @param maximum  the maximum number of tiles to format. If there is more tiles, a message will be
+     *                 formatted below the table. A reasonable value like 5000 is recommended because
+     *                 attempt to format millions of tiles leads to {@link OutOfMemoryError}.
+     * @throws IOException if an error occurred while writing to the given writer.
+     */
+    public static void writeTable(final Collection<Tile> tiles, final Writer out, final int maximum) throws IOException {
+        int remaining = maximum;
+        final TableAppender table = new TableAppender(out);
+        table.setMultiLinesCells(false);
+        table.nextLine('═');
+        table.append("Format\tInput\tindex\tx\ty\twidth\theight\tdx\tdy\n");
+        table.nextLine('─');
+        table.setMultiLinesCells(true);
+        for (final Tile tile : tiles) {
+            if (--remaining < 0) {
+                break;
+            }
+            table.setCellAlignment(TableAppender.ALIGN_LEFT);
+            tile.getName(false).ifPresent(table::append); table.nextColumn();
+            tile.getName(true) .ifPresent(table::append); table.nextColumn();
+            table.setCellAlignment(TableAppender.ALIGN_RIGHT);
+            table.append(String.valueOf(tile.getImageIndex()));
+            table.nextColumn();
+            /*
+             * Extract now the tile information that we are going to format. Those information may
+             * be replaced by the information provided by getter methods (they should be the same,
+             * unless a subclass override those methods).
+             */
+            int x            = tile.x;
+            int y            = tile.y;
+            int width        = tile.width;
+            int height       = tile.height;
+            int xSubsampling = tile.xSubsampling;
+            int ySubsampling = tile.ySubsampling;
+            try {
+                final Dimension subsampling = tile.getSubsampling();
+                xSubsampling = subsampling.width;
+                ySubsampling = subsampling.height;
+                try {
+                    final Rectangle region = tile.getRegion();
+                    x      = region.x;
+                    y      = region.y;
+                    width  = region.width;
+                    height = region.height;
+                } catch (IOException e) {
+                    /*
+                     * The (x,y) are likely to be correct since only (width,height) are read
+                     * from the image file. So set only (width,height) to "unknown" and keep
+                     * the remaining, with (x,y) obtained from direct access to Tile fields.
+                     */
+                    width  = 0;
+                    height = 0;
+                }
+            } catch (IllegalStateException e) {
+                // Ignore. Format using the information read from the fields as a fallback.
+            }
+            table.append(String.valueOf(x));
+            table.nextColumn();
+            table.append(String.valueOf(y));
+            if ((width | height) != 0) {
+                table.nextColumn();
+                table.append(String.valueOf(width));
+                table.nextColumn();
+                table.append(String.valueOf(height));
+            } else {
+                table.nextColumn();
+                table.nextColumn();
+            }
+            if ((xSubsampling | ySubsampling) != 0) {
+                table.nextColumn();
+                table.append(String.valueOf(xSubsampling));
+                table.nextColumn();
+                table.append(String.valueOf(ySubsampling));
+            }
+            table.nextLine();
+        }
+        table.nextLine('═');
+        /*
+         * Table completed. Flushs to the writer and appends additional text if we have
+         * not formatted every tiles. IOException may be thrown starting from this point
+         * (the above code is not expected to thrown any IOException).
+         */
+        table.flush();
+        if (remaining < 0) {
+            out.write(Vocabulary.getResources((Locale) null).getString(Vocabulary.Keys.More_1, tiles.size() - maximum));
+            out.write(System.lineSeparator());
+        }
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/TileOrganizer.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/TileOrganizer.java
new file mode 100644
index 0000000..8be2c96
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/TileOrganizer.java
@@ -0,0 +1,515 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing.j2d;
+
+import java.util.List;
+import java.util.Arrays;
+import java.util.ArrayList;
+import java.util.Map;
+import java.util.HashMap;
+import java.util.IdentityHashMap;
+import java.util.Collection;
+import java.util.Comparator;
+import java.util.Iterator;
+import java.io.IOException;
+import java.awt.Point;
+import java.awt.Dimension;
+import java.awt.Rectangle;
+import java.awt.geom.AffineTransform;
+import java.awt.geom.NoninvertibleTransformException;
+import java.awt.geom.Rectangle2D;
+import org.apache.sis.referencing.operation.matrix.AffineTransforms2D;
+
+
+/**
+ * Creates a collection of {@link Tile}s from their <cite>grid to CRS</cite> affine transforms.
+ * When the {@link Rectangle} that describe the destination region is known for each tiles,
+ * the {@link Tile#Tile(Rectangle, Dimension)} constructor should be invoked directly.
+ * But in some cases the destination rectangle is not known directly. Instead we have a set of tiles,
+ * all of them with an upper-left corner located at (0,0), but different <cite>grid to CRS</cite>
+ * affine transforms read from <a href="https://en.wikipedia.org/wiki/World_file">World Files</a>.
+ * This {@code TileOrganizer} class infers the destination regions automatically
+ * from the set of affine transforms.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public class TileOrganizer {
+    /**
+     * Small number for floating point comparisons.
+     */
+    private static final double EPS = 1E-10;
+
+    /**
+     * The location of the final bounding box (the one including every tiles).
+     * Tiles will be translated as needed in order to fit this location. This
+     * is usually zero, but not necessarily.
+     */
+    private final int xLocation, yLocation;
+
+    /**
+     * Tiles for which we should compute the bounding box after we have them all.
+     * Their bounding box (region) will need to be adjusted for the affine transform.
+     */
+    private final Map<AffineTransform,Tile> tiles;
+
+    /**
+     * Creates an initially empty tile collection with the given location.
+     *
+     * @param  location  the location, or {@code null} for (0,0).
+     */
+    public TileOrganizer(final Point location) {
+        if (location != null) {
+            xLocation = location.x;
+            yLocation = location.y;
+        } else {
+            xLocation = yLocation = 0;
+        }
+        /*
+         * We really need an IdentityHashMap, not an ordinary HashMap, because we will
+         * put many AffineTransforms that are equal in the sense of Object.equals, but
+         * we still want to associate them to different Tile instances.
+         */
+        tiles = new IdentityHashMap<>();
+    }
+
+    /**
+     * Returns the location of the tile collections to be created. The location is often
+     * (0,0) as expected in {@link java.awt.image.BufferedImage}, but does not have to.
+     *
+     * @return origin of the tile collections to be created.
+     */
+    public Point getLocation() {
+        return new Point(xLocation, yLocation);
+    }
+
+    /**
+     * Adds a tile to the collection of tiles to process.
+     * Each tile can be added only once.
+     *
+     * @param  tile  the tile to add.
+     * @return {@code true} if the tile has been successfully added, or
+     *         {@code false} if the tile does not need to be processed by this class.
+     */
+    public boolean add(final Tile tile) {
+        final AffineTransform gridToCRS = tile.getPendingGridToCRS();
+        if (gridToCRS == null) {
+            return false;
+        }
+        if (tiles.putIfAbsent(gridToCRS, tile) != null) {
+            throw new IllegalStateException();              // Tile already present.
+        }
+        return true;
+    }
+
+    /**
+     * Returns the tiles. Keys are pyramid geometry (containing mosaic bounds and <cite>grid to CRS</cite>
+     * transforms) and values are the tiles in that pyramid. This method usually returns a singleton map,
+     * but more entries may be present if this method was not able to build a single pyramid using all
+     * provided tiles.
+     *
+     * <p><strong>Invoking this method clear the collection</strong>. On return, this instance is empty.
+     * This is because current implementation modify its workspace directly for efficiency.</p>
+     *
+     * @return all tiles added to this {@code TileOrganizer}, grouped by pyramids.
+     * @throws IOException if a call to {@link Tile#getSize()} or {@link Tile#getRegion()} failed,
+     *         and {@link #unavailableSize(Tile, IOException)} did not consumed the exception.
+     */
+    public Map<Tile,Tile[]> tiles() throws IOException {
+        final Map<Tile,Tile[]> results = new HashMap<>(4);
+        for (final Map<AffineTransform,Dimension> tilesAT : computePyramidLevels(tiles.keySet())) {
+            /*
+             * Picks an affine transform to be used as the reference one. We need the finest one.
+             * If more than one have the finest resolution, we pickup the one that will lead to a
+             * (0,0) translation at the end of this method. This is because while the final result
+             * is expected to have integer translation terms, the intermediates results before the
+             * final translation may have fractional terms. Since those intermediate results are
+             * stored as integers in Tile fields, it can leads to errors.
+             */
+            AffineTransform reference = null;
+            double xMin  = Double.POSITIVE_INFINITY;
+            double xLead = Double.POSITIVE_INFINITY;            // Minimum on the first row only.
+            double yMin  = Double.POSITIVE_INFINITY;
+            double scale = Double.POSITIVE_INFINITY;
+            for (final AffineTransform tr : tilesAT.keySet()) {
+                final double s = AffineTransforms2D.getScale(tr);
+                double y = tr.getTranslateY(); if (tr.getScaleY() < 0 || tr.getShearY() < 0) y = -y;
+                double x = tr.getTranslateX(); if (tr.getScaleX() < 0 || tr.getShearX() < 0) x = -x;
+                if (!(Math.abs(s - scale) <= EPS)) {
+                    if (!(s < scale)) continue;                 // '!' is for catching NaN.
+                    scale = s;                                  // Found a smaller scale.
+                    yMin = y;
+                    xMin = x;
+                } else {                                        // Found a transform with the same scale.
+                    if (x < xMin) xMin = x;
+                    if (!(Math.abs(y - yMin) <= EPS)) {
+                        if (!(y < yMin)) continue;
+                        yMin = y;                               // Found a smaller y.
+                    } else if (!(x < xLead)) continue;
+                }
+                xLead = x;
+                reference = tr;
+            }
+            /*
+             * If there is missing tiles at the beginning of the first row, then the x location
+             * of the first tile is greater than the "true" minimum. We will need to adjust.
+             */
+            if (reference == null) {
+                continue;
+            }
+            xLead -= xMin;
+            if (xLead > EPS) {
+                final double[] matrix = new double[6];
+                reference.getMatrix(matrix);
+                matrix[4] = xMin;
+                reference = new AffineTransform(matrix);
+            } else {
+                reference = new AffineTransform(reference);     // Protects from upcoming changes.
+            }
+            /*
+             * Transform the image bounding box from its own space to the reference space.
+             * If `computePyramidLevels(…)` did its job correctly, the transform should contain only a
+             * scale and a translation - no shear (we do not put assertions because of rounding errors).
+             * In such particular case, transforming a Rectangle2D is accurate. We round (we do not clip
+             * as in the default Rectangle implementation) because we really expect integer results.
+             */
+            final AffineTransform toGrid;
+            try {
+                toGrid = reference.createInverse();
+            } catch (NoninvertibleTransformException e) {
+                throw new IllegalStateException(e);
+            }
+            int index = 0;
+            Rectangle groupBounds = null;
+            final Rectangle2D.Double envelope = new Rectangle2D.Double();
+            final Tile[] tilesArray = new Tile[tilesAT.size()];
+            for (final Map.Entry<AffineTransform,Dimension> entry : tilesAT.entrySet()) {
+                final AffineTransform tr = entry.getKey();
+                Tile tile = tiles.remove(tr);                   // Should never be null.
+                tr.preConcatenate(toGrid);
+                /*
+                 * Compute the transformed bounds. If we fail to obtain it, there is probably something wrong
+                 * with the tile but this is not fatal to this method. In such case we will transform only the
+                 * location instead of the full box, which sometime implies a lost of accuracy but not always.
+                 */
+                Rectangle bounds;
+                synchronized (tile) {
+                    tile.setSubsampling(entry.getValue());
+                    try {
+                        bounds = tile.getRegion();
+                    } catch (IOException exception) {
+                        if (!unavailableSize(tile, exception)) {
+                            throw exception;
+                        }
+                        bounds = null;
+                    }
+                    if (bounds != null) {
+                        AffineTransforms2D.transform(tr, bounds, envelope);
+                        bounds.x      = Math.toIntExact(Math.round(envelope.x));
+                        bounds.y      = Math.toIntExact(Math.round(envelope.y));
+                        bounds.width  = Math.toIntExact(Math.round(envelope.width));
+                        bounds.height = Math.toIntExact(Math.round(envelope.height));
+                    } else {
+                        final Point location = tile.getLocation();
+                        tr.transform(location, location);
+                        bounds = new Rectangle(location.x, location.y, 0, 0);
+                    }
+                    tile.setAbsoluteRegion(bounds);
+                }
+                if (groupBounds == null) {
+                    groupBounds = bounds;
+                } else {
+                    groupBounds.add(bounds);
+                }
+                tilesArray[index++] = tile;
+            }
+            tilesAT.clear();                                            // Lets GC do its work.
+            /*
+             * Translate the tiles in such a way that the upper-left corner has the coordinates
+             * specified by (xLocation, yLocation). Adjust the tile affine transform consequently.
+             * After this block, tiles having the same subsampling will share the same immutable
+             * affine transform instance.
+             */
+            if (groupBounds != null) {
+                final int dx = xLocation - groupBounds.x;
+                final int dy = yLocation - groupBounds.y;
+                if (dx != 0 || dy != 0) {
+                    reference.translate(-dx, -dy);
+                    groupBounds.translate(dx, dy);
+                }
+                reference = new AffineTransform2D(reference);               // Make immutable.
+                final Map<Dimension,TileTranslation> pool = new HashMap<>();
+                for (final Tile tile : tilesArray) {
+                    final Dimension subsampling = tile.getSubsampling();
+                    TileTranslation translated = pool.get(subsampling);
+                    if (translated == null) {
+                        translated = new TileTranslation(subsampling, reference, dx, dy);
+                        pool.put(subsampling, translated);
+                    }
+                    translated.applyTo(tile);
+                }
+                results.put(new Tile(reference, groupBounds), tilesArray);
+            }
+        }
+        return results;
+    }
+
+    /**
+     * Sorts affine transform by increasing X scales in absolute value.
+     * For {@link #computePyramidLevels(Collection)} internal working only.
+     */
+    private static final Comparator<AffineTransform> X_COMPARATOR = new Comparator<AffineTransform>() {
+        @Override public int compare(final AffineTransform tr1, final AffineTransform tr2) {
+            return Double.compare(AffineTransforms2D.getScaleX0(tr1), AffineTransforms2D.getScaleX0(tr2));
+        }
+    };
+
+    /**
+     * Sorts affine transform by increasing Y scales in absolute value.
+     * For {@link #computePyramidLevels(Collection)} internal working only.
+     */
+    private static final Comparator<AffineTransform> Y_COMPARATOR = new Comparator<AffineTransform>() {
+        @Override public int compare(final AffineTransform tr1, final AffineTransform tr2) {
+            return Double.compare(AffineTransforms2D.getScaleY0(tr1), AffineTransforms2D.getScaleY0(tr2));
+        }
+    };
+
+    /**
+     * From a set of arbitrary affine transforms, computes pyramid levels that can be given to
+     * {@link Tile} constructors. This method tries to locate the affine transform with finest resolution.
+     * This is typically (but not always, depending on rotation or axis flip) the transform with smallest
+     * {@linkplain AffineTransform#getScaleX scale X} and {@linkplain AffineTransform#getScaleY scale Y}
+     * coefficients in absolute value. That transform is given a "resolution" of (1,1) and stored in an
+     * {@link IdentityHashMap}. Other transforms are stored in the same map with their resolution relative
+     * to the first one, or discarded if the relative resolution is not an integer. In the later case, the
+     * transforms that were discarded from the first pass will be put in a new map to be added as the second
+     * element in the returned list. A new pass is run, discarded transforms from the second pass are put in
+     * the third element of the list, <i>etc</i>.
+     *
+     * @param  gridToCRS  the <cite>grid to CRS</cite> affine transforms computed from the image to use in a pyramid.
+     *         The collection and the transform elements are not modified by this method (they may be modified by the
+     *         caller however).
+     * @return a subset of the given transforms with their relative resolution. This method typically returns one map,
+     *         but more could be returned if the scale ratio is not an integer for every transforms.
+     */
+    private static List<Map<AffineTransform,Dimension>> computePyramidLevels(final Collection<AffineTransform> gridToCRS) {
+        final List<Map<AffineTransform,Dimension>> results = new ArrayList<>(2);
+        /*
+         * First, compute the pyramid levels along the X axis. Transforms that we were unable
+         * to classify will be discarded from the first run and put in a subsequent run.
+         */
+        AffineTransform[] transforms = gridToCRS.toArray(new AffineTransform[gridToCRS.size()]);
+        Arrays.sort(transforms, X_COMPARATOR);
+        int length = transforms.length;
+        while (length != 0) {
+            final Map<AffineTransform,Dimension> result = new IdentityHashMap<>();
+            if (length <= (length = computePyramidLevels(transforms, length, result, false))) {
+                throw new AssertionError(length);               // Should always be decreasing.
+            }
+            results.add(result);
+        }
+        /*
+         * Next, compute the pyramid levels along the Y axis. If we fail to compute the
+         * pyramid level for some AffineTransform, they will be removed from the map.
+         * If a map became empty because of that, the whole map will be removed.
+         */
+        final Iterator<Map<AffineTransform,Dimension>> iterator = results.iterator();
+        while (iterator.hasNext()) {
+            final Map<AffineTransform,Dimension> result = iterator.next();
+            length = result.size();
+            transforms = result.keySet().toArray(transforms);
+            Arrays.sort(transforms, 0, length, Y_COMPARATOR);
+            length = computePyramidLevels(transforms, length, result, true);
+            while (--length >= 0) {
+                if (result.remove(transforms[length]) == null) {
+                    throw new AssertionError(length);
+                }
+            }
+            if (result.isEmpty()) {
+                iterator.remove();
+            }
+        }
+        return results;
+    }
+
+    /**
+     * Computes the pyramid level for the given affine transforms along the X or Y axis,
+     * and stores the result in the given map.
+     *
+     * @param  gridToCRS  the AffineTransform to analyze. This array <strong>must</strong>
+     *                    be sorted along the dimension specified by {@code isY}.
+     * @param  length     the number of valid entries in the {@code gridToCRS} array.
+     * @param  result     an initially empty map in which to store the results.
+     * @param  isY        {@code false} for analyzing the X axis, or {@code true} for the Y axis.
+     * @return the number of entries remaining in {@code gridToCRS}.
+     */
+    private static int computePyramidLevels(final AffineTransform[] gridToCRS, final int length,
+            final Map<AffineTransform,Dimension> result, final boolean isY)
+    {
+        int processing = 0;             // Index of the AffineTransform under process.
+        int remaining  = 0;             // Count of AffineTransforms that this method did not processed.
+        AffineTransform base;
+        double scale, shear;
+        boolean scaleIsNull, shearIsNull;
+        for (;;) {
+            if (processing >= length) {
+                return remaining;
+            }
+            base = gridToCRS[processing++];
+            if (isY) {
+                scale = base.getScaleY();
+                shear = base.getShearY();
+            } else {
+                scale = base.getScaleX();
+                shear = base.getShearX();
+            }
+            scaleIsNull = Math.abs(scale) < EPS;
+            shearIsNull = Math.abs(shear) < EPS;
+            if (!(scaleIsNull & shearIsNull)) break;
+            result.remove(base);
+        }
+        if (isY) {
+            // If we get a NullPointerException here, it would be a bug in the algorithm.
+            result.get(base).height = 1;
+        } else {
+            assert result.isEmpty() : result;
+            result.put(base, new Dimension(1,0));
+        }
+        /*
+         * From this point, consider 'base', 'scale', 'shear', 'scaleIsNull', 'shearIsNull' as final.
+         * They describe the AffineTransform with finest resolution along one axis (X or Y).
+         */
+        while (processing < length) {
+            final AffineTransform candidate = gridToCRS[processing++];
+            final double scale2, shear2;
+            if (isY) {
+                scale2 = candidate.getScaleY();
+                shear2 = candidate.getShearY();
+            } else {
+                scale2 = candidate.getScaleX();
+                shear2 = candidate.getShearX();
+            }
+            final int level;
+            if (scaleIsNull) {
+                if (!(Math.abs(scale2) < EPS)) {
+                    // Expected a null scale but was not.
+                    gridToCRS[remaining++] = candidate;
+                    continue;
+                }
+                level = level(shear2 / shear);
+            } else {
+                level = level(scale2 / scale);
+                if (shearIsNull ? !(Math.abs(shear2) < EPS) : (level(shear2 / shear) != level)) {
+                    // Expected (a null shear) : (the same pyramid level), but was not.
+                    gridToCRS[remaining++] = candidate;
+                    continue;
+                }
+            }
+            if (level == 0) {
+                // Not a pyramid level (the ratio is not an integer).
+                gridToCRS[remaining++] = candidate;
+                continue;
+            }
+            /*
+             * Stores the pyramid level either as the width or as the height, depending on the `isY` value.
+             * The map is assumed initially empty for the X values, and is assumed containing every required
+             * entries for the Y values.
+             */
+            if (isY) {
+                // If we get a NullPointerException here, it would be a bug in the algorithm.
+                result.get(candidate).height = level;
+            } else {
+                if (result.put(candidate, new Dimension(level,0)) != null) {
+                    throw new AssertionError(candidate);                                // Should never happen.
+                }
+            }
+        }
+        Arrays.fill(gridToCRS, remaining, length, null);
+        return remaining;
+    }
+
+    /**
+     * Computes the pyramid level from the ratio between two affine transform coefficients.
+     * If the ratio has been computed from {@code entry2.scaleX / entry1.scaleX}, then a return value of:
+     *
+     * <ul>
+     *   <li>1 means that both entries are at the same level.</li>
+     *   <li>2 means that the second entry has pixels twice as large as first entry.</li>
+     *   <li>3 means that the second entry has pixels three time larger than first entry.</li>
+     *   <li><i>etc...</i></li>
+     *   <li>A negative number means that the second entry has pixels smaller than first entry.</li>
+     *   <li>0 means that the ratio between entries is not an integer number.</li>
+     * </ul>
+     *
+     * @param  ratio  the ratio between affine transform coefficients.
+     * @return the pixel size (actually subsampling) relative to the smallest pixel, or 0 if it can not be computed.
+     *         If the ratio is between 0 and 1, then this method returns a negative number.
+     */
+    private static int level(double ratio) {
+        if (ratio > 0 && ratio < Double.POSITIVE_INFINITY) {
+            /*
+             * The 0.75 threshold could be anything between 0.5 and 1.
+             * We take a middle value for being safe regarding rounding errors.
+             */
+            final boolean inverse = (ratio < 0.75);
+            if (inverse) {
+                ratio = 1 / ratio;
+            }
+            final double integer = Math.rint(ratio);
+            if (integer < Integer.MAX_VALUE && Math.abs(ratio - integer) < EPS) {
+                /*
+                 * Found an integer ratio.
+                 * Inverse the sign (just as a matter of convention) if smaller than 1.
+                 */
+                int level = (int) integer;
+                if (inverse) {
+                    level = -level;
+                }
+                return level;
+            }
+        }
+        return 0;
+    }
+
+    /**
+     * Invoked when an I/O error occurred in {@link Tile#getSize()} or {@link Tile#getRegion()}.
+     * This error is non-fatal since {@code TileOrganizer} can fallback on calculation based
+     * on tile location only (without size).
+     *
+     * <p>The default implementation returns {@code false}, which instructs the caller to let
+     * the exception propagate.</p>
+     *
+     * @param  tile       the tile on which an error occurred.
+     * @param  exception  the error that occurred.
+     * @return {@code true} if the exception has been consumed, or {@code false} for re-throwing it.
+     */
+    protected boolean unavailableSize(final Tile tile, final IOException exception) {
+        return false;
+    }
+
+    /**
+     * Returns a string representation of the tiles contained in this object. Since this method is
+     * for debugging purpose, only the first tiles may be formatted in order to avoid consuming to
+     * much space in the debugger.
+     */
+    @Override
+    public String toString() {
+        return Tile.toString(tiles.values(), 400);
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/TileTranslation.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/TileTranslation.java
new file mode 100644
index 0000000..095ad78
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/TileTranslation.java
@@ -0,0 +1,75 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing.j2d;
+
+import java.awt.Dimension;
+import java.awt.geom.AffineTransform;
+
+
+/**
+ * An affine transform which is translated relative to an original transform.
+ * The translation terms are stored separately without modifying the transform.
+ * This class if for internal use by {@link TileOrganizer} only.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+final class TileTranslation {
+    /**
+     * The translated "grid to real world" transform, as an immutable instance.
+     */
+    private final AffineTransform gridToCRS;
+
+    /**
+     * The translation in "absolute units". This is the same units than for tiles at subsampling (1,1).
+     */
+    private final int dx, dy;
+
+    /**
+     * Creates a new translated transform. The translation is specified in "absolute units",
+     * i.e. in the same units than for tiles at subsampling (1,1).
+     *
+     * @param  subsampling  the {@linkplain Tile#getSubsampling() tile subsampling}.
+     * @param  reference    the "grid to real world" transform at subsampling (1,1).
+     * @param  dx           the translation along <var>x</var> axis in "absolute units".
+     * @param  dy           the translation along <var>y</var> axis in "absolute units".
+     */
+    TileTranslation(final Dimension subsampling, AffineTransform reference, int dx, int dy) {
+        this.dx = dx / subsampling.width;                           // It is okay to round toward zero.
+        this.dy = dy / subsampling.height;
+        dx %= subsampling.width;
+        dy %= subsampling.height;
+        reference = new AffineTransform(reference);
+        reference.scale(subsampling.width, subsampling.height);
+        reference.translate(dx, dy);                                // Correction for non-integer division of (dx,dy).
+        gridToCRS = new ImmutableAffineTransform(reference);
+    }
+
+    /**
+     * Applies the translation and the new "grid to CRS" transform on the given tile.
+     *
+     * @param  tile  the tile on which to apply the translation.
+     */
+    final void applyTo(final Tile tile) {
+        synchronized (tile) {
+            tile.translate(dx, dy);
+            tile.setGridToCRS(gridToCRS);
+        }
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/package-info.java
index 1df2079..3d4904d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/package-info.java
@@ -26,7 +26,7 @@
  * may change in incompatible ways in any future version without notice.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see org.apache.sis.internal.feature.j2d
  *
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AbstractProvider.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AbstractProvider.java
index d3b067f..aa77cd5 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AbstractProvider.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/AbstractProvider.java
@@ -38,13 +38,15 @@ import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactor
 import org.apache.sis.util.resources.Vocabulary;
 import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.Workaround;
+import org.apache.sis.util.logging.Logging;
+import org.apache.sis.internal.system.Loggers;
 
 
 /**
  * Base class for all providers defined in this package.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.1
  * @since   0.6
  * @module
  */
@@ -262,4 +264,16 @@ public abstract class AbstractProvider extends DefaultOperationMethod implements
     public boolean isInvertible() {
         return false;
     }
+
+    /**
+     * Convenience method for reporting a non-fatal error at transform construction time.
+     * This method assumes that the error occurred (indirectly) during execution of
+     * {@link #createMathTransform(MathTransformFactory, ParameterValueGroup)}.
+     *
+     * @param  caller  the provider class in which the error occurred.
+     * @param  e       the error that occurred.
+     */
+    static void recoverableException(final Class<? extends AbstractProvider> caller, Exception e) {
+        Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION), caller, "createMathTransform", e);
+    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java
index 3527244..139760e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java
@@ -19,6 +19,7 @@ package org.apache.sis.internal.referencing.provider;
 import java.util.Arrays;
 import javax.measure.Quantity;
 import org.apache.sis.math.DecimalFunctions;
+import org.apache.sis.internal.util.Numerics;
 
 
 /**
@@ -28,7 +29,7 @@ import org.apache.sis.math.DecimalFunctions;
  * 5 digits in base 10 in ASCII files.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @param <C>  dimension of the coordinate unit (usually {@link javax.measure.quantity.Angle}).
  * @param <T>  dimension of the translation unit (usually {@link javax.measure.quantity.Angle}
@@ -44,13 +45,15 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
     private static final long serialVersionUID = 4847888093457104917L;
 
     /**
-     * Maximal grid index along the <var>y</var> axis.
+     * Maximal grid index along the <var>y</var> axis, inclusive.
      * This is the number of grid cells minus 2.
      */
     private final int ymax;
 
     /**
      * An "average" value for the offset in each dimension.
+     *
+     * @see #getCellMean(int)
      */
     private final double[] averages;
 
@@ -122,6 +125,9 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
 
     /**
      * Returns a new grid with the same geometry than this grid but different data arrays.
+     * This method is invoked by {@link #useSharedData()} when it detects that a newly created
+     * grid uses the same data than an existing grid. The {@code other} object is the old grid,
+     * so we can share existing data.
      */
     @Override
     protected final DatumShiftGridFile<C,T> setData(final Object[] other) {
@@ -253,8 +259,7 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
     public boolean equals(final Object other) {
         if (super.equals(other)) {
             final DatumShiftGridCompressed<?,?> that = (DatumShiftGridCompressed<?,?>) other;
-            return Double.doubleToLongBits(scale) == Double.doubleToLongBits(that.scale)
-                   && Arrays.equals(averages, that.averages);
+            return Numerics.equals(scale, that.scale) && Arrays.equals(averages, that.averages);
         }
         return false;
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
index ebc74ee..f77bba7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
@@ -17,19 +17,38 @@
 package org.apache.sis.internal.referencing.provider;
 
 import java.util.Arrays;
+import java.util.Collection;
+import java.util.LinkedHashMap;
+import java.util.Locale;
+import java.util.Map;
+import java.util.logging.Level;
 import java.lang.reflect.Array;
 import java.nio.file.Path;
 import javax.measure.Unit;
 import javax.measure.Quantity;
+import javax.measure.quantity.Angle;
+import org.opengis.geometry.Envelope;
+import org.opengis.util.FactoryException;
 import org.opengis.parameter.ParameterDescriptor;
 import org.opengis.parameter.ParameterDescriptorGroup;
 import org.opengis.parameter.GeneralParameterDescriptor;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.measure.Units;
 import org.apache.sis.math.DecimalFunctions;
 import org.apache.sis.util.collection.Cache;
+import org.apache.sis.util.collection.TreeTable;
+import org.apache.sis.util.collection.TableColumn;
+import org.apache.sis.util.collection.DefaultTreeTable;
+import org.apache.sis.util.collection.Containers;
+import org.apache.sis.util.resources.Errors;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.referencing.datum.DatumShiftGrid;
 import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
 
 
 /**
@@ -41,7 +60,7 @@ import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
  * sharing data and for {@link #equals(Object)} and {@link #hashCode()} implementations.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @param <C>  dimension of the coordinate unit (usually {@link javax.measure.quantity.Angle}).
  * @param <T>  dimension of the translation unit (usually {@link javax.measure.quantity.Angle}
@@ -66,8 +85,12 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     static final Cache<Object, DatumShiftGridFile<?,?>> CACHE = new Cache<Object, DatumShiftGridFile<?,?>>(4, 32*1024, true) {
         @Override protected int cost(final DatumShiftGridFile<?,?> grid) {
             int p = 1;
-            for (final Object array : grid.getData()) {
-                p *= Array.getLength(array);
+            for (final Object data : grid.getData()) {
+                if (data instanceof DatumShiftGridFile<?,?>) {
+                    p += cost((DatumShiftGridFile<?,?>) data);          // When `grid` is a DatumShiftGridGroup.
+                } else {
+                    p *= Array.getLength(data);                         // short[], float[] or double[].
+                }
             }
             return p;
         }
@@ -95,14 +118,32 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
      * The best translation accuracy that we can expect from this file.
      * The unit of measurement depends on {@link #isCellValueRatio()}.
      *
-     * <p>This field is initialized to {@link Double#NaN}. It is loader responsibility
-     * to assign a value to this field after {@code DatumShiftGridFile} construction.</p>
+     * <p>This field is initialized to zero. It is loader responsibility to assign
+     * a value to this field after {@code DatumShiftGridFile} construction.</p>
      *
      * @see #getCellPrecision()
      */
     double accuracy;
 
     /**
+     * The sub-grids, or {@code null} if none. The domain of validity of each sub-grid should be contained
+     * in the domain of validity of this grid. Children do not change the way this {@code DatumShiftGrid}
+     * performs its calculation; this list is used only at the time of building {@link MathTransform} tree.
+     *
+     * <div class="note"><b>Design note:</b>
+     * we do not provide sub-grids functionality in the {@link DatumShiftGrid} parent class because
+     * the {@link MathTransform} tree will depend on assumptions about {@link #getCoordinateToGrid()},
+     * in particular that it contains only translations and scales (no rotation, no shear).
+     * Those assumptions are enforced by the {@link DatumShiftGridFile} constructor.</div>
+     *
+     * This field has protected access for usage by {@link DatumShiftGridGroup} subclass only.
+     * No access to this field should be done except by subclasses.
+     *
+     * @see #setSubGrids(Collection)
+     */
+    protected DatumShiftGridFile<C,T>[] subgrids;
+
+    /**
      * Creates a new datum shift grid for the given grid geometry.
      * The actual offset values need to be provided by subclasses.
      *
@@ -132,11 +173,11 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
         this.descriptor = descriptor;
         this.files      = files;
         this.nx         = nx;
-        this.accuracy   = Double.NaN;
     }
 
     /**
      * Creates a new datum shift grid with the same grid geometry than the given grid.
+     * This is used by {@link DatumShiftGridCompressed} for replacing a grid by another one.
      *
      * @param  other  the other datum shift grid from which to copy the grid geometry.
      */
@@ -146,6 +187,82 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
         files      = other.files;
         nx         = other.nx;
         accuracy   = other.accuracy;
+        subgrids   = other.subgrids;
+    }
+
+    /**
+     * Creates a new datum shift grid with the same configuration than the given grid,
+     * except the size and transform which are set to the given values.
+     * The {@link #accuracy} is initialized to zero and should be updated by the caller.
+     *
+     * @param other             the other datum shift grid from which to copy parameters.
+     * @param coordinateToGrid  conversion from the "real world" coordinates to grid indices including fractional parts.
+     * @param nx                number of cells along the <var>x</var> axis in the grid.
+     * @param ny                number of cells along the <var>y</var> axis in the grid.
+     */
+    DatumShiftGridFile(final DatumShiftGridFile<C,T> other,
+                       final AffineTransform2D coordinateToGrid, final int nx, final int ny)
+    {
+        super(other.getCoordinateUnit(), coordinateToGrid, new int[] {nx, ny},
+              other.isCellValueRatio(), other.getTranslationUnit());
+        descriptor = other.descriptor;
+        files      = other.files;
+        this.nx    = nx;
+        // Accuracy to be set by caller. Initial value needs to be zero.
+    }
+
+    /**
+     * Sets the sub-grids that are direct children of this grid.
+     * This method can be invoked only once.
+     */
+    @SuppressWarnings({"rawtypes", "unchecked"})
+    final void setSubGrids(final Collection<DatumShiftGridFile<C,T>> children) {
+        if (subgrids != null) throw new IllegalStateException();
+        subgrids = children.toArray(new DatumShiftGridFile[children.size()]);
+    }
+
+    /**
+     * Returns the number of grids, including this grid and all sub-grids counted recursively.
+     * This is used for information purpose only.
+     *
+     * @see #toTree(TreeTable.Node)
+     */
+    private int getGridCount() {
+        int n = 1;
+        if (subgrids != null) {
+            for (final DatumShiftGridFile<C,T> subgrid : subgrids) {
+                n += subgrid.getGridCount();
+            }
+        }
+        return n;
+    }
+
+    /**
+     * Returns a string representation of this grid for debugging purpose.
+     * If this grid has children, then it will be formatted as a tree.
+     */
+    @Override
+    public final String toString() {
+        if (subgrids == null) {
+            return super.toString();
+        }
+        final TreeTable tree = new DefaultTreeTable(TableColumn.NAME);
+        toTree(tree.getRoot());
+        return tree.toString();
+    }
+
+    /**
+     * Formats this grid as a tree with its children.
+     */
+    private void toTree(final TreeTable.Node branch) {
+        String label = super.toString();
+        if (subgrids != null) {
+            label = label + " (" + getGridCount() + " grids)";
+            for (final DatumShiftGridFile<C,T> subgrid : subgrids) {
+                subgrid.toTree(branch.newChild());
+            }
+        }
+        branch.setValue(TableColumn.NAME, label);
     }
 
     /**
@@ -203,6 +320,37 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     protected abstract Object[] getData();
 
     /**
+     * Returns {@code true} if the given object is a grid containing the same data than this grid.
+     * This method compares the data provided by {@link #getData()}.
+     *
+     * @param  other  the other object to compare with this datum shift grid.
+     * @return {@code true} if the given object is non-null, of the same class than this {@code DatumShiftGrid}
+     *         and contains the same data.
+     */
+    @Override
+    public boolean equals(final Object other) {
+        if (other == this) {                        // Optimization for a common case.
+            return true;
+        }
+        if (super.equals(other)) {
+            final DatumShiftGridFile<?,?> that = (DatumShiftGridFile<?,?>) other;
+            return Arrays.equals(files, that.files) && Arrays.deepEquals(getData(), that.getData());
+        }
+        return false;
+    }
+
+    /**
+     * Returns a hash code value for this datum shift grid. The hash code is based on metadata
+     * such as filename, but not on {@link #getData()} for performance reason.
+     *
+     * @return a hash code based on metadata.
+     */
+    @Override
+    public int hashCode() {
+        return super.hashCode() + Arrays.hashCode(files);
+    }
+
+    /**
      * Suggests a precision for the translation values in this grid.
      * This information is used for deciding when to stop iterations in inverse transformations.
      * The default implementation returns the {@linkplain #accuracy} divided by an arbitrary value.
@@ -245,32 +393,38 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     }
 
     /**
-     * Returns {@code true} if the given object is a grid containing the same data than this grid.
+     * Creates a transformation between two geodetic CRS, including the sub-grid transforms.
+     * If the given grid has no sub-grid, then this method is equivalent to a direct call to
+     * {@link InterpolatedTransform#createGeodeticTransformation(MathTransformFactory, DatumShiftGrid)}.
      *
-     * @param  other  the other object to compare with this datum shift grid.
-     * @return {@code true} if the given object is non-null, of the same class than this {@code DatumShiftGrid}
-     *         and contains the same data.
+     * @param  provider  the provider which is creating a transform.
+     * @param  factory   the factory to use for creating the transform.
+     * @param  grid      the grid of datum shifts from source to target datum.
+     * @return the transformation between geodetic coordinates.
+     * @throws FactoryException if an error occurred while creating a transform.
+     *
+     * @see InterpolatedTransform#createGeodeticTransformation(MathTransformFactory, DatumShiftGrid)
      */
-    @Override
-    public boolean equals(final Object other) {
-        if (other == this) {                        // Optimization for a common case.
-            return true;
+    public static MathTransform createGeodeticTransformation(final Class<? extends AbstractProvider> provider,
+            final MathTransformFactory factory, final DatumShiftGridFile<Angle,Angle> grid) throws FactoryException
+    {
+        MathTransform global = InterpolatedTransform.createGeodeticTransformation(factory, grid);
+        final DatumShiftGridFile<Angle,Angle>[] subgrids = grid.subgrids;
+        if (subgrids == null) {
+            return global;
         }
-        if (super.equals(other)) {
-            final DatumShiftGridFile<?,?> that = (DatumShiftGridFile<?,?>) other;
-            return Arrays.equals(files, that.files) && Arrays.deepEquals(getData(), that.getData());
+        final Map<Envelope,MathTransform> specializations = new LinkedHashMap<>(Containers.hashMapCapacity(subgrids.length));
+        for (final DatumShiftGridFile<Angle,Angle> sg : subgrids) try {
+            final Envelope domain = sg.getDomainOfValidity(Units.DEGREE);
+            final MathTransform st = createGeodeticTransformation(provider, factory, sg);
+            if (specializations.putIfAbsent(domain, st) != null) {
+                DatumShiftGridLoader.log(provider, Errors.getResources((Locale) null)
+                        .getLogRecord(Level.FINE, Errors.Keys.DuplicatedElement_1, domain));
+            }
+        } catch (TransformException e) {
+            throw new FactoryException(e);
         }
-        return false;
-    }
-
-    /**
-     * Returns a hash code value for this datum shift grid.
-     *
-     * @return {@inheritDoc}
-     */
-    @Override
-    public int hashCode() {
-        return super.hashCode() + Arrays.hashCode(files);
+        return MathTransforms.specialize(global, specializations);
     }
 
 
@@ -308,6 +462,8 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
         /**
          * Creates a new datum shift grid with the given grid geometry, filename and number of shift dimensions.
          * All {@code double} values given to this constructor will be converted from degrees to radians.
+         *
+         * @param  dim  number of dimensions of translation vectors.
          */
         Float(final int dim,
               final Unit<C> coordinateUnit,
@@ -320,11 +476,7 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
               final Path... files) throws NoninvertibleTransformException
         {
             super(coordinateUnit, translationUnit, isCellValueRatio, x0, y0, Δx, Δy, nx, ny, descriptor, files);
-            offsets = new float[dim][];
-            final int size = Math.multiplyExact(nx, ny);
-            for (int i=0; i<dim; i++) {
-                Arrays.fill(offsets[i] = new float[size], java.lang.Float.NaN);
-            }
+            offsets = new float[dim][Math.multiplyExact(nx, ny)];
         }
 
         /**
@@ -337,6 +489,9 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
 
         /**
          * Returns a new grid with the same geometry than this grid but different data arrays.
+         * This method is invoked by {@link #useSharedData()} when it detects that a newly created
+         * grid uses the same data than an existing grid. The {@code other} object is the old grid,
+         * so we can share existing data.
          */
         @Override
         protected final DatumShiftGridFile<C,T> setData(final Object[] other) {
@@ -377,5 +532,108 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
         public final double getCellValue(final int dim, final int gridX, final int gridY) {
             return DecimalFunctions.floatToDouble(offsets[dim][gridX + gridY*nx]);
         }
+
+        /**
+         * Returns the average translation parameters from source to target.
+         * There is no need to use double-double arithmetic here since all data have only single precision.
+         *
+         * @param  dim  the dimension for which to get an average value.
+         * @return a value close to the average for the given dimension.
+         */
+        @Override
+        public double getCellMean(final int dim) {
+            final float[] data = offsets[dim];
+            double sum = 0;
+            for (final float value : data) {
+                sum += value;
+            }
+            return sum / data.length;
+        }
+    }
+
+
+
+
+    /**
+     * An implementation of {@link DatumShiftGridFile} which stores the offset values in {@code double[]} arrays.
+     * See {@link DatumShiftGridFile.Float} for more information (most comments apply to this class as well).
+     *
+     * @author  Martin Desruisseaux (Geomatys)
+     * @version 1.1
+     * @since   1.1
+     * @module
+     */
+    static final class Double<C extends Quantity<C>, T extends Quantity<T>> extends DatumShiftGridFile<C,T> {
+        /**
+         * Serial number for inter-operability with different versions.
+         */
+        private static final long serialVersionUID = 3999271636016362364L;
+
+        /**
+         * The translation values. See {@link DatumShiftGridFile.Float#offsets} for more documentation.
+         */
+        final double[][] offsets;
+
+        /**
+         * Creates a new datum shift grid with the given grid geometry, filename and number of shift dimensions.
+         * All {@code double} values given to this constructor will be converted from degrees to radians.
+         */
+        Double(final int dim,
+               final Unit<C> coordinateUnit,
+               final Unit<T> translationUnit,
+               final boolean isCellValueRatio,
+               final double x0, final double y0,
+               final double Δx, final double Δy,
+               final int    nx, final int    ny,
+               final ParameterDescriptorGroup descriptor,
+               final Path... files) throws NoninvertibleTransformException
+        {
+            super(coordinateUnit, translationUnit, isCellValueRatio, x0, y0, Δx, Δy, nx, ny, descriptor, files);
+            offsets = new double[dim][Math.multiplyExact(nx, ny)];
+        }
+
+        /**
+         * Creates a new grid of the same geometry than the given grid but using a different data array.
+         */
+        private Double(final DatumShiftGridFile<C,T> grid, final double[][] offsets) {
+            super(grid);
+            this.offsets = offsets;
+        }
+
+        /**
+         * Returns a new grid with the same geometry than this grid but different data arrays.
+         * See {@link DatumShiftGridFile.Float#setData(Object[])} for more documentation.
+         */
+        @Override
+        protected final DatumShiftGridFile<C,T> setData(final Object[] other) {
+            return new Double<>(this, (double[][]) other);
+        }
+
+        /**
+         * Returns direct references (not cloned) to the data arrays.
+         * See {@link DatumShiftGridFile.Float#getData()} for more documentation.
+         */
+        @Override
+        @SuppressWarnings("ReturnOfCollectionOrArrayField")
+        protected final Object[] getData() {
+            return offsets;
+        }
+
+        /**
+         * Returns the number of shift dimensions.
+         */
+        @Override
+        public final int getTranslationDimensions() {
+            return offsets.length;
+        }
+
+        /**
+         * Returns the cell value at the given dimension and grid index.
+         * See {@link DatumShiftGridFile.Float#getCellValue(int, int, int)} for more documentation.
+         */
+        @Override
+        public final double getCellValue(final int dim, final int gridX, final int gridY) {
+            return offsets[dim][gridX + gridY*nx];
+        }
     }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
new file mode 100644
index 0000000..df2cefa
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
@@ -0,0 +1,306 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing.provider;
+
+import java.util.Map;
+import java.util.List;
+import java.nio.file.Path;
+import java.io.IOException;
+import java.awt.Dimension;
+import java.awt.Rectangle;
+import java.awt.geom.AffineTransform;
+import java.util.LinkedHashMap;
+import javax.measure.Quantity;
+import org.opengis.util.FactoryException;
+import org.opengis.referencing.operation.NoninvertibleTransformException;
+import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
+import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+import org.apache.sis.internal.referencing.j2d.IntervalRectangle;
+import org.apache.sis.internal.referencing.j2d.TileOrganizer;
+import org.apache.sis.internal.referencing.j2d.Tile;
+import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.internal.util.CollectionsExt;
+import org.apache.sis.util.collection.Containers;
+
+
+/**
+ * A group of datum shift grids. This is used when a NTv2 file contains more than one grid with no common parent.
+ * This class creates a synthetic parent which always delegates its work to a child (as opposed to more classical
+ * transform trees where the parent can do some work if no child can). Coordinate transformations will be applied
+ * as below:
+ *
+ * <ol>
+ *   <li>{@link org.apache.sis.referencing.operation.transform.SpecializableTransform} will try to locate the
+ *       most appropriate grid for given coordinates. This is the class where to put our optimization efforts,
+ *       for example by checking the last used grid before to check all other grids.</li>
+ *   <li>Only if {@code SpecializableTransform} did not found a better transform, it will fallback on a transform
+ *       backed by this {@code DatumShiftGridGroup}. In such case, {@link InterpolatedTransform} will perform its
+ *       calculation by invoking {@link #interpolateInCell(double, double, double[])}. That method tries again to
+ *       locate the best grid, but performance is less important there since that method is only a fallback.</li>
+ *   <li>The default {@link DatumShiftGridFile#interpolateInCell(double, double, double[])} implementation invokes
+ *       {@link #getCellValue(int, int, int)}. We provide that method for consistency, but it should not be invoked
+ *       since we overrode {@link #interpolateInCell(double, double, double[])}.</li>
+ * </ol>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+final class DatumShiftGridGroup<C extends Quantity<C>, T extends Quantity<T>> extends DatumShiftGridFile<C,T> {
+    /**
+     * For cross-version compatibility.
+     */
+    private static final long serialVersionUID = -1602724619897451422L;
+
+    /**
+     * The bounds of a sub-grid, together with the subsampling level compared to the grid having the finest resolution.
+     * All values in this class are integers, but nevertheless stored as {@code double} for avoiding to cast them every
+     * time {@link DatumShiftGridGroup#interpolateInCell(double, double, double[])} is executed.
+     */
+    @SuppressWarnings("CloneableImplementsClone")
+    private static final class Region extends IntervalRectangle {
+        /** Subsampling compared to the grid having finest resolution. */
+        private final double sx, sy;
+
+        /** Creates a new instance from the given {@link TileOrganizer} result. */
+        Region(final Tile tile) throws IOException {
+            final Rectangle r = tile.getAbsoluteRegion();       // In units of the grid having finest resolution.
+            final Dimension s = tile.getSubsampling();
+            xmin = r.getMinX();
+            xmax = r.getMaxX();
+            ymin = r.getMinY();
+            ymax = r.getMaxY();
+            sx   = s.width;
+            sy   = s.height;
+        }
+
+        /** Converts a coordinate from the parent grid to this grid. */
+        final double x(final double p) {return (p - xmin) / sx;}
+        final double y(final double p) {return (p - ymin) / sy;}
+
+        /** Returns the subsampling (compared to the grid having finest resolution) in the specified dimension. */
+        final double relativeCellSize(final int dim) {
+            switch (dim) {
+                case 0:  return sx;
+                case 1:  return sy;
+                default: throw new IndexOutOfBoundsException();
+            }
+        }
+    }
+
+    /**
+     * For each {@code subgrids[i]}, {@code regions[i]} is the range of indices valid for that grid.
+     * This array will be used only as a fallback if {@code SpecializableTransform} has not been able
+     * to find the sub-grid itself. Since it should be rarely used, we do not bother using a R-Tree.
+     */
+    private final Region[] regions;
+
+    /**
+     * Creates a new group for the given list of sub-grids. That list shall contain at least 2 elements.
+     * The first sub-grid is taken as a template for setting parameter values such as filename (all list
+     * elements should declare the same filename parameters, so the selected element should not matter).
+     *
+     * @param  tiles      the tiles computed by {@link TileOrganizer}.
+     * @param  grids      sub-grids associated to tiles computed by {@link TileOrganizer}.
+     * @param  gridToCRS  conversion from grid indices to "real world" coordinates.
+     * @param  gridSize   number of cells along the <var>x</var> and <var>y</var> axes in the grid.
+     * @throws IOException declared because {@link Tile#getRegion()} declares it, but should not happen.
+     */
+    @SuppressWarnings({"rawtypes", "unchecked"})                        // For generic array creation.
+    private DatumShiftGridGroup(final Tile[] tiles,
+                                final Map<Tile,DatumShiftGridFile<C,T>> grids,
+                                final AffineTransform2D gridToCRS,
+                                final Dimension gridSize)
+            throws IOException, NoninvertibleTransformException
+    {
+        super(grids.get(tiles[0]), gridToCRS.inverse(), gridSize.width, gridSize.height);
+        final int n = grids.size();
+        regions  = new Region[n];
+        subgrids = new DatumShiftGridFile[n];
+        for (int i=0; i<n; i++) {
+            final Tile tile = tiles[i];
+            final DatumShiftGridFile<C,T> grid = grids.get(tile);
+            regions [i] = new Region(tile);
+            subgrids[i] = grid;
+            if (grid.accuracy > accuracy) {
+                accuracy = grid.accuracy;           // Conservatively set accuracy to the largest value.
+            }
+        }
+    }
+
+    /**
+     * Puts the given sub-grid in a group. This method infers itself what would be the size
+     * of a grid containing all given sub-grids.
+     *
+     * @param  file      filename to report in case of error.
+     * @param  subgrids  the sub-grids to put under a common root.
+     * @throws FactoryException if the sub-grid can not be combined in a single mosaic or pyramid.
+     * @throws IOException declared because {@link Tile#getRegion()} declares it, but should not happen.
+     */
+    static <C extends Quantity<C>, T extends Quantity<T>> DatumShiftGridGroup<C,T> create(
+            final Path file, final List<DatumShiftGridFile<C,T>> subgrids)
+            throws IOException, FactoryException, NoninvertibleTransformException
+    {
+        final TileOrganizer mosaic = new TileOrganizer(null);
+        final Map<Tile,DatumShiftGridFile<C,T>> grids = new LinkedHashMap<>(Containers.hashMapCapacity(subgrids.size()));
+        for (final DatumShiftGridFile<C,T> grid : subgrids) {
+            final int[] size = grid.getGridSize();
+            final Tile  tile = new Tile(new Rectangle(size[0], size[1]),
+                    (AffineTransform) grid.getCoordinateToGrid().inverse());
+            /*
+             * Assertions below would fail if the tile has already been processed by TileOrganizer,
+             * or if it duplicates another tile. Since we created that tile just above, a failure
+             * would be a bug in Tile or TileOrganizer.
+             */
+            if (!mosaic.add(tile) || grids.put(tile, grid) != null) {
+                throw new AssertionError(tile);
+            }
+        }
+        /*
+         * After processing by TileOrganizer, we should have only one group of tiles. If we have more groups,
+         * it would mean that the cell size of the grid having larger cells is not a multiple of cell size of
+         * the grid having smallest cells, or that cell indices in some grids, when expressed in units of the
+         * smallest cells, would be fractional numbers. It should not happen in a NTv2 compliant file.
+         */
+        final Map.Entry<Tile,Tile[]> result = CollectionsExt.singletonOrNull(mosaic.tiles().entrySet());
+        if (result == null) {
+            throw new FactoryException(Resources.format(Resources.Keys.MisalignedDatumShiftGrid_1, file));
+        }
+        final Tile global = result.getKey();
+        return new DatumShiftGridGroup<>(result.getValue(), grids, global.getGridToCRS(), global.getSize());
+    }
+
+    /**
+     * Creates a new grid sharing the same data than an existing grid.
+     * This constructor is for {@link #setData(Object[])} usage only.
+     */
+    private DatumShiftGridGroup(final DatumShiftGridGroup<C,T> other, final DatumShiftGridFile<C,T>[] data) {
+        super(other);
+        subgrids = data;
+        regions  = other.regions;
+    }
+
+    /**
+     * Returns a new grid with the same geometry than this grid but different data arrays.
+     * This method is invoked by {@link #useSharedData()} when it detects that a newly created
+     * grid uses the same data than an existing grid. The {@code other} object is the old grid,
+     * so we can share existing data.
+     */
+    @Override
+    @SuppressWarnings("unchecked")
+    protected final DatumShiftGridFile<C,T> setData(final Object[] other) {
+        return new DatumShiftGridGroup<>(this, (DatumShiftGridFile<C,T>[]) other);
+    }
+
+    /**
+     * Returns direct references (not cloned) to the data arrays. This method is for cache management,
+     * {@link #equals(Object)} and {@link #hashCode()} implementations only and should not be invoked
+     * in other context.
+     */
+    @Override
+    @SuppressWarnings("ReturnOfCollectionOrArrayField")
+    protected Object[] getData() {
+        return subgrids;
+    }
+
+    /**
+     * Returns the number of dimensions of the translation vectors interpolated by this datum shift grid.
+     * This implementation takes the first sub-grid as a template. The choice of the grid does not matter
+     * since all grids have the same number of target dimensions.
+     */
+    @Override
+    public int getTranslationDimensions() {
+        return subgrids[0].getTranslationDimensions();
+    }
+
+    /**
+     * Returns the translation stored at the given two-dimensional grid indices for the given dimension.
+     * This method is defined for consistency with {@link #interpolateInCell(double, double, double[])}
+     * but should never be invoked. The {@link InterpolatedTransform} class will rather invoke the
+     * {@code interpolateInCell(…)} method for efficiency.
+     *
+     * <p>Caller must ensure that all arguments given to this method are in their expected ranges.
+     * The behavior of this method is undefined if any argument value is out-of-range.</p>
+     *
+     * @param  dim    the dimension of the translation vector component to get.
+     * @param  gridX  the grid index on the <var>x</var> axis, from 0 inclusive to {@code gridSize[0]} exclusive.
+     * @param  gridY  the grid index on the <var>y</var> axis, from 0 inclusive to {@code gridSize[1]} exclusive.
+     * @return the translation for the given dimension in the grid cell at the given index.
+     */
+    @Override
+    public double getCellValue(final int dim, final int gridX, final int gridY) {
+        for (int i=0; i<regions.length; i++) {
+            final Region r = regions[i];
+            if (r.containsInclusive(gridX, gridY)) {
+                double shift = subgrids[i].getCellValue(dim,
+                        Math.toIntExact(Math.round(r.x(gridX))),
+                        Math.toIntExact(Math.round(r.y(gridY))));
+                /*
+                 * If the translations have been divided by the cell size, we may need to compensate.
+                 * The size of the cells of the grid used below may be bigger than the cells of this
+                 * pseudo-grid.
+                 */
+                if (isCellValueRatio()) {
+                    shift *= r.relativeCellSize(dim);
+                }
+                return shift;
+            }
+        }
+        /*
+         * May be in the valid range of this DatumShiftGridGroup but not in the range of a subgrid.
+         * This situation may happen if there is holes in the data coverage provided by subgrids.
+         */
+        throw new IllegalArgumentException();
+    }
+
+    /**
+     * Interpolates the translation to apply for the given two-dimensional grid indices.
+     * During forward coordinate transformations, this method is invoked only if {@code SpecializableTransform}
+     * has been unable to use directly one of the child transforms — so performance is not the priority in that
+     * situation. During inverse transformations, this method is invoked for estimating an initial position before
+     * iterative refinements. The given point may be outside all sub-grids (otherwise {@code SpecializableTransform}
+     * would have done the work itself at least in the forward transformation case). Consequently searching a sub-grid
+     * containing the given point is not sufficient; we need to search for the nearest grid even if the point is outside.
+     *
+     * @param  gridX   first grid coordinate of the point for which to get the translation.
+     * @param  gridY   second grid coordinate of the point for which to get the translation.
+     * @param  vector  a pre-allocated array where to write the translation vector.
+     */
+    @Override
+    public void interpolateInCell(final double gridX, final double gridY, final double[] vector) {
+        int ni = 0;
+        Region nearest = regions[ni];
+        double dmin = nearest.distanceSquared(gridX, gridY);
+        for (int i=1; i<regions.length; i++) {
+            final Region r = regions[i];
+            final double d = r.distanceSquared(gridX, gridY);
+            if (d < dmin) {
+                dmin     = d;
+                nearest  = r;
+                ni       = i;
+                if (d == 0) break;
+            }
+        }
+        subgrids[ni].interpolateInCell(nearest.x(gridX), nearest.y(gridY), vector);
+        if (isCellValueRatio()) {
+            for (int dim=0; dim < INTERPOLATED_DIMENSIONS; dim++) {
+                vector[dim] *= nearest.relativeCellSize(dim);
+            }
+        }
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridLoader.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridLoader.java
index f701d51..ba35251 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridLoader.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridLoader.java
@@ -143,8 +143,18 @@ class DatumShiftGridLoader {
      *                 the source method will be set to {@code "createMathTransform"}.
      * @param  file    the grid file, as a {@link String} or a {@link Path}.
      */
-    static void log(final Class<?> caller, final Object file) {
-        final LogRecord record = Resources.forLocale(null).getLogRecord(Level.FINE, Resources.Keys.LoadingDatumShiftFile_1, file);
+    static void startLoading(final Class<?> caller, final Object file) {
+        log(caller, Resources.forLocale(null).getLogRecord(Level.FINE, Resources.Keys.LoadingDatumShiftFile_1, file));
+    }
+
+    /**
+     * Logs the given record.
+     *
+     * @param  caller  the provider to logs as the source class.
+     *                 the source method will be set to {@code "createMathTransform"}.
+     * @param record   the record to log.
+     */
+    static void log(final Class<?> caller, final LogRecord record) {
         record.setLoggerName(Loggers.COORDINATE_OPERATION);
         Logging.log(caller, "createMathTransform", record);
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolation.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolation.java
index 453fa3c..4e4d795 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolation.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolation.java
@@ -16,6 +16,7 @@
  */
 package org.apache.sis.internal.referencing.provider;
 
+import java.util.Arrays;
 import java.util.Collections;
 import java.util.Locale;
 import java.util.NoSuchElementException;
@@ -340,7 +341,7 @@ public class FranceGeocentricInterpolation extends GeodeticOperation {
                 grid = handler.peek();
                 if (grid == null) {
                     try (BufferedReader in = Files.newBufferedReader(resolved)) {
-                        DatumShiftGridLoader.log(FranceGeocentricInterpolation.class, file);
+                        DatumShiftGridLoader.startLoading(FranceGeocentricInterpolation.class, file);
                         final DatumShiftGridFile.Float<Angle,Length> g = load(in, file);
                         grid = DatumShiftGridCompressed.compress(g, averages, scale);
                     } catch (IOException | NoninvertibleTransformException | RuntimeException e) {
@@ -426,6 +427,10 @@ public class FranceGeocentricInterpolation extends GeodeticOperation {
                             grid = new DatumShiftGridFile.Float<>(3,
                                     Units.DEGREE, Units.METRE, false,
                                     x0, y0, Δx, Δy, nx, ny, PARAMETERS, file);
+                            grid.accuracy = Double.NaN;
+                            for (final float[] data : grid.offsets) {
+                                Arrays.fill(data, Float.NaN);
+                            }
                         }
                         break;
                     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/GeocentricAffine.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/GeocentricAffine.java
index 501236a..e2b3c06 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/GeocentricAffine.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/GeocentricAffine.java
@@ -341,7 +341,7 @@ public abstract class GeocentricAffine extends GeodeticOperation {
         if (datumShift != null) try {
             parameters.setPositionVectorTransformation(datumShift, BURSAWOLF_TOLERANCE);
         } catch (IllegalArgumentException e) {
-            log(Loggers.COORDINATE_OPERATION, "createParameters", e);
+            recoverableException(GeocentricAffine.class, e);
             return null;
         } else {
             /*
@@ -435,7 +435,7 @@ public abstract class GeocentricAffine extends GeodeticOperation {
                          * Should not occur, except sometime on inverse transform of relatively complex datum shifts
                          * (more than just translation terms). We can fallback on formatting the full matrix.
                          */
-                        log(Loggers.WKT, "asDatumShift", e);
+                        Logging.recoverableException(Logging.getLogger(Loggers.WKT), GeocentricAffine.class, "asDatumShift", e);
                         continue;
                     }
                     final boolean isTranslation = parameters.isTranslation();
@@ -459,11 +459,4 @@ public abstract class GeocentricAffine extends GeodeticOperation {
         return (actual instanceof Parameterized) &&
                IdentifiedObjects.isHeuristicMatchForName(((Parameterized) actual).getParameterDescriptors(), expected);
     }
-
-    /**
-     * Logs a warning about a failure to compute the Bursa-Wolf parameters.
-     */
-    private static void log(final String logger, final String method, final Exception e) {
-        Logging.recoverableException(Logging.getLogger(logger), GeocentricAffine.class, method, e);
-    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Molodensky.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Molodensky.java
index 7bb3083..6ae4795 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Molodensky.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Molodensky.java
@@ -35,11 +35,9 @@ import org.apache.sis.referencing.datum.DefaultEllipsoid;
 import org.apache.sis.referencing.operation.transform.MolodenskyTransform;
 import org.apache.sis.internal.referencing.NilReferencingObject;
 import org.apache.sis.internal.referencing.Formulas;
-import org.apache.sis.internal.system.Loggers;
 import org.apache.sis.internal.util.Constants;
 import org.apache.sis.measure.Units;
 import org.apache.sis.util.resources.Errors;
-import org.apache.sis.util.logging.Logging;
 import org.apache.sis.util.Debug;
 
 
@@ -329,8 +327,7 @@ public final class Molodensky extends GeocentricAffineBetweenGeographic {
                 values.getOrCreate(parameter).setValue(value, unit);
             } catch (ParameterNotFoundException | InvalidParameterValueException e) {
                 // Nonn-fatal since this attempt was only for information purpose.
-                Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION),
-                        Molodensky.class, "createMathTransform", e);
+                recoverableException(Molodensky.class, e);
             }
         }
 
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NADCON.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NADCON.java
index 67e9823..d6594ae 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NADCON.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NADCON.java
@@ -36,7 +36,6 @@ import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.Transformation;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
-import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
 import org.apache.sis.parameter.ParameterBuilder;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.util.CharSequences;
@@ -153,7 +152,7 @@ public final class NADCON extends AbstractProvider {
             throws ParameterNotFoundException, FactoryException
     {
         final Parameters pg  = Parameters.castOrWrap(values);
-        return InterpolatedTransform.createGeodeticTransformation(factory,
+        return DatumShiftGridFile.createGeodeticTransformation(NADCON.class, factory,
                 getOrLoad(pg.getMandatoryValue(LATITUDE), pg.getMandatoryValue(LONGITUDE)));
     }
 
@@ -184,7 +183,7 @@ public final class NADCON extends AbstractProvider {
                         final ByteBuffer buffer = ByteBuffer.allocate(4096).order(ByteOrder.LITTLE_ENDIAN);
                         final FloatBuffer fb = buffer.asFloatBuffer();
                         try (ReadableByteChannel in = Files.newByteChannel(rlat)) {
-                            DatumShiftGridLoader.log(NADCON.class, CharSequences.commonPrefix(
+                            DatumShiftGridLoader.startLoading(NADCON.class, CharSequences.commonPrefix(
                                     latitudeShifts.toString(), longitudeShifts.toString()).toString() + '…');
                             loader = new Loader(in, buffer, file);
                             loader.readGrid(fb, null, longitudeShifts);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv1.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv1.java
new file mode 100644
index 0000000..6243465
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv1.java
@@ -0,0 +1,90 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.internal.referencing.provider;
+
+import javax.xml.bind.annotation.XmlTransient;
+import org.opengis.util.FactoryException;
+import org.opengis.parameter.ParameterValueGroup;
+import org.opengis.parameter.ParameterDescriptorGroup;
+import org.opengis.parameter.ParameterNotFoundException;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.MathTransformFactory;
+import org.opengis.referencing.operation.Transformation;
+import org.apache.sis.parameter.ParameterBuilder;
+
+
+/**
+ * The provider for <cite>"National Transformation version 1"</cite> (EPSG:9614).
+ * This transform requires data that are not bundled by default with Apache SIS.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+@XmlTransient
+public final class NTv1 extends AbstractProvider {
+    /**
+     * Serial number for inter-operability with different versions.
+     */
+    private static final long serialVersionUID = 3485687719315248009L;
+
+    /**
+     * The group of all parameters expected by this coordinate operation.
+     */
+    public static final ParameterDescriptorGroup PARAMETERS;
+    static {
+        final ParameterBuilder builder = builder();
+        PARAMETERS = builder
+                .addIdentifier("9614")
+                .addName("NTv1")
+                .createGroup(NTv2.FILE);
+    }
+
+    /**
+     * Creates a new provider.
+     */
+    public NTv1() {
+        super(2, 2, PARAMETERS);
+    }
+
+    /**
+     * Returns the base interface of the {@code CoordinateOperation} instances that use this method.
+     *
+     * @return fixed to {@link Transformation}.
+     */
+    @Override
+    public Class<Transformation> getOperationType() {
+        return Transformation.class;
+    }
+
+    /**
+     * Creates a transform from the specified group of parameter values.
+     *
+     * @param  factory  the factory to use if this constructor needs to create other math transforms.
+     * @param  values   the group of parameter values.
+     * @return the created math transform.
+     * @throws ParameterNotFoundException if a required parameter was not found.
+     * @throws FactoryException if an error occurred while loading the grid.
+     */
+    @Override
+    public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup values)
+            throws ParameterNotFoundException, FactoryException
+    {
+        return NTv2.createMathTransform(NTv1.class, factory, values, 1);
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv2.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv2.java
index c65afec..9742bf0 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv2.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/NTv2.java
@@ -19,10 +19,11 @@ package org.apache.sis.internal.referencing.provider;
 import java.util.Map;
 import java.util.HashMap;
 import java.util.LinkedHashMap;
+import java.util.List;
+import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Locale;
 import java.util.logging.Level;
-import java.util.logging.LogRecord;
 import java.io.IOException;
 import java.nio.ByteOrder;
 import java.nio.ByteBuffer;
@@ -42,14 +43,14 @@ import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.Transformation;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
-import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
-import org.apache.sis.internal.system.Loggers;
 import org.apache.sis.internal.system.DataDirectory;
 import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.internal.util.Strings;
 import org.apache.sis.parameter.ParameterBuilder;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.util.collection.Cache;
-import org.apache.sis.util.logging.Logging;
+import org.apache.sis.util.collection.Containers;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.resources.Messages;
 import org.apache.sis.measure.Units;
@@ -61,7 +62,7 @@ import org.apache.sis.measure.Units;
  *
  * @author  Simon Reynard (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.7
  * @module
  */
@@ -86,7 +87,7 @@ public final class NTv2 extends AbstractProvider {
      *   <li>No default value</li>
      * </ul>
      */
-    private static final ParameterDescriptor<Path> FILE;
+    static final ParameterDescriptor<Path> FILE;
 
     /**
      * The group of all parameters expected by this coordinate operation.
@@ -134,18 +135,41 @@ public final class NTv2 extends AbstractProvider {
     public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup values)
             throws ParameterNotFoundException, FactoryException
     {
+        return createMathTransform(NTv2.class, factory, values, 2);
+    }
+
+    /**
+     * Creates a transform from the specified group of parameter values.
+     *
+     * @param  provider  the provider which is creating a transform: {@link NTv2} or {@link NTv1}.
+     * @param  factory   the factory to use if this constructor needs to create other math transforms.
+     * @param  values    the group of parameter values.
+     * @param  version   the expected version (1 or 2).
+     * @return the created math transform.
+     * @throws ParameterNotFoundException if a required parameter was not found.
+     * @throws FactoryException if an error occurred while loading the grid.
+     */
+    static MathTransform createMathTransform(final Class<? extends AbstractProvider> provider,
+            final MathTransformFactory factory, final ParameterValueGroup values, final int version)
+            throws ParameterNotFoundException, FactoryException
+    {
         final Parameters pg = Parameters.castOrWrap(values);
-        return InterpolatedTransform.createGeodeticTransformation(factory, getOrLoad(pg.getMandatoryValue(FILE)));
+        final DatumShiftGridFile<Angle,Angle> grid = getOrLoad(provider, pg.getMandatoryValue(FILE), version);
+        return DatumShiftGridFile.createGeodeticTransformation(provider, factory, grid);
     }
 
     /**
      * Returns the grid of the given name. This method returns the cached instance if it still exists,
      * or load the grid otherwise.
      *
-     * @param  file  name of the datum shift grid file to load.
+     * @param  provider  the provider which is creating a transform.
+     * @param  file      name of the datum shift grid file to load.
+     * @param  version   the expected version (1 or 2).
      */
     @SuppressWarnings("null")
-    static DatumShiftGridFile<Angle,Angle> getOrLoad(final Path file) throws FactoryException {
+    static DatumShiftGridFile<Angle,Angle> getOrLoad(final Class<? extends AbstractProvider> provider,
+            final Path file, final int version) throws FactoryException
+    {
         final Path resolved = DataDirectory.DATUM_CHANGES.resolve(file).toAbsolutePath();
         DatumShiftGridFile<?,?> grid = DatumShiftGridFile.CACHE.peek(resolved);
         if (grid == null) {
@@ -154,12 +178,12 @@ public final class NTv2 extends AbstractProvider {
                 grid = handler.peek();
                 if (grid == null) {
                     try (ReadableByteChannel in = Files.newByteChannel(resolved)) {
-                        DatumShiftGridLoader.log(NTv2.class, file);
-                        final Loader loader = new Loader(in, file);
-                        grid = loader.readGrid();
-                        loader.reportWarnings();
+                        DatumShiftGridLoader.startLoading(provider, file);
+                        final Loader loader = new Loader(in, file, version);
+                        grid = loader.readAllGrids();
+                        loader.report(provider);
                     } catch (IOException | NoninvertibleTransformException | RuntimeException e) {
-                        throw DatumShiftGridLoader.canNotLoad("NTv2", file, e);
+                        throw DatumShiftGridLoader.canNotLoad(provider.getSimpleName(), file, e);
                     }
                     grid = grid.useSharedData();
                 }
@@ -175,10 +199,21 @@ public final class NTv2 extends AbstractProvider {
 
     /**
      * Loaders of NTv2 data. Instances of this class exist only at loading time.
+     * More information on that file format can be found with
+     * <a href="https://github.com/Esri/ntv2-file-routines">ESRI NTv2 routines</a>.
+     *
+     * <p>A NTv2 file contains an arbitrary number of sub-files, where each sub-file is a grid.
+     * There is at least one grid (the parent), and potentially many sub-grids of higher density.
+     * At the beginning is an overview header block of information that is common to all sub-files.
+     * Then there is other headers specific to each sub-files.</p>
+     *
+     * <p>While this loader is primarily targeted at loading NTv2 files, it can also opportunistically
+     * read NTv1 files. The two file formats differ by some header records having different names (but
+     * same meanings), the possibility to have sub-grids and the presence of accuracy information.</p>
      *
      * @author  Simon Reynard (Geomatys)
      * @author  Martin Desruisseaux (Geomatys)
-     * @version 1.0
+     * @version 1.1
      * @since   0.7
      * @module
      */
@@ -190,79 +225,115 @@ public final class NTv2 extends AbstractProvider {
         private static final int RECORD_LENGTH = 16;
 
         /**
-         * Maximum number of characters of a key in a header record.
+         * Maximum number of characters for a key in a header record.
+         * Expected keys are listed in the {@link #TYPES} map.
          */
         private static final int KEY_LENGTH = 8;
 
         /**
-         * Type of data allowed in header records.
+         * Type of data allowed in header records. Each record header identified by a key contains a value
+         * of a type hard-coded by the NTv2 specification; the type is not specified in the file itself.
          */
-        private static final int STRING_TYPE = 0, INTEGER_TYPE = 1, DOUBLE_TYPE = 2;
+        private enum DataType {STRING, INTEGER, DOUBLE};
 
         /**
-         * Some known keywords that may appear in NTv2 header records.
+         * Some known keywords that may appear in NTv2 header records, associated the the expected type of values.
+         * The type is not encoded in a NTv2 file; it has to be hard-coded in this table. The first 11 entries in
+         * this map (ignoring entries marked by "NTv1") are typically found in overview header, and the remaining
+         * entries in the sub-grid headers.
          */
-        private static final Map<String,Integer> TYPES;
+        private static final Map<String,DataType> TYPES;
         static {
-            final Map<String,Integer> types = new HashMap<>(32);
-            final Integer string  = STRING_TYPE;    // Autoboxing
-            final Integer integer = INTEGER_TYPE;
-            final Integer real    = DOUBLE_TYPE;
-            types.put("NUM_OREC", integer);         // Number of records in the header - usually 11
-            types.put("NUM_SREC", integer);         // Number of records in the header of sub-grids - usually 11
-            types.put("NUM_FILE", integer);         // Number of sub-grids
-            types.put("GS_TYPE",  string);          // Units: "SECONDS", "MINUTES" or "DEGREES"
-            types.put("VERSION",  string);          // Grid version
-            types.put("SYSTEM_F", string);          // Source CRS
-            types.put("SYSTEM_T", string);          // Target CRS
-            types.put("MAJOR_F",  real);            // Semi-major axis of source ellipsoid (in metres)
-            types.put("MINOR_F",  real);            // Semi-minor axis of source ellipsoid (in metres)
-            types.put("MAJOR_T",  real);            // Semi-major axis of target ellipsoid (in metres)
-            types.put("MINOR_T",  real);            // Semi-minor axis of target ellipsoid (in metres)
-            types.put("SUB_NAME", string);          // Sub-grid identifier
-            types.put("PARENT",   string);          // Parent grid
-            types.put("CREATED",  string);          // Creation time
-            types.put("UPDATED",  string);          // Update time
-            types.put("S_LAT",    real);            // Southmost φ value
-            types.put("N_LAT",    real);            // Northmost φ value
-            types.put("E_LONG",   real);            // Eastmost λ value - west is positive, east is negative
-            types.put("W_LONG",   real);            // Westmost λ value - west is positive, east is negative
-            types.put("LAT_INC",  real);            // Increment on φ axis
-            types.put("LONG_INC", real);            // Increment on λ axis - positive toward west
-            types.put("GS_COUNT", integer);         // Number of sub-grid records following
+            final Map<String,DataType> types = new HashMap<>(38);
+/* NTv1 */  types.put("HEADER",   DataType.INTEGER);        // Number of header records (replaced by NUM_OREC)
+            types.put("NUM_OREC", DataType.INTEGER);        // Number of records in the header - usually 11
+            types.put("NUM_SREC", DataType.INTEGER);        // Number of records in the header of sub-grids - usually 11
+            types.put("NUM_FILE", DataType.INTEGER);        // Number of sub-grids
+/* NTv1 */  types.put("TYPE",     DataType.STRING);         // Grid shift data type (replaced by GS_TYPE)
+            types.put("GS_TYPE",  DataType.STRING);         // Units: "SECONDS", "MINUTES" or "DEGREES"
+            types.put("VERSION",  DataType.STRING);         // Grid version
+/* NTv1 */  types.put("FROM",     DataType.STRING);         // Source CRS (replaced by SYSTEM_F)
+/* NTv1 */  types.put("TO",       DataType.STRING);         // Target CRS (replaced by SYSTEM_T)
+            types.put("SYSTEM_F", DataType.STRING);         // Source CRS
+            types.put("SYSTEM_T", DataType.STRING);         // Target CRS
+            types.put("DATUM_F",  DataType.STRING);         // Source datum (some time replace SYSTEM_F)
+            types.put("DATUM_T",  DataType.STRING);         // Target datum (some time replace SYSTEM_T)
+            types.put("MAJOR_F",  DataType.DOUBLE);         // Semi-major axis of source ellipsoid (in metres)
+            types.put("MINOR_F",  DataType.DOUBLE);         // Semi-minor axis of source ellipsoid (in metres)
+            types.put("MAJOR_T",  DataType.DOUBLE);         // Semi-major axis of target ellipsoid (in metres)
+            types.put("MINOR_T",  DataType.DOUBLE);         // Semi-minor axis of target ellipsoid (in metres)
+            types.put("SUB_NAME", DataType.STRING);         // Sub-grid identifier
+            types.put("PARENT",   DataType.STRING);         // Parent grid
+            types.put("CREATED",  DataType.STRING);         // Creation time
+            types.put("UPDATED",  DataType.STRING);         // Update time
+            types.put("S_LAT",    DataType.DOUBLE);         // Southmost φ value
+            types.put("N_LAT",    DataType.DOUBLE);         // Northmost φ value
+            types.put("E_LONG",   DataType.DOUBLE);         // Eastmost λ value - west is positive, east is negative
+            types.put("W_LONG",   DataType.DOUBLE);         // Westmost λ value - west is positive, east is negative
+/* NTv1 */  types.put("N_GRID",   DataType.DOUBLE);         // Latitude grid interval (replaced by LAT_INC)
+/* NTv1 */  types.put("W_GRID",   DataType.DOUBLE);         // Longitude grid interval (replaced by LONG_INC)
+            types.put("LAT_INC",  DataType.DOUBLE);         // Increment on φ axis
+            types.put("LONG_INC", DataType.DOUBLE);         // Increment on λ axis - positive toward west
+            types.put("GS_COUNT", DataType.INTEGER);        // Number of sub-grid records following
             TYPES = types;
+            /*
+             * NTv1 as two last unnamed records of DataType.DOUBLE: "Semi_Major_Axis_From"
+             * and "Semi_Major_Axis_To". Those records are currently ignored.
+             */
         }
 
         /**
-         * The header content. Keys are strings like {@code "VERSION"}, {@code "SYSTEM_F"},
-         * <var>etc.</var>. Values are {@link String}, {@link Integer} or {@link Double}.
-         * If some keys are unrecognized, they will be put in this map with the {@code null} value
-         * and the {@link #hasUnrecognized} field will be set to {@code true}.
+         * The headers content, as the union of the overview header and the header in process of being read.
+         * Keys are strings like {@code "VERSION"}, {@code "SYSTEM_F"}, {@code "LONG_INC"}, <i>etc.</i>.
+         * Values are {@link String}, {@link Integer} or {@link Double}. If some keys are unrecognized,
+         * they will be put in this map with the {@code null} value and the {@link #hasUnrecognized} flag
+         * will be set to {@code true}.
          */
         private final Map<String,Object> header;
 
         /**
+         * Keys of {@link #header} for entries that were declared in the overview header.
+         * This is used after {@link #readGrid(Map, Map)} execution for discarding all
+         * entries specific to sub-grids, for avoiding to mix entries from two sub-grids.
+         */
+        private final String[] overviewKeys;
+
+        /**
+         * {@code true} if we are reading a NTv2 file, or {@code false} if we are reading a NTv1 file.
+         */
+        private final boolean isV2;
+
+        /**
          * {@code true} if the {@code header} map contains at least one key associated to a null value.
          */
         private boolean hasUnrecognized;
 
         /**
-         * Number of grids remaining in the file. This value is set in the constructor,
-         * then decremented at every call to {@link #readGrid()}.
+         * Number of grids expected in the file.
+         */
+        private final int numGrids;
+
+        /**
+         * Dates at which the grid has been created or updated, or {@code null} if unknown.
+         * Used for information purpose only.
          */
-        private int remainingGrids;
+        private String created, updated;
 
         /**
          * Creates a new reader for the given channel.
          * This constructor parses the header immediately, but does not read any grid.
+         * A hint about expected NTv2 version is given, but this constructor may override
+         * that hint with information found in the file.
          *
          * @param  channel  where to read data from.
-         * @param  file     path to the longitude and latitude difference file. Used for parameter declaration and error reporting.
+         * @param  file     path to the longitude and latitude difference file.
+         *                  Used for parameter declaration and error reporting.
+         * @param  version  the expected version (1 or 2).
          * @throws FactoryException if a data record can not be parsed.
          */
-        Loader(final ReadableByteChannel channel, final Path file) throws IOException, FactoryException {
+        Loader(final ReadableByteChannel channel, final Path file, int version) throws IOException, FactoryException {
             super(channel, ByteBuffer.allocate(4096), file);
-            this.header = new LinkedHashMap<>();
+            header = new LinkedHashMap<>();
             ensureBufferContains(RECORD_LENGTH);
             if (isLittleEndian(buffer.getInt(KEY_LENGTH))) {
                 buffer.order(ByteOrder.LITTLE_ENDIAN);
@@ -272,10 +343,160 @@ public final class NTv2 extends AbstractProvider {
              * NUM_OREC, NUM_SREC, NUM_FILE, GS_TYPE, VERSION, SYSTEM_F, SYSTEM_T, MAJOR_F, MINOR_F, MAJOR_T,
              * MINOR_T.
              */
-            readHeader(11, "NUM_OREC");
-            remainingGrids = (Integer) get("NUM_FILE");
-            if (remainingGrids < 1) {
-                throw new FactoryException(Errors.format(Errors.Keys.UnexpectedValueInElement_2, "NUM_FILE", remainingGrids));
+            readHeader(version >= 2 ? 11 : 12, "NUM_OREC");
+            /*
+             * The version number is a string like "NTv2.0". If there is no version number, it is probably NTv1
+             * since the "VERSION" record was introduced only in version 2. In such case the `version` parameter
+             * should have been 1; in case of doubt we do not modify the provided value.
+             */
+            final String vs = (String) get("VERSION", false);
+            if (vs != null) {
+                for (int i=0; i<vs.length(); i++) {
+                    final char c = vs.charAt(i);
+                    if (c >= '0' && c <= '9') {
+                        version = c - '0';
+                        break;
+                    }
+                }
+            }
+            /*
+             * Subgrids are NTv2 features which did not existed in NTv1. If we expect a NTv2 file,
+             * the record is mandatory. If we expect a NTv1 file, the record should not be present
+             * but we nevertheless check in case we have been misleaded by a missing "VERSION" record.
+             */
+            final Integer n = (Integer) get("NUM_FILE", (vs != null) && version >= 2);
+            isV2 = (n != null);
+            if (isV2) {
+                numGrids = n;
+                if (numGrids < 1) {
+                    throw new FactoryException(Errors.format(Errors.Keys.UnexpectedValueInElement_2, "NUM_FILE", n));
+                }
+            } else {
+                numGrids = 1;
+            }
+            overviewKeys = header.keySet().toArray(new String[header.size()]);
+        }
+
+        /**
+         * Returns {@code true} if the given value seems to be stored in little endian order.
+         * The strategy is to read an integer that we expect to be small (the HEADER or NUM_OREC
+         * value which should be 12 or 11) and to check which order gives the smallest value.
+         */
+        private static boolean isLittleEndian(final int n) {
+            return Integer.compareUnsigned(n, Integer.reverseBytes(n)) > 0;
+        }
+
+        /**
+         * Reads a string at the current buffer position, assuming ASCII encoding.
+         * After this method call, the buffer position will be the first byte after
+         * the string. The buffer content is unmodified.
+         *
+         * @param  length  number of bytes to read.
+         */
+        private String readString(int length) {
+            final byte[] array = buffer.array();
+            final int position = buffer.position();
+            buffer.position(position + length);     // Update before we modify `length`.
+            while (length > position && array[position + length - 1] <= ' ') length--;
+            return new String(array, position, length, StandardCharsets.US_ASCII).trim();
+        }
+
+        /**
+         * Reads all records found in the header, starting from the current buffer position.
+         * The header may be the overview header (in which case we expect a number of records
+         * given by {@code HEADER} or {@code NUM_OREC} value) or a sub-grid header (in which
+         * case we expect {@code NUM_SREC} records).
+         *
+         * <p>The {@code numRecords} given in argument is a default value.
+         * It will be updated as soon as the {@code numKey} record is found.</p>
+         *
+         * @param  numRecords  default number of expected records (usually 11).
+         * @param  numkey      key of the record giving the number of records: {@code "NUM_OREC"} or {@code "NUM_SREC"}.
+         */
+        private void readHeader(int numRecords, final String numkey) throws IOException, FactoryException {
+            for (int i=0; i < numRecords; i++) {
+                ensureBufferContains(RECORD_LENGTH);
+                final String key = readString(KEY_LENGTH).toUpperCase(Locale.US).replace(' ', '_');
+                final DataType type = TYPES.get(key);
+                final Comparable<?> value;
+                if (type == null) {
+                    value = null;
+                    hasUnrecognized = true;
+                } else switch (type) {                              // TODO: check if we can simplify in JDK14.
+                    default: throw new AssertionError(type);
+                    case STRING: value = readString(RECORD_LENGTH - KEY_LENGTH); break;
+                    case DOUBLE: value = buffer.getDouble(); break;
+                    case INTEGER: {
+                        final int n = buffer.getInt();
+                        buffer.position(buffer.position() + Integer.BYTES);
+                        if (key.equals(numkey) || key.equals("HEADER")) {
+                            /*
+                             * HEADER (NTv1), NUM_OREC (NTv2) or NUM_SREC specify the number of records expected
+                             * in the header, which may the the header that we are reading right now. If value
+                             * applies to the reader we are reading, we need to update `numRecords` on the fly.
+                             */
+                            numRecords = n;
+                        }
+                        value = n;
+                        break;
+                    }
+                }
+                final Object old = header.put(key, value);
+                if (old != null && !old.equals(value)) {
+                    throw new FactoryException(Errors.format(Errors.Keys.KeyCollision_1, key));
+                }
+            }
+            if (created == null) created = Strings.trimOrNull((String) get("CREATED", false));
+            if (updated == null) updated = Strings.trimOrNull((String) get("UPDATED", false));
+        }
+
+        /**
+         * Reads all grids and returns the root grid. After reading all grids, this method rearrange
+         * them in a child-parent relationship. The result is a tree with a single root containing
+         * sub-grids (if any) as children.
+         */
+        final DatumShiftGridFile<Angle,Angle> readAllGrids() throws IOException, FactoryException, NoninvertibleTransformException {
+            final Map<String,      DatumShiftGridFile<Angle,Angle>>  grids    = new HashMap<>(Containers.hashMapCapacity(numGrids));
+            final Map<String, List<DatumShiftGridFile<Angle,Angle>>> children = new LinkedHashMap<>();   // Should have few entries.
+            while (grids.size() < numGrids) {
+                readGrid(grids, children);
+            }
+            /*
+             * Assign the sub-grids to their parent only after we finished to read all grids.
+             * Doing this work last is more robust to cases where grids are in random order.
+             *
+             * Notes: if the parent-child graph contains cycles (deeper than a child declaring itself as its parent),
+             *        the grids in cycles will be lost. This is because we need a grid without parent for getting the
+             *        graph added in the roots list. There is currently no mechanism for detecting those problems.
+             */
+            final List<DatumShiftGridFile<Angle,Angle>> roots = new ArrayList<>();
+            for (final Map.Entry<String, List<DatumShiftGridFile<Angle,Angle>>> entry : children.entrySet()) {
+                final DatumShiftGridFile<Angle,Angle> parent = grids.get(entry.getKey());
+                final List<DatumShiftGridFile<Angle,Angle>> subgrids = entry.getValue();
+                if (parent != null) {
+                    /*
+                     * Verify that the children does not declare themselves as their parent.
+                     * It may happen if SUB_GRID and PARENT have the same value, typically a
+                     * null or empty value if those records were actually unspecified.
+                     */
+                    for (int i=subgrids.size(); --i >= 0;) {
+                        if (subgrids.get(i) == parent) {      // Want identity check, no need for equals(Object).
+                            subgrids.remove(i);
+                            roots.add(parent);
+                            break;
+                        }
+                    }
+                    if (!subgrids.isEmpty()) {
+                        parent.setSubGrids(subgrids);
+                    }
+                } else {
+                    roots.addAll(subgrids);
+                }
+            }
+            switch (roots.size()) {
+                case 0:  throw new FactoryException(Errors.format(Errors.Keys.CanNotRead_1, file));
+                case 1:  return roots.get(0);
+                default: return DatumShiftGridGroup.create(file, roots);
             }
         }
 
@@ -285,47 +506,48 @@ public final class NTv2 extends AbstractProvider {
          * The first grid can cover a large area with a coarse resolution, and next grids cover smaller
          * areas overlapping the first grid but with finer resolution.
          *
-         * Current SIS implementation does not yet handle the above-cited hierarchy of grids.
-         * For now we just take the first one.
-         *
          * <p>NTv2 grids contain also information about shifts accuracy. This is not yet handled by SIS,
          * except for determining an approximate grid cell resolution.</p>
+         *
+         * @param  addTo     the map where to add the grid with the grid name as the key.
+         * @param  children  the map where to add children with the parent name as the key.
          */
-        final DatumShiftGridFile<Angle,Angle> readGrid() throws IOException, FactoryException, NoninvertibleTransformException {
-            if (--remainingGrids < 0) {
-                throw new FactoryException(Errors.format(Errors.Keys.CanNotRead_1, file));
+        private void readGrid(final Map<String, DatumShiftGridFile<Angle,Angle>> addTo,
+                final Map<String, List<DatumShiftGridFile<Angle,Angle>>> children)
+                throws IOException, FactoryException, NoninvertibleTransformException
+        {
+            if (isV2) {
+                readHeader((Integer) get("NUM_SREC", null, null), "NUM_SREC");
             }
-            final Object[] overviewKeys = header.keySet().toArray();
-            readHeader((Integer) get("NUM_SREC"), "NUM_SREC");
             /*
              * Extract the geographic bounding box and cell size. While different units are allowed,
-             * in practice we usually have seconds of angle. This units has the advantage of allowing
+             * in practice we usually have seconds of angle. This unit has the advantage of allowing
              * all floating-point values to be integers.
              *
              * Note that the longitude values in NTv2 files are positive WEST.
              */
             final Unit<Angle> unit;
             final double precision;
-            final String name = (String) get("GS_TYPE");
-            if (name.equalsIgnoreCase("SECONDS")) {                 // Most common value
+            final String type = (String) get("GS_TYPE", "TYPE", null);
+            if (type.equalsIgnoreCase("SECONDS")) {                 // Most common value
                 unit = Units.ARC_SECOND;
                 precision = SECOND_PRECISION;                       // Used only as a hint; will not hurt if wrong.
-            } else if (name.equalsIgnoreCase("MINUTES")) {
+            } else if (type.equalsIgnoreCase("MINUTES")) {
                 unit = Units.ARC_MINUTE;
                 precision = SECOND_PRECISION / 60;                  // Used only as a hint; will not hurt if wrong.
-            } else if (name.equalsIgnoreCase("DEGREES")) {
+            } else if (type.equalsIgnoreCase("DEGREES")) {
                 unit = Units.DEGREE;
                 precision = SECOND_PRECISION / DEGREES_TO_SECONDS;  // Used only as a hint; will not hurt if wrong.
             } else {
-                throw new FactoryException(Errors.format(Errors.Keys.UnexpectedValueInElement_2, "GS_TYPE", name));
+                throw new FactoryException(Errors.format(Errors.Keys.UnexpectedValueInElement_2, "GS_TYPE", type));
             }
-            final double  ymin     = (Double)  get("S_LAT");
-            final double  ymax     = (Double)  get("N_LAT");
-            final double  xmin     = (Double)  get("E_LONG");       // Sign reversed compared to usual convention.
-            final double  xmax     = (Double)  get("W_LONG");       // Idem.
-            final double  dy       = (Double)  get("LAT_INC");
-            final double  dx       = (Double)  get("LONG_INC");     // Positive toward west.
-            final Integer declared = (Integer) header.get("GS_COUNT");
+            final double  ymin     = (Double)  get("S_LAT",    null,     null);
+            final double  ymax     = (Double)  get("N_LAT",    null,     null);
+            final double  xmin     = (Double)  get("E_LONG",   null,     null);   // Sign reversed compared to usual convention.
+            final double  xmax     = (Double)  get("W_LONG",   null,     null);   // Idem.
+            final double  dy       = (Double)  get("LAT_INC",  "N_GRID", null);
+            final double  dx       = (Double)  get("LONG_INC", "W_GRID", null);   // Positive toward west.
+            final Integer declared = (Integer) get("GS_COUNT", false);
             final int     width    = Math.toIntExact(Math.round((xmax - xmin) / dx + 1));
             final int     height   = Math.toIntExact(Math.round((ymax - ymin) / dy + 1));
             final int     count    = Math.multiplyExact(width, height);
@@ -340,17 +562,38 @@ public final class NTv2 extends AbstractProvider {
              * free us from reversing the sign of longitude translations in the code below; instead, this reversal
              * will be handled by grid.coordinateToGrid MathTransform and its inverse.
              */
-            final DatumShiftGridFile.Float<Angle,Angle> grid = new DatumShiftGridFile.Float<>(2,
-                    unit, unit, true, -xmin, ymin, -dx, dy, width, height, PARAMETERS, file);
-            @SuppressWarnings("MismatchedReadAndWriteOfArray") final float[] tx = grid.offsets[0];
-            @SuppressWarnings("MismatchedReadAndWriteOfArray") final float[] ty = grid.offsets[1];
-            for (int i=0; i<count; i++) {
-                ensureBufferContains(4 * Float.BYTES);
-                ty[i] = (float) (buffer.getFloat() / dy);   // Division by dx and dy because isCellValueRatio = true.
-                tx[i] = (float) (buffer.getFloat() / dx);
-                final double accuracy = Math.min(buffer.getFloat() / dy, buffer.getFloat() / dx);
-                if (accuracy > 0 && !(accuracy >= grid.accuracy)) {   // Use '!' for replacing the initial NaN.
-                    grid.accuracy = accuracy;                         // Smallest non-zero accuracy.
+            final double size = Math.max(dx, dy);
+            final DatumShiftGridFile<Angle,Angle> grid;
+            if (isV2) {
+                final DatumShiftGridFile.Float<Angle,Angle> data;
+                data = new DatumShiftGridFile.Float<>(2, unit, unit, true,
+                        -xmin, ymin, -dx, dy, width, height, PARAMETERS, file);
+                @SuppressWarnings("MismatchedReadAndWriteOfArray") final float[] tx = data.offsets[0];
+                @SuppressWarnings("MismatchedReadAndWriteOfArray") final float[] ty = data.offsets[1];
+                data.accuracy = Double.NaN;
+                for (int i=0; i<count; i++) {
+                    ensureBufferContains(4 * Float.BYTES);
+                    ty[i] = (float) (buffer.getFloat() / dy);   // Division by dx and dy because isCellValueRatio = true.
+                    tx[i] = (float) (buffer.getFloat() / dx);
+                    final double accuracy = Math.min(buffer.getFloat() / dy, buffer.getFloat() / dx);
+                    if (accuracy > 0 && !(accuracy >= data.accuracy)) {     // Use '!' for replacing the initial NaN.
+                        data.accuracy = accuracy;                           // Smallest non-zero accuracy.
+                    }
+                }
+                grid = DatumShiftGridCompressed.compress(data, null, precision / size);
+            } else {
+                /*
+                 * NTv1: same as NTv2 but using double precision and without accuracy information.
+                 */
+                final DatumShiftGridFile.Double<Angle,Angle> data;
+                grid = data = new DatumShiftGridFile.Double<>(2, unit, unit, true,
+                        -xmin, ymin, -dx, dy, width, height, PARAMETERS, file);
+                @SuppressWarnings("MismatchedReadAndWriteOfArray") final double[] tx = data.offsets[0];
+                @SuppressWarnings("MismatchedReadAndWriteOfArray") final double[] ty = data.offsets[1];
+                for (int i=0; i<count; i++) {
+                    ensureBufferContains(2 * Double.BYTES);
+                    ty[i] = buffer.getDouble() / dy;
+                    tx[i] = buffer.getDouble() / dx;
                 }
             }
             /*
@@ -358,91 +601,83 @@ public final class NTv2 extends AbstractProvider {
              * during inverse transformations. If we did not found that information in the file, compute
              * an arbitrary default accuracy.
              */
-            final double size = Math.max(dx, dy);
-            if (Double.isNaN(grid.accuracy)) {
+            if (!(grid.accuracy > 0)) {                 // Use ! for catching NaN values (paranoiac check).
                 grid.accuracy = Units.DEGREE.getConverterTo(unit).convert(Formulas.ANGULAR_TOLERANCE) / size;
             }
-            header.keySet().retainAll(Arrays.asList(overviewKeys));   // Keep only overview records.
-            return DatumShiftGridCompressed.compress(grid, null, precision / size);
-        }
-
-        /**
-         * Returns {@code true} if the given value seems to be stored in little endian order.
-         */
-        private static boolean isLittleEndian(final int n) {
-            return Integer.compareUnsigned(n, Integer.reverseBytes(n)) > 0;
-        }
-
-        /**
-         * Reads a string at the given position in the buffer.
-         */
-        private String readString(final int position, int length) {
-            final byte[] array = buffer.array();
-            while (length > position && array[position + length - 1] <= ' ') length--;
-            return new String(array, position, length, StandardCharsets.US_ASCII).trim();
+            /*
+             * Add the grid to two collection. The first collection associates this grid to its name, and the
+             * second collection associates the grid to its parent. We do not try to resolve the child-parent
+             * relationship here; we will do that after all sub-grids have been read.
+             */
+            final String name = (String) get("SUB_NAME", numGrids > 1);
+            if (addTo.put(name, grid) != null) {
+                throw new FactoryException(Errors.format(Errors.Keys.DuplicatedIdentifier_1, name));
+            }
+            children.computeIfAbsent((String) get("PARENT", numGrids > 1), (k) -> new ArrayList<>()).add(grid);
+            /*
+             * End of grid parsing. Remove all header entries that are specific to this sub-grid.
+             * After this operation, `header` will contain only overview records.
+             */
+            header.keySet().retainAll(Arrays.asList(overviewKeys));
         }
 
         /**
-         * Reads all records found in the header, starting from the current buffer position.
-         * It may be the overview header (in which case we expect {@code NUM_OREC} records)
-         * or a sub-grid header (in which case we expect {@code NUM_SREC} records).
+         * Gets the value for the given key. If the value is absent, then this method throws an exception
+         * if {@code mandatory} is {@code true} or returns {@code null} otherwise.
          *
-         * @param  numRecords  default number of expected records (usually 11).
-         * @param  numkey      key of the record giving the number of records: {@code "NUM_OREC"} or {@code "NUM_SREC"}.
+         * @param  key        key of the value to search.
+         * @param  mandatory  whether to throw an exception if the value is not found.
+         * @return value associated to the given key, or {@code null} if none and not mandatory.
          */
-        private void readHeader(int numRecords, final String numkey) throws IOException, FactoryException {
-            int position = buffer.position();
-            for (int i=0; i < numRecords; i++) {
-                ensureBufferContains(RECORD_LENGTH);
-                final String key = readString(position, KEY_LENGTH).toUpperCase(Locale.US);
-                position += KEY_LENGTH;
-                final Integer type = TYPES.get(key);
-                final Comparable<?> value;
-                if (type == null) {
-                    value = null;
-                    hasUnrecognized = true;
-                } else switch (type) {
-                    case STRING_TYPE: {
-                        value = readString(position, RECORD_LENGTH - KEY_LENGTH);
-                        break;
-                    }
-                    case INTEGER_TYPE: {
-                        final int n = buffer.getInt(position);
-                        if (key.equals(numkey)) {
-                            numRecords = n;
-                        }
-                        value = n;
-                        break;
-                    }
-                    case DOUBLE_TYPE: {
-                        value = buffer.getDouble(position);
-                        break;
-                    }
-                    default: throw new AssertionError(type);
-                }
-                final Object old = header.put(key, value);
-                if (old != null && !old.equals(value)) {
-                    throw new FactoryException(Errors.format(Errors.Keys.KeyCollision_1, key));
-                }
-                buffer.position(position += RECORD_LENGTH - KEY_LENGTH);
+        private Object get(final String key, final boolean mandatory) throws FactoryException {
+            final Object value = header.get(key);
+            if (value != null || !mandatory) {
+                return value;
             }
+            throw new FactoryException(Errors.format(Errors.Keys.PropertyNotFound_2, file, key));
         }
 
         /**
          * Returns the value for the given key, or thrown an exception if the value is not found.
+         * Before to fail if the key is not found, this method searches for a value associated to
+         * an alternative name. That alternative should be the name used in legacy NTv1.
+         *
+         * @param  key  key of the value to search.
+         * @param  alt  alternative key name, or name used in NTv1, or {@code null} if none.
+         * @param  kv1  name used in NTv1, or {@code null} if none.
+         * @return value associated to the given key (never {@code null}).
          */
-        private Object get(final String key) throws FactoryException {
-            final Object value = header.get(key);
-            if (value != null) {
-                return value;
+        private Object get(final String key, final String alt, final String kv1) throws FactoryException {
+            Object value = header.get(key);
+            if (value == null) {
+                value = header.get(alt);
+                if (value == null) {
+                    value = header.get(kv1);
+                    if (value == null) {
+                        throw new FactoryException(Errors.format(Errors.Keys.PropertyNotFound_2, file, key));
+                    }
+                }
             }
-            throw new FactoryException(Errors.format(Errors.Keys.PropertyNotFound_2, file, key));
+            return value;
         }
 
         /**
          * If we had any warnings during the loading process, report them now.
+         *
+         * @param  caller  the provider which created this loader.
          */
-        void reportWarnings() {
+        void report(final Class<? extends AbstractProvider> caller) {
+            try {
+                final String source = (String) get("SYSTEM_F", "DATUM_F", "FROM");
+                final String target = (String) get("SYSTEM_T", "DATUM_T", "TO");
+                log(caller, Resources.forLocale(null).getLogRecord(Level.FINE,
+                            Resources.Keys.UsingDatumShiftGrid_4, source, target,
+                            (created != null) ? created : "?",
+                            (updated != null) ? updated : "?"));
+            } catch (FactoryException e) {
+                recoverableException(caller, e);
+                // Ignore since above code is only for information purpose.
+            }
             if (hasUnrecognized) {
                 final StringBuilder keywords = new StringBuilder();
                 for (final Map.Entry<String,Object> entry : header.entrySet()) {
@@ -453,10 +688,8 @@ public final class NTv2 extends AbstractProvider {
                         keywords.append(entry.getKey());
                     }
                 }
-                final LogRecord record = Messages.getResources(null).getLogRecord(Level.WARNING,
-                        Messages.Keys.UnknownKeywordInRecord_2, file, keywords.toString());
-                record.setLoggerName(Loggers.COORDINATE_OPERATION);
-                Logging.log(NTv2.class, "createMathTransform", record);
+                log(caller, Messages.getResources(null).getLogRecord(Level.WARNING,
+                        Messages.Keys.UnknownKeywordInRecord_2, file, keywords.toString()));
             }
         }
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
index ac74508..3313b97 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
@@ -21,6 +21,7 @@ import java.util.Objects;
 import java.io.Serializable;
 import javax.measure.Unit;
 import javax.measure.Quantity;
+import javax.measure.UnitConverter;
 import org.opengis.geometry.Envelope;
 import org.opengis.parameter.ParameterDescriptorGroup;
 import org.opengis.referencing.operation.Matrix;
@@ -117,6 +118,15 @@ import org.apache.sis.measure.Units;
  * in more than two dimensions. See the above <cite>datum shift by geocentric translations</cite> use case for
  * an example.
  *
+ * <h2>Sub-grids</h2>
+ * Some datum shift grid files provide a grid valid on a wide region, refined with denser sub-grids in smaller regions.
+ * For each point to transform, the {@link org.opengis.referencing.operation.MathTransform} should search and use the
+ * densest sub-grid containing the point. This functionality is not supported directly by {@code DatumShiftGrid},
+ * but can be achieved by organizing many transforms in a tree. The first step is to create an instance of
+ * {@link org.apache.sis.referencing.operation.transform.InterpolatedTransform} for each {@code DatumShiftGrid}.
+ * Then, those transforms with their domain of validity can be given to
+ * {@link org.apache.sis.referencing.operation.transform.MathTransforms#specialize MathTransforms.specialize(…)}.
+ *
  * <h2>Serialization</h2>
  * Serialized objects of this class are not guaranteed to be compatible with future Apache SIS releases.
  * Serialization support is appropriate for short term storage or RMI between applications running the
@@ -124,7 +134,7 @@ import org.apache.sis.measure.Units;
  * NTv2 should be preferred.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @param <C>  dimension of the coordinate unit (usually {@link javax.measure.quantity.Angle}).
  * @param <T>  dimension of the translation unit (usually {@link javax.measure.quantity.Angle}
@@ -266,11 +276,14 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
 
     /**
      * Returns the domain of validity of input coordinates that can be specified to the
-     * {@link #interpolateAt interpolateAt(…)} method. Coordinates outside that domain of
-     * validity will still be accepted, but the extrapolated results may be very wrong.
+     * {@link #interpolateAt interpolateAt(…)} method. Coordinates outside that domain
+     * will still be accepted, but results will be extrapolations possibly far from reality.
      *
-     * <p>The unit of measurement for the coordinate values in the returned envelope is
-     * given by {@link #getCoordinateUnit()}. The envelope CRS is undefined.</p>
+     * <p>The envelope coordinates are computed at cell centers; the envelope does not contain
+     * the margin of 0.5 cell between cell center and cell border at the edges of the envelope.
+     * The unit of measurement for the coordinate values in the returned envelope is given by
+     * {@link #getCoordinateUnit()}. The envelope CRS is not set, but its value is implicitly
+     * the CRS of grid input coordinates.</p>
      *
      * @return the domain covered by this grid.
      * @throws TransformException if an error occurred while computing the envelope.
@@ -278,12 +291,47 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
     public Envelope getDomainOfValidity() throws TransformException {
         final GeneralEnvelope env = new GeneralEnvelope(gridSize.length);
         for (int i=0; i<gridSize.length; i++) {
-            env.setRange(i, -0.5, gridSize[i] - 0.5);
+            /*
+             * Note: a previous version was using the following code in an attempt to encompass
+             * fully all cells (keeping in mind that the `coordinatetoGrid` maps cell centers):
+             *
+             *    env.setRange(i, -0.5, gridSize[i] - 0.5);
+             *
+             * However it was causing spurious overlaps when two grids are side-by-side
+             * (no overlapping) but one grid has larger cells than the other other grid.
+             * The 0.5 cell expansion caused the grid with larger cells to overlap the
+             * grid with smaller cells. This case happens with NTv2 datum shift grid.
+             */
+            env.setRange(i, 0, gridSize[i] - 1);
         }
         return Envelopes.transform(getCoordinateToGrid().inverse(), env);
     }
 
     /**
+     * Returns the domain of validity converted to the specified unit of measurement.
+     * A common use case for this method is for converting the domain of a NADCON or
+     * NTv2 datum shift grid file, which are expressed in {@link Units#ARC_SECOND},
+     * to {@link Units#DEGREE}.
+     *
+     * @param  unit  the desired unit of measurement.
+     * @return the domain covered by this grid, converted to the given unit of measurement.
+     * @throws TransformException if an error occurred while computing the envelope.
+     *
+     * @since 1.1
+     */
+    public Envelope getDomainOfValidity(final Unit<C> unit) throws TransformException {
+        final UnitConverter uc = getCoordinateUnit().getConverterTo(unit);
+        if (uc.isIdentity()) {
+            return getDomainOfValidity();
+        }
+        final GeneralEnvelope domain = GeneralEnvelope.castOrCopy(getDomainOfValidity());
+        for (int i=domain.getDimension(); --i >= 0;) {
+            domain.setRange(i, uc.convert(domain.getLower(i)), uc.convert(domain.getUpper(i)));
+        }
+        return domain;
+    }
+
+    /**
      * Returns the unit of measurement of input values, before conversion to grid indices.
      * The coordinate unit is usually {@link Units#DEGREE}, but other units are allowed.
      *
@@ -726,22 +774,26 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
     public abstract void getParameterValues(Parameters parameters);
 
     /**
-     * Returns a string representation of this {@code DatumShiftGrid}. The default implementation
-     * formats the {@linkplain #getParameterValues(Parameters) parameter values}.
+     * Returns a string representation of this {@code DatumShiftGrid} for debugging purposes.
      *
-     * @return a string representation of the grid parameters.
+     * @return a string representation of this datum shift grid.
      *
      * @since 1.0
      */
     @Override
     public String toString() {
-        final ParameterDescriptorGroup d = getParameterDescriptors();
-        if (d != null) {
-            final Parameters p = Parameters.castOrWrap(d.createValue());
-            getParameterValues(p);
-            return p.toString();
+        final StringBuffer buffer = new StringBuffer("DatumShift[");
+        for (int i=0; i<gridSize.length; i++) {
+            if (i != 0) buffer.append(" × ");
+            buffer.append(gridSize[i]);
+        }
+        String s = String.valueOf(coordinateUnit);  if (s.isEmpty()) s = "1";
+        String t = String.valueOf(translationUnit); if (t.isEmpty()) t = "1";
+        buffer.append(" cells; units = ").append(s).append(" → ").append(t);
+        if (isCellValueRatio) {
+            buffer.append("∕cellSize");
         }
-        return super.toString();
+        return buffer.append(']').toString();
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
index 1abb8f0..412f7d7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
@@ -532,9 +532,10 @@ search: for (int j=numPoints; --j >= 0;) {
     }
 
     /**
-     * Sets all matching control point pairs, overwriting any previous setting. The source positions are the keys in
-     * the given map, and the target positions are the associated values in the map. The map should not contain two
-     * entries with the same source position. Coordinate reference systems are ignored.
+     * Sets all control point (source, target) pairs, overwriting any previous setting.
+     * The source positions are the keys in given map, and the target positions are the associated values.
+     * The map should not contain two entries with the same source position.
+     * Coordinate reference systems are ignored.
      * Null positions are silently ignored.
      * Positions with NaN or infinite coordinates cause an exception to be thrown.
      *
@@ -575,7 +576,7 @@ search: for (int j=numPoints; --j >= 0;) {
             final DirectPosition tgt = position(entry.getValue()); if (tgt == null) continue;
             /*
              * The first time that we get a non-null source and target coordinate, allocate the arrays.
-             * The sources arrays are allocated only if the source coordiantes are randomly distributed.
+             * The sources arrays are allocated only if the source coordinates are randomly distributed.
              */
             if (targets == null) {
                 tgtDim = tgt.getDimension();
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/AffineTransforms2D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/AffineTransforms2D.java
index a3441d4..2a40411 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/AffineTransforms2D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/AffineTransforms2D.java
@@ -38,7 +38,7 @@ import static java.awt.geom.AffineTransform.*;
  * Those {@code AffineTransform} instances can be viewed as 3×3 matrices.
  *
  * @author  Martin Desruisseaux (IRD, Geomatys)
- * @version 0.4
+ * @version 1.1
  * @since   0.4
  * @module
  */
@@ -414,4 +414,16 @@ public final class AffineTransforms2D extends Static {
         if (scale == 0) return abs(shear);                  // Not as common as above, but still common enough.
         return hypot(scale, shear);
     }
+
+    /**
+     * Returns a global scale factor for the specified affine transform. This scale factor combines
+     * {@link #getScaleX0 getScaleX0(tr)} and {@link #getScaleY0 getScaleY0(tr)}. The way to compute
+     * such a "global" scale is somewhat arbitrary and may change in any future version.
+     *
+     * @param  tr  the affine transform to inspect.
+     * @return a "global" scale factor.
+     */
+    public static double getScale(final AffineTransform tr) {
+        return 0.5 * (AffineTransforms2D.getScaleX0(tr) + AffineTransforms2D.getScaleY0(tr));
+    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
index 1a23ff2..4e62d0a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
@@ -311,13 +311,16 @@ public final class MathTransforms extends Static {
     public static MathTransform specialize(final MathTransform global, final Map<Envelope,MathTransform> specializations) {
         ArgumentChecks.ensureNonNull("generic", global);
         ArgumentChecks.ensureNonNull("specializations", specializations);
+        final SpecializableTransform tr;
         if (specializations.isEmpty()) {
             return global;
         } else if (global.getSourceDimensions() == 2 && global.getTargetDimensions() == 2) {
-            return new SpecializableTransform2D(global, specializations);
+            tr = new SpecializableTransform2D(global, specializations);
         } else {
-            return new SpecializableTransform(global, specializations);
+            tr = new SpecializableTransform(global, specializations);
         }
+        final MathTransform substitute = tr.getSubstitute();
+        return (substitute != null) ? substitute : tr;
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java
index 99110de..1d9e306 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java
@@ -18,24 +18,22 @@ package org.apache.sis.referencing.operation.transform;
 
 import java.util.Map;
 import java.util.List;
-import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.Objects;
+import java.util.Collections;
 import java.io.Serializable;
 import org.opengis.geometry.Envelope;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.geometry.MismatchedDimensionException;
-import org.apache.sis.geometry.MismatchedReferenceSystemException;
 import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
-import org.opengis.referencing.crs.CoordinateReferenceSystem;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.apache.sis.internal.referencing.DirectPositionView;
 import org.apache.sis.internal.referencing.Resources;
-import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.internal.referencing.RTreeNode;
 import org.apache.sis.io.wkt.Formatter;
-import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.ComparisonMode;
 import org.apache.sis.util.Utilities;
 
@@ -46,7 +44,7 @@ import org.apache.sis.util.Utilities;
  * The lower and upper values of given envelopes are inclusive.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see MathTransforms#specialize(MathTransform, Map)
  *
@@ -69,12 +67,12 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
      * Contains also a chain of {@code SubArea}s fully included in this area.
      * Shall be unmodified after {@link SpecializableTransform} construction.
      */
-    @SuppressWarnings("CloneableClassWithoutClone")                             // We will not use clone().
-    private static final class SubArea extends GeneralEnvelope {
+    @SuppressWarnings("CloneableImplementsClone")                               // We will not use clone().
+    private static final class SubArea extends RTreeNode {
         /**
          * For cross-version compatibility.
          */
-        private static final long serialVersionUID = 4197316795428796526L;
+        private static final long serialVersionUID = -7668993675003269862L;
 
         /**
          * The transform to apply in this area.
@@ -85,22 +83,11 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
          * The inverse of the transform, computed when first needed.
          * Synchronization for multi-threading is done (indirectly) in {@link SpecializableTransform#inverse()}.
          *
-         * @see #createInverse(SubArea)
+         * @see #createInverseTransform()
          */
         MathTransform inverse;
 
         /**
-         * Specialization, or {@code null} if none. If non-null, that sub-area shall be fully included
-         * in this {@code SubArea}. The specialization may itself contain another specialization. This
-         * form a chain from this wider area to smallest area, where each step is a smaller area.
-         *
-         * <p>Note that this is note a substitute for an R-Tree. This is an optimization for the common
-         * case where coordinates are close to each other (e.g. when iterating in a geometry), in which
-         * case we can check immediately for inclusion in smaller areas before to check the wider area.</p>
-         */
-        private SubArea specialization;
-
-        /**
          * Creates a new area where a transform is valid.
          */
         SubArea(final Envelope area, final MathTransform transform) {
@@ -109,97 +96,23 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
         }
 
         /**
-         * Tries to add a nested sub-area (a specialization of a specialization).
-         *
-         * @return whether the given area has been added.
-         */
-        boolean addSpecialization(final SubArea candidate) {
-            if (specialization == null) {
-                if (!contains(candidate)) {
-                    return false;
-                }
-            } else if (!candidate.addSpecialization(specialization)) {
-                return specialization.addSpecialization(candidate);
-            }
-            specialization = candidate;
-            return true;
-        }
-
-        /**
-         * Sets the CRS of all given ares to a common value. An exception is thrown if incompatible CRS are found.
-         * This method does not verify the number of dimensions; this check should have been done by the caller.
-         */
-        static void uniformize(final SubArea[] domains) {
-            CoordinateReferenceSystem common = null;
-            for (SubArea area : domains) {
-                do {
-                    final CoordinateReferenceSystem crs = area.getCoordinateReferenceSystem();
-                    if (common == null) {
-                        common = crs;
-                    } else if (crs != null && !Utilities.equalsIgnoreMetadata(common, crs)) {
-                        throw new MismatchedReferenceSystemException(Errors.format(Errors.Keys.MismatchedCRS));
-                    }
-                } while ((area = area.specialization) != null);
-            }
-            for (SubArea area : domains) {
-                do area.setCoordinateReferenceSystem(common);
-                while ((area = area.specialization) != null);
-            }
-        }
-
-        /**
          * Creates the inverse transforms. This method should be invoked only once when first needed
          * in a block synchronized (indirectly) by {@link SpecializableTransform#inverse()}.
          */
-        static void createInverse(SubArea area) throws NoninvertibleTransformException {
-            do {
-                area.inverse = area.transform.inverse();
-                area = area.specialization;
-            } while (area != null);
-        }
-
-        /**
-         * Returns the area that contains the given position, or {@code null} if none.
-         * This method may be replaced by an R-Tree in a future Apache SIS version.
-         */
-        static SubArea find(final SubArea[] domains, final DirectPosition pos) {
-            for (SubArea area : domains) {
-                if (area.contains(pos)) {
-                    SubArea next = area.specialization;
-                    while (next != null && next.contains(pos)) {
-                        area = next;
-                        next = next.specialization;
-                    }
-                    return area;
-                }
-            }
-            return null;
-        }
-
-        /**
-         * Returns the area that contains the given position, looking only in the given area or its specializations.
-         * Returns {@code null} if no area has been found.
-         */
-        static SubArea find(SubArea area, final DirectPosition pos) {
-            SubArea found = null;
-            while (area.contains(pos)) {
-                found = area;
-                area = area.specialization;
-                if (area == null) break;
+        final void createInverseTransform() throws NoninvertibleTransformException {
+            inverse = transform.inverse();
+            for (final RTreeNode node : getChildren()) {
+                ((SubArea) node).createInverseTransform();
             }
-            return found;
         }
 
         /**
-         * Formats the given area and its transform as a pseudo-WKT.
+         * Formats this area and its transform as a pseudo-WKT.
          * For {@link SpecializableTransform#formatTo(Formatter)} implementation only.
          */
-        static void format(SubArea area, final Formatter formatter) {
-            while (area != null) {
-                formatter.newLine(); formatter.append(area);
-                formatter.newLine(); formatter.append(area.transform);
-                area = area.specialization;
-            }
+        final void format(final Formatter formatter) {
+            formatter.newLine(); formatter.append(this);
+            formatter.newLine(); formatter.append(transform);
         }
 
         /**
@@ -207,11 +120,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
          */
         @Override
         public int hashCode() {
-            int code = super.hashCode() ^ transform.hashCode();
-            if (specialization != null) {
-                code += 37 * specialization.hashCode();
-            }
-            return code;
+            return super.hashCode() ^ transform.hashCode();
         }
 
         /**
@@ -219,34 +128,18 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
          */
         @Override
         public boolean equals(final Object obj) {
-            if (super.equals(obj)) {
-                final SubArea other = (SubArea) obj;
-                return transform.equals(other.transform) && Objects.equals(specialization, other.specialization);
-            }
-            return false;
-        }
-
-        /**
-         * Formats this envelope as a "{@code DOMAIN}" element (non-standard).
-         * This is used by {@link SpecializableTransform#formatTo(Formatter)}.
-         */
-        @Override
-        protected String formatTo(final Formatter formatter) {
-            super.formatTo(formatter);
-            return "Domain";
+            return super.equals(obj) && transform.equals(((SubArea) obj).transform);
         }
     }
 
     /**
-     * Domains where specialized transforms are valid. The array should be very small.
-     * In current implementation, elements in this array shall not overlap.
-     * This array may be replaced by an R-Tree in a future Apache SIS version.
+     * Domains where specialized transforms are valid. This is the root of an R-Tree.
      */
-    private final SubArea[] domains;
+    private final RTreeNode domains;
 
     /**
      * The inverse of this transform, computed when first needed.
-     * Part of serialization for avoiding rounding error issues.
+     * This object is included in serialization for avoiding rounding error issues.
      *
      * @see #inverse()
      */
@@ -260,36 +153,42 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
      */
     SpecializableTransform(final MathTransform global, final Map<Envelope,MathTransform> specializations) {
         this.global = global;
+        SubArea root = null;
         final int sourceDim = global.getSourceDimensions();
         final int targetDim = global.getTargetDimensions();
-        final List<SubArea> areas = new ArrayList<>(specializations.size());
         for (final Map.Entry<Envelope,MathTransform> entry : specializations.entrySet()) {
             MathTransform tr = entry.getValue();
             ensureDimensionMatches(0, sourceDim, tr.getSourceDimensions());
             ensureDimensionMatches(1, targetDim, tr.getTargetDimensions());
-            SubArea[] inherited = null;
+            /*
+             * If the given MathTransform is another SpecializableTransform, then instead of storing nested
+             * SpecializableTransforms we will store directly the specializations that it contains. It will
+             * reduce the amountof steps when transforming coordinates.
+             */
+            List<SubArea> inherited = Collections.emptyList();
             if (tr instanceof SpecializableTransform) {
-                inherited = ((SpecializableTransform) tr).domains;
-                tr = ((SpecializableTransform) tr).global;
+                inherited = ((SpecializableTransform) tr).roots();
+                tr        = ((SpecializableTransform) tr).global;
             }
             final SubArea area = new SubArea(entry.getKey(), tr);
-            addSpecialization(area, areas, sourceDim);
-            /*
-             * At this point we are usually done for the current SubArea. But if the given MathTransform
-             * is another SpecializableTransform, then instead of storing nested SpecializableTransforms
-             * we will store directly the specializations that it contains.  This will reduce the amount
-             * of steps when transforming coordinates.
-             */
-            if (inherited != null) {
+            ArgumentChecks.ensureDimensionMatches("envelope", sourceDim, area);
+            if (!area.isEmpty()) {
+                if (root == null) root = area;
+                else root.addNode(area);
+                /*
+                 * If the transform was another SpecializableTransform, copies the nested RTreeNode.
+                 * A copy is necessary: we shall not modify the nodes of the given transform.
+                 */
                 for (final SubArea other : inherited) {
                     final SubArea e = new SubArea(other, other.transform);
                     e.intersect(area);
-                    addSpecialization(e, areas, sourceDim);
+                    if (!e.isEmpty()) {
+                        root.addNode(e);
+                    }
                 }
             }
         }
-        domains = areas.toArray(new SubArea[areas.size()]);
-        SubArea.uniformize(domains);
+        domains = (root != null) ? root.finish() : null;
     }
 
     /**
@@ -305,35 +204,35 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
     }
 
     /**
-     * Verifies if the given {@code area} has the expected number of dimensions,
-     * then adds it to {@code domains} list (eventually as a child of an existing node).
-     *
-     * @param  area     the new sub-area to add.
-     * @param  domains  where to add the sub-area (not necessarily directly; maybe as a child of an existing node).
-     * @param  dim      expected number of dimensions, for verification purpose.
+     * Returns the {@link SubArea} instances at the root of this class. This is the {@link #domains} node,
+     * unless that node is a synthetic node created by {@link RTreeNode} when it needs to contain more than
+     * one children.
      */
-    private static void addSpecialization(final SubArea area, final List<SubArea> domains, final int dim) {
-        if (!area.isEmpty()) {
-            if (area.getDimension() != dim) {
-                throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_3,
-                            "envelope", dim, area.getDimension()));
-            }
-            for (final SubArea previous : domains) {
-                if (previous.addSpecialization(area)) {
-                    return;
-                }
-            }
-            for (final SubArea previous : domains) {
-                if (area.intersects(previous)) {
-                    // Pending implementation of R-Tree in Apache SIS.
-                    throw new IllegalArgumentException("Current implementation does not accept overlapping envelopes.");
-                }
-            }
-            domains.add(area);
+    @SuppressWarnings("unchecked")
+    private List<SubArea> roots() {
+        if (domains == null) {
+            return Collections.emptyList();
+        } else if (domains instanceof SubArea) {
+            return Collections.singletonList((SubArea) domains);
+        } else {
+            /*
+             * We are cheating here since we have a `List<RTreeNode>`. But this `SpecializableTransform`
+             * class adds only `SubArea` instance in this list, so a ClassCastException in the code that
+             * use this list would be a bug in our conditions of `RTreeNode` use.
+             */
+            return (List) domains.getChildren();
         }
     }
 
     /**
+     * If this transform has no children, then returns the transform that we should use instead.
+     * Otherwise returns {@code null}.
+     */
+    final MathTransform getSubstitute() {
+        return (domains == null) ? global : null;
+    }
+
+    /**
      * Gets the dimension of input points.
      */
     @Override
@@ -350,10 +249,20 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
     }
 
     /**
-     * Returns the transform to use for the given domain.
+     * Returns the node that contains the given position, or {@code null} if none.
+     * This method searches from the root of the tree. It should be invoked only when
+     * we do not know what was the last search result.
      */
-    private MathTransform forDomain(final SubArea domain) {
-        return (domain != null) ? domain.transform : global;
+    private SubArea locate(final DirectPosition pos) {
+        return (SubArea) RTreeNode.locate(domains, pos);
+    }
+
+    /**
+     * Returns the transform to use for the given position, or {@link #global} if none.
+     */
+    private MathTransform forDomain(final DirectPosition pos) {
+        final RTreeNode domain = RTreeNode.locate(domains, pos);
+        return (domain != null) ? ((SubArea) domain).transform : global;
     }
 
     /**
@@ -362,7 +271,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
      */
     @Override
     public final DirectPosition transform(final DirectPosition ptSrc, DirectPosition ptDst) throws TransformException {
-        return forDomain(SubArea.find(domains, ptSrc)).transform(ptSrc, ptDst);
+        return forDomain(ptSrc).transform(ptSrc, ptDst);
     }
 
     /**
@@ -371,7 +280,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
      */
     @Override
     public final Matrix derivative(final DirectPosition point) throws TransformException {
-        return forDomain(SubArea.find(domains, point)).derivative(point);
+        return forDomain(point).derivative(point);
     }
 
     /**
@@ -384,7 +293,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
                                   boolean derivate) throws TransformException
     {
         final DirectPositionView pos = new DirectPositionView.Double(srcPts, srcOff, global.getSourceDimensions());
-        final MathTransform tr = forDomain(SubArea.find(domains, pos));
+        final MathTransform tr = forDomain(pos);
         if (tr instanceof AbstractMathTransform) {
             return ((AbstractMathTransform) tr).transform(srcPts, srcOff, dstPts, dstOff, derivate);
         } else {
@@ -422,28 +331,26 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
             int dstOff, int srcInc, int dstInc, int numPts) throws TransformException
     {
         final boolean downard = (srcInc < 0);
-        SubArea domain = SubArea.find(domains, src);
+        SubArea domain = locate(src);
         while (numPts > 0) {
             int srcOff = src.offset;
             final MathTransform tr;
             if (domain == null) {
-                tr = global;                               // The transform to apply when no specialization is found.
+                tr = global;                                // The transform to apply when no specialization is found.
                 do {                                        // Count how many points will use that transform.
                     src.offset += srcInc;
                     if (--numPts <= 0) break;
-                    domain = SubArea.find(domains, src);    // More expansive check than the case where domain is non-null.
+                    domain = locate(src);                   // More expansive check than the case where domain is non-null.
                 } while (domain == null);
             } else {
-                final SubArea previous = domain;
+                RTreeNode next = domain;
                 tr = domain.transform;                      // The specialized transform to apply.
                 do {                                        // Count how many points will use that transform.
                     src.offset += srcInc;
                     if (--numPts <= 0) break;
-                    domain = SubArea.find(domain, src);     // Cheaper check compared to the case where domain is null.
-                } while (domain == previous);
-                if (domain == null) {
-                    domain = SubArea.find(domains, src);    // Need to update with the more expansive check.
-                }
+                    next = SubArea.locate(domain, src);     // Cheaper check compared to the case where domain is null.
+                } while (next == domain);
+                domain = (SubArea) next;
             }
             final int num = (src.offset - srcOff) / srcInc;
             int dstLow = dstOff;
@@ -560,7 +467,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
      */
     @Override
     protected final int computeHashCode() {
-        return super.computeHashCode() + 7*global.hashCode() ^ Arrays.hashCode(domains);
+        return super.computeHashCode() + 7*global.hashCode() ^ domains.hashCode();
     }
 
     /**
@@ -571,7 +478,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
         if (super.equals(object, mode)) {
             final SpecializableTransform other = (SpecializableTransform) object;
             return Utilities.deepEquals(global, other.global, mode) &&
-                   Utilities.deepEquals(domains, other.domains, mode);
+                   Objects.equals(domains, other.domains);
         }
         return false;
     }
@@ -589,9 +496,11 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
     protected final String formatTo(final Formatter formatter) {
         formatter.newLine();
         formatter.append(global);
-        for (SubArea domain : domains) {
-            SubArea.format(domain, formatter);
-        }
+        RTreeNode.walk(domains, (node) -> {
+            if (node instanceof SubArea) {
+                ((SubArea) node).format(formatter);
+            }
+        });
         formatter.setInvalidWKT(SpecializableTransform.class, null);
         return "Specializable_MT";
     }
@@ -640,8 +549,8 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
         Inverse(final SpecializableTransform forward) throws NoninvertibleTransformException {
             this.forward = forward;
             this.global = forward.global.inverse();
-            for (final SubArea domain : forward.domains) {
-                SubArea.createInverse(domain);
+            for (final SubArea domain : forward.roots()) {
+                domain.createInverseTransform();
             }
         }
 
@@ -660,7 +569,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
         public final DirectPosition transform(final DirectPosition ptSrc, DirectPosition ptDst) throws TransformException {
             final double[] source = ptSrc.getCoordinate();      // Needs to be first in case ptDst overwrites ptSrc.
             ptDst = global.transform(ptSrc, ptDst);
-            final SubArea domain = SubArea.find(forward.domains, ptDst);
+            final SubArea domain = forward.locate(ptDst);
             if (domain != null) {
                 ptDst = domain.inverse.transform(new DirectPositionView.Double(source), ptDst);
             }
@@ -709,7 +618,7 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
                     derivative = derivate ? tr.derivative(new DirectPositionView.Double(srcPts, srcOff, srcInc)) : null;
                 }
                 if (secondTry) break;
-                final SubArea domain = SubArea.find(forward.domains, new DirectPositionView.Double(dstPts, dstOff, dstInc));
+                final SubArea domain = forward.locate(new DirectPositionView.Double(dstPts, dstOff, dstInc));
                 if (domain != null) {
                     tr = domain.inverse;
                     secondTry = true;
@@ -726,18 +635,17 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
         private void transform(final TransformCall transform, final double[] dstPts,
                 int srcOff, int dstOff, int srcInc, int dstInc, int numPts) throws TransformException
         {
-            final SubArea[] domains = forward.domains;
             transform.apply(global, srcOff, dstOff, numPts);
             final DirectPositionView dst = new DirectPositionView.Double(dstPts, dstOff, dstInc);
             while (numPts > 0) {
-                SubArea domain = SubArea.find(domains, dst);
+                SubArea domain = forward.locate(dst);
                 if (domain == null) {
                     dst.offset += dstInc;
                     numPts--;
                     continue;
                 }
                 do {
-                    final SubArea specialized = domain;             // Contains the specialized transform to use.
+                    RTreeNode next = domain;                        // Contains the specialized transform to use.
                     int num = (dst.offset - dstOff) / dstInc;       // Number of points that are not retransformeD.
                     srcOff += num * srcInc;                         // Skip the source coordinates that are not retransformed.
                     dstOff = dst.offset;                            // Destination index of the first coordinate to retransform.
@@ -747,10 +655,11 @@ class SpecializableTransform extends AbstractMathTransform implements Serializab
                             domain = null;
                             break;
                         }
-                        domain = SubArea.find(domain, dst);
-                    } while (domain == specialized);
+                        next = RTreeNode.locate(domain, dst);
+                    } while (next == domain);
                     num = (dst.offset - dstOff) / dstInc;           // Number of points to retransform.
-                    transform.apply(specialized.inverse, srcOff, dstOff, num);
+                    transform.apply(domain.inverse, srcOff, dstOff, num);
+                    domain = (SubArea) next;
                     srcOff += srcInc * num;
                     dstOff = dst.offset;
                 } while (domain != null);
diff --git a/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod b/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
index 5aa20e9..ae5ddfe 100644
--- a/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
+++ b/core/sis-referencing/src/main/resources/META-INF/services/org.opengis.referencing.operation.OperationMethod
@@ -60,6 +60,7 @@ org.apache.sis.internal.referencing.provider.Sinusoidal
 org.apache.sis.internal.referencing.provider.Polyconic
 org.apache.sis.internal.referencing.provider.Mollweide
 org.apache.sis.internal.referencing.provider.NTv2
+org.apache.sis.internal.referencing.provider.NTv1
 org.apache.sis.internal.referencing.provider.NADCON
 org.apache.sis.internal.referencing.provider.FranceGeocentricInterpolation
 org.apache.sis.internal.referencing.provider.Interpolation1D
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolationTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolationTest.java
index 145c6c5..895c263 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolationTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/FranceGeocentricInterpolationTest.java
@@ -182,10 +182,10 @@ public final strictfp class FranceGeocentricInterpolationTest extends DatumShift
      */
     private static void verifyGrid(final DatumShiftGridFile<Angle,Length> grid) throws TransformException {
         final Envelope envelope = grid.getDomainOfValidity();
-        assertEquals("xmin",  2.2 - 0.05, envelope.getMinimum(0), 1E-12);
-        assertEquals("xmax",  2.5 + 0.05, envelope.getMaximum(0), 1E-12);
-        assertEquals("ymin", 48.5 - 0.05, envelope.getMinimum(1), 1E-12);
-        assertEquals("ymax", 49.0 + 0.05, envelope.getMaximum(1), 1E-12);
+        assertEquals("xmin",  2.2, envelope.getMinimum(0), 1E-12);
+        assertEquals("xmax",  2.5, envelope.getMaximum(0), 1E-12);
+        assertEquals("ymin", 48.5, envelope.getMinimum(1), 1E-12);
+        assertEquals("ymax", 49.0, envelope.getMaximum(1), 1E-12);
         /*
          * The values in the NTG_88 document are:
          *
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NADCONTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NADCONTest.java
index 0bec535..a9b34d2 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NADCONTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NADCONTest.java
@@ -93,12 +93,11 @@ public final strictfp class NADCONTest extends DatumShiftTestCase {
      * The point used for this test is given by {@link #samplePoint(int)}.
      *
      * @throws URISyntaxException if the URL to the test file can not be converted to a path.
-     * @throws IOException if an error occurred while loading the grid.
-     * @throws FactoryException if an error occurred while computing the grid.
+     * @throws FactoryException if an error occurred while loading or computing the grid.
      * @throws TransformException if an error occurred while computing the envelope or testing the point.
      */
     @Test
-    public void testLoader() throws URISyntaxException, IOException, FactoryException, TransformException {
+    public void testLoader() throws URISyntaxException, FactoryException, TransformException {
         testNADCON(getResource(TEST_FILE + ".laa"),     // Latitude shifts
                    getResource(TEST_FILE + ".loa"),     // Longitude shifts
                    -99.75, -98.0, 37.5, 39.75);
@@ -112,12 +111,11 @@ public final strictfp class NADCONTest extends DatumShiftTestCase {
      *
      * @param  latitudeShifts   path to the official {@code "conus.las"} file.
      * @param  longitudeShifts  path to the official {@code "conus.los"} file.
-     * @throws IOException if an error occurred while loading the grid.
-     * @throws FactoryException if an error occurred while computing the grid.
+     * @throws FactoryException if an error occurred while loading or computing the grid.
      * @throws TransformException if an error occurred while computing the envelope or testing the point.
      */
     public static void testNADCON(final Path latitudeShifts, final Path longitudeShifts)
-            throws IOException, FactoryException, TransformException
+            throws FactoryException, TransformException
     {
         testNADCON(latitudeShifts, longitudeShifts, -131, -63, 20, 50);
     }
@@ -132,7 +130,7 @@ public final strictfp class NADCONTest extends DatumShiftTestCase {
      */
     private static void testNADCON(final Path latitudeShifts, final Path longitudeShifts,
             final double xmin, final double xmax, final double ymin, final double ymax)
-            throws IOException, FactoryException, TransformException
+            throws FactoryException, TransformException
     {
         final DatumShiftGridFile<Angle,Angle> grid = NADCON.getOrLoad(latitudeShifts, longitudeShifts);
         assertInstanceOf("Should not be compressed.", DatumShiftGridFile.Float.class, grid);
@@ -147,10 +145,10 @@ public final strictfp class NADCONTest extends DatumShiftTestCase {
          */
         final double cellSize = 0.25;
         final Envelope envelope = grid.getDomainOfValidity();
-        assertEquals("xmin", xmin - cellSize/2, envelope.getMinimum(0), STRICT);
-        assertEquals("xmax", xmax + cellSize/2, envelope.getMaximum(0), STRICT);
-        assertEquals("ymin", ymin - cellSize/2, envelope.getMinimum(1), STRICT);
-        assertEquals("ymax", ymax + cellSize/2, envelope.getMaximum(1), STRICT);
+        assertEquals("xmin", xmin, envelope.getMinimum(0), STRICT);
+        assertEquals("xmax", xmax, envelope.getMaximum(0), STRICT);
+        assertEquals("ymin", ymin, envelope.getMinimum(1), STRICT);
+        assertEquals("ymax", ymax, envelope.getMaximum(1), STRICT);
         assertMatrixEquals("coordinateToGrid",
                 new Matrix3(cellSize,  0,  xmin,
                             0,  cellSize,  ymin,
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
index c87ef93..53f6954 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
@@ -21,6 +21,7 @@ import java.io.IOException;
 import java.nio.ByteBuffer;
 import java.nio.ByteOrder;
 import java.nio.file.Path;
+import java.nio.file.Paths;
 import java.nio.file.Files;
 import java.nio.file.StandardOpenOption;
 import java.nio.channels.WritableByteChannel;
@@ -33,19 +34,29 @@ import org.apache.sis.referencing.operation.matrix.Matrix3;
 import org.apache.sis.geometry.Envelope2D;
 import org.apache.sis.geometry.Envelopes;
 import org.apache.sis.measure.Units;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.internal.system.DataDirectory;
+import org.apache.sis.test.DependsOn;
 import org.junit.Test;
 
+import static org.junit.Assume.assumeTrue;
 import static org.apache.sis.test.Assert.*;
+import static org.apache.sis.internal.referencing.provider.DatumShiftGridLoader.DEGREES_TO_SECONDS;
 
 
 /**
  * Tests the {@link NTv2} grid loader.
+ * It will also indirectly tests {@link DatumShiftGridGroup} class.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.1
+ *
+ * @see org.apache.sis.referencing.operation.transform.MolodenskyTransformTest#testFranceGeocentricInterpolationPoint()
+ *
  * @since 0.7
  * @module
  */
+@DependsOn(DatumShiftGridFileTest.class)
 public final strictfp class NTv2Test extends DatumShiftTestCase {
     /**
      * Name of the file containing a small extract of the "{@code NTF_R93.gsb}" file.
@@ -54,6 +65,13 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
     public static final String TEST_FILE = "NTF_R93-extract.gsb";
 
     /**
+     * Name of the file to load for testing the multi-grids support.
+     * This file should be present in the {@code $SIS_DATA/DatumChanges} directory.
+     * The test will be skipped if this file is absent.
+     */
+    private static final String MULTIGRID_TEST_FILE = "NTv2_0.gsb";
+
+    /**
      * Best accuracy found in the "{@code NTF_R93.gsb}" file.
      */
     private static final float ACCURACY = 0.001618f;
@@ -83,11 +101,10 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
      * explicitly if they can provide a path to the {@code "NTF_R93.gsb"} file.
      *
      * @param  file  path to the official {@code "NTF_R93.gsb"} file.
-     * @throws IOException if an error occurred while loading the grid.
-     * @throws FactoryException if an error occurred while computing the grid.
+     * @throws FactoryException if an error occurred while loading or computing the grid.
      * @throws TransformException if an error occurred while computing the envelope or testing the point.
      */
-    public static void testRGF93(final Path file) throws IOException, FactoryException, TransformException {
+    public static void testRGF93(final Path file) throws FactoryException, TransformException {
         testRGF93(file, -19800, 36000, 147600, 187200);
     }
 
@@ -100,25 +117,25 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
      * @param  ymax  value of the {@code "N_LAT"} record.
      */
     private static void testRGF93(final Path file, final double xmin, final double xmax,
-            final double ymin, final double ymax) throws IOException, FactoryException, TransformException
+            final double ymin, final double ymax) throws FactoryException, TransformException
     {
         final double cellSize = 360;
-        final DatumShiftGridFile<Angle,Angle> grid = NTv2.getOrLoad(file);
+        final DatumShiftGridFile<Angle,Angle> grid = NTv2.getOrLoad(NTv2.class, file, 2);
         assertInstanceOf("Should not be compressed.", DatumShiftGridFile.Float.class, grid);
         assertEquals("coordinateUnit",  Units.ARC_SECOND, grid.getCoordinateUnit());
         assertEquals("translationUnit", Units.ARC_SECOND, grid.getTranslationUnit());
-        assertEquals("translationDimensions", 2, grid.getTranslationDimensions());
-        assertTrue  ("isCellValueRatio", grid.isCellValueRatio());
+        assertEquals("translationDimensions", 2,          grid.getTranslationDimensions());
+        assertTrue  ("isCellValueRatio",                  grid.isCellValueRatio());
         assertEquals("cellPrecision", (ACCURACY / 10) / cellSize, grid.getCellPrecision(), 0.5E-6 / cellSize);
         /*
          * Verify the envelope and the conversion between geographic coordinates and grid indices.
          * The cells are expected to have the same size (360″ or 0.1°) in longitudes and latitudes.
          */
         final Envelope envelope = grid.getDomainOfValidity();
-        assertEquals("xmin", xmin - cellSize/2, envelope.getMinimum(0), STRICT);
-        assertEquals("xmax", xmax + cellSize/2, envelope.getMaximum(0), STRICT);
-        assertEquals("ymin", ymin - cellSize/2, envelope.getMinimum(1), STRICT);
-        assertEquals("ymax", ymax + cellSize/2, envelope.getMaximum(1), STRICT);
+        assertEquals("xmin", xmin, envelope.getMinimum(0), STRICT);
+        assertEquals("xmax", xmax, envelope.getMaximum(0), STRICT);
+        assertEquals("ymin", ymin, envelope.getMinimum(1), STRICT);
+        assertEquals("ymax", ymax, envelope.getMaximum(1), STRICT);
         assertMatrixEquals("coordinateToGrid",
                 new Matrix3(-cellSize,  0,  xmax,
                             0,  +cellSize,  ymin,
@@ -138,8 +155,8 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
         final double[] indices  = new double[position.length];
         final double[] vector   = new double[2];
         for (int i=0; i<expected.length; i++) {
-            position[i] *= DatumShiftGridLoader.DEGREES_TO_SECONDS;
-            expected[i] *= DatumShiftGridLoader.DEGREES_TO_SECONDS;
+            position[i] *= DEGREES_TO_SECONDS;
+            expected[i] *= DEGREES_TO_SECONDS;
             expected[i] -= position[i];  // We will test the interpolated shifts rather than final coordinates.
         }
         grid.getCoordinateToGrid().transform(position, 0, indices, 0, 1);
@@ -147,12 +164,81 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
         vector[0] *= -cellSize;   // Was positive toward west.
         vector[1] *= +cellSize;
         assertArrayEquals("interpolateInCell", expected, vector,
-                FranceGeocentricInterpolationTest.ANGULAR_TOLERANCE * DatumShiftGridLoader.DEGREES_TO_SECONDS);
+                FranceGeocentricInterpolationTest.ANGULAR_TOLERANCE * DEGREES_TO_SECONDS);
 
         // Same test than above, but let DatumShiftGrid do the conversions for us.
         assertArrayEquals("interpolateAt", expected, grid.interpolateAt(position),
-                FranceGeocentricInterpolationTest.ANGULAR_TOLERANCE * DatumShiftGridLoader.DEGREES_TO_SECONDS);
-        assertSame("Grid should be cached.", grid, NTv2.getOrLoad(file));
+                FranceGeocentricInterpolationTest.ANGULAR_TOLERANCE * DEGREES_TO_SECONDS);
+        assertSame("Grid should be cached.", grid, NTv2.getOrLoad(NTv2.class, file, 2));
+    }
+
+    /**
+     * Tests using a file containing many grids. This tests depends on the {@value #MULTIGRID_TEST_FILE}
+     * to be present in the {@code $SIS_DATA/DatumChanges} directory. This test is executed only if the
+     * {@link #RUN_EXTENSIVE_TESTS} flag is set.
+     *
+     * @throws FactoryException if an error occurred while loading or computing the grid.
+     * @throws TransformException if an error occurred while computing the envelope or testing the point.
+     */
+    @Test
+    public void testMultiGrids() throws FactoryException, TransformException {
+        assumeTrue(RUN_EXTENSIVE_TESTS);
+        final Path file = DataDirectory.DATUM_CHANGES.resolve(Paths.get(MULTIGRID_TEST_FILE));
+        assumeTrue(Files.exists(file));
+        final DatumShiftGridFile<Angle,Angle> grid = NTv2.getOrLoad(NTv2.class, file, 2);
+        assertInstanceOf("Should contain many grids.", DatumShiftGridGroup.class, grid);
+        assertEquals("coordinateUnit",  Units.ARC_SECOND, grid.getCoordinateUnit());
+        assertEquals("translationUnit", Units.ARC_SECOND, grid.getTranslationUnit());
+        assertEquals("translationDimensions", 2,          grid.getTranslationDimensions());
+        assertTrue  ("isCellValueRatio",                  grid.isCellValueRatio());
+        /*
+         * Area of use declared in EPSG database for coordinate operation EPSG::1693:
+         *
+         *     40.04°N  to  83.17°N    and    141.01°W  to  47.74°W.
+         *
+         * In the assertions below, the `cellSize` value has been verified in the NTv2
+         * file header but the envelope bounds have been determined empirically.
+         */
+        final double cellSize = 300;                                    // Number of arc-seconds in a cell.
+        final Envelope envelope = grid.getDomainOfValidity();
+        assertEquals("xmin", -142.25 * DEGREES_TO_SECONDS, envelope.getMinimum(0), 1E-10);
+        assertEquals("xmax",  -44.00 * DEGREES_TO_SECONDS, envelope.getMaximum(0), 1E-10);
+        assertEquals("ymin",   40.00 * DEGREES_TO_SECONDS, envelope.getMinimum(1), 1E-10);
+        assertEquals("ymax",   84.00 * DEGREES_TO_SECONDS, envelope.getMaximum(1), 1E-10);
+        /*
+         * Test a point. This point is located on the 3th grid in the NTv2 file.
+         * Consequently if the NTv2 implementation just pickups the first grid,
+         * then this test would fail with an error around 100 metres.
+         */
+        final double[] position = {-134.998106062 * DEGREES_TO_SECONDS, 61.000285047 * DEGREES_TO_SECONDS};
+        final double[] expected = {-135.0         * DEGREES_TO_SECONDS, 61.0         * DEGREES_TO_SECONDS};
+        final double[] indices  = new double[position.length];
+        grid.getCoordinateToGrid().transform(position, 0, indices, 0, 1);
+        final int gridX = Math.toIntExact(Math.round(indices[0]));
+        final int gridY = Math.toIntExact(Math.round(indices[1]));
+        assertEquals("gridX", 1092, gridX);                                 // Value determined empirically.
+        assertEquals("gridY",  252, gridY);
+        /*
+         * First check the value computed by `getCellValue(…)`. This method is only a fallback and
+         * should not be invoked in normal usage, so a direct invocation is the only way to test it.
+         */
+        final double[] result = new double[] {
+            position[0] - grid.getCellValue(0, gridX, gridY) * cellSize,    // Positive translation is toward west.
+            position[1] + grid.getCellValue(1, gridX, gridY) * cellSize
+        };
+        assertArrayEquals("getCellValue", expected, result, 0.001);
+        /*
+         * Check `interpolateInCell(…)`, which is the method invoked by `InterpolatedTransform`
+         * when `SpecializableTransform` has not been able to find the most appropriate grid.
+         */
+        grid.interpolateInCell(indices[0], indices[1], result);
+        result[0] = position[0] - result[0] * cellSize;                     // Positive translation is toward west.
+        result[1] = position[1] + result[1] * cellSize;
+        assertArrayEquals("interpolateInCell", expected, result, Formulas.ANGULAR_TOLERANCE * DEGREES_TO_SECONDS);
+        /*
+         * Verify that the caching mechanism works for DatumShiftGridGroup too.
+         */
+        assertSame("Grid should be cached.", grid, NTv2.getOrLoad(NTv2.class, file, 2));
     }
 
 
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
index 25e3b8b..12002e6 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/ProvidersTest.java
@@ -111,6 +111,7 @@ public final strictfp class ProvidersTest extends TestCase {
             Polyconic.class,
             Mollweide.class,
             NTv2.class,
+            NTv1.class,
             NADCON.class,
             FranceGeocentricInterpolation.class,
             MolodenskyInterpolation.class,
diff --git a/core/sis-utility/src/main/java/org/apache/sis/internal/util/Strings.java b/core/sis-utility/src/main/java/org/apache/sis/internal/util/Strings.java
index 9558156..1167d1f 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/internal/util/Strings.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/internal/util/Strings.java
@@ -30,7 +30,7 @@ import org.apache.sis.util.CharSequences;
  * Most of those methods are for {@link Object#toString()} implementations.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.3
  * @module
  */
@@ -54,6 +54,36 @@ public final class Strings extends Static {
     }
 
     /**
+     * Trims the leading and trailing spaces of the given string.
+     * If the string is null, empty or contains only spaces, then this method returns {@code null}.
+     *
+     * <p>Note that this method strips white spaces too, including no-break spaces.
+     * In some cases this is not wanted, for example if the text is a programmatic identifier
+     * (maybe the developer really wanted no-break spaces). To preserve no-break spaces, the
+     * following can be used instead:</p>
+     *
+     * {@preformat java
+     *     if (text != null && !(text = text.trim()).isEmpty()) {
+     *         // Use text here.
+     *     }
+     * }
+     *
+     * @param  text  the text to trim, or {@code null}.
+     * @return the trimmed text, or {@code null} if the given text was null or blank.
+     *
+     * @since 1.1
+     */
+    public static String trimOrNull(String text) {
+        if (text != null) {
+            text = CharSequences.trimWhitespaces(text.trim());
+            if (text.isEmpty()) {
+                return null;
+            }
+        }
+        return text;
+    }
+
+    /**
      * Appends to the given buffer only the characters that are valid for a Unicode identifier.
      * The given separator character is append before the given {@code text} only if the buffer
      * is not empty and at least one {@code text} character is valid.
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 351b915..36f4bf6 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
@@ -577,6 +577,11 @@ public final class Vocabulary extends IndexedResourceBundle {
         public static final short ModifiedJulian = 71;
 
         /**
+         * … {0} more…
+         */
+        public static final short More_1 = 197;
+
+        /**
          * Multiplicity
          */
         public static final short Multiplicity = 147;
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 15af0b8..1a67b2c 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
@@ -118,6 +118,7 @@ MissingValue            = Missing value
 Measures                = Measures
 Methods                 = Methods
 ModifiedJulian          = Modified Julian
+More_1                  = \u2026 {0} more\u2026
 Multiplicity            = Multiplicity
 Name                    = Name
 Nodata                  = No data
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 cb9f195..054f01d 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
@@ -125,6 +125,7 @@ MissingValue            = Valeur manquante
 Measures                = Mesures
 Methods                 = M\u00e9thodes
 ModifiedJulian          = Julien modifi\u00e9
+More_1                  = \u2026 {0} de plus\u2026
 Multiplicity            = Multiplicit\u00e9
 Name                    = Nom
 Nodata                  = Absence de donn\u00e9es
diff --git a/core/sis-utility/src/test/java/org/apache/sis/test/Assume.java b/core/sis-utility/src/test/java/org/apache/sis/test/Assume.java
index b92957f..7573214 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/test/Assume.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/test/Assume.java
@@ -54,7 +54,7 @@ public final strictfp class Assume extends org.junit.Assume {
         Path path = type.getDirectory();
         assumeNotNull("$SIS_DATA/" + type + " directory not found.", path);
         path = path.resolve(file);
-        assumeTrue("Specified directory not found.", Files.isDirectory(path));
+        assumeTrue("Specified file or directory not found.", Files.exists(path));
         assumeTrue("Specified directory not readable.", Files.isReadable(path));
         return path;
     }
diff --git a/storage/sis-netcdf/src/test/java/org/apache/sis/internal/netcdf/SatelliteGroundTrackTest.java b/storage/sis-netcdf/src/test/java/org/apache/sis/internal/netcdf/SatelliteGroundTrackTest.java
index c2fdad6..fd39239 100644
--- a/storage/sis-netcdf/src/test/java/org/apache/sis/internal/netcdf/SatelliteGroundTrackTest.java
+++ b/storage/sis-netcdf/src/test/java/org/apache/sis/internal/netcdf/SatelliteGroundTrackTest.java
@@ -47,7 +47,7 @@ public final strictfp class SatelliteGroundTrackTest extends MathTransformTestCa
      */
     private void createTransform() throws TransformException, FactoryException {
         final LocalizationGridBuilder grid = new LocalizationGridBuilder(WIDTH, HEIGHT);
-        final AffineTransform tr = AffineTransform.getRotateInstance(StrictMath.random()/2 + 0.25);     // Between 14 and 43°.
+        final AffineTransform tr = AffineTransform.getRotateInstance(StrictMath.toRadians(31));
         tr.translate(-WIDTH / 2, -HEIGHT / 2);
         final double[] point = new double[2];
         for (int y=0; y<HEIGHT; y++) {


Mime
View raw message