sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 05/05: Read GeoTIFF sample values from a potentially n-dimensional data cube. Whether the file has 3 or more dimensions is assumed already known; this commit does not yet include a convention for specifying that.
Date Tue, 06 Jul 2021 15:39:57 GMT
This is an automated email from the ASF dual-hosted git repository.

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

commit bf4cb0ccce600bb449fdae8dcb49993f54fa77d5
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Tue May 4 20:55:49 2021 +0200

    Read GeoTIFF sample values from a potentially n-dimensional data cube.
    Whether the file has 3 or more dimensions is assumed already known;
    this commit does not yet include a convention for specifying that.
---
 .../java/org/apache/sis/image/BandSelectImage.java |   2 +-
 .../java/org/apache/sis/image/PlanarImage.java     |   8 +-
 .../internal/coverage/j2d/ColorModelFactory.java   |  82 ++-
 .../sis/internal/coverage/j2d/Colorizer.java       |   2 +-
 .../sis/internal/coverage/j2d/TiledImage.java      |  37 +-
 .../org/apache/sis/internal/geotiff/Resources.java |  15 +
 .../sis/internal/geotiff/Resources.properties      |   3 +
 .../sis/internal/geotiff/Resources_fr.properties   |   3 +
 .../org/apache/sis/storage/geotiff/CRSBuilder.java |   8 +-
 .../org/apache/sis/storage/geotiff/DataCube.java   | 151 +++++
 .../org/apache/sis/storage/geotiff/DataSubset.java | 297 +++++++++
 .../org/apache/sis/storage/geotiff/GeoTIFF.java    |  10 +-
 .../apache/sis/storage/geotiff/GeoTiffStore.java   |  11 +
 .../sis/storage/geotiff/GridGeometryBuilder.java   |   7 +-
 .../sis/storage/geotiff/ImageFileDirectory.java    | 395 +++++++++---
 .../org/apache/sis/storage/geotiff/Reader.java     |  16 +-
 .../sis/internal/storage/AbstractGridResource.java |  18 +-
 .../org/apache/sis/internal/storage/Resources.java |   5 +
 .../sis/internal/storage/Resources.properties      |   1 +
 .../sis/internal/storage/Resources_fr.properties   |   1 +
 .../sis/internal/storage/TiledGridCoverage.java    | 678 +++++++++++++++++++++
 .../sis/internal/storage/TiledGridResource.java    | 199 ++++++
 .../internal/storage/io/HyperRectangleReader.java  |  24 +-
 .../org/apache/sis/internal/storage/io/Region.java |  13 +-
 .../sis/test/storage/CoverageReadConsistency.java  | 301 +++++++++
 .../apache/sis/test/storage/SubsampledImage.java   | 247 ++++++++
 .../org/apache/sis/test/storage/package-info.java  |  34 ++
 27 files changed, 2442 insertions(+), 126 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/BandSelectImage.java b/core/sis-feature/src/main/java/org/apache/sis/image/BandSelectImage.java
index f209b49..097d425 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/image/BandSelectImage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/image/BandSelectImage.java
@@ -82,7 +82,7 @@ final class BandSelectImage extends SourceAlignedImage {
             return source;
         }
         ArgumentChecks.ensureNonEmpty("bands", bands, 0, numBands - 1, false);
-        final ColorModel cm = ColorModelFactory.createSubsetColorModel(source.getColorModel(), bands);
+        final ColorModel cm = ColorModelFactory.createSubset(source.getColorModel(), bands);
         /*
          * If the image is an instance of `BufferedImage`, create the subset immediately
          * (reminder: this operation will not copy pixel data). It allows us to return a
diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java b/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java
index 397dec2..0dd5392 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java
@@ -551,17 +551,17 @@ public abstract class PlanarImage implements RenderedImage {
             if (sm.getHeight() < tileHeight) return "tileHeight";
         }
         long size = getWidth();
-        long remainder = size - JDK9.multiplyFull(getNumXTiles(), tileWidth);
+        long remainder = JDK9.multiplyFull(getNumXTiles(), tileWidth) - size;
         if (remainder != 0) {
             return (remainder >= 0 && remainder < size) ? "width" : "numXTiles";
         }
         size = getHeight();
-        remainder = size - JDK9.multiplyFull(getNumYTiles(), tileHeight);
+        remainder = JDK9.multiplyFull(getNumYTiles(), tileHeight) - size;
         if (remainder != 0) {
             return (remainder >= 0 && remainder < size) ? "height" : "numYTiles";
         }
-        if (JDK9.multiplyFull(getMinTileX(),  tileWidth)  + getTileGridXOffset() != getMinX()) return "tileX";
-        if (JDK9.multiplyFull(getMinTileY(),  tileHeight) + getTileGridYOffset() != getMinY()) return "tileY";
+        if (JDK9.multiplyFull(getMinTileX(), tileWidth)  + getTileGridXOffset() != getMinX()) return "tileX";
+        if (JDK9.multiplyFull(getMinTileY(), tileHeight) + getTileGridYOffset() != getMinY()) return "tileY";
         return null;
     }
 
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/ColorModelFactory.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/ColorModelFactory.java
index eab9478..474724b 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/ColorModelFactory.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/ColorModelFactory.java
@@ -26,10 +26,12 @@ import java.awt.color.ColorSpace;
 import java.awt.image.ColorModel;
 import java.awt.image.IndexColorModel;
 import java.awt.image.PackedColorModel;
+import java.awt.image.DirectColorModel;
 import java.awt.image.ComponentColorModel;
 import java.awt.image.DataBuffer;
 import org.apache.sis.measure.NumberRange;
 import org.apache.sis.util.ArraysExt;
+import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.collection.WeakHashSet;
 import org.apache.sis.util.collection.WeakValueHashMap;
 import org.apache.sis.util.Debug;
@@ -310,10 +312,10 @@ public final class ColorModelFactory {
      *
      * @see Colorizer
      */
-    public static ColorModel createColorModel(final int dataType, final int numBands, final int visibleBand,
-                                              final Collection<Map.Entry<NumberRange<?>,Color[]>> colors)
+    public static ColorModel createPiecewise(final int dataType, final int numBands, final int visibleBand,
+                                             final Collection<Map.Entry<NumberRange<?>,Color[]>> colors)
     {
-        return createColorModel(dataType, numBands, visibleBand, ColorsForRange.list(colors));
+        return createPiecewise(dataType, numBands, visibleBand, ColorsForRange.list(colors));
     }
 
     /**
@@ -335,8 +337,8 @@ public final class ColorModelFactory {
      * @param  colors       the colors associated to ranges of sample values.
      * @return a color model suitable for {@link java.awt.image.RenderedImage} objects with values in the given ranges.
      */
-    static ColorModel createColorModel(final int dataType, final int numBands, final int visibleBand,
-                                       final ColorsForRange[] colors)
+    static ColorModel createPiecewise(final int dataType, final int numBands, final int visibleBand,
+                                      final ColorsForRange[] colors)
     {
         final ColorModelFactory key = new ColorModelFactory(dataType, numBands, visibleBand, colors);
         synchronized (PIECEWISES) {
@@ -357,7 +359,7 @@ public final class ColorModelFactory {
      * @param  transparent  the transparent pixel, or -1 for auto-detection.
      * @return An index color model for the specified array.
      */
-    public static IndexColorModel createIndexColorModel(final int numBands, final int visibleBand, final int[] ARGB, int transparent) {
+    public static IndexColorModel createIndexColorModel(final int numBands, final int visibleBand, final int[] ARGB, final int transparent) {
         /*
          * No need to scan the ARGB values in search of a transparent pixel;
          * the IndexColorModel constructor does that for us.
@@ -376,6 +378,26 @@ public final class ColorModelFactory {
     }
 
     /**
+     * Returns a color model interpolated for the given range of values. This is a convenience method for
+     * {@link #createColorModel(int, int, int, Collection)} when the collection contains only one element.
+     *
+     * @param  dataType     the color model type.
+     * @param  numBands     the number of bands for the color model (usually 1).
+     * @param  visibleBand  the band to be made visible (usually 0). All other bands (if any) will be ignored.
+     * @param  minimum      the minimum value, inclusive.
+     * @param  maximum      the maximum value, inclusive.
+     * @param  colors       the colors to use for the range of sample values.
+     * @return a color model suitable for {@link java.awt.image.RenderedImage} objects with values in the given ranges.
+     */
+    public static ColorModel createColorScale(final int dataType, final int numBands, final int visibleBand,
+                                              final double minimum, final double maximum, final Color... colors)
+    {
+        return createPiecewise(dataType, numBands, visibleBand, new ColorsForRange[] {
+            new ColorsForRange(null, new NumberRange<>(Double.class, minimum, true, maximum, true), colors)
+        });
+    }
+
+    /**
      * Creates a color model for opaque images storing pixels as real numbers.
      * The color model can have an arbitrary number of bands, but in current implementation only one band is used.
      *
@@ -453,21 +475,55 @@ public final class ColorModelFactory {
     }
 
     /**
+     * Creates a RGB color model. The {@code packed} argument should be
+     * {@code true}  for color model used with {@link java.awt.image.SinglePixelPackedSampleModel}, and
+     * {@code false} for color model used with {@link java.awt.image.BandedSampleModel}.
+     *
+     * @param  bitsPerSample  number of bits per sample, between 1 and 8 inclusive.
+     * @param  packed         whether sample values are packed in a single element.
+     * @param  hasAlpha       whether the color model should have an alpha channel.
+     * @return the color model.
+     */
+    public static ColorModel createRGB(final int bitsPerSample, final boolean packed, final boolean hasAlpha) {
+        if ((hasAlpha & packed) && bitsPerSample == Byte.SIZE) {
+            return ColorModel.getRGBdefault();
+        }
+        ArgumentChecks.ensureBetween("bitsPerSample", 1, Byte.SIZE, bitsPerSample);
+        final int mask = (1 << bitsPerSample) - 1;
+        final ColorModel cm;
+        if (packed) {
+            cm = new DirectColorModel((hasAlpha ? 4 : 3) * bitsPerSample,
+                    mask << (bitsPerSample * 2),        // Red
+                    mask <<  bitsPerSample,             // Green
+                    mask,                               // Blue
+                    hasAlpha ? mask << (bitsPerSample * 3) : 0);
+        } else {
+            final int[] numBits = new int[hasAlpha ? 4 : 3];
+            Arrays.fill(numBits, bitsPerSample);
+            cm = new ComponentColorModel(ColorSpace.getInstance(ColorSpace.CS_sRGB), numBits, hasAlpha, false,
+                            hasAlpha ? Transparency.TRANSLUCENT : Transparency.OPAQUE, DataBuffer.TYPE_BYTE);
+        }
+        return unique(cm);
+    }
+
+    /**
      * Creates a color model with only a subset of the bands of the given color model.
      *
      * @param  cm     the color model, or {@code null}.
      * @param  bands  the bands to select.
      * @return the subset color model, or {@code null} if it can not be created.
      */
-    public static ColorModel createSubsetColorModel(final ColorModel cm, final int[] bands) {
+    public static ColorModel createSubset(final ColorModel cm, final int[] bands) {
+        final ColorModel subset;
         if (cm instanceof MultiBandsIndexColorModel) {
-            return unique(((MultiBandsIndexColorModel) cm).createSubsetColorModel(bands));
-        }
-        if (cm instanceof ScaledColorModel) {
-            return unique(((ScaledColorModel) cm).createSubsetColorModel(bands));
+            subset = ((MultiBandsIndexColorModel) cm).createSubsetColorModel(bands);
+        } else if (cm instanceof ScaledColorModel) {
+            subset = ((ScaledColorModel) cm).createSubsetColorModel(bands);
+        } else {
+            // TODO: handle other color models.
+            return Colorizer.NULL_COLOR_MODEL;
         }
-        // TODO: handle other color models.
-        return Colorizer.NULL_COLOR_MODEL;
+        return unique(subset);
     }
 
     /**
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/Colorizer.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/Colorizer.java
index da2eacb..5259934 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/Colorizer.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/Colorizer.java
@@ -474,7 +474,7 @@ reuse:  if (source != null) {
         checkInitializationStatus(true);
         ArgumentChecks.ensureStrictlyPositive("numBands", numBands);
         ArgumentChecks.ensureBetween("visibleBand", 0, numBands - 1, visibleBand);
-        return ColorModelFactory.createColorModel(dataType, numBands, visibleBand, entries);
+        return ColorModelFactory.createPiecewise(dataType, numBands, visibleBand, entries);
     }
 
     /**
diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/TiledImage.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/TiledImage.java
index 8cce878..abecce8 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/TiledImage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/j2d/TiledImage.java
@@ -80,8 +80,7 @@ public class TiledImage extends PlanarImage {
      * @param height      number of pixels along Y axis in the whole rendered image.
      * @param minTileX    minimum tile index in the X direction.
      * @param minTileY    minimum tile index in the Y direction.
-     * @param tiles       the tiles. Must contains at least one element.
-     *                    This array is not cloned.
+     * @param tiles       the tiles. Must contains at least one element. This array is not cloned.
      */
     public TiledImage(final Map<String,Object> properties, final ColorModel colorModel,
                       final int width, final int height, final int minTileX, final int minTileY,
@@ -100,6 +99,40 @@ public class TiledImage extends PlanarImage {
     }
 
     /**
+     * Verifies whether image layout information and tile coordinates are consistent.
+     * This method verifies the size and minimum pixel coordinates of all tiles,
+     * in addition to the verifications documented in the super-class.
+     *
+     * @return {@code null} if image layout information are consistent,
+     *         or the name of inconsistent attribute if a problem is found.
+     */
+    @Override
+    public String verify() {
+        final int minX       = getMinX();
+        final int minY       = getMinY();
+        final int numXTiles  = getNumXTiles();
+        final int tileWidth  = getTileWidth();
+        final int tileHeight = getTileHeight();
+        for (int i=0; i < tiles.length; i++) {
+            final Raster tile = tiles[i];
+            final int tx = i % numXTiles;
+            final int ty = i / numXTiles;
+            if (tile.getWidth()  != tileWidth)              return property(tx, ty, "width");
+            if (tile.getHeight() != tileHeight)             return property(tx, ty, "height");
+            if (tile.getMinX()   != tileWidth  * tx + minX) return property(tx, ty, "x");
+            if (tile.getMinY()   != tileHeight * ty + minY) return property(tx, ty, "y");
+        }
+        return super.verify();
+    }
+
+    /**
+     * Label returned by {@link #verify()} for identifying an error in a specified tile.
+     */
+    private static String property(final int tx, final int ty, final String name) {
+        return "tiles[" + tx + ", " + ty + "]." + name;
+    }
+
+    /**
      * Returns the color model, or {@code null} if none.
      */
     @Override
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.java b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.java
index 0986f36..9adc70c 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.java
@@ -169,11 +169,21 @@ public final class Resources extends IndexedResourceBundle {
         public static final short UnexpectedTileCount_3 = 18;
 
         /**
+         * Tile has an unexpected length: {1} bytes while {0} were expected.
+         */
+        public static final short UnexpectedTileLength_2 = 28;
+
+        /**
          * TIFF file “{0}” uses an unknown coordinate reference system.
          */
         public static final short UnknownCRS_1 = 22;
 
         /**
+         * Compression method “{0}” is unsupported.
+         */
+        public static final short UnsupportedCompressionMethod_1 = 27;
+
+        /**
          * Coordinate system kind {0} is unsupported.
          */
         public static final short UnsupportedCoordinateSystemKind_1 = 19;
@@ -192,6 +202,11 @@ public final class Resources extends IndexedResourceBundle {
          * TIFF file “{0}” uses an unsupported map projection.
          */
         public static final short UnsupportedProjectionMethod_1 = 23;
+
+        /**
+         * Unsupported value “{1}” for TIFF tag “{0}”.
+         */
+        public static final short UnsupportedTagValue_2 = 29;
     }
 
     /**
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.properties b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.properties
index b88f5aa..9713ab5 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.properties
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources.properties
@@ -40,8 +40,11 @@ UndefinedDataFormat_1             = The \u201c{0}\u201d GeoTIFF file does not sp
 UnexpectedListOfValues_2          = A single value was expected for the \u201c{0}\u201d key but {1} values have been found.
 UnexpectedParameter_2             = The \u201c{1}\u201d parameter was not expected for the \u201c{0}\u201d projection method.
 UnexpectedTileCount_3             = Found {2} tiles or strips in the \u201c{0}\u201d file while {1} were expected.
+UnexpectedTileLength_2            = Tile has an unexpected length: {1} bytes while {0} were expected.
 UnknownCRS_1                      = TIFF file \u201c{0}\u201d uses an unknown coordinate reference system.
+UnsupportedCompressionMethod_1    = Compression method \u201c{0}\u201d is unsupported.
 UnsupportedCoordinateSystemKind_1 = Coordinate system kind {0} is unsupported.
 UnsupportedGeoKeyDirectory_1      = Version {0}\u00a0of GeoTIFF key directory is not supported.
 UnsupportedGeoKeyStorage_1        = Unsupported storage location for the \u201c{0}\u201d GeoTIFF value.
 UnsupportedProjectionMethod_1     = TIFF file \u201c{0}\u201d uses an unsupported map projection.
+UnsupportedTagValue_2             = Unsupported value \u201c{1}\u201d for TIFF tag \u201c{0}\u201d.
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources_fr.properties b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources_fr.properties
index cbab49f..1747215 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources_fr.properties
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/internal/geotiff/Resources_fr.properties
@@ -45,8 +45,11 @@ UndefinedDataFormat_1             = Le fichier GeoTIFF \u00ab\u202f{0}\u202f\u00
 UnexpectedListOfValues_2          = Une seule valeur \u00e9tait attendue pour la cl\u00e9 \u00ab\u202f{0}\u202f\u00bb, mais on en a trouv\u00e9es {1}.
 UnexpectedParameter_2             = Le param\u00e8tre \u00ab\u202f{1}\u202f\u00bb est inattendu pour la m\u00e9thode de projection \u00ab\u202f{0}\u202f\u00bb.
 UnexpectedTileCount_3             = {2} tuiles ont \u00e9t\u00e9 trouv\u00e9es dans le fichier \u00ab\u202f{0}\u202f\u00bb alors qu\u2019on en attendait {1}.
+UnexpectedTileLength_2            = Une tuile a une longueur inattendue\u00a0: {1} octets alors qu\u2019on en attendait {0}.
 UnknownCRS_1                      = Le fichier TIFF \u00ab\u202f{0}\u202f\u00bb utilise un syst\u00e8me de r\u00e9f\u00e9rence des coordonn\u00e9es inconnu.
+UnsupportedCompressionMethod_1    = La m\u00e9thode de compression  \u00ab\u202f{0}\u202f\u00bb n\u2019est pas support\u00e9e.
 UnsupportedCoordinateSystemKind_1 = Le type de syst\u00e8me de coordonn\u00e9es {0} n\u2019est pas support\u00e9.
 UnsupportedGeoKeyDirectory_1      = La version {0} du r\u00e9pertoire de cl\u00e9s GeoTIFF n\u2019est pas support\u00e9e.
 UnsupportedGeoKeyStorage_1        = La valeur GeoTIFF \u00ab\u202f{0}\u202f\u00bb utilise un mode de stockage non-support\u00e9.
 UnsupportedProjectionMethod_1     = Le fichier TIFF \u00ab\u202f{0}\u202f\u00bb utilise une projection cartographique non-support\u00e9e.
+UnsupportedTagValue_2             = La valeur \u201c{1}\u201d n\u2019est pas support\u00e9e pour le tag TIFF \u201c{0}\u201d.
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CRSBuilder.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CRSBuilder.java
index 82d1dd6..abd0d11 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CRSBuilder.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CRSBuilder.java
@@ -244,7 +244,7 @@ final class CRSBuilder extends ReferencingFactoryContainer {
      */
     private void warning(final short key, final Object... args) {
         final LogRecord r = reader.resources().getLogRecord(Level.WARNING, key, args);
-        reader.owner.warning(r);
+        reader.store.warning(r);
     }
 
     /**
@@ -463,7 +463,7 @@ final class CRSBuilder extends ReferencingFactoryContainer {
                 try {
                     expected = Integer.parseInt(id.getCode());
                 } catch (NumberFormatException e) {
-                    reader.owner.warning(null, e);                  // Should not happen.
+                    reader.store.warning(null, e);                  // Should not happen.
                     return;
                 }
                 if (code != expected) {
@@ -720,7 +720,7 @@ final class CRSBuilder extends ReferencingFactoryContainer {
         if (epsg != null) try {
             return getCSAuthorityFactory().createCartesianCS(epsg.toString());
         } catch (NoSuchAuthorityCodeException e) {
-            reader.owner.warning(null, e);
+            reader.store.warning(null, e);
         }
         return (CartesianCS) CoordinateSystems.replaceLinearUnit(cs, unit);
     }
@@ -738,7 +738,7 @@ final class CRSBuilder extends ReferencingFactoryContainer {
         if (epsg != null) try {
             return getCSAuthorityFactory().createEllipsoidalCS(epsg.toString());
         } catch (NoSuchAuthorityCodeException e) {
-            reader.owner.warning(null, e);
+            reader.store.warning(null, e);
         }
         return (EllipsoidalCS) CoordinateSystems.replaceAngularUnit(cs, unit);
     }
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java
new file mode 100644
index 0000000..5089a02
--- /dev/null
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataCube.java
@@ -0,0 +1,151 @@
+/*
+ * 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.storage.geotiff;
+
+import java.nio.file.Path;
+import java.awt.image.ColorModel;
+import java.awt.image.SampleModel;
+import org.apache.sis.storage.DataStoreException;
+import org.apache.sis.storage.DataStoreContentException;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.internal.geotiff.Resources;
+import org.apache.sis.internal.storage.TiledGridResource;
+import org.apache.sis.internal.storage.ResourceOnFileSystem;
+import org.apache.sis.util.resources.Errors;
+import org.apache.sis.math.Vector;
+
+
+/**
+ * One or many GeoTIFF images packaged as a single resource.
+ * This is typically a single two-dimensional image represented as a {@link ImageFileDirectory}.
+ * But it can also be a stack of images organized in a <var>n</var>-dimensional data cube,
+ * or a pyramid of images with their overviews used when low resolution images is requested.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+abstract class DataCube extends TiledGridResource implements ResourceOnFileSystem {
+    /**
+     * The GeoTIFF reader which contain this {@code DataCube}.
+     * Used for fetching information like the input channel and where to report warnings.
+     */
+    final Reader reader;
+
+    /**
+     * Creates a new data cube.
+     *
+     * @param  reader  information about the input stream to read, the metadata and the character encoding.
+     */
+    DataCube(final Reader reader) {
+        super(reader.store.listeners());
+        this.reader = reader;
+    }
+
+    /**
+     * Shortcut for a frequently requested information.
+     */
+    final String filename() {
+        return reader.input.filename;
+    }
+
+    /**
+     * Gets the paths to files used by this resource, or an empty array if unknown.
+     */
+    @Override
+    public final Path[] getComponentFiles() {
+        final Path location = reader.store.path;
+        return (location != null) ? new Path[] {location} : new Path[0];
+    }
+
+    /**
+     * Returns the Java2D sample model describing pixel type and layout.
+     * The raster size is the tile size as stored in the GeoTIFF file.
+     *
+     * <h4>Multi-dimensional data cube</h4>
+     * If this resource has more than 2 dimensions, then this model is for the two first ones (usually horizontal).
+     * The images for all levels in additional dimensions shall use the same sample model.
+     *
+     * @throws DataStoreContentException if the type is not recognized.
+     */
+    protected abstract SampleModel getSampleModel() throws DataStoreException;
+
+    /**
+     * Returns the Java2D color model for rendering images.
+     *
+     * @throws DataStoreContentException if the type is not recognized.
+     */
+    protected abstract ColorModel getColorModel() throws DataStoreException;
+
+    /**
+     * Returns the value to use for filling empty spaces in rasters,
+     * or {@code null} if none, not different than zero or not valid for the target data type.
+     * This value is used if a tile contains less pixels than expected.
+     * The zero value is excluded because tiles are already initialized to zero by default.
+     */
+    protected abstract Number fillValue();
+
+    /**
+     * Gets the stream position and the length in bytes of compressed tile arrays in the GeoTIFF file.
+     * Values in the returned vector are {@code long} primitive type.
+     *
+     * @return stream position (relative to file beginning) and length of compressed tile arrays, in bytes.
+     */
+    abstract Vector[] getTileArrayInfo();
+
+    /**
+     * Returns the compression method, or {@code null} if unspecified.
+     */
+    abstract Compression getCompression();
+
+    /**
+     * Creates a {@link GridCoverage} which will load pixel data in the given domain.
+     *
+     * @param  domain  desired grid extent and resolution, or {@code null} for reading the whole domain.
+     * @param  range   0-based index of sample dimensions to read, or an empty sequence for reading all ranges.
+     * @return the grid coverage for the specified domain and range.
+     * @throws DataStoreException if an error occurred while reading the grid coverage data.
+     */
+    @Override
+    public final GridCoverage read(final GridGeometry domain, final int... range) throws DataStoreException {
+        final long startTime = System.nanoTime();
+        final GridCoverage coverage;
+        try {
+            synchronized (reader.store) {
+                final Subset subset = new Subset(this, domain, range);
+                final Compression compression = getCompression();
+                if (compression == null) {
+                    throw new DataStoreContentException(reader.resources().getString(
+                            Resources.Keys.MissingValue_2, Tags.name(Tags.Compression)));
+                }
+                switch (compression) {
+                    case NONE: coverage = new DataSubset(this, subset); break;
+                    default: {
+                        throw new DataStoreContentException(reader.resources().getString(
+                                Resources.Keys.UnsupportedCompressionMethod_1, compression));
+                    }
+                }
+            }
+        } catch (RuntimeException e) {
+            throw new DataStoreException(reader.errors().getString(Errors.Keys.CanNotRead_1, filename()), e);
+        }
+        logReadOperation(reader.store.path, coverage.getGridGeometry(), startTime);
+        return coverage;
+    }
+}
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java
new file mode 100644
index 0000000..7306a1d
--- /dev/null
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/DataSubset.java
@@ -0,0 +1,297 @@
+/*
+ * 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.storage.geotiff;
+
+import java.util.Arrays;
+import java.util.Locale;
+import java.nio.Buffer;
+import java.io.IOException;
+import java.awt.Point;
+import java.awt.image.BandedSampleModel;
+import java.awt.image.DataBuffer;
+import java.awt.image.WritableRaster;
+import org.apache.sis.image.DataType;
+import org.apache.sis.storage.DataStoreException;
+import org.apache.sis.storage.DataStoreContentException;
+import org.apache.sis.internal.storage.io.Region;
+import org.apache.sis.internal.storage.io.HyperRectangleReader;
+import org.apache.sis.internal.storage.TiledGridCoverage;
+import org.apache.sis.internal.storage.TiledGridResource;
+import org.apache.sis.internal.coverage.j2d.ImageUtilities;
+import org.apache.sis.internal.coverage.j2d.RasterFactory;
+import org.apache.sis.internal.geotiff.Resources;
+import org.apache.sis.util.Localized;
+import org.apache.sis.math.Vector;
+
+import static java.lang.Math.addExact;
+import static java.lang.Math.subtractExact;
+import static java.lang.Math.multiplyExact;
+import static org.apache.sis.internal.jdk9.JDK9.multiplyFull;
+
+
+/**
+ * Raster data obtained from a GeoTIFF file in the domain requested by user. The number of dimensions is 2
+ * for standard TIFF files, but this class accepts higher number of dimensions if 3- or 4-dimensional data
+ * are stored in a GeoTIFF file using some convention. This base class transfers uncompressed data.
+ * Compressed data are handled by specialized subclasses.
+ *
+ * <h2>Cell Coordinates</h2>
+ * When there is no subsampling, {@code DataSubset} uses the same cell coordinates than {@link DataCube}.
+ * When there is a subsampling, cell coordinates in this subset are divided by the subsampling factors.
+ * Conversion is done by {@link #toFullResolution(long, int)}.
+ *
+ * <h2>Tile Matrix Coordinates</h2>
+ * In each {@code DataSubset}, indices of tiles starts at (0, 0, …). This class does not use
+ * the same tile indices than {@link DataCube} in order to avoid integer overflow.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+class DataSubset extends TiledGridCoverage implements Localized {
+    /**
+     * The resource which contain this {@code DataSubset}.
+     * Used for fetching information like the input channel and where to report warnings.
+     */
+    private final DataCube source;
+
+    /**
+     * For each tile, the byte offset of that tile as compressed and stored on disk.
+     * Tile X index varies fastest, followed by tile Y index, then tile Z index if any.
+     * The first tile included in this {@code DataSubset} is at index {@link #indexOfFirstTile}.
+     *
+     * @see ImageFileDirectory#tileOffsets
+     * @see #indexOfFirstTile
+     * @see #tileStrides
+     */
+    private final Vector tileOffsets;
+
+    /**
+     * For each tile, the number of (compressed) bytes in that tile.
+     * Elements are in the same order than {@link #tileOffsets}.
+     *
+     * @see ImageFileDirectory#tileByteCounts
+     * @see #indexOfFirstTile
+     * @see #tileStrides
+     */
+    private final Vector tileByteCounts;
+
+    /**
+     * Creates a new data subset. All parameters should have been validated
+     * by {@link ImageFileDirectory#validateMandatoryTags()} before this call.
+     * This constructor should be invoked inside a synchronized block.
+     *
+     * @param  source   the resource which contain this {@code DataSubset}.
+     * @param  subset   description of the {@code owner} subset to cover.
+     * @param  rasters  potentially shared cache of rasters read by this {@code DataSubset}.
+     * @throws ArithmeticException if the number of tiles overflows 32 bits integer arithmetic.
+     */
+    DataSubset(final DataCube source, final TiledGridResource.Subset subset) throws DataStoreException {
+        super(subset, source.getSampleModel(), source.getColorModel(), source.fillValue());
+        final Vector[] tileArrayInfo = source.getTileArrayInfo();
+        this.source         = source;
+        this.tileOffsets    = tileArrayInfo[0];
+        this.tileByteCounts = tileArrayInfo[1];
+    }
+
+    /**
+     * Returns the locale for warning or error messages, or {@code null} if unspecified.
+     */
+    @Override
+    public final Locale getLocale() {
+        return source.getLocale();
+    }
+
+    /**
+     * Returns an human-readable identification of this coverage for error messages.
+     */
+    @Override
+    protected final String getDisplayName() {
+        return source.filename();
+    }
+
+    /**
+     * Information about a tile to be read. A list of {@code Tile} is created and sorted by increasing offsets
+     * before the read operation begins, in order to read tiles in the order they are written in the TIFF file.
+     */
+    private static final class Tile extends Snapshot implements Comparable<Tile> {
+        /**
+         * Value of {@link DataSubset#tileOffsets} at index {@link #indexInTileVector}.
+         * This is the value that we want in increasing order.
+         *
+         * @see #compareTo(Tile)
+         */
+        final long byteOffset;
+
+        /**
+         * Stores information about a tile to be loaded.
+         *
+         * @param iterator    the iterator for which to create a snapshot of its current position.
+         * @param byteOffset  value of {@code tileOffsets.longValue(indexInTileVector)}.
+         */
+        Tile(final AOI domain, final long byteOffset) {
+            super(domain);
+            this.byteOffset = byteOffset;
+        }
+
+        /**
+         * Compares this tile with the specified tile for order in which to perform read operations.
+         */
+        @Override public int compareTo(final Tile other) {
+            return Long.compare(byteOffset, other.byteOffset);
+        }
+    }
+
+    /**
+     * Returns all tiles in the given area of interest. Tile indices are relative to this {@code DataSubset}:
+     * (0,0) is the tile in the upper-left corner of this {@code DataSubset} (not necessarily the upper-left
+     * corner of the image stored in the TIFF file).
+     *
+     * The {@link WritableRaster#getMinX()} and {@code getMinY()} coordinates of returned rasters
+     * will start at the given {@code offsetAOI} values.
+     *
+     * <p>This method is thread-safe.</p>
+     *
+     * @param  iterator  an iterator over the tiles that intersect the Area Of Interest specified by user.
+     * @return tiles decoded from the TIFF file.
+     * @throws IOException if an I/O error occurred.
+     * @throws DataStoreException if a logical error occurred.
+     * @throws RuntimeException if the Java2D image can not be created for another reason
+     *         (too many exception types to list them all).
+     */
+    @Override
+    protected final WritableRaster[] readTiles(final AOI iterator) throws IOException, DataStoreContentException {
+        /*
+         * Prepare an array for all tiles to be returned. Tiles that are already in memory will be stored
+         * in this array directly. Other tiles will be declared in the `missings` array and loaded later.
+         */
+        final WritableRaster[] result = new WritableRaster[iterator.tileCountInQuery];
+        final Tile[] missings = new Tile[iterator.tileCountInQuery];
+        int numMissings = 0;
+        synchronized (source.reader.store) {
+            do {
+                final WritableRaster tile = iterator.getCachedTile();
+                if (tile != null) {
+                    result[iterator.getIndexInResultArray()] = tile;
+                } else {
+                    // Tile not yet loaded. Add to a queue of tiles to load later.
+                    missings[numMissings++] = new Tile(iterator, tileOffsets.longValue(iterator.getIndexInTileVector()));
+                }
+            } while (iterator.next());
+            Arrays.sort(missings, 0, numMissings);
+            /*
+             * At this point we finished to list all tiles inside the Area Of Interest (AOI), both the ones that
+             * were already in memory and the new ones. The loop below processes only the new tiles, by reading
+             * them in the order they appear in the TIFF file.
+             *
+             * TODO: Use `tile.byteCount` for checking if two tiles are consecutive in the TIFF file.
+             * In such case we should send only one HTTP range request.
+             */
+            final long[] lower       = new long[BIDIMENSIONAL + 1];   // Coordinates of the first pixel to read relative to the tile.
+            final long[] upper       = new long[BIDIMENSIONAL + 1];   // Coordinates after the last pixel to read relative to the tile.
+            final int[]  subsampling = new int [BIDIMENSIONAL + 1];
+            final Point  origin      = new Point();
+            for (int i=0; i<numMissings; i++) {
+                final Tile tile = missings[i];
+                origin.x = tile.originX;
+                origin.y = tile.originY;
+                tile.getRegionInsideTile(lower, upper, subsampling, BIDIMENSIONAL);
+                result[tile.indexInResultArray] = tile.cache(readSlice(
+                        addExact(source.reader.origin, tile.byteOffset),
+                        tileByteCounts.longValue(tile.indexInTileVector),
+                        lower, upper, subsampling, origin));
+            }
+        }
+        return result;
+    }
+
+    /**
+     * Reads a two-dimensional slice of the data cube from the given input channel. This method is usually
+     * invoked for reading the tile in full, in which case the {@code lower} argument is (0,0,0) and the
+     * {@code upper} argument is the tile size. But those arguments may identify a smaller region if the
+     * {@link DataSubset} contains only one (potentially large) tile.
+     *
+     * <p>The length of {@code lower}, {@code upper} and {@code subsampling} arrays shall be 3,
+     * with the third dimension reserved for bands (this is not the <var>z</var> dimension).
+     * This method may modify those arrays in-place.</p>
+     *
+     * <p>The default implementation in this base class assumes uncompressed data.
+     * Subclasses must override for handling decompression.</p>
+     *
+     * @param  offset       position in the channel where tile data begins.
+     * @param  byteCount    number of bytes for the compressed tile data.
+     * @param  lower        (<var>x</var>, <var>y</var>, 0) coordinates of the first pixel to read relative to the tile.
+     * @param  upper        (<var>x</var>, <var>y</var>, 0) coordinates after the last pixel to read relative to the tile.
+     * @param  subsampling  (<var>sx</var>, <var>sy</var>, 1) subsampling factors.
+     * @param  location     pixel coordinates in the upper-left corner of the tile to return.
+     * @return image decoded from the GeoTIFF file.
+     * @throws IOException if an I/O error occurred.
+     * @throws DataStoreContentException if the sample model is not supported.
+     * @throws RuntimeException if the Java2D image can not be created for another reason
+     *         (too many exception types to list them all).
+     */
+    WritableRaster readSlice(final long offset, final long byteCount, final long[] lower, final long[] upper,
+                             final int[] subsampling, final Point location) throws IOException, DataStoreContentException
+    {
+        final int  type   = model.getDataType();
+        final long width  = subtractExact(upper[0], lower[0]);
+        final long height = subtractExact(upper[1], lower[1]);
+        final long length = multiplyExact(multiplyExact(width, height), DataBuffer.getDataTypeSize(type) / Byte.SIZE);
+        if (length > byteCount) {
+            // We may have less bytes to read if only a subregion is read.
+            throw new DataStoreContentException(source.reader.resources().getString(
+                    Resources.Keys.UnexpectedTileLength_2, length, byteCount));
+        }
+        final HyperRectangleReader hr = new HyperRectangleReader(ImageUtilities.toNumberEnum(type), source.reader.input, offset);
+        /*
+         * "Banks" (in `java.awt.image.DataBuffer` sense) are synonymous to "bands" for planar image only.
+         * Otherwise there is only one bank not matter the amount of bands. Each bank is read separately.
+         */
+        int numBanks = 1;
+        int numInterleaved = model.getNumBands();
+        if (model instanceof BandedSampleModel) {
+            numBanks = numInterleaved;
+            numInterleaved = 1;
+        }
+        final Buffer[] banks    = new Buffer[numBanks];
+        final long[]   size     = new long[] {multiplyFull(numInterleaved, getTileSize(0)), getTileSize(1), numBanks};
+        final int      capacity = multiplyExact(model.getWidth(), model.getHeight());
+        for (int b=0; b<numBanks; b++) {
+            lower[BIDIMENSIONAL] = b;
+            upper[BIDIMENSIONAL] = b + 1;
+            subsampling[BIDIMENSIONAL] = 1;
+            final Region region = new Region(size, lower, upper, subsampling);
+            final Buffer buffer = hr.readAsBuffer(region, capacity);
+            /*
+             * The buffer returned by `readAsBuffer(…)` may have less data than the buffer capacity
+             * if the current tile is smaller than the expected tile size (e.g. truncated last tile).
+             * Following code applies the fill value if it is different than the default value (zero).
+             */
+            if (fillValue != null) {
+                final int end = buffer.limit();
+                if (end != capacity) {
+                    Vector.create(buffer.limit(capacity), ImageUtilities.isUnsignedType(model))
+                          .fill(end, capacity, fillValue);
+                }
+            }
+            banks[b] = buffer.limit(capacity);
+        }
+        final DataBuffer buffer = RasterFactory.wrap(DataType.forDataBufferType(type), banks);
+        return WritableRaster.createWritableRaster(model, buffer, location);
+    }
+}
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTIFF.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTIFF.java
index 1ff3a40..3f9d7da 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTIFF.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTIFF.java
@@ -64,7 +64,7 @@ abstract class GeoTIFF implements Closeable {
     /**
      * The store which created this reader or writer.
      */
-    final GeoTiffStore owner;
+    final GeoTiffStore store;
 
     /**
      * The object to use for parsing and formatting dates. Created when first needed.
@@ -74,22 +74,22 @@ abstract class GeoTIFF implements Closeable {
     /**
      * For subclass constructors.
      */
-    GeoTIFF(final GeoTiffStore owner) {
-        this.owner = owner;
+    GeoTIFF(final GeoTiffStore store) {
+        this.store = store;
     }
 
     /**
      * Returns the resources to use for formatting error messages.
      */
     final Errors errors() {
-        return Errors.getResources(owner.getLocale());
+        return Errors.getResources(store.getLocale());
     }
 
     /**
      * Returns the GeoTIFF-specific resource for error messages and warnings.
      */
     final Resources resources() {
-        return Resources.forLocale(owner.getLocale());
+        return Resources.forLocale(store.getLocale());
     }
 
     /**
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTiffStore.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTiffStore.java
index d8175b6..aba07c6 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTiffStore.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GeoTiffStore.java
@@ -24,6 +24,7 @@ import java.net.URI;
 import java.io.IOException;
 import java.nio.charset.Charset;
 import java.nio.charset.StandardCharsets;
+import java.nio.file.Path;
 import java.nio.file.StandardOpenOption;
 import org.opengis.util.NameSpace;
 import org.opengis.util.NameFactory;
@@ -90,6 +91,15 @@ public class GeoTiffStore extends DataStore implements Aggregate {
     private final URI location;
 
     /**
+     * Same value than {@link #location} but as a path, or {@code null} if none.
+     * Stored separately because conversion from path to URI back to path is not
+     * looseness (relative paths become absolutes).
+     *
+     * @todo May become an array later if we want to handle TFW and PRJ file heres.
+     */
+    final Path path;
+
+    /**
      * The data store identifier created from the filename, or {@code null} if none.
      * Defined as a namespace for use as the scope of children resources (the images).
      */
@@ -128,6 +138,7 @@ public class GeoTiffStore extends DataStore implements Aggregate {
                     connector.getStorage(), connector.getOption(OptionKey.OPEN_OPTIONS));
         }
         location = connector.getStorageAs(URI.class);
+        path = connector.getStorageAs(Path.class);
         connector.closeAllExcept(input);
         try {
             reader = new Reader(this, input);
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
index 3a6973f..3ca072d 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
@@ -182,6 +182,7 @@ final class GridGeometryBuilder {
 
     /**
      * The grid geometry to be created by {@link #build(long, long)}.
+     * It has 2 or 3 dimensions, depending on whether the CRS declares a vertical axis or not.
      */
     public GridGeometry gridGeometry;
 
@@ -295,7 +296,7 @@ final class GridGeometryBuilder {
                 if (e instanceof NoSuchAuthorityCodeException) {
                     key = Resources.Keys.UnknownCRS_1;
                 }
-                reader.owner.warning(reader.resources().getString(key, reader.owner.getDisplayName()), e);
+                reader.store.warning(reader.resources().getString(key, reader.store.getDisplayName()), e);
             } catch (IllegalArgumentException | NoSuchElementException | ClassCastException e) {
                 if (!helper.alreadyReported) {
                     canNotCreate(e);
@@ -304,7 +305,7 @@ final class GridGeometryBuilder {
         }
         /*
          * If the CRS is non-null, then it is either two- or three-dimensional.
-         * The 'affine' matrix may be for a greater number of dimensions, so it
+         * The `affine` matrix may be for a greater number of dimensions, so it
          * may need to be reduced.
          */
         int n = (crs != null) ? crs.getCoordinateSystem().getDimension() : 2;
@@ -392,6 +393,6 @@ final class GridGeometryBuilder {
      * Logs a warning telling that we can not create a grid geometry for the given reason.
      */
     private void canNotCreate(final Exception e) {
-        reader.owner.warning(reader.resources().getString(Resources.Keys.CanNotComputeGridGeometry_1, reader.input.filename), e);
+        reader.store.warning(reader.resources().getString(Resources.Keys.CanNotComputeGridGeometry_1, reader.input.filename), e);
     }
 }
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java
index 5456d32..305a591 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/ImageFileDirectory.java
@@ -25,6 +25,13 @@ import java.util.Optional;
 import java.util.logging.Level;
 import java.util.logging.LogRecord;
 import java.nio.charset.Charset;
+import java.awt.Color;
+import java.awt.image.ColorModel;
+import java.awt.image.SampleModel;
+import java.awt.image.BandedSampleModel;
+import java.awt.image.MultiPixelPackedSampleModel;
+import java.awt.image.PixelInterleavedSampleModel;
+import java.awt.image.SinglePixelPackedSampleModel;
 import javax.measure.Unit;
 import javax.measure.quantity.Length;
 import org.opengis.metadata.citation.DateType;
@@ -33,36 +40,43 @@ import org.opengis.util.GenericName;
 import org.opengis.util.InternationalString;
 import org.apache.sis.internal.geotiff.Resources;
 import org.apache.sis.internal.storage.MetadataBuilder;
-import org.apache.sis.internal.storage.AbstractGridResource;
 import org.apache.sis.internal.storage.io.ChannelDataInput;
+import org.apache.sis.internal.coverage.j2d.ColorModelFactory;
 import org.apache.sis.internal.util.UnmodifiableArrayList;
 import org.apache.sis.storage.DataStoreException;
 import org.apache.sis.storage.DataStoreContentException;
-import org.apache.sis.coverage.grid.GridCoverage;
 import org.apache.sis.coverage.grid.GridGeometry;
 import org.apache.sis.coverage.grid.GridExtent;
 import org.apache.sis.coverage.SampleDimension;
 import org.apache.sis.util.resources.Vocabulary;
+import org.apache.sis.util.resources.Errors;
+import org.apache.sis.util.ArraysExt;
+import org.apache.sis.util.Numbers;
 import org.apache.sis.math.Vector;
 import org.apache.sis.measure.Units;
+import org.apache.sis.image.DataType;
 
 
 /**
  * An Image File Directory (FID) in a TIFF image.
  *
+ * <h2>Thread-safety</h2>
+ * Public methods should be synchronized because they can be invoked directly by users.
+ * Package-private methods are not synchronized; synchronization is caller's responsibility.
+ *
  * @author  Rémi Maréchal (Geomatys)
  * @author  Alexis Manin (Geomatys)
  * @author  Johann Sorel (Geomatys)
  * @author  Thi Phuong Hao Nguyen (VNSC)
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see <a href="http://www.awaresystems.be/imaging/tiff/tifftags.html">TIFF Tag Reference</a>
  *
  * @since 0.8
  * @module
  */
-final class ImageFileDirectory extends AbstractGridResource {
+final class ImageFileDirectory extends DataCube {
     /**
      * Possible value for the {@link #tileTagFamily} field. That field tells whether image tiling
      * was specified using the {@code Tile*} family of TIFF tags or the {@code Strip*} family.
@@ -77,12 +91,6 @@ final class ImageFileDirectory extends AbstractGridResource {
     private static final byte SIGNED = 1, UNSIGNED = 0, FLOAT = 3;
 
     /**
-     * The GeoTIFF reader which contain this {@code ImageFileDirectory}.
-     * Used for fetching information like the input channel and where to report warnings.
-     */
-    private final Reader reader;
-
-    /**
      * The identifier as a sequence number in the namespace of the {@link GeoTiffStore}.
      * The first image has the sequence number "1".
      *
@@ -140,7 +148,7 @@ final class ImageFileDirectory extends AbstractGridResource {
      * The offset is specified with respect to the beginning of the TIFF file.
      * Each tile has a location independent of the locations of other tiles
      *
-     * <p>Offsets are ordered left-to-right and top-to-bottom. if {@link #isPlanar} is {@code true}
+     * <p>Offsets are ordered left-to-right and top-to-bottom. If {@link #isPlanar} is {@code true}
      * (i.e. components are stored in separate “component planes”), then the offsets for the first
      * component plane are stored first, followed by all the offsets for the second component plane,
      * and so on.</p>
@@ -319,8 +327,17 @@ final class ImageFileDirectory extends AbstractGridResource {
     private Unit<Length> resolutionUnit = Units.INCH;
 
     /**
+     * The "no data" or background pixel value, or NaN if undefined.
+     *
+     * @see #fillValue()
+     */
+    private double noData = Double.NaN;
+
+    /**
      * The compression method, or {@code null} if unknown. If the compression method is unknown
      * or unsupported we can not read the image, but we still can read the metadata.
+     *
+     * @see #getCompression()
      */
     private Compression compression;
 
@@ -340,13 +357,6 @@ final class ImageFileDirectory extends AbstractGridResource {
     private GridGeometryBuilder referencing;
 
     /**
-     * The sample dimensions, or {@code null} if not yet created.
-     *
-     * @see #getSampleDimensions()
-     */
-    private List<SampleDimension> sampleDimensions;
-
-    /**
      * Returns {@link #referencing}, created when first needed. We delay its creation since
      * this object is not needed for ordinary TIFF files (i.e. without the GeoTIFF extension).
      */
@@ -358,15 +368,36 @@ final class ImageFileDirectory extends AbstractGridResource {
     }
 
     /**
+     * The sample dimensions, or {@code null} if not yet created.
+     *
+     * @see #getSampleDimensions()
+     */
+    private List<SampleDimension> sampleDimensions;
+
+    /**
+     * The image sample model, created when first needed. The raster size is the tile size.
+     * Sample models with different size and number of bands can be derived from this model.
+     *
+     * @see #getSampleModel()
+     */
+    private SampleModel sampleModel;
+
+    /**
+     * The image color model, created when first needed.
+     *
+     * @see #getColorModel()
+     */
+    private ColorModel colorModel;
+
+    /**
      * Creates a new image file directory.
      *
      * @param reader  information about the input stream to read, the metadata and the character encoding.
      * @param index   the image index as a sequence number starting with 0 for the first image.
      */
     ImageFileDirectory(final Reader reader, final int index) {
-        super(reader.owner.listeners());
-        this.reader = reader;
-        identifier = reader.nameFactory.createLocalName(reader.owner.identifier, String.valueOf(index + 1));
+        super(reader);
+        identifier = reader.nameFactory.createLocalName(reader.store.identifier, String.valueOf(index + 1));
     }
 
     /**
@@ -379,15 +410,8 @@ final class ImageFileDirectory extends AbstractGridResource {
     /**
      * Shortcut for a frequently requested information.
      */
-    private String filename() {
-        return input().filename;
-    }
-
-    /**
-     * Shortcut for a frequently requested information.
-     */
     private Charset encoding() {
-        return reader.owner.encoding;
+        return reader.store.encoding;
     }
 
     /**
@@ -573,9 +597,9 @@ final class ImageFileDirectory extends AbstractGridResource {
             case Tags.BitsPerSample: {
                 final Vector values = type.readVector(input(), count);
                 /*
-                 * The current implementation requires that all 'bitsPerSample' elements have the same value.
+                 * The current implementation requires that all `bitsPerSample` elements have the same value.
                  * This restriction may be revisited in future Apache SIS versions.
-                 * Note: 'count' is never zero when this method is invoked, so we do not need to check bounds.
+                 * Note: `count` is never zero when this method is invoked, so we do not need to check bounds.
                  */
                 bitsPerSample = values.shortValue(0);
                 final int length = values.size();
@@ -619,7 +643,7 @@ final class ImageFileDirectory extends AbstractGridResource {
              * 1 = BlackIsZero. For bilevel and grayscale images: 0 is imaged as black.
              * 2 = RGB. RGB value of (0,0,0) represents black, and (65535,65535,65535) represents white.
              * 3 = Palette color. The value of the component is used as an index into the RGB values of the ColorMap.
-             * 4 = Transparency Mask the defines an irregularly shaped region of another image in the same TIFF file.
+             * 4 = Transparency Mask. Defines an irregularly shaped region of another image in the same TIFF file.
              */
             case Tags.PhotometricInterpretation: {
                 final short value = type.readShort(input(), count);
@@ -706,7 +730,7 @@ final class ImageFileDirectory extends AbstractGridResource {
                 break;
             }
             /*
-             * Stores all of the 'double' valued GeoKeys, referenced by the GeoKeyDirectory.
+             * Stores all of the `double` valued GeoKeys, referenced by the GeoKeyDirectory.
              */
             case Tags.GeoDoubleParams: {
                 referencing().numericParameters = type.readVector(input(), count);
@@ -777,7 +801,7 @@ final class ImageFileDirectory extends AbstractGridResource {
             ////                                                                                        ////
             ////    Metadata for discovery purposes, conditions of use, etc.                            ////
             ////    Those metadata are not "critical" information for reading the image.                ////
-            ////    Should not write anything under 'metadata/contentInfo' node.                        ////
+            ////    Should not write anything under `metadata/contentInfo` node.                        ////
             ////                                                                                        ////
             ////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -939,6 +963,17 @@ final class ImageFileDirectory extends AbstractGridResource {
                 warning(Level.FINE, Resources.Keys.IgnoredTag_1, Tags.name(tag));
                 break;
             }
+
+            ////////////////////////////////////////////////////////////////////////////////////////////////
+            ////                                                                                        ////
+            ////    Extensions defined by GDAL.                                                         ////
+            ////                                                                                        ////
+            ////////////////////////////////////////////////////////////////////////////////////////////////
+
+            case Tags.GDAL_NODATA: {
+                noData = type.readDouble(input(), count);
+                break;
+            }
         }
         return null;
     }
@@ -973,7 +1008,7 @@ final class ImageFileDirectory extends AbstractGridResource {
                 i = s;
                 final Vector t = a; a = b; b = t;
             }
-            while (--i >= 0) {                      // At this point, 'b' shall be the longest vector.
+            while (--i >= 0) {                      // At this point, `b` shall be the longest vector.
                 final double va = a.doubleValue(i);
                 final double vb = b.doubleValue(i);
                 if (Double.isNaN(vb) || (max ? va > vb : va < vb)) {
@@ -1060,11 +1095,11 @@ final class ImageFileDirectory extends AbstractGridResource {
         }
         if (samplesPerPixel == 0) {
             samplesPerPixel = 1;
-            missingTag(Tags.SamplesPerPixel, 1, false);
+            missingTag(Tags.SamplesPerPixel, 1, false, false);
         }
         if (bitsPerSample == 0) {
             bitsPerSample = 1;
-            missingTag(Tags.BitsPerSample, 1, false);
+            missingTag(Tags.BitsPerSample, 1, false, false);
         }
         if (colorMap != null) {
             ensureSameLength(Tags.ColorMap, Tags.BitsPerSample, colorMap.size(),  3 * (1 << bitsPerSample));
@@ -1091,7 +1126,7 @@ final class ImageFileDirectory extends AbstractGridResource {
         /*
          * All of tile width, height and length information should be provided. But if only one of them is missing,
          * we can compute it provided that the file does not use any compression method. If there is a compression,
-         * then we set a bit for preventing the 'switch' block to perform a calculation but we let the code performs
+         * then we set a bit for preventing the `switch` block to perform a calculation but we let the code performs
          * the other checks in order to get an exception to be thrown with a good message.
          */
         int missing = !isPlanar && compression.equals(Compression.NONE) ? 0 : 0b1000;
@@ -1105,20 +1140,20 @@ final class ImageFileDirectory extends AbstractGridResource {
             }
             case 0b0001: {          // Compute missing tile width.
                 tileWidth = computeTileSize(tileHeight);
-                missingTag(Tags.TileWidth, tileWidth, true);
+                missingTag(Tags.TileWidth, tileWidth, true, true);
                 break;
             }
             case 0b0010: {          // Compute missing tile height.
                 tileHeight = computeTileSize(tileWidth);
-                missingTag(Tags.TileLength, tileHeight, true);
+                missingTag(Tags.TileLength, tileHeight, true, true);
                 break;
             }
-            case 0b0100: {          // Compute missing tile byte count.
+            case 0b0100: {          // Compute missing tile byte count in uncompressed case.
                 final long tileByteCount = pixelToByteCount(Math.multiplyExact(tileWidth, tileHeight));
                 final long[] tileByteCountArray = new long[tileOffsets.size()];
                 Arrays.fill(tileByteCountArray, tileByteCount);
                 tileByteCounts = Vector.create(tileByteCountArray, true);
-                missingTag(byteCountsTag, tileByteCount, true);
+                missingTag(byteCountsTag, tileByteCount, true, true);
                 break;
             }
             default: {
@@ -1224,6 +1259,12 @@ final class ImageFileDirectory extends AbstractGridResource {
         }
     }
 
+    /**
+     * Invoked the first time that {@link #getMetadata()} is invoked.
+     *
+     * @param  metadata  the builder where to set metadata properties.
+     * @throws DataStoreException if an error occurred while reading metadata from the data store.
+     */
     @Override
     protected void createMetadata(final MetadataBuilder metadata) throws DataStoreException {
         super.createMetadata(metadata);
@@ -1242,54 +1283,249 @@ final class ImageFileDirectory extends AbstractGridResource {
 
     /**
      * Returns an object containing the image size, the CRS and the conversion from pixel indices to CRS coordinates.
+     * The grid geometry has 2 or 3 dimensions, depending on whether the CRS declares a vertical axis or not.
+     *
+     * <h4>Thread-safety</h4>
+     * This method is thread-safe because it can be invoked directly by user.
+     *
+     * @see #getTileSize()
      */
     @Override
     public GridGeometry getGridGeometry() throws DataStoreContentException {
-        if (referencing != null) {
-            GridGeometry gridGeometry = referencing.gridGeometry;
-            if (gridGeometry == null) try {
-                gridGeometry = referencing.build(imageWidth, imageHeight);
-            } catch (FactoryException e) {
-                throw new DataStoreContentException(reader.resources().getString(Resources.Keys.CanNotComputeGridGeometry_1, filename()), e);
-            }
-            return gridGeometry;
-        } else {
-            return new GridGeometry(new GridExtent(imageWidth, imageHeight), null, null);
+        synchronized (reader.store) {
+            if (referencing != null) {
+                GridGeometry gridGeometry = referencing.gridGeometry;
+                if (gridGeometry == null) try {
+                    gridGeometry = referencing.build(imageWidth, imageHeight);
+                } catch (FactoryException e) {
+                    throw new DataStoreContentException(reader.resources().getString(Resources.Keys.CanNotComputeGridGeometry_1, filename()), e);
+                }
+                return gridGeometry;
+            } else {
+                return new GridGeometry(new GridExtent(imageWidth, imageHeight), null, null);
+            }
         }
     }
 
     /**
      * Returns the ranges of sample values together with the conversion from samples to real values.
+     *
+     * <h4>Thread-safety</h4>
+     * This method is thread-safe because it can be invoked directly by user.
      */
     @Override
     @SuppressWarnings("ReturnOfCollectionOrArrayField")
     public List<SampleDimension> getSampleDimensions() throws DataStoreContentException {
-        if (sampleDimensions == null) {
-            final SampleDimension[] dimensions = new SampleDimension[samplesPerPixel];
-            final SampleDimension.Builder builder = new SampleDimension.Builder();
-            final InternationalString name = Vocabulary.formatInternational(Vocabulary.Keys.Value);
-            for (int band = 0; band < samplesPerPixel;) {
-                builder.addQualitative(name, minValues.get(Math.min(band, minValues.size()-1)),
-                                             maxValues.get(Math.min(band, maxValues.size()-1)));
-                dimensions[band] = builder.setName(++band).build();
-                builder.clear();
-            }
-            sampleDimensions = UnmodifiableArrayList.wrap(dimensions);
+        synchronized (reader.store) {
+            if (sampleDimensions == null) {
+                final SampleDimension[] dimensions = new SampleDimension[samplesPerPixel];
+                final SampleDimension.Builder builder = new SampleDimension.Builder();
+                final InternationalString name = Vocabulary.formatInternational(Vocabulary.Keys.Value);
+                for (int band = 0; band < samplesPerPixel;) {
+                    builder.addQualitative(name, minValues.get(Math.min(band, minValues.size()-1)),
+                                                 maxValues.get(Math.min(band, maxValues.size()-1)));
+                    dimensions[band] = builder.setName(++band).build();
+                    builder.clear();
+                }
+                sampleDimensions = UnmodifiableArrayList.wrap(dimensions);
+            }
+            return sampleDimensions;        // Safe because unmodifiable.
         }
-        return sampleDimensions;        // Safe because unmodifiable.
     }
 
     /**
-     * Loads a subset of the grid coverage represented by this resource.
+     * Returns the Java2D sample model describing pixel type and layout.
+     * The raster size is the tile size and the number of bands is {@link #samplesPerPixel}.
+     * A sample model of different size or number of bands can be derived after construction
+     * by call to one of {@code SampleModel.create…} methods.
+     *
+     * @throws DataStoreContentException if the data type is not supported.
      *
-     * @param  domain  desired grid extent and resolution, or {@code null} for reading the whole domain.
-     * @param  range   0-based index of sample dimensions to read, or an empty sequence for reading all ranges.
-     * @return the grid coverage for the specified domain and range.
-     * @throws DataStoreException if an error occurred while reading the grid coverage data.
+     * @see SampleModel#createCompatibleSampleModel(int, int)
+     * @see SampleModel#createSubsetSampleModel(int[])
+     * @see #getColorModel()
+     * @see #getTileSize()
      */
     @Override
-    public GridCoverage read(final GridGeometry domain, final int... range) throws DataStoreException {
-        throw new DataStoreException("Not yet implemented.");   // TODO
+    protected SampleModel getSampleModel() throws DataStoreContentException {
+        if (sampleModel == null) {
+            final DataType type = getDataType();
+            try {
+                final int bt = type.ordinal();
+                if (bitsPerSample != type.size()) {
+                    /*
+                     * Supported types for both sample models are TYPE_BYTE, TYPE_USHORT, and TYPE_INT.
+                     * Note that if the TIFF data are signed bytes, then we have TYPE_SHORT which will
+                     * cause an exception to be thrown be sample model constructor.
+                     */
+                    if (samplesPerPixel == 1) {
+                        sampleModel = new MultiPixelPackedSampleModel(bt, tileWidth, tileHeight, bitsPerSample);
+                    } else if (!isPlanar) {
+                        final int[] bitMasks = new int[samplesPerPixel];
+                        bitMasks[0] = (1 << bitsPerSample) - 1;
+                        for (int i=1; i<bitMasks.length; i++) {
+                            bitMasks[i] = bitMasks[i-1] << bitsPerSample;
+                        }
+                        sampleModel = new SinglePixelPackedSampleModel(bt, tileWidth, tileHeight, bitMasks);
+                    } else {
+                        // TODO: we can support that with a little bit more work.
+                        throw new DataStoreContentException(Errors.format(Errors.Keys.UnsupportedType_1, type));
+                    }
+                } else if (isPlanar) {
+                    sampleModel = new BandedSampleModel(bt, tileWidth, tileHeight, samplesPerPixel);
+                } else {
+                    sampleModel = new PixelInterleavedSampleModel(bt, tileWidth, tileHeight, samplesPerPixel,
+                            Math.multiplyExact(samplesPerPixel, tileWidth), ArraysExt.range(0, samplesPerPixel));
+                }
+            } catch (IllegalArgumentException e) {
+                throw new DataStoreContentException(Errors.format(Errors.Keys.UnsupportedType_1, type), e);
+            }
+        }
+        return sampleModel;
+    }
+
+    /**
+     * Returns the size of tiles. This is also the size of the image sample model.
+     * The number of dimensions is always 2 for {@code ImageFileDirectory}.
+     *
+     * @see #getSampleModel()
+     */
+    @Override
+    protected int[] getTileSize() {
+        return new int[] {tileWidth, tileHeight};
+    }
+
+    /**
+     * Returns the type of raster data. The enumeration values are restricted to types compatible with Java2D,
+     * at the cost of using more bits than {@link #bitsPerSample} if there is no exact match.
+     *
+     * @throws DataStoreContentException if the type is not recognized.
+     */
+    private DataType getDataType() throws DataStoreContentException {
+        final String format;
+        switch (sampleFormat) {
+            case SIGNED: {
+                if (bitsPerSample <  Byte   .SIZE) return DataType.BYTE;
+                if (bitsPerSample <= Short  .SIZE) return DataType.SHORT;
+                if (bitsPerSample <= Integer.SIZE) return DataType.INT;
+                format = "int";
+                break;
+            }
+            case UNSIGNED: {
+                if (bitsPerSample <= Byte   .SIZE) return DataType.BYTE;
+                if (bitsPerSample <= Short  .SIZE) return DataType.USHORT;
+                if (bitsPerSample <= Integer.SIZE) return DataType.INT;
+                format = "unsigned";
+                break;
+            }
+            case FLOAT: {
+                if (bitsPerSample == Float  .SIZE) return DataType.FLOAT;
+                if (bitsPerSample == Double .SIZE) return DataType.DOUBLE;
+                format = "float";
+                break;
+            }
+            default: {
+                format = "?";
+                break;
+            }
+        }
+        throw new DataStoreContentException(Errors.format(
+                Errors.Keys.UnsupportedType_1, format + ' ' + bitsPerSample + " bits"));
+    }
+
+    /**
+     * Returns the Java2D color model.
+     *
+     * @throws DataStoreContentException if the data type is not supported.
+     *
+     * @see #getSampleModel()
+     */
+    @Override
+    protected ColorModel getColorModel() throws DataStoreContentException {
+        if (colorModel == null) {
+            final SampleModel sm  = getSampleModel();
+            final int dataType    = sm.getDataType();
+            final int visibleBand = 0;      // May be configurable in a future version.
+            short missing = 0;              // Non-zero if there is a warning about missing information.
+            switch (photometricInterpretation) {
+                default: {                  // For any unrecognized code, fallback on grayscale with 0 as black.
+                    unsupportedTagValue(Tags.PhotometricInterpretation, photometricInterpretation);
+                    break;
+                }
+                case -1: missing = Tags.PhotometricInterpretation; break;
+                case  0: break;             // WhiteIsZero: 0 is imaged as white.
+                case  1: break;             // BlackIsZero: 0 is imaged as black.
+                case  3: {                  // PaletteColor
+                    if (colorMap == null) {
+                        missing = Tags.ColorMap;
+                        break;
+                    }
+                    final int[] ARGB = new int[colorMap.size() / 3];
+                    for (int i=0, j=0; i < ARGB.length; i++) {
+                        ARGB[i] = ((colorMap.intValue(j++) & 0xFF00) << Byte.SIZE)
+                                | ((colorMap.intValue(j++) & 0xFF00))
+                                | ((colorMap.intValue(j++) & 0xFF00) >>> Byte.SIZE);
+                    }
+                    return colorModel = ColorModelFactory.createIndexColorModel(samplesPerPixel, visibleBand, ARGB,
+                                                          Double.isFinite(noData) ? (int) Math.round(noData) : -1);
+                }
+                case 2: {                   // RGB: (0,0,0) is black and (255,255,255) is white.
+                    colorModel = ColorModelFactory.createRGB(bitsPerSample, sm instanceof SinglePixelPackedSampleModel, false);
+                    break;
+                }
+            }
+            if (missing != 0) {
+                missingTag(missing, "GrayScale", false, true);
+            }
+            final Color[] colors = {Color.BLACK, Color.WHITE};
+            if (photometricInterpretation == 0) {
+                ArraysExt.swap(colors, 0, 1);
+            }
+            colorModel = ColorModelFactory.createColorScale(dataType, samplesPerPixel, visibleBand,
+                    Math.max(minValues.doubleValue(visibleBand), 0),
+                    Math.min(maxValues.doubleValue(visibleBand), (1 << bitsPerSample) - 1), colors);
+        }
+        return colorModel;
+    }
+
+    /**
+     * Returns the value to use for filling empty spaces in the raster, or {@code null} if none,
+     * not different than zero or not valid for the target data type.
+     */
+    @Override
+    protected Number fillValue() {
+        if (Double.isFinite(noData) && noData != 0) {
+            final long min, max;
+            switch (sampleFormat) {
+                case UNSIGNED: max = 1L << (bitsPerSample    ); min =    0; break;
+                case SIGNED:   max = 1L << (bitsPerSample - 1); min = ~max; break;
+                default: return noData;
+            }
+            final long value = Math.round(noData);
+            if (value >= min && value <= max && value != 0) {
+                return Numbers.narrowestNumber(value);
+            }
+        }
+        return null;
+    }
+
+    /**
+     * Gets the stream position or the length in bytes of compressed tile arrays in the GeoTIFF file.
+     * Values in the returned vector are {@code long} primitive type.
+     *
+     * @return stream position (relative to file beginning) and length of compressed tile arrays, in bytes.
+     */
+    @Override
+    Vector[] getTileArrayInfo() {
+        return new Vector[] {tileOffsets, tileByteCounts};
+    }
+
+    /**
+     * Returns the compression method, or {@code null} if unspecified.
+     */
+    @Override
+    Compression getCompression() {
+        return compression;
     }
 
     /**
@@ -1301,7 +1537,7 @@ final class ImageFileDirectory extends AbstractGridResource {
      */
     private void warning(final Level level, final short key, final Object... parameters) {
         final LogRecord r = reader.resources().getLogRecord(level, key, parameters);
-        reader.owner.warning(r);
+        reader.store.warning(r);
     }
 
     /**
@@ -1319,14 +1555,25 @@ final class ImageFileDirectory extends AbstractGridResource {
      * @param  missing   the numerical value of the missing tag.
      * @param  value     the default value or the computed value.
      * @param  computed  whether the default value has been computed.
+     * @param  warning   whether the problem should be reported as a warning.
      */
-    private void missingTag(final short missing, final long value, final boolean computed) {
-        warning(computed ? Level.WARNING : Level.FINE,
+    private void missingTag(final short missing, final Object value, final boolean computed, final boolean warning) {
+        warning(warning ? Level.WARNING : Level.FINE,
                 computed ? Resources.Keys.ComputedValueForAttribute_2 : Resources.Keys.DefaultValueForAttribute_2,
                 Tags.name(missing), value);
     }
 
     /**
+     * Reports a warning for an unsupported TIFF tag value.
+     *
+     * @param  tag    the numerical value of the tag.
+     * @param  value  the unsupported value.
+     */
+    private void unsupportedTagValue(final short tag, final Object value) {
+        warning(Level.WARNING, Resources.Keys.UnsupportedTagValue_2, Tags.name(tag), value);
+    }
+
+    /**
      * Builds an exception for a missing TIFF tag for which no default value can be computed.
      *
      * @param  missing  the numerical value of the missing tag.
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/Reader.java b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/Reader.java
index 8302217..e6f90be 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/Reader.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/Reader.java
@@ -50,7 +50,7 @@ import org.apache.sis.util.resources.Errors;
  * @author  Alexis Manin (Geomatys)
  * @author  Johann Sorel (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.8
  * @module
  */
@@ -63,7 +63,7 @@ final class Reader extends GeoTIFF {
     /**
      * Stream position of the first byte of the GeoTIFF file. This is usually zero.
      */
-    private final long origin;
+    final long origin;
 
     /**
      * A multiplication factor for the size of pointers, expressed as a power of 2.
@@ -137,8 +137,8 @@ final class Reader extends GeoTIFF {
      * @throws IOException if an error occurred while reading first bytes from the stream.
      * @throws DataStoreException if the file is not encoded in the TIFF or BigTIFF format.
      */
-    Reader(final GeoTiffStore owner, final ChannelDataInput input) throws IOException, DataStoreException {
-        super(owner);
+    Reader(final GeoTiffStore store, final ChannelDataInput input) throws IOException, DataStoreException {
+        super(store);
         this.input       = input;
         this.origin      = input.getStreamPosition();
         this.metadata    = new MetadataBuilder();
@@ -186,7 +186,7 @@ final class Reader extends GeoTIFF {
             }
         }
         // Do not invoke this.errors() yet because GeoTiffStore construction may not be finished. Owner.error() is okay.
-        throw new DataStoreContentException(owner.errors().getString(Errors.Keys.UnexpectedFileFormat_2, "TIFF", input.filename));
+        throw new DataStoreContentException(store.errors().getString(Errors.Keys.UnexpectedFileFormat_2, "TIFF", input.filename));
     }
 
     /**
@@ -215,7 +215,7 @@ final class Reader extends GeoTIFF {
         if (pointer >= 0) {
             return pointer;
         }
-        throw new DataStoreContentException(owner.getLocale(), "TIFF", input.filename, null);
+        throw new DataStoreContentException(store.getLocale(), "TIFF", input.filename, null);
     }
 
     /**
@@ -232,7 +232,7 @@ final class Reader extends GeoTIFF {
         if (entry >= 0) {
             return entry;
         }
-        throw new DataStoreContentException(owner.getLocale(), "TIFF", input.filename, null);
+        throw new DataStoreContentException(store.getLocale(), "TIFF", input.filename, null);
     }
 
     /**
@@ -388,7 +388,7 @@ final class Reader extends GeoTIFF {
             exception = null;
         }
         args[0] = Tags.name(tag);
-        owner.warning(errors().getString(key, args), exception);
+        store.warning(errors().getString(key, args), exception);
     }
 
     /**
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java
index 43f798e..31d0f68 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/AbstractGridResource.java
@@ -50,6 +50,21 @@ import org.apache.sis.internal.jdk9.JDK9;
 
 /**
  * Base class for implementations of {@link GridCoverageResource}.
+ * This class provides default implementations of {@code Resource} methods.
+ * Those default implementations get data from the following abstract methods:
+ *
+ * <ul>
+ *   <li>{@link #getGridGeometry()}</li>
+ *   <li>{@link #getSampleDimensions()}</li>
+ * </ul>
+ *
+ * This class also provides the following helper methods for implementation
+ * of the {@link #read(GridGeometry, int...) read(…)} method in subclasses:
+ *
+ * <ul>
+ *   <li>{@link #validateRangeArgument(int, int[])} for validation of the {@code range} argument.</li>
+ *   <li>{@link #logReadOperation(Object, GridGeometry, long)} for logging a notice about a read operation.</li>
+ * </ul>
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @version 1.1
@@ -105,7 +120,8 @@ public abstract class AbstractGridResource extends AbstractResource implements G
      * This method verifies that all indices are between 0 and {@code numSampleDimensions}
      * and that there is no duplicated index.
      *
-     * @param  numSampleDimensions  number of sample dimensions.
+     * @param  numSampleDimensions  number of sample dimensions in the resource.
+     *         Equal to <code>{@linkplain #getSampleDimensions()}.size()</code>.
      * @param  range  the {@code range} argument given by the user. May be null or empty.
      * @return the {@code range} argument encapsulated with a set of convenience tools.
      * @throws IllegalArgumentException if a range index is invalid.
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.java
index b053137..bb4515d 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.java
@@ -115,6 +115,11 @@ public final class Resources extends IndexedResourceBundle {
         public static final short CanNotRemoveResource_2 = 49;
 
         /**
+         * Can not render an image for the “{0}” coverage.
+         */
+        public static final short CanNotRenderImage_1 = 61;
+
+        /**
          * Can not save resources of type ‘{1}’ in a “{0}” store.
          */
         public static final short CanNotStoreResourceType_2 = 41;
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.properties b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.properties
index a19aac8..92372c0 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.properties
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources.properties
@@ -30,6 +30,7 @@ CanNotReadFile_2                  = Can not read \u201c{1}\u201d as a file in th
 CanNotReadFile_3                  = Can not read line {2} of \u201c{1}\u201d as part of a file in the {0} format.
 CanNotReadFile_4                  = Can not read after column {3} of line {2} of \u201c{1}\u201d as part of a file in the {0} format.
 CanNotRemoveResource_2            = Can not remove resource \u201c{1}\u201d from aggregate \u201c{0}\u201d.
+CanNotRenderImage_1               = Can not render an image for the \u201c{0}\u201d coverage.
 CanNotStoreResourceType_2         = Can not save resources of type \u2018{1}\u2019 in a \u201c{0}\u201d store.
 ClosedStorageConnector            = This storage connector is closed.
 ClosedReader_1                    = This {0} reader is closed.
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources_fr.properties b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources_fr.properties
index 9f9d58d..862f835 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources_fr.properties
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/Resources_fr.properties
@@ -35,6 +35,7 @@ CanNotReadFile_2                  = Ne peut pas lire \u00ab\u202f{1}\u202f\u00bb
 CanNotReadFile_3                  = Ne peut pas lire la ligne {2} de \u00ab\u202f{1}\u202f\u00bb comme une partie d\u2019un fichier au format {0}.
 CanNotReadFile_4                  = Ne peut pas lire apr\u00e8s la colonne {3} de la ligne {2} de \u00ab\u202f{1}\u202f\u00bb comme une partie d\u2019un fichier au format {0}.
 CanNotRemoveResource_2            = Ne peut pas supprimer la ressource \u00ab\u202f{1}\u202f\u00bb de l\u2019agr\u00e9gat \u00ab\u202f{0}\u202f\u00bb.
+CanNotRenderImage_1               = Ne peut pas produire une image pour la couverture de donn\u00e9es \u00ab\u202f{0}\u202f\u00bb.
 CanNotStoreResourceType_2         = Ne peut pas enregistrer des ressources de type \u2018{1}\u2019 dans un entrep\u00f4t de donn\u00e9es \u00ab\u202f{0}\u202f\u00bb.
 ClosedStorageConnector            = Ce connecteur est ferm\u00e9.
 ClosedReader_1                    = Ce lecteur {0} est ferm\u00e9.
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.java
new file mode 100644
index 0000000..af68b33
--- /dev/null
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridCoverage.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.storage;
+
+import java.util.Map;
+import java.util.Locale;
+import java.io.IOException;
+import java.awt.image.ColorModel;
+import java.awt.image.SampleModel;
+import java.awt.image.RenderedImage;
+import java.awt.image.WritableRaster;
+import org.opengis.coverage.CannotEvaluateException;
+import org.opengis.geometry.MismatchedDimensionException;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.coverage.grid.DisjointExtentException;
+import org.apache.sis.internal.coverage.j2d.TiledImage;
+import org.apache.sis.storage.DataStoreException;
+import org.apache.sis.util.resources.Errors;
+
+import static java.lang.Math.addExact;
+import static java.lang.Math.subtractExact;
+import static java.lang.Math.multiplyExact;
+import static java.lang.Math.incrementExact;
+import static java.lang.Math.decrementExact;
+import static java.lang.Math.toIntExact;
+import static java.lang.Math.floorDiv;
+import static org.apache.sis.internal.util.Numerics.ceilDiv;
+
+
+/**
+ * Base class of grid coverage read from a resource where data are stored in tiles.
+ * This grid coverage may represent only a subset of the coverage resource.
+ * Tiles are read from the storage only when first needed.
+ *
+ * <h2>Cell Coordinates</h2>
+ * When there is no subsampling, this coverage uses the same cell coordinates than the originating resource.
+ * When there is a subsampling, cell coordinates in this coverage are divided by the subsampling factors.
+ * Conversions are done by {@link #toFullResolution(long, int)}.
+ *
+ * <h2>Tile coordinate matrix</h2>
+ * In each {@code TiledGridCoverage}, indices of tiles starts at (0, 0, …).
+ * This class does not use the same tile indices than the coverage resource
+ * in order to avoid integer overflow.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public abstract class TiledGridCoverage extends GridCoverage {
+    /**
+     * Number of dimensions in a rendered image.
+     * Used for identifying codes where a two-dimensional slice is assumed.
+     */
+    protected static final int BIDIMENSIONAL = 2;
+
+    /**
+     * The area to read in unit of the full coverage (without subsampling).
+     * This is the intersection between user-specified domain and the source
+     * {@link TiledGridResource} domain, expanded to an integer number of tiles.
+     */
+    private final GridExtent readExtent;
+
+    /**
+     * Size of all tiles in the domain of this {@code TiledGridCoverage}, without clipping and subsampling.
+     * All coverages created from the same {@link TiledGridResource} have the same tile size values.
+     * The length of this array is the number of dimensions in the source {@link GridExtent}.
+     * This is often {@value #BIDIMENSIONAL} but can also be more.
+     */
+    private final int[] tileSize;
+
+    /**
+     * Values by which to multiply each tile coordinates for obtaining the index in the tile vector.
+     * The length of this array is the same as {@link #tileSize}. All coverages created from the same
+     * {@link TiledGridResource} have the same stride values.
+     */
+    private final int[] tileStrides;
+
+    /**
+     * Index of the first {@code TiledGridCoverage} tile in a row-major array of tiles.
+     * This is the value to add to the index computed with {@link #tileStrides} before to access vector elements.
+     */
+    private final int indexOfFirstTile;
+
+    /**
+     * The Tile Matrix Coordinates (TMC) of the first tile.
+     * This is the value to subtract from tile indices computed from pixel coordinates.
+     *
+     * @see #indexOfFirstTile
+     */
+    private final long[] tmcOfFirstTile;
+
+    /**
+     * Conversion from pixel coordinates in this (potentially subsampled) coverage
+     * to pixel coordinates in the resource coverage at full resolution.
+     * The conversion from (<var>x</var>, <var>y</var>) to (<var>x′</var>, <var>y′</var>) is as below,
+     * where <var>s</var> are subsampling factors and <var>t</var> are subsampling offsets:
+     *
+     * <ul>
+     *   <li><var>x′</var> = s₀⋅<var>x</var> + t₀</li>
+     *   <li><var>y′</var> = s₁⋅<var>y</var> + t₁</li>
+     * </ul>
+     *
+     * This transform maps {@linkplain org.opengis.referencing.datum.PixelInCell#CELL_CORNER pixel corners}.
+     *
+     * @see #toFullResolution(long, int)
+     */
+    private final int[] subsampling, subsamplingOffsets;
+
+    /**
+     * Cache of rasters read by this {@code TiledGridCoverage}. This cache may be shared with other coverages
+     * created for the same {@link TiledGridResource} resource. For each value, the raster {@code minX} and
+     * {@code minY} values can be anything, depending which {@code TiledGridCoverage} was first to load the tile.
+     *
+     * @see TiledGridResource#rasters
+     * @see AOI#getCachedTile()
+     */
+    private final Map<Integer, WritableRaster> rasters;
+
+    /**
+     * The sample model for all rasters. The size of this sample model is the values of
+     * the two first elements of {@link #tileSize} divided by subsampling after clipping.
+     */
+    protected final SampleModel model;
+
+    /**
+     * The Java2D color model for images rendered from this coverage.
+     */
+    protected final ColorModel colors;
+
+    /**
+     * The value to use for filling empty spaces in rasters, or {@code null} if zero.
+     */
+    protected final Number fillValue;
+
+    /**
+     * Creates a new tiled grid coverage. All parameters should have been validated before this call.
+     *
+     * @param  subset     description of the {@link TiledGridResource} subset to cover.
+     * @param  model      description of image layout before clipping and subsampling.
+     * @param  colors     Java2D color model for images rendered from this coverage.
+     * @param  fillValue  the value to use for filling empty spaces in rasters, of {@code null}.
+     * @throws ArithmeticException if the number of tiles overflows 32 bits integer arithmetic.
+     */
+    protected TiledGridCoverage(final TiledGridResource.Subset subset, SampleModel model,
+                                final ColorModel colors, final Number fillValue)
+    {
+        super(subset.domain,  subset.ranges);
+        final GridExtent extent = subset.domain.getExtent();
+        final int dimension = subset.sourceExtent.getDimension();
+        readExtent          = subset.readExtent;
+        subsampling         = subset.subsampling;
+        subsamplingOffsets  = subset.subsamplingOffsets;
+        rasters             = subset.cache;
+        tileSize            = subset.tileSize;
+        tmcOfFirstTile      = new long[dimension];
+        tileStrides         = new int [dimension];
+        final int[] subSize = new int [dimension];
+        int  tileStride       = 1;
+        long indexOfFirstTile = 0;
+        for (int i=0; i<dimension; i++) {
+            final int ts        = tileSize[i];
+            tmcOfFirstTile[i]   = floorDiv(readExtent.getLow(i), ts);
+            tileStrides   [i]   = tileStride;
+            subSize       [i]   = (int) Math.min(((ts-1) / subsampling[i]) + 1, extent.getSize(i));
+            indexOfFirstTile    = addExact(indexOfFirstTile, multiplyExact(tmcOfFirstTile[i], tileStride));
+            final int tileCount = toIntExact(ceilDiv(subset.sourceExtent.getSize(i), ts));
+            tileStride          = multiplyExact(tileCount, tileStride);
+        }
+        this.indexOfFirstTile = toIntExact(indexOfFirstTile);
+        /*
+         * At this point, `tileStride` is the total number of tiles in source.
+         * This value is not stored but its computation is still useful because
+         * we want `ArithmeticException` to be thrown if the value is too high.
+         *
+         * Note: the compatible sample model will require only the memory needed for the new size.
+         * But the subset sample model will use a `DataBuffer` that still consume all the memory of
+         * the sample model for all bands (contrarily to the sample model for a smaller raster size).
+         * It is okay for `PixelInterleavedSampleModel` because we will need to read all bands anyway.
+         * It will be more complicated for the `BandedSampleModel` case.
+         */
+        if (model.getWidth() != subSize[0] || model.getHeight() != subSize[1]) {
+            model = model.createCompatibleSampleModel(subSize[0], subSize[1]);
+        }
+        if (subset.selectedBands != null) {
+            model = model.createSubsetSampleModel(subset.selectedBands);
+        }
+        this.model     = model;
+        this.colors    = colors;
+        this.fillValue = fillValue;
+    }
+
+    /**
+     * Returns a name that identifies this coverage. It may be the filename.
+     * This is used for error messages.
+     *
+     * @return an human-readable identification of this coverage.
+     */
+    protected abstract String getDisplayName();
+
+    /**
+     * Returns the locale for error messages, or {@code null} for the default.
+     *
+     * @return the locale for warning or error messages, or {@code null} if unspecified.
+     */
+    protected Locale getLocale() {
+        return null;
+    }
+
+    /**
+     * Returns the size of all tiles in the domain of this {@code TiledGridCoverage}, without clipping and subsampling.
+     *
+     * @param  dimension  dimension for which to get tile size.
+     * @return tile size in the given dimension, without clipping and subsampling.
+     */
+    protected final int getTileSize(final int dimension) {
+        return tileSize[dimension];
+    }
+
+    /**
+     * Converts a cell coordinate from this {@code TiledGridCoverage} coordinate space to full resolution.
+     * This method removes the subsampling effect. Note that since this {@code TiledGridCoverage} uses the
+     * same coordinate space than {@link TiledGridResource}, the converted coordinates should be valid in
+     * the full resource as well.
+     *
+     * @param  coordinate  coordinate in this {@code TiledGridCoverage} domain.
+     * @param  dimension   dimension of the coordinate.
+     * @return coordinate in this {@code TiledGridCoverage} as if no subsampling was applied.
+     * @throws ArithmeticException if the coordinate can not be represented as a long integer.
+     */
+    private long toFullResolution(final long coordinate, final int dimension) {
+        return addExact(multiplyExact(coordinate, subsampling[dimension]), subsamplingOffsets[dimension]);
+    }
+
+    /**
+     * Converts a cell coordinate from {@link TiledGridResource} space to {@code TiledGridCoverage} coordinate.
+     * This is the converse of {@link #toFullResolution(long, int)}. Note that there is a possible accuracy lost.
+     *
+     * @param  coordinate  coordinate in the {@code TiledGridResource} domain.
+     * @param  dimension   dimension of the coordinate.
+     * @return coordinates in this subsampled {@code TiledGridCoverage} domain.
+     * @throws ArithmeticException if the coordinate can not be represented as a long integer.
+     */
+    private long toSubsampledPixel(final long coordinate, final int dimension) {
+        return floorDiv(subtractExact(coordinate, subsamplingOffsets[dimension]), subsampling[dimension]);
+    }
+
+    /**
+     * Converts a cell coordinate from this {@code TiledGridCoverage} coordinate space
+     * to the Tile Matrix Coordinate (TMC) of the tile which contains that cell.
+     * The TMC is relative to the full {@link TiledGridResource},
+     * i.e. without subtraction of {@link #tmcOfFirstTile}.
+     *
+     * @param  coordinate  coordinates in this {@code TiledGridCoverage} domain.
+     * @param  dimension   dimension of the coordinate.
+     * @return Tile Matrix Coordinate (TMC) of the tile which contains the specified cell.
+     * @throws ArithmeticException if the coordinate can not be represented as an integer.
+     */
+    private long toTileMatrixCoordinate(final long coordinate, final int dimension) {
+        return floorDiv(toFullResolution(coordinate, dimension), tileSize[dimension]);
+    }
+
+    /**
+     * Returns a two-dimensional slice of grid data as a rendered image.
+     *
+     * @param  sliceExtent  a subspace of this grid coverage extent, or {@code null} for the whole image.
+     * @return the grid slice as a rendered image. Image location is relative to {@code sliceExtent}.
+     */
+    @Override
+    public RenderedImage render(GridExtent sliceExtent) {
+        final GridExtent available = getGridGeometry().getExtent();
+        final int dimension = available.getDimension();
+        if (sliceExtent == null) {
+            sliceExtent = available;
+        } else {
+            final int sd = sliceExtent.getDimension();
+            if (sd != dimension) {
+                throw new MismatchedDimensionException(Errors.format(
+                        Errors.Keys.MismatchedDimension_3, "sliceExtent", dimension, sd));
+            }
+        }
+        final int[] sd = sliceExtent.getSubspaceDimensions(BIDIMENSIONAL);
+        if (sd[1] != 1) {
+            // TODO
+            throw new UnsupportedOperationException("Non-horizontal slices not yet implemented.");
+        }
+        final TiledImage image;
+        try {
+            final int[] tileLower = new int[dimension];         // Indices of first tile to read, inclusive.
+            final int[] tileUpper = new int[dimension];         // Indices of last tile to read, exclusive.
+            final int[] offsetAOI = new int[dimension];         // Pixel offset compared to Area Of Interest.
+            final int[] imageSize = new int[dimension];         // Subsampled image size.
+            for (int i=0; i<dimension; i++) {
+                final long min    = available  .getLow (i);     // Lowest valid coordinate in subsampled image.
+                final long max    = available  .getHigh(i);     // Highest valid coordinate, inclusive.
+                final long aoiMin = sliceExtent.getLow (i);     // Requested coordinate in subsampled image.
+                final long aoiMax = sliceExtent.getHigh(i);
+                final long tileUp = incrementExact(toTileMatrixCoordinate(Math.min(aoiMax, max), i));
+                final long tileLo =                toTileMatrixCoordinate(Math.max(aoiMin, min), i);
+                if (tileUp <= tileLo) {
+                    final String message = Errors.getResources(getLocale())
+                            .getString(Errors.Keys.IllegalRange_2, aoiMin, aoiMax);
+                    if (aoiMin > aoiMax) {
+                        throw new IllegalArgumentException(message);
+                    } else {
+                        throw new DisjointExtentException(message);
+                    }
+                }
+                final long lower =                Math.max(toSubsampledPixel(               multiplyExact(tileLo, tileSize[i]),  i), min);
+                final long upper = incrementExact(Math.min(toSubsampledPixel(decrementExact(multiplyExact(tileUp, tileSize[i])), i), max));
+                imageSize[i] = toIntExact(subtractExact(upper, lower));
+                offsetAOI[i] = toIntExact(subtractExact(lower, aoiMin));
+                tileLower[i] = toIntExact(subtractExact(tileLo, tmcOfFirstTile[i]));
+                tileUpper[i] = toIntExact(subtractExact(tileUp, tmcOfFirstTile[i]));
+            }
+            /*
+             * Get all tiles in the specified region. I/O operations, if needed, happen here.
+             */
+            final WritableRaster[] result = readTiles(new AOI(tileLower, tileUpper, offsetAOI, dimension));
+            image = new TiledImage(null, colors, imageSize[0], imageSize[1], tileLower[0], tileLower[1], result);
+        } catch (Exception e) {     // Too many exception types for listing them all.
+            throw new CannotEvaluateException(Resources.forLocale(getLocale()).getString(
+                    Resources.Keys.CanNotRenderImage_1, getDisplayName()), e);
+        }
+        assert validate(image);
+        return image;
+    }
+
+    /**
+     * The Area Of Interest specified by user in a call to {@link #render(GridExtent)}.
+     * This class is also an iterator over tiles in the region of interest.
+     */
+    protected final class AOI {
+        /**
+         * Total number of tiles in the AOI, from {@link #tileLower} inclusive to {@link #tileUpper} exclusive.
+         * This is the length of the array to be returned by {@link #readTiles(AOI)}.
+         */
+        public final int tileCountInQuery;
+
+        /**
+         * Indices (relative to enclosing {@code TiledGridCoverage}) of the upper-left tile to read.
+         * Tile indices starts at (0, 0, …), not at the indices of the corresponding tile in resource,
+         * for avoiding integer overflow.
+         */
+        private final int[] tileLower;
+
+        /**
+         * Indices (relative to enclosing {@code TiledGridCoverage}) after the bottom-right tile to read.
+         */
+        private final int[] tileUpper;
+
+        /**
+         * Pixel coordinates to assign to the upper-left corner of the region to read, ignoring subsampling.
+         */
+        private final int[] offsetAOI;
+
+        /**
+         * Tile Matrix Coordinates (TMC) of current iterator position relative to enclosing {@code TiledGridCoverage}.
+         * Initial position is a clone of {@link #tileLower}. This array is modified by calls to {@link #next()}.
+         */
+        private final int[] tmcInSubset;
+
+        /**
+         * Pixel coordinates of current iterator position relative to the Area Of Interest specified by user.
+         * Initial position is a clone of {@link #offsetAOI}. This array is modified by calls to {@link #next()}.
+         */
+        private final int[] tileOffsetAOI;
+
+        /**
+         * Current iterator position as an index in the array of tiles to be returned by {@link #readTiles(AOI)}.
+         * The initial position is zero. This field is incremented by calls to {@link #next()}.
+         */
+        private int indexInResultArray;
+
+        /**
+         * Current iterator position as an index in the vector of tiles in the {@link TiledGridResource}.
+         * Tiles are assumed stored in a row-major fashion. This field is incremented by calls to {@link #next()}.
+         */
+        private int indexInTileVector;
+
+        /**
+         * Creates a new Area Of Interest for the given tile indices.
+         *
+         * @param  tileLower  indices (relative to enclosing {@code TiledGridCoverage}) of the upper-left tile to read.
+         * @param  tileUpper  indices (relative to enclosing {@code TiledGridCoverage}) after the bottom-right tile to read.
+         * @param  offsetAOI  pixel coordinates to assign to the upper-left corner of the region to read, ignoring subsampling.
+         * @param  dimension  number of dimension of the {@code TiledGridCoverage} grid extent.
+         */
+        AOI(final int[] tileLower, final int[] tileUpper, final int[] offsetAOI, final int dimension) {
+            this.tileLower = tileLower;
+            this.tileUpper = tileUpper;
+            this.offsetAOI = offsetAOI;
+            /*
+             * Initialize variables to values for the first tile to read. The loop does arguments validation and
+             * converts the `tileLower` coordinates to index in the `tileOffsets` and `tileByteCounts` vectors.
+             */
+            indexInTileVector = indexOfFirstTile;
+            int tileCountInQuery = 1;
+            for (int i=0; i<dimension; i++) {
+                final int lower   = tileLower[i];
+                final int count   = subtractExact(tileUpper[i], lower);
+                indexInTileVector = addExact(indexInTileVector, multiplyExact(tileStrides[i], lower));
+                tileCountInQuery  = multiplyExact(tileCountInQuery, count);
+                /*
+                 * Following is the pixel coordinate after the last pixel in current dimension.
+                 * This is not stored; the intent is to get a potential `ArithmeticException`
+                 * now instead than in a call to `next()` during iteration. A negative value
+                 * would mean that the AOI does not intersect the region requested by user.
+                 */
+                final int max = addExact(offsetAOI[i], multiplyExact(tileSize[i], count));
+                assert max > Math.max(offsetAOI[i], 0) : max;
+            }
+            this.tileCountInQuery = tileCountInQuery;
+            this.tmcInSubset      = tileLower.clone();
+            this.tileOffsetAOI    = offsetAOI.clone();
+        }
+
+        /**
+         * Returns the enclosing coverage.
+         */
+        final TiledGridCoverage getCoverage() {
+            return TiledGridCoverage.this;
+        }
+
+        /**
+         * Returns the current iterator position as an index in the array of tiles to be returned
+         * by {@link #readTiles(AOI)}. The initial position is zero.
+         * The position is incremented by 1 in each call to {@link #next()}.
+         *
+         * @return current iterator position in result array.
+         */
+        public final int getIndexInResultArray() {
+            return indexInResultArray;
+        }
+
+        /**
+         * Returns the current iterator position as an index in the vector of tiles in the {@link TiledGridResource}.
+         * Tiles are assumed stored in a row-major fashion. This position is incremented by calls to {@link #next()}.
+         *
+         * @return current iterator position in the vector of all tiles in the resource.
+         */
+        public final int getIndexInTileVector() {
+            return indexInTileVector;
+        }
+
+        /**
+         * Returns the cached tile for current iterator position.
+         *
+         * @return cached tile at current iterator position, or {@code null} if none.
+         *
+         * @see Snapshot#cache(WritableRaster)
+         */
+        public WritableRaster getCachedTile() {
+            WritableRaster tile = rasters.get(indexInTileVector);
+            if (tile != null) {
+                // Found a tile. Make sure that the raster starts at the expected coordinates.
+                final int x = getTileOrigin(0);
+                final int y = getTileOrigin(1);
+                if (tile.getMinX() != x || tile.getMinY() != y) {
+                    tile = tile.createWritableTranslatedChild(x, y);
+                }
+            }
+            return tile;
+        }
+
+        /**
+         * Returns the origin to assign to the tile at current iterator position.
+         * Note that if there is more than one tile, then the subsampling should be
+         * a divisor of tile size, otherwise a drift in pixel coordinates will appear.
+         */
+        final int getTileOrigin(final int dimension) {
+            return floorDiv(tileOffsetAOI[dimension], subsampling[dimension]);
+        }
+
+        /**
+         * Moves the iterator position to next tile. This method should be invoked in a loop as below:
+         *
+         * {@preformat java
+         *     do {
+         *         // Process current tile.
+         *     } while (domain.next());
+         * }
+         *
+         * @return {@code true} on success, or {@code false} if the iteration is finished.
+         */
+        public boolean next() {
+            if (++indexInResultArray >= tileCountInQuery) {
+                return false;
+            }
+            /*
+             * Iterates over all tiles in the region specified to this method by maintaining 4 indices:
+             *
+             *   - `indexInResultArray` is the index in the `rasters` array to be returned.
+             *   - `indexInTileVector`  is the corresponding index in the `tileOffsets` vector.
+             *   - `tmcInSubset[]`      contains the Tile Matrix Coordinates (TMC) relative to this `DataSubset`.
+             *   - `tileOffsetAOI[]`    contains the pixel coordinates relative to the user-specified AOI.
+             *
+             * We do not check for integer overflow in this method because if an overflow is possible,
+             * then `ArithmeticException` should have occurred in `TiledGridCoverage` constructor.
+             */
+            for (int i=0; i<tmcInSubset.length; i++) {
+                indexInTileVector += tileStrides[i];
+                if (++tmcInSubset[i] < tileUpper[i]) {
+                    tileOffsetAOI[i] += tileSize[i];
+                    break;
+                }
+                // Rewind to index for tileLower[i].
+                indexInTileVector -= (tmcInSubset[i] - tileLower[i]) * tileStrides[i];
+                tmcInSubset  [i]   = tileLower[i];
+                tileOffsetAOI[i]   = offsetAOI[i];
+            }
+            return true;
+        }
+    }
+
+    /**
+     * Snapshot of a {@link AOI} iterator position. Those snapshots can be created during an iteration
+     * for processing a tile later. For example a {@link #readTiles(AOI)} method implementation may want
+     * to create a list of all tiles to load before to start the actual reading process in order to read
+     * the tiles in some optimal order, or for combining multiple read operations in a single operation.
+     */
+    protected static class Snapshot {
+        /**
+         * The source coverage.
+         */
+        private final TiledGridCoverage coverage;
+
+        /**
+         * Tile Matrix Coordinates (TMC) relative to the enclosing {@link DataSubset}.
+         * Tile (0,0) is the tile in the upper-left corner of this {@link DataSubset},
+         * not necessarily the tile in the upper-left corner of the image in the resource.
+         */
+        private final int[] tmcInSubset;
+
+        /**
+         * Index of this tile in the array of tiles returned by {@link #readTiles(AOI)}.
+         */
+        public final int indexInResultArray;
+
+        /**
+         * Index of this tile in the {@link DataCube} resource. This is the index of the tile in the
+         * {@link ImageFileDirectory#tileOffsets} and {@link ImageFileDirectory#tileByteCounts} vectors.
+         * This index is also used as key in the {@link TiledGridCoverage#rasters} map.
+         */
+        public final int indexInTileVector;
+
+        /**
+         * Pixel coordinates of the upper-left corner of the tile.
+         */
+        public final int originX, originY;
+
+        /**
+         * Stores information about a tile to be loaded.
+         *
+         * @param iterator  the iterator for which to create a snapshot of its current position.
+         */
+        public Snapshot(final AOI iterator) {
+            coverage = iterator.getCoverage();
+            this.tmcInSubset        = iterator.tmcInSubset.clone();
+            this.indexInResultArray = iterator.indexInResultArray;
+            this.indexInTileVector  = iterator.indexInTileVector;
+            this.originX            = iterator.getTileOrigin(0);
+            this.originY            = iterator.getTileOrigin(1);
+        }
+
+        /**
+         * Returns the coordinate of the pixel to read <em>inside</em> the tile, ignoring subsampling.
+         * The tile upper-left corner is assumed (0,0). Consequently the lower coordinates are usually
+         * (0,0) and the upper coordinates are usually the tile size, but those value may be different
+         * if the enclosing {@link TiledGridCoverage} contains only one (potentially big) tile.
+         * In that case, the reading process is more like untiled image reading.
+         *
+         * @param  lower        a pre-allocated array where to store relative coordinates of the first pixel.
+         * @param  upper        a pre-allocated array where to store relative coordinates after the last pixel.
+         * @param  subsampling  a pre-allocated array where to store subsampling.
+         * @param  dimension    number of elements to write in the {@code lower} and {@code upper} arrays.
+         */
+        public void getRegionInsideTile(final long[] lower, final long[] upper, final int[] subsampling, int dimension) {
+            System.arraycopy(coverage.subsampling, 0, subsampling, 0, dimension);
+            while (--dimension >= 0) {
+                final int  tileSize  = coverage.tileSize[dimension];
+                final long tileIndex = addExact(coverage.tmcOfFirstTile[dimension], tmcInSubset[dimension]);
+                final long tileBase  = multiplyExact(tileIndex, tileSize);
+                /*
+                 * The `offset` value is usually zero or negative because the tile to read should be inside the AOI,
+                 * e.g. at the right of the AOI left border. It may be positive if the `TiledGridCoverage` contains
+                 * only one (potentially big) tile, so the tile reading process become a reading of untiled data.
+                 */
+                long offset = subtractExact(coverage.readExtent.getLow(dimension), tileBase);
+                upper[dimension] = Math.min(addExact(offset, coverage.readExtent.getSize(dimension)), tileSize);
+                if (offset < 0) {
+                    /*
+                     * Example: for `tileSize` = 10 pixels and `subsampling` = 3,
+                     * the pixels to read are represented by black small squares below:
+                     *
+                     *  -10          0         10         20         30
+                     *    ┼──────────╫──────────┼──────────┼──────────╫
+                     *    │▪▫▫▪▫▫▪▫▫▪║▫▫▪▫▫▪▫▫▪▫│▫▪▫▫▪▫▫▪▫▫│▪▫▫▪▫▫▪▫▫▪║
+                     *    ┼──────────╫──────────┼──────────┼──────────╫
+                     *
+                     * If reading the second tile, then `tileBase` = 10 and `offset` = -10.
+                     * The first pixel to read in the second tile has a subsampling offset.
+                     * We usually try to avoid this situation because it causes a variable
+                     * number of white squares in tiles (4,3,3,4 in above example), except
+                     * when there is only 1 tile to read in which case offset is tolerated.
+                     */
+                    final int s = coverage.subsampling[dimension];
+                    offset %= s;
+                    if (offset != 0) {
+                        offset += s;
+                    }
+                }
+                lower[dimension] = offset;
+            }
+        }
+
+        /**
+         * Stores the given raster in the cache.
+         *
+         * @param  raster  the raster to cache.
+         * @return the cached raster. Should be the given {@code raster} instance,
+         *         but this method check for concurrent caching as a paranoiac check.
+         *
+         * @see AOI#getCachedTile()
+         */
+        public WritableRaster cache(final WritableRaster raster) {
+            final WritableRaster existing = coverage.rasters.putIfAbsent(indexInTileVector, raster);
+            return (existing != null) ? existing : raster;
+        }
+    }
+
+    /**
+     * Returns all tiles in the given area of interest. Tile indices are relative to this {@code TiledGridCoverage}:
+     * (0,0) is the tile in the upper-left corner of this {@code TiledGridCoverage} (not necessarily the upper-left
+     * corner of the image in the {@link TiledGridResource}).
+     *
+     * The {@link WritableRaster#getMinX()} and {@code getMinY()} coordinates of returned rasters
+     * shall start at the given {@code offsetAOI} values.
+     *
+     * <p>This method must be thread-safe.</p>
+     *
+     * @param  iterator  an iterator over the tiles that intersect the Area Of Interest specified by user.
+     * @return tiles decoded from the {@link TiledGridResource}.
+     * @throws IOException if an I/O error occurred.
+     * @throws DataStoreException if a logical error occurred.
+     * @throws RuntimeException if the Java2D image can not be created for another reason
+     *         (too many exception types to list them all).
+     */
+    protected abstract WritableRaster[] readTiles(AOI iterator) throws IOException, DataStoreException;
+
+    /**
+     * Verifies if the given image is consistent. This method tolerates smaller width or height compared to
+     * what we would expect from an integer amount of tiles, because this may get smaller height from images
+     * encoded as strips.
+     */
+    private static boolean validate(final TiledImage image) {
+        final String error = image.verify();
+        if (error == null || error.equals("width") || error.equals("height")) {
+            return true;
+        }
+        throw new AssertionError(error);
+    }
+}
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java
new file mode 100644
index 0000000..fd0c6ee
--- /dev/null
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/TiledGridResource.java
@@ -0,0 +1,199 @@
+/*
+ * 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.storage;
+
+import java.util.List;
+import java.util.Arrays;
+import java.awt.image.WritableRaster;
+import org.apache.sis.coverage.SampleDimension;
+import org.apache.sis.coverage.grid.GridDerivation;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.GridRoundingMode;
+import org.apache.sis.storage.DataStoreException;
+import org.apache.sis.storage.event.StoreListeners;
+import org.apache.sis.util.collection.WeakValueHashMap;
+
+
+/**
+ * Base class of grid coverage resource storing data in tiles.
+ * Word "tile" is used for simplicity but can be understood as "chunk"
+ * in a <var>n</var>-dimensional generalization.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public abstract class TiledGridResource extends AbstractGridResource {
+    /**
+     * All tiles loaded by any {@link TiledGridCoverage} created from this resource.
+     * Keys are tile indices in a row-major array of tiles.
+     * For each value, the {@link WritableRaster#getMinX()} and {@code minY} values
+     * can be anything, depending which {@link TiledGridResource} was first to load the tile.
+     *
+     * @see TiledGridCoverage#rasters
+     * @see TiledGridCoverage.AOI#getCachedTile()
+     */
+    private final WeakValueHashMap<Integer, WritableRaster> rasters;
+
+    /**
+     * Creates a new resource.
+     *
+     * @param  parent  listeners of the parent resource, or {@code null} if none.
+     */
+    protected TiledGridResource(final StoreListeners parent) {
+        super(parent);
+        rasters = new WeakValueHashMap<>(Integer.class);
+    }
+
+    /**
+     * Returns the size of tiles in this resource.
+     * The length of the returned array is the number of dimensions.
+     *
+     * @return the size of tiles in this resource.
+     */
+    protected abstract int[] getTileSize();
+
+    /**
+     * Parameters that describe the resource subset to be accepted by the {@link TiledGridCoverage} constructor.
+     * This is a temporary class used only for transferring information from {@link TiledGridResource}.
+     * This class does not perform I/O operations.
+     */
+    public static final class Subset {
+        /**
+         * The full size of the coverage in the enclosing {@link TiledGridResource}.
+         */
+        final GridExtent sourceExtent;
+
+        /**
+         * The area to read in unit of the full coverage (without subsampling).
+         * This is the intersection between user-specified domain and enclosing
+         * {@link TiledGridResource} domain, expanded to an integer number of tiles.
+         */
+        final GridExtent readExtent;
+
+        /**
+         * The sub-region extent, CRS and conversion from cell indices to CRS.
+         * This is the domain of the grid coverage to create.
+         */
+        final GridGeometry domain;
+
+        /**
+         * Sample dimensions for each image band.
+         * This is the range of the grid coverage to create.
+         */
+        final List<? extends SampleDimension> ranges;
+
+        /**
+         * Indices of {@link TiledGridResource} bands which have been retained
+         * for inclusion in the {@link TiledGridCoverage} to construct.
+         * This is {@code null} if all bands shall be included.
+         */
+        final int[] selectedBands;
+
+        /**
+         * Coordinate conversion from subsampled grid to the grid at full resolution.
+         * This array contains the factors by which to divide {@link TiledGridResource}
+         * cell coordinates in order to obtain {@link TiledGridCoverage} cell coordinates.
+         */
+        final int[] subsampling;
+
+        /**
+         * Remainder of the divisions of {@link TiledGridResource} cell coordinates by subsampling factors.
+         */
+        final int[] subsamplingOffsets;
+
+        /**
+         * Size of tiles in each dimension.
+         */
+        final int[] tileSize;
+
+        /**
+         * Cache to use for tiles loaded by the {@link TiledGridCoverage}.
+         * It is a reference to {@link TiledGridResource#rasters} if shareable.
+         */
+        final WeakValueHashMap<Integer, WritableRaster> cache;
+
+        /**
+         * Creates parameters for the given domain and range.
+         *
+         * @param  caller  the resource for which a subset is created.
+         * @param  domain  the domain argument specified by user in a call to {@code GridCoverageResource.read(…)}.
+         * @param  range   the range argument specified by user in a call to {@code GridCoverageResource.read(…)}.
+         * @throws ArithmeticException if pixel indices exceed 64 bits integer capacity.
+         * @throws DataStoreException if a call to {@link TiledGridResource} method failed.
+         */
+        public Subset(final TiledGridResource caller, GridGeometry domain, final int[] range) throws DataStoreException {
+            List<SampleDimension> bands        = caller.getSampleDimensions();
+            final RangeArgument   rangeIndices = caller.validateRangeArgument(bands.size(), range);
+            final GridGeometry    gridGeometry = caller.getGridGeometry();
+            sourceExtent = gridGeometry.getExtent();
+            tileSize = caller.getTileSize();
+            boolean sharedCache = true;
+            if (domain == null) {
+                domain             = gridGeometry;
+                readExtent         = sourceExtent;
+                subsamplingOffsets = new int[gridGeometry.getDimension()];
+                subsampling        = new int[subsamplingOffsets.length];
+                Arrays.fill(subsampling, 1);
+            } else {
+                /*
+                 * If an area of interest has been specified, we may need to expand it to an integer amount of tiles.
+                 * But we do not need to do that if the image is untiled; it is okay to read only a sub-region of the
+                 * single tile. We disable the "integer amount of tiles" restriction by setting the tile size to 1.
+                 * Note that it is possible to disable this restriction in a single dimension, typically the X one
+                 * when reading a TIFF image using strips instead of tiles.
+                 */
+                int tileWidth  = tileSize[0];
+                int tileHeight = tileSize[1];
+                if (tileWidth  >= sourceExtent.getSize(0)) {tileWidth  = 1; sharedCache = false;}
+                if (tileHeight >= sourceExtent.getSize(1)) {tileHeight = 1; sharedCache = false;}
+                final GridDerivation target = gridGeometry.derive().chunkSize(tileWidth, tileHeight)
+                                              .rounding(GridRoundingMode.ENCLOSING).subgrid(domain);
+                domain             = target.build();
+                readExtent         = target.getIntersection();
+                subsampling        = target.getSubsampling();
+                subsamplingOffsets = target.getSubsamplingOffsets();
+                if (sharedCache) {
+                    for (final int s : subsampling) {
+                        if (s != 1) {
+                            sharedCache = false;
+                            break;
+                        }
+                    }
+                }
+            }
+            int[] selectedBands = null;
+            if (!rangeIndices.isIdentity()) {
+                bands = Arrays.asList(rangeIndices.select(bands));
+                selectedBands = new int[rangeIndices.getNumBands()];
+                for (int i=0; i<selectedBands.length; i++) {
+                    selectedBands[rangeIndices.getTargetIndex(i)] = rangeIndices.getSourceIndex(i);
+                }
+            }
+            this.domain        = domain;
+            this.ranges        = bands;
+            this.selectedBands = selectedBands;
+            /*
+             * All `TiledGridCoverage` instances can share the same cache if they read all tiles fully.
+             * If they read only sub-regions or apply subsampling, then they will need their own cache.
+             */
+            cache = sharedCache ? caller.rasters : new WeakValueHashMap<>(Integer.class);
+        }
+    }
+}
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java
index 7cb3c93..585c931 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/HyperRectangleReader.java
@@ -123,26 +123,35 @@ public final class HyperRectangleReader {
      * @throws ArithmeticException if the region to read is too large or too far from origin.
      */
     public Object read(final Region region) throws IOException {
-        return read(region, false);
+        return read(region, 0, false);
     }
 
     /**
      * Reads data in the given region as a buffer. This method performs the same work
      * than {@link #read(Region)} except that the array is wrapped in a heap buffer.
+     * The {@code length} parameter is the minimal length of the array to allocate.
+     * The actual length of data read will be the {@linkplain Buffer#limit() limit}
+     * of the returned buffer.
      *
-     * @param  region  the sub-area to read and the subsampling to use.
+     * @param  region    the sub-area to read and the subsampling to use.
+     * @param  capacity  minimal length of the array to allocate, or 0 for automatic.
      * @return the data in a buffer backed by an array on the heap.
      * @throws IOException if an error occurred while transferring data from the channel.
      * @throws ArithmeticException if the region to read is too large or too far from origin.
      */
-    public Buffer readAsBuffer(final Region region) throws IOException {
-        return (Buffer) read(region, true);
+    public Buffer readAsBuffer(final Region region, final int capacity) throws IOException {
+        return (Buffer) read(region, capacity, true);
     }
 
     /**
      * Implementation of {@link #read(Region)} and {@link #readAsBuffer(Region)}.
+     *
+     * @param  region    the sub-area to read and the subsampling to use.
+     * @param  capacity  minimal length of the array to allocate, or 0 for automatic.
+     * @param  asBuffer  {@code true} for wrapping the array in a {@link Buffer}.
+     * @return the data as an array or wrapped in a buffer, depending on {@code asBuffer} value.
      */
-    private Object read(final Region region, final boolean asBuffer) throws IOException {
+    private Object read(final Region region, final int capacity, final boolean asBuffer) throws IOException {
         final int contiguousDataDimension = region.contiguousDataDimension();
         final int contiguousDataLength = region.targetLength(contiguousDataDimension);
         final long[] strides = new long[region.getDimension() - contiguousDataDimension];
@@ -154,8 +163,9 @@ public final class HyperRectangleReader {
             strides[i] = region.stride(i + contiguousDataDimension, contiguousDataLength, sampleSize);
             assert (strides[i] > 0) : i;
         }
+        final int limit = region.targetLength(region.getDimension());
         try {
-            reader.createDataArray(region.targetLength(region.getDimension()));
+            reader.createDataArray(Math.max(capacity, limit));
             final Buffer view = reader.view();
 loop:       do {
                 reader.seek(streamPosition);
@@ -178,7 +188,7 @@ loop:       do {
                 }
                 break;
             } while (true);
-            return asBuffer ? reader.dataArrayAsBuffer() : reader.dataArray();
+            return asBuffer ? reader.dataArrayAsBuffer().limit(limit) : reader.dataArray();
         } finally {
             reader.setDest(null);
         }
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/Region.java b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/Region.java
index 0bb9c80..fbbd802 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/Region.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/internal/storage/io/Region.java
@@ -82,6 +82,12 @@ public final class Region {
     private long[] skipBytes;
 
     /**
+     * Total length of the region.
+     * This is the product of all values in the {@code size} argument given to the constructor.
+     */
+    public final long length;
+
+    /**
      * Creates a new region. It is caller's responsibility to ensure that:
      * <ul>
      *   <li>all arrays have the same length</li>
@@ -122,6 +128,7 @@ public final class Region {
             skips[++i]    = skip;
         }
         startAt = position;
+        length  = stride;
     }
 
     /**
@@ -220,11 +227,11 @@ public final class Region {
      * This method takes in account only the given number of dimensions.
      */
     final int targetLength(final int dimension) {
-        long length = 1;
+        long size = 1;
         for (int i=0; i<dimension; i++) {
-            length *= targetSize[i];
+            size *= targetSize[i];
         }
-        return toIntExact(length);
+        return toIntExact(size);
     }
 
     /**
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java
new file mode 100644
index 0000000..ea1de7d
--- /dev/null
+++ b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/CoverageReadConsistency.java
@@ -0,0 +1,301 @@
+/*
+ * 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.test.storage;
+
+import java.util.Arrays;
+import java.util.Random;
+import java.awt.Point;
+import java.awt.Rectangle;
+import java.awt.image.RenderedImage;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridDerivation;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.storage.DataStoreException;
+import org.apache.sis.storage.GridCoverageResource;
+import org.apache.sis.image.PixelIterator;
+import org.apache.sis.test.DependsOnMethod;
+import org.apache.sis.test.TestUtilities;
+import org.apache.sis.test.TestCase;
+import org.junit.Test;
+
+import static org.junit.Assert.*;
+
+
+/**
+ * Base class for testing the consistency of grid coverage read operations.
+ * The test reads the grid coverage in full, then reads random sub-regions.
+ * The sub-regions pixels are compared with the original image.
+ *
+ * <h2>Assumptions</h2>
+ * Assuming that the code reading the full extent is correct, this class can detect some bugs
+ * in the code reading sub-regions or applying sub-sampling. This assumption is reasonable if
+ * we consider that the code reading the full extent is usually simpler than the code reading
+ * a subset of data.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+public strictfp class CoverageReadConsistency extends TestCase {
+    /**
+     * A constant for identifying the codes working on two dimensional slices.
+     */
+    private static final int BIDIMENSIONAL = 2;
+
+    /**
+     * The resource to test.
+     */
+    private final GridCoverageResource resource;
+
+    /**
+     * The coverage at full extent, full resolution and with all bands.
+     * This coverage will be used as a reference for verifying values read in sub-domains.
+     */
+    private final GridCoverage full;
+
+    /**
+     * The random number generator to use.
+     */
+    private final Random random;
+
+    /**
+     * Creates a new tester. This constructor reads immediately the coverage at full extent and full resolution.
+     * That full coverage will be used as a reference for verifying values read in sub-domains.
+     *
+     * @param  resource  the resource to test.
+     * @throws DataStoreException if the full coverage can not be read.
+     */
+    public CoverageReadConsistency(final GridCoverageResource resource) throws DataStoreException {
+        this.resource = resource;
+        full = resource.read(null, null);
+        random = TestUtilities.createRandomNumberGenerator();
+    }
+
+    /**
+     * Tests reading in random sub-regions starting at coordinates (0,0).
+     * Data are read at full resolution (no-subsampling) and all bands are read. This is the simplest test.
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    public void testSubRegionAtOrigin() throws DataStoreException {
+        readAndCompareRandomRegions(10, false, false);
+    }
+
+    /**
+     * Tests reading in random sub-regions starting at random offsets.
+     * Data are read at full resolution (no-subsampling) and all bands are read.
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    @DependsOnMethod("testSubRegionAtOrigin")
+    public void testSubRegionsAnywhere() throws DataStoreException {
+        readAndCompareRandomRegions(10, true, false);
+    }
+
+    /**
+     * Tests reading in random sub-regions starting at coordinates (0,0) with subsampling applied.
+     * All bands are read.
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    @DependsOnMethod("testSubRegionAtOrigin")
+    public void testSubsamplingAtOrigin() throws DataStoreException {
+        readAndCompareRandomRegions(10, false, true);
+    }
+
+    /**
+     * Tests reading in random sub-regions starting at random offsets with subsampling applied.
+     * All bands are read.
+     *
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    @Test
+    @DependsOnMethod({"testSubsamplingAtOrigin", "testSubRegionsAnywhere"})
+    public void testSubsamplingAnywhere() throws DataStoreException {
+        readAndCompareRandomRegions(100, true, true);
+    }
+
+    /**
+     * Implementation of methods testing reading in random sub-regions with random sub-samplings.
+     *
+     * @param  numIterations      number of random sub-regions to test.
+     * @param  allowOffsets       whether to allow sub-regions to start elsewhere than (0,0).
+     * @param  allowSubsamplings  whether to use random subsamplings.
+     * @throws DataStoreException if an error occurred while using the resource.
+     */
+    private void readAndCompareRandomRegions(int numIterations, final boolean allowOffsets, final boolean allowSubsamplings)
+            throws DataStoreException
+    {
+        final GridGeometry gg = resource.getGridGeometry();
+        final GridExtent fullExtent = gg.getExtent();
+        final int    dimension   = fullExtent.getDimension();
+        final long[] low         = new long[dimension];
+        final long[] high        = new long[dimension];
+        final int [] subsampling = new int [dimension];
+        final int [] subOffsets  = new int [dimension];
+        while (--numIterations >= 0) {
+            /*
+             * Create a random domain to be used as a query on the `GridCoverageResource`.
+             */
+            for (int d=0; d<dimension; d++) {
+                final int span = StrictMath.toIntExact(fullExtent.getSize(d));
+                final int rs = random.nextInt(span);                    // Span of the sub-region - 1.
+                if (allowOffsets) {
+                    low[d] = random.nextInt(span - rs);                 // Note: (span - rs) > 0.
+                }
+                high[d] = low[d] + rs;
+                subsampling[d] = 1;
+                if (allowSubsamplings) {
+                    subsampling[d] += random.nextInt(StrictMath.max(rs / 16, 1));
+                }
+            }
+            final GridGeometry domain = gg.derive()
+                    .subgrid(new GridExtent(null, low, high, true), subsampling).build();
+            /*
+             * Read a coverage containing the requested sub-domain. Note that the reader is free to read
+             * more data than requested. The extent actually read is `actualReadExtent`. It shall contain
+             * fully the requested `domain`.
+             */
+            final GridCoverage subset = resource.read(domain, null);
+            final GridExtent actualReadExtent = subset.getGridGeometry().getExtent();
+            assertEquals("Unexpected number of dimensions.", dimension, actualReadExtent.getDimension());
+            for (int d=0; d<dimension; d++) {
+                if (subsampling[d] == 1) {
+                    assertTrue("Actual extent is too small.", actualReadExtent.getLow (d) <= low [d]);
+                    assertTrue("Actual extent is too small.", actualReadExtent.getHigh(d) >= high[d]);
+                }
+            }
+            /*
+             * If subsampling was enabled, the factors selected by the reader may be different than
+             * the subsampling factors that we specified. The following block updates those values.
+             */
+            if (allowSubsamplings) {
+                final GridDerivation change = full.getGridGeometry().derive().subgrid(subset.getGridGeometry());
+                System.arraycopy(change.getSubsampling(),        0, subsampling, 0, dimension);
+                System.arraycopy(change.getSubsamplingOffsets(), 0, subOffsets,  0, dimension);
+            }
+            /*
+             * Iterate over all dimensions greater than 2. In the common case where we are reading a
+             * two-dimensional image, the following loop will be executed only once. If reading a 3D
+             * or 4D image, the loop is executed for all possible two-dimensional slices in the cube.
+             */
+            final long[] sliceMin = actualReadExtent.getLow() .getCoordinateValues();
+            final long[] sliceMax = actualReadExtent.getHigh().getCoordinateValues();
+nextSlice:  for (;;) {
+                System.arraycopy(sliceMin, BIDIMENSIONAL, sliceMax, BIDIMENSIONAL, dimension - BIDIMENSIONAL);
+                final PixelIterator itr = iterator(full,   sliceMin, sliceMax, subsampling, subOffsets, allowSubsamplings);
+                final PixelIterator itc = iterator(subset, sliceMin, sliceMax, subsampling, subOffsets, false);
+                if (itr != null) {
+                    assertEquals(itr.getDomain().getSize(), itc.getDomain().getSize());
+                    double[] expected = null, actual = null;
+                    while (itr.next()) {
+                        assertTrue(itc.next());
+                        expected = itr.getPixel(expected);
+                        actual   = itc.getPixel(actual);
+                        if (!Arrays.equals(expected, actual)) {
+                            final Point pr = itr.getPosition();
+                            final Point pc = itc.getPosition();
+                            assertArrayEquals("Mismatch at position (" + pr.x + ", " + pr.y + ") in full image " +
+                                              "and (" + pc.x + ", " + pc.y + ") in tested sub-image",
+                                              expected, actual, STRICT);
+                        }
+                    }
+                    assertFalse(itc.next());
+                } else {
+                    // Unable to create a reference image. Just check that no exception is thrown.
+                    double[] actual = null;
+                    while (itc.next()) {
+                        actual = itc.getPixel(actual);
+                    }
+                }
+                /*
+                 * Move to the next two-dimensional slice and read again.
+                 * We stop the loop after we have read all 2D slices.
+                 */
+                for (int d=dimension; --d >= BIDIMENSIONAL;) {
+                    if (sliceMin[d]++ <= actualReadExtent.getHigh(d)) continue nextSlice;
+                    sliceMin[d] = actualReadExtent.getLow(d);
+                }
+                break;
+            }
+        }
+    }
+
+    /**
+     * Creates a pixel iterator for a sub-region in a slice of the specified coverage.
+     * All coordinates given to this method are in the coordinate space of subsampled coverage subset.
+     * This method returns {@code null} if the arguments are valid but the image can not be created
+     * because of a restriction in {@code PixelInterleavedSampleModel} constructor.
+     *
+     * @param  coverage     the coverage from which to get the iterator.
+     * @param  sliceMin     lower bounds of the <var>n</var>-dimensional region of the coverage for which to get an iterator.
+     * @param  sliceMax     upper bounds of the <var>n</var>-dimensional region of the coverage for which to get an iterator.
+     * @param  subsampling  subsampling factors to apply on the image.
+     * @param  subOffsets   offsets to add after multiplication by subsampling factors.
+     * @return pixel iterator over requested area, or {@code null} if unavailable.
+     */
+    private static PixelIterator iterator(final GridCoverage coverage, long[] sliceMin, long[] sliceMax,
+            final int[] subsampling, final int[] subOffsets, final boolean allowSubsamplings)
+    {
+        /*
+         * Same extent than `areaOfInterest` but in two dimensions and with (0,0) origin.
+         * We use that for clipping iteration to the area that we requested even if the
+         * coverage gave us a larger area.
+         */
+        final Rectangle sliceAOI = new Rectangle(StrictMath.toIntExact(sliceMax[0] - sliceMin[0] + 1),
+                                                 StrictMath.toIntExact(sliceMax[1] - sliceMin[1] + 1));
+        /*
+         * If the given coordinates were in a subsampled space while the coverage is at full resolution,
+         * convert the coordinates to full resolution.
+         */
+        if (allowSubsamplings) {
+            sliceMin = sliceMin.clone();
+            sliceMax = sliceMax.clone();
+            for (int i=0; i<sliceMin.length; i++) {
+                sliceMin[i] = sliceMin[i] * subsampling[i] + subOffsets[i];
+                sliceMax[i] = sliceMax[i] * subsampling[i] + subOffsets[i];
+            }
+        }
+        RenderedImage image = coverage.render(new GridExtent(null, sliceMin, sliceMax, true));
+        /*
+         * The subsampling offsets were included in the extent given to above `render` method call, so in principle
+         * they should not be given again to `SubsampledImage` constructor.  However the `render` method is free to
+         * return an image with a larger extent, which may result in different offsets. The result can be "too much"
+         * offset. We want to compensate by subtracting the surplus. But because we can not have negative offsets,
+         * we shift the whole `sliceAOI` (which is equivalent to subtracting `subX|Y` in full resolution coordinates)
+         * and set the offset to the complement.
+         */
+        if (allowSubsamplings) {
+            final int subX = subsampling[0];
+            final int subY = subsampling[1];
+            int offX = StrictMath.floorMod(image.getMinX(), subX);
+            int offY = StrictMath.floorMod(image.getMinY(), subY);
+            if (offX != 0) {sliceAOI.x--; offX = subX - offX;}
+            if (offY != 0) {sliceAOI.y--; offY = subY - offY;}
+            image = SubsampledImage.create(image, subX, subY, offX, offY);
+            if (image == null) {
+                return null;
+            }
+        }
+        return new PixelIterator.Builder().setRegionOfInterest(sliceAOI).create(image);
+    }
+}
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/test/storage/SubsampledImage.java b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/SubsampledImage.java
new file mode 100644
index 0000000..93dbb07
--- /dev/null
+++ b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/SubsampledImage.java
@@ -0,0 +1,247 @@
+/*
+ * 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.test.storage;
+
+import java.awt.Point;
+import java.util.Arrays;
+import java.util.Vector;
+import java.awt.image.Raster;
+import java.awt.image.ColorModel;
+import java.awt.image.PixelInterleavedSampleModel;
+import java.awt.image.RenderedImage;
+import java.awt.image.SampleModel;
+import org.apache.sis.image.PlanarImage;
+
+import static java.lang.StrictMath.floorDiv;
+import static org.junit.Assert.*;
+
+
+/**
+ * An image which is a subsampling of another image.
+ * This image is a view; sample values are not copied.
+ *
+ * <h2>Limitations</h2>
+ * <ul>
+ *   <li>Sample model must be an instance of {@link PixelInterleavedSampleModel}.</li>
+ *   <li>Subsampling must be a divisor of tile size, except in dimensions having only one tile.</li>
+ *   <li>Conversion from source coordinates to target coordinates is a division only, without subsampling offsets.</li>
+ * </ul>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+final class SubsampledImage extends PlanarImage {
+    /**
+     * The image at full resolution.
+     */
+    private final RenderedImage source;
+
+    /**
+     * The subsampling to apply.
+     */
+    private final int subX, subY;
+
+    /**
+     * The subsampled model.
+     */
+    private final SampleModel model;
+
+    /**
+     * Creates a new image as a subsampling of the given image.
+     * The offsets given to this constructor apply to the {@link java.awt.image.DataBuffer},
+     * <strong>not</strong> to the image coordinate system. An offset needs to be added for
+     * accesses to the data buffers because those buffers are shared with the source image,
+     * and have a fixed origin regardless the origin of image coordinate systems.
+     *
+     * Current version does not take a subsampling offset for the image coordinate system
+     * because in the context of {@link CoverageReadConsistency} test, coordinates (0,0) of
+     * {@linkplain #source} image is the first pixel in the Area Of Interest specified by user,
+     * so the subsampling offsets on image coordinates are already applied at this stage.
+     *
+     * @param  source  the image at full resolution.
+     * @param  subX    the subsampling factor along X axis.
+     * @param  subY    the subsampling factor along Y axis.
+     * @param  offX    the offset along X axis in the shared data buffer.
+     * @param  offY    the offset along Y axis in the shared data buffer.
+     */
+    private SubsampledImage(final RenderedImage source, final int subX, final int subY, final int offX, final int offY) {
+        this.source = source;
+        this.subX   = subX;
+        this.subY   = subY;
+        final PixelInterleavedSampleModel sm = (PixelInterleavedSampleModel) source.getSampleModel();
+        final int   pixelStride    = sm.getPixelStride();
+        final int   scanlineStride = sm.getScanlineStride();
+        final int   offset         = pixelStride*offX + scanlineStride*offY;
+        final int[] offsets        = sm.getBandOffsets();
+        /*
+         * Conversion from subsampled coordinate x′ to full resolution x is:
+         *
+         *    x = (x′ × subsampling) + offset
+         *
+         * We simulate the offset addition by adding the value in the offset bands.
+         * PixelInterleavedSampleModel uses that value for computing array index as below:
+         *
+         *    y*scanlineStride + x*pixelStride + bandOffsets[b]
+         */
+        for (int i=0; i<offsets.length; i++) {
+            offsets[i] += offset;
+        }
+        model = new PixelInterleavedSampleModel(sm.getDataType(),
+                divExclusive(sm.getWidth(),  subX),
+                divExclusive(sm.getHeight(), subY),
+                pixelStride*subX, scanlineStride*subY, offsets);
+        /*
+         * Conditions documented in class javadoc.
+         */
+        if (getNumXTiles() > 1) assertEquals(0, sm.getWidth()  % subX);
+        if (getNumYTiles() > 1) assertEquals(0, sm.getHeight() % subY);
+    }
+
+    /**
+     * Returns an image as a subsampling of the given image.
+     * This method returns {@code null} if the arguments are valid but the image can not be
+     * created because of a restriction in {@link PixelInterleavedSampleModel} constructor.
+     *
+     * @param  source  the image at full resolution.
+     * @param  subX    the subsampling factor along X axis.
+     * @param  subY    the subsampling factor along Y axis.
+     * @param  offX    the subsampling offset along X axis.
+     * @param  offY    the subsampling offset along Y axis.
+     */
+    static RenderedImage create(final RenderedImage source, final int subX, final int subY, final int offX, final int offY) {
+        if (subX == 1 && subY == 1) {
+            return source;
+        } else {
+            final SubsampledImage image;
+            try {
+                image = new SubsampledImage(source, subX, subY, offX, offY);
+            } catch (IllegalArgumentException e) {
+                /*
+                 * PixelInterleavedSampleModel constructor has the following argument check:
+                 *
+                 *     if (pixelStride * width > scanlineStride) {
+                 *         throw new IllegalArgumentException("Pixel stride times width must be less than or equal to the scanline stride");
+                 *     }
+                 *
+                 * However this check rejects some valid layouts. Consider an image of size 16 × 3 pixels
+                 * with a single band and subsamplig factors (5,1). In the illustration below, "X" and "-"
+                 * are pixels from the source images and "X" are pixels retained in the subsampled image.
+                 * Note that the last column of the source image is included in the subsampled image.
+                 *
+                 *     X----X----X----X
+                 *     X----X----X----X
+                 *     X----X----X----X
+                 *
+                 * With a a pixelStride = 5, a width = 4 and scanlineStride = 16 we get!
+                 *
+                 *     pixelStride*w > scanlineStride
+                 *               5*4 > 16
+                 *                20 > 16           →  true  →  IllegalArgumentException is thrown.
+                 *
+                 * Above condition is equivalent to requiring image to be
+                 *
+                 *     like that:   X----X----X----X----
+                 *     instead of:  X----X----X----X
+                 *
+                 * The amended check below checks if there is enough room for storing the last sample value of a row,
+                 * ignoring the remaining of pixel stride that are just skipped.
+                 *
+                 *     pixelStride*(w-1) + maxBandOff >= scanlineStride
+                 *               5*(4-1) + 0          >= 16
+                 *                   15               >= 16     →  false  →  no exception thrown.
+                 */
+                final PixelInterleavedSampleModel sm = (PixelInterleavedSampleModel) source.getSampleModel();
+                final int pixelStride    = sm.getPixelStride() * subX;
+                final int scanlineStride = sm.getScanlineStride() * subY;
+                final int width = divExclusive(sm.getWidth(), subX);
+                if (pixelStride * width > scanlineStride) {
+                    final int maxBandOff = Arrays.stream(sm.getBandOffsets()).max().getAsInt();
+                    if (pixelStride * (width - 1) + maxBandOff < scanlineStride) {
+                        return null;
+                    }
+                }
+                throw e;
+            }
+            final String warning = image.verify();
+            if (warning != null && (source instanceof PlanarImage)) {
+                assertEquals(((PlanarImage) source).verify(), warning);
+            }
+            return image;
+        }
+    }
+
+    /**
+     * Returns the image at full resolution.
+     */
+    @Override
+    @SuppressWarnings("UseOfObsoleteCollectionType")
+    public Vector<RenderedImage> getSources() {
+        final Vector<RenderedImage> sources = new Vector<>(1);
+        sources.add(source);
+        return sources;
+    }
+
+    @Override public SampleModel getSampleModel() {return model;}
+    @Override public ColorModel  getColorModel()  {return source.getColorModel();}
+    @Override public int         getNumXTiles()   {return source.getNumXTiles();}
+    @Override public int         getNumYTiles()   {return source.getNumYTiles();}
+    @Override public int         getMinTileX()    {return source.getMinTileX();}
+    @Override public int         getMinTileY()    {return source.getMinTileY();}
+    @Override public int         getTileWidth()   {return divExclusive(source.getTileWidth(),  subX);}
+    @Override public int         getTileHeight()  {return divExclusive(source.getTileHeight(), subY);}
+    @Override public int         getWidth()       {return divExclusive(source.getWidth(),  subX);}
+    @Override public int         getHeight()      {return divExclusive(source.getHeight(), subY);}
+    @Override public int         getMinX()        {return divInclusive(source.getMinX(), subX);}
+    @Override public int         getMinY()        {return divInclusive(source.getMinY(), subY);}
+
+    /**
+     * Computes {@code (coordinate - offset) / subsampling} rounded toward 0.
+     * The subsampling offset is assumed 0 in current version.
+     *
+     * <div class="note"><b>Implementation note:</b>
+     * in principle we should subtract the <var>subsampling offset</var>. However that offset is
+     * zero in the context of {@link CoverageReadConsistency} test, because coordinates (0,0) of
+     * {@linkplain #source} image is the first pixel in the Area Of Interest specified by user,
+     * so there is no more offset at this stage. Note that we are talking about offset in image
+     * coordinate system, not to be confused with offset relative to the data bank
+     * (given to the {@link SampleModel} at construction time).</div>
+     */
+    private static int divInclusive(final int coordinate, final int subsampling) {
+        return floorDiv(coordinate, subsampling);
+    }
+
+    /**
+     * Computes {@code (coordinate - offset) / subsampling}, but handling the given {@code coordinate} as exclusive.
+     * The subsampling offset is assumed 0 in current version (see {@link #divInclusive(int, int)} for explanation).
+     */
+    private static int divExclusive(final int coordinate, final int subsampling) {
+        return floorDiv(coordinate - 1, subsampling) + 1;
+    }
+
+    /**
+     * Returns the tile at the given index.
+     */
+    @Override
+    public Raster getTile(final int tileX, final int tileY) {
+        final Raster tile = source.getTile(tileX, tileY);
+        return Raster.createRaster(model, tile.getDataBuffer(),
+                new Point(divInclusive(tile.getMinX(), subX),
+                          divInclusive(tile.getMinY(), subY)));
+    }
+}
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/test/storage/package-info.java b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/package-info.java
new file mode 100644
index 0000000..a674acd
--- /dev/null
+++ b/storage/sis-storage/src/test/java/org/apache/sis/test/storage/package-info.java
@@ -0,0 +1,34 @@
+/*
+ * 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.
+ */
+
+/**
+ * Helper classes for testing {@link org.apache.sis.storage.DataStore} implementations.
+ * The {@link org.apache.sis.test.storage.CoverageReadConsistency} class reads a coverage fully,
+ * then requests various sub-regions. Sub-regions are than compared to the corresponding regions
+ * in the full image. It is not a proof that data values are correct, but it shows at least that
+ * read operations are consistent for the test data set.
+ *
+ * <p>The classes in this package can be used by other modules (netCDF, GeoTIFF, <i>etc.</i>).
+ * Each module will need to provide a test data file in the format to be tested.
+ * That test data shall be small enough for fitting in memory.</p>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+package org.apache.sis.test.storage;

Mime
View raw message