diff --git a/CMakeLists.txt b/CMakeLists.txt
index bc24a2fee9..2b9363006c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -4,7 +4,7 @@
 
 cmake_minimum_required (VERSION 3.18.2...4.0)
 
-set (OpenImageIO_VERSION "3.1.2.0")
+set (OpenImageIO_VERSION "3.1.3.0")
 set (OpenImageIO_VERSION_OVERRIDE "" CACHE STRING
      "Version override (use with caution)!")
 mark_as_advanced (OpenImageIO_VERSION_OVERRIDE)
diff --git a/src/doc/imageoutput.rst b/src/doc/imageoutput.rst
index 8e3c9b2ecb..e6d2e7d20b 100644
--- a/src/doc/imageoutput.rst
+++ b/src/doc/imageoutput.rst
@@ -246,26 +246,20 @@ individual plugin.
 
 .. tabs::
 
-   .. code-tab:: c++
-
-      unsigned char tile[tilesize*tilesize*channels];
-      int z = 0;   // Always zero for 2D images
-      for (int y = 0;  y < yres;  y += tilesize) {
-          for (int x = 0;  x < xres;  x += tilesize) {
-              ... generate data in tile[] ..
-              out->write_tile (x, y, z, TypeDesc::UINT8, tile);
-          }
-      }
-      out->close ();
+   .. tab:: C++
+      .. literalinclude:: ../../testsuite/docs-examples-cpp/src/docs-examples-imageoutput.cpp
+          :language: c++
+          :start-after: BEGIN-imageoutput-tiles
+          :end-before: END-imageoutput-tiles
+          :dedent: 4
 
-   .. code-tab:: py
+   .. tab:: Python
 
-      z = 0  # Always zero for 2D images
-      for y in range(0, yres, tilesize) :
-          for x in range(0, xres, tilesize) :
-              # ... generate data in tile[][][] ..
-              out.write_tile (x, y, z, tile)
-      out.close ()
+      .. literalinclude:: ../../testsuite/docs-examples-python/src/docs-examples-imageoutput.py
+          :language: py
+          :start-after: BEGIN-imageoutput-tiles
+          :end-before: END-imageoutput-tiles
+          :dedent: 8
 
 The first three arguments to ``write_tile()`` specify which tile is being
 written by the pixel coordinates of any pixel contained in the tile: *x*
diff --git a/src/include/OpenImageIO/imageio.h b/src/include/OpenImageIO/imageio.h
index 1255021f50..c717fca241 100644
--- a/src/include/OpenImageIO/imageio.h
+++ b/src/include/OpenImageIO/imageio.h
@@ -60,6 +60,10 @@ using stride_t = int64_t;
 inline constexpr stride_t AutoStride = std::numeric_limits<stride_t>::min();
 #endif
 
+// Signal that this version of ImageBuf has constructors from spans
+#define OIIO_IMAGEINPUT_IMAGE_SPAN_SUPPORT 1
+#define OIIO_IMAGEOUTPUT_IMAGE_SPAN_SUPPORT 1
+
 
 
 /// Pointer to a function called periodically by read_image and
@@ -1234,8 +1238,343 @@ class OIIO_API ImageInput {
     }
 #endif
 
+    // clang-format on
+
+    /// @name Reading pixels ("safe" methods with bounded spans)
+    ///
+    /// Common features of all the `read` methods:
+    ///
+    /// * There is a base case that takes a `image_span<byte>` describing
+    ///   untyped memory layout and a `TypeDesc` describing the data type
+    ///   that the values should be converted to (or TypeUnknown to keep
+    ///   the data in its "native" file types with no conversion).
+    ///
+    /// * The type-aware versions that accept an `image_span<T>` (for
+    ///   potentially non-contiguous data) or `span<T>` (for contiguous data)
+    ///   and understand to convert the data into the given `T` type.
+    ///
+    /// * The image_span (in either case) includes the memory bounds and
+    ///   stride lengths (in bytes) between channels, pixels, scanlines, and
+    ///   volumetric slices.
+    ///
+    /// * Any *range* parameters (such as `ybegin` and `yend`) describe a
+    ///   "half open interval", meaning that `begin` is the first item and
+    ///   `end` is *one past the last item*. That means that the number of
+    ///   items is `end - begin`.
+    ///
+    /// * For ordinary 2D (non-volumetric) images, any `z` or `zbegin`
+    ///   coordinates should be 0 and any `zend` should be 1, indicating
+    ///   that only a single image "plane" exists.
+    ///
+    /// * Some read methods take a channel range [chbegin,chend) to allow
+    ///   reading of a contiguous subset of channels (chbegin=0,
+    ///   chend=spec.nchannels reads all channels).
+    ///
+    /// * ImageInput readers are expected to give the appearance of random
+    ///   access -- in other words, if it can't randomly seek to the given
+    ///   scanline or tile, it should transparently close, reopen, and
+    ///   sequentially read through prior scanlines.
+    ///
+    /// * All read functions return `true` for success, `false` for failure
+    ///   (after which a call to `geterror()` may retrieve a specific error
+    ///   message).
+    ///
+
+    /// Read the entire image of `spec.width x spec.height x spec.depth`
+    /// pixels into a buffer with the given strides and in the desired
+    /// data format.
+    ///
+    /// Depending on the spec, this will read either all tiles or all
+    /// scanlines. Assume that data points to a layout in row-major order.
+    ///
+    /// This version of read_image, because it passes explicit subimage and
+    /// miplevel, does not require a separate call to seek_subimage, and is
+    /// guaranteed to be thread-safe against other concurrent calls to any
+    /// of the read_* methods that take an explicit subimage/miplevel (but
+    /// not against any other ImageInput methods).
+    ///
+    /// Added in OIIO 3.1, this is the "safe" preferred alternative to
+    /// the version of read_image that takes raw pointers.
+    ///
+    /// @param  subimage    The subimage to read from (starting with 0).
+    /// @param  miplevel    The MIP level to read (0 is the highest
+    ///                     resolution level).
+    /// @param  chbegin/chend
+    ///                     The channel range to read. If chend is -1, it
+    ///                     will automatically be set to spec.nchannels.
+    /// @param  format      A TypeDesc describing the type of the pixel data
+    ///                     that `data`'s memory contains. Use `TypeUnknown`
+    ///                     to indicate that you want the native data format
+    ///                     with no type conversion.
+    /// @param  data        An `image_span<byte>` describing the memory
+    ///                     extent of the data buffer and including the sizes
+    ///                     and byte strides for each dimension (channel, x,
+    ///                     y, and z).
+    /// @returns            `true` upon success, or `false` upon failure.
+    ///
+    virtual bool read_image(int subimage, int miplevel, int chbegin, int chend,
+                            TypeDesc format, const image_span<std::byte>& data);
+
+    /// A version of `read_image()` taking an `image_span<T>`, where the type
+    /// of the underlying data is `T`.  This is a convenience wrapper around
+    /// the `read_image()` that takes an `image_span<std::byte>`.
+    template<typename T>
+    bool read_image(int subimage, int miplevel, int chbegin, int chend,
+                    const image_span<T>& data)
+    {
+        static_assert(!std::is_const_v<T>,
+                      "read_image() does not accept image_span<const T>");
+        return read_image(subimage, miplevel, chbegin, chend,
+                          TypeDescFromC<T>::value(),
+                          as_image_span_writable_bytes(data));
+    }
+
+    /// A version of `read_image()` taking a `span<T>`, which assumes
+    /// contiguous strides in all dimensions. This is a convenience wrapper
+    /// around the `read_image()` that takes an `image_span<T>`.
+    template<typename T>
+    bool read_image(int subimage, int miplevel, int chbegin, int chend,
+                    span<T> data)
+    {
+        static_assert(!std::is_const_v<T>,
+                      "read_image() does not accept span<const T>");
+        ImageSpec spec = spec_dimensions(subimage, miplevel);
+        if (chend < 0 || chend > spec.nchannels)
+            chend = spec.nchannels;
+        auto ispan = image_span<T>(data.data(), chend - chbegin, spec.width,
+                                   spec.height, spec.depth);
+        OIIO_DASSERT(data.size_bytes() == ispan.size_bytes()
+                     && ispan.is_contiguous());
+        return read_image(subimage, miplevel, chbegin, chend, ispan);
+    }
+
+    /// Read multiple scanlines that include pixels (*,y) for all ybegin <= y
+    /// < yend in the specified subimage and mip level, into `data`, using the
+    /// strides given and converting to the requested data `format`
+    /// (TypeUnknown indicates no conversion, just copy native data types).
+    /// Only channels [chbegin,chend) will be read/copied (chbegin=0,
+    /// chend=nchannels reads all channels, and chend=-1 also means the same
+    /// thing).
+    ///
+    /// This method, does not require a separate call to seek_subimage, and is
+    /// guaranteed to be thread-safe against other concurrent calls to any of
+    /// the read_* methods that take an explicit subimage/miplevel (but not
+    /// against any other non-thread-safe ImageInput methods).
+    ///
+    /// @param  subimage    The subimage to read from (starting with 0).
+    /// @param  miplevel    The MIP level to read (0 is the highest
+    ///                     resolution level).
+    /// @param  ybegin/yend The y range of the scanlines to read.
+    /// @param  chbegin/chend
+    ///                     The channel range to read. If chend is -1, it
+    ///                     will automatically be set to spec.nchannels.
+    /// @param  format      A TypeDesc describing the type of the pixel data
+    ///                     that `data`'s memory contains. Use `TypeUnknown`
+    ///                     to indicate that you want the native data format
+    ///                     with no type conversion.
+    /// @param  data        An `image_span<byte>` describing the memory
+    ///                     extent of the data buffer and including the sizes
+    ///                     and byte strides for each dimension (channel, x,
+    ///                     y, and z).
+    /// @returns            `true` upon success, or `false` upon failure.
+    ///
+    /// Added in OIIO 3.1, this is the "safe" preferred alternative to
+    /// the version of read_scanlines that takes raw pointers.
+    ///
+    virtual bool read_scanlines(int subimage, int miplevel, int ybegin,
+                                int yend, int chbegin, int chend,
+                                TypeDesc format,
+                                const image_span<std::byte>& data);
+
+    /// A version of `read_scanlines()` taking an `image_span<T>`, where the
+    /// type of the underlying data is `T`.  This is a convenience wrapper
+    /// around the `read_scanlines()` that takes an `image_span<std::byte>`.
+    template<typename T>
+    bool read_scanlines(int subimage, int miplevel, int ybegin, int yend,
+                        int chbegin, int chend, const image_span<T>& data)
+    {
+        static_assert(!std::is_const_v<T>,
+                      "read_scanlines() does not accept span<const T>");
+        // reduce to type + image_span<byte>
+        return read_scanlines(subimage, miplevel, ybegin, yend, chbegin, chend,
+                              TypeDescFromC<T>::value(),
+                              as_image_span_writable_bytes(data));
+    }
+
+    /// A version of `read_scanlines()` taking a `span<T>`, which assumes
+    /// contiguous strides in all dimensions. This is a convenience wrapper
+    /// around the `read_scanlines()` that takes an `image_span<T>`.
+    template<typename T>
+    bool read_scanlines(int subimage, int miplevel, int ybegin, int yend,
+                        int chbegin, int chend, span<T> data)
+    {
+        static_assert(!std::is_const_v<T>,
+                      "read_scanlines() does not accept span<const T>");
+        // reduce to type + image_span<byte>
+        ImageSpec spec = spec_dimensions(subimage, miplevel);
+        if (chend < 0 || chend > spec.nchannels)
+            chend = spec.nchannels;
+        auto ispan = image_span<T>(data.data(), chend - chbegin, spec.width,
+                                   spec.height, 1);
+        return read_scanlines(subimage, miplevel, ybegin, yend, chbegin, chend,
+                              ispan);
+    }
+
+    /// Read the block of multiple tiles that include all pixels in
+    ///
+    ///     [xbegin,xend) X [ybegin,yend) X [zbegin,zend)
+    ///
+    /// The begin/end pairs must correctly delineate tile boundaries, with
+    /// the exception that it may also be the end of the image data if the
+    /// image resolution is not a whole multiple of the tile size.
+    ///
+    /// This version of read_tiles, because it passes explicit subimage and
+    /// miplevel, does not require a separate call to seek_subimage, and is
+    /// guaranteed to be thread-safe against other concurrent calls to any
+    /// of the read_* methods that take an explicit subimage/miplevel (but
+    /// not against any other ImageInput methods).
+    ///
+    /// @param  subimage    The subimage to read from (starting with 0).
+    /// @param  miplevel    The MIP level to read (0 is the highest
+    ///                     resolution level).
+    /// @param  xbegin/xend The x range of the pixels covered by the group
+    ///                     of tiles being read.
+    /// @param  ybegin/yend The y range of the pixels covered by the tiles.
+    /// @param  zbegin/zend The z range of the pixels covered by the tiles
+    ///                     (for a 2D image, zbegin=0 and zend=1).
+    /// @param  chbegin/chend
+    ///                     The channel range to read.
+    /// @param  format      A TypeDesc describing the type of the pixel data
+    ///                     that `data`'s memory contains. Use `TypeUnknown`
+    ///                     to indicate that you want the native data format
+    ///                     with no type conversion.
+    /// @param  data        An `image_span<byte>` describing the memory
+    ///                     extent of the data buffer and including the sizes
+    ///                     and byte strides for each dimension (channel, x,
+    ///                     y, and z).
+    /// @returns            `true` upon success, or `false` upon failure.
+    ///
+    /// Added in OIIO 3.1, this is the "safe" preferred alternative to
+    /// the version of read_tiles that takes raw pointers.
+    ///
+    /// @note The call will fail if the image is not tiled, or if the pixel
+    /// ranges do not fall along tile (or image) boundaries, or if it is not
+    /// a valid tile range.
+    virtual bool read_tiles(int subimage, int miplevel, int xbegin, int xend,
+                            int ybegin, int yend, int zbegin, int zend,
+                            int chbegin, int chend, TypeDesc format,
+                            const image_span<std::byte>& data);
+
+    /// A version of `read_tiles()` taking an `image_span<T>`, where the type
+    /// of the underlying data is `T`.  This is a convenience wrapper around
+    /// the `read_tiles()` that takes an `image_span<const std::byte>`.
+    template<typename T>
+    bool read_tiles(int subimage, int miplevel, int xbegin, int xend,
+                    int ybegin, int yend, int zbegin, int zend, int chbegin,
+                    int chend, const image_span<T>& data)
+    {
+        return read_tiles(subimage, miplevel, xbegin, xend, ybegin, yend,
+                          zbegin, zend, chbegin, chend,
+                          TypeDescFromC<T>::value(),
+                          as_image_span_writable_bytes(data));
+    }
+
+    /// A version of `read_tiles()` taking a `cspan<T>`, which assumes
+    /// contiguous strides in all dimensions. This is a convenience wrapper
+    /// around the `read_tiles()` that takes an `image_span<const T>`.
+    template<typename T>
+    bool read_tiles(int subimage, int miplevel, int xbegin, int xend,
+                    int ybegin, int yend, int zbegin, int zend, int chbegin,
+                    int chend, span<T> data)
+    {
+        static_assert(!std::is_const_v<T>,
+                      "read_tiles() does not accept span<const T>");
+        // reduce to type + image_span<byte>
+        ImageSpec spec = spec_dimensions(subimage, miplevel);
+        if (chend < 0 || chend > spec.nchannels)
+            chend = spec.nchannels;
+        auto ispan = image_span<T>(data.data(), chend - chbegin, xend - xbegin,
+                                   yend - ybegin, zend - zbegin);
+        OIIO_DASSERT(data.size_bytes() >= ispan.size_bytes()
+                     && ispan.is_contiguous());
+        return read_tiles(subimage, miplevel, xbegin, xend, ybegin, yend,
+                          zbegin, zend, chbegin, chend, ispan);
+    }
+
+    // clang-format off
+
+    /// Read deep scanlines containing pixels (*,y,z), for all y in the
+    /// range [ybegin,yend) into `deepdata`. This will fail if it is not a
+    /// deep file.
+    ///
+    /// @param  subimage    The subimage to read from (starting with 0).
+    /// @param  miplevel    The MIP level to read (0 is the highest
+    ///                     resolution level).
+    /// @param  chbegin/chend
+    ///                     The channel range to read.
+    /// @param  ybegin/yend The y range of the scanlines being passed.
+    /// @param  z           The z coordinate of the scanline.
+    /// @param  deepdata    A `DeepData` object into which the data for
+    ///                     these scanlines will be placed.
+    /// @returns            `true` upon success, or `false` upon failure.
+    virtual bool read_native_deep_scanlines (int subimage, int miplevel,
+                                             int ybegin, int yend, int z,
+                                             int chbegin, int chend,
+                                             DeepData &deepdata);
+
+    /// Read into `deepdata` the block of native deep data tiles that
+    /// include all pixels and channels specified by pixel range.
+    ///
+    /// @param  subimage    The subimage to read from (starting with 0).
+    /// @param  miplevel    The MIP level to read (0 is the highest
+    ///                     resolution level).
+    /// @param  xbegin/xend The x range of the pixels covered by the group
+    ///                     of tiles being read.
+    /// @param  ybegin/yend The y range of the pixels covered by the tiles.
+    /// @param  zbegin/zend The z range of the pixels covered by the tiles
+    ///                     (for a 2D image, zbegin=0 and zend=1).
+    /// @param  chbegin/chend
+    ///                     The channel range to read.
+    /// @param  deepdata    A `DeepData` object into which the data for
+    ///                     these tiles will be placed.
+    /// @returns            `true` upon success, or `false` upon failure.
+    ///
+    /// @note The call will fail if the image is not tiled, or if the pixel
+    /// ranges do not fall along tile (or image) boundaries, or if it is not
+    /// a valid tile range.
+    virtual bool read_native_deep_tiles (int subimage, int miplevel,
+                                         int xbegin, int xend,
+                                         int ybegin, int yend,
+                                         int zbegin, int zend,
+                                         int chbegin, int chend,
+                                         DeepData &deepdata);
+
+    /// Read the entire deep data image of spec.width x spec.height x
+    /// spec.depth pixels, all channels, into `deepdata`.
+    ///
+    /// @param  subimage    The subimage to read from (starting with 0).
+    /// @param  miplevel    The MIP level to read (0 is the highest
+    ///                     resolution level).
+    /// @param  deepdata    A `DeepData` object into which the data for
+    ///                     the image will be placed.
+    /// @returns            `true` upon success, or `false` upon failure.
+    virtual bool read_native_deep_image (int subimage, int miplevel,
+                                         DeepData &deepdata);
+
+    /// @}
+
     /// @{
-    /// @name Reading pixels
+    /// @name Reading pixels (unsafe methods with pointers and strides)
+    ///
+    /// These methods are the "unsafe" versions of the `read` methods, which
+    /// take raw pointers and strides.  They are provided for backwards
+    /// compatibility with older code, but the preferred interface is to use
+    /// the `read` methods that take `image_span<T>` or `span<T>` arguments,
+    /// which are bounds-safe and type-aware.
+    ///
+    /// These pointer-based versions are considered "soft-deprecated" in
+    /// OpenImageIO 3.1, will be marked/warned as deprecated in 3.2, and will
+    /// be removed in 4.0.
     ///
     /// Common features of all the `read` methods:
     ///
@@ -1457,7 +1796,7 @@ class OIIO_API ImageInput {
     ///                     resolution level).
     /// @param  chbegin/chend
     ///                     The channel range to read. If chend is -1, it
-    ///                     will be set to spec.nchannels.
+    ///                     will automatically be set to spec.nchannels.
     /// @param  format      A TypeDesc describing the type of `data`.
     /// @param  data        Pointer to the pixel data.
     /// @param  xstride/ystride/zstride
@@ -1475,64 +1814,6 @@ class OIIO_API ImageInput {
                              ProgressCallback progress_callback=NULL,
                              void *progress_callback_data=NULL);
 
-    /// Read deep scanlines containing pixels (*,y,z), for all y in the
-    /// range [ybegin,yend) into `deepdata`. This will fail if it is not a
-    /// deep file.
-    ///
-    /// @param  subimage    The subimage to read from (starting with 0).
-    /// @param  miplevel    The MIP level to read (0 is the highest
-    ///                     resolution level).
-    /// @param  chbegin/chend
-    ///                     The channel range to read.
-    /// @param  ybegin/yend The y range of the scanlines being passed.
-    /// @param  z           The z coordinate of the scanline.
-    /// @param  deepdata    A `DeepData` object into which the data for
-    ///                     these scanlines will be placed.
-    /// @returns            `true` upon success, or `false` upon failure.
-    virtual bool read_native_deep_scanlines (int subimage, int miplevel,
-                                             int ybegin, int yend, int z,
-                                             int chbegin, int chend,
-                                             DeepData &deepdata);
-
-    /// Read into `deepdata` the block of native deep data tiles that
-    /// include all pixels and channels specified by pixel range.
-    ///
-    /// @param  subimage    The subimage to read from (starting with 0).
-    /// @param  miplevel    The MIP level to read (0 is the highest
-    ///                     resolution level).
-    /// @param  xbegin/xend The x range of the pixels covered by the group
-    ///                     of tiles being read.
-    /// @param  ybegin/yend The y range of the pixels covered by the tiles.
-    /// @param  zbegin/zend The z range of the pixels covered by the tiles
-    ///                     (for a 2D image, zbegin=0 and zend=1).
-    /// @param  chbegin/chend
-    ///                     The channel range to read.
-    /// @param  deepdata    A `DeepData` object into which the data for
-    ///                     these tiles will be placed.
-    /// @returns            `true` upon success, or `false` upon failure.
-    ///
-    /// @note The call will fail if the image is not tiled, or if the pixel
-    /// ranges do not fall along tile (or image) boundaries, or if it is not
-    /// a valid tile range.
-    virtual bool read_native_deep_tiles (int subimage, int miplevel,
-                                         int xbegin, int xend,
-                                         int ybegin, int yend,
-                                         int zbegin, int zend,
-                                         int chbegin, int chend,
-                                         DeepData &deepdata);
-
-    /// Read the entire deep data image of spec.width x spec.height x
-    /// spec.depth pixels, all channels, into `deepdata`.
-    ///
-    /// @param  subimage    The subimage to read from (starting with 0).
-    /// @param  miplevel    The MIP level to read (0 is the highest
-    ///                     resolution level).
-    /// @param  deepdata    A `DeepData` object into which the data for
-    ///                     the image will be placed.
-    /// @returns            `true` upon success, or `false` upon failure.
-    virtual bool read_native_deep_image (int subimage, int miplevel,
-                                         DeepData &deepdata);
-
     /// @}
 
     /// @{
@@ -2251,7 +2532,7 @@ class OIIO_API ImageOutput {
     ///   no conversion needed).
     ///
     /// * The type-aware versions that accept an `image_span<T>` (for
-    ///   optionally non-contiguous data) or `span<T>` (for contiguous data)
+    ///   potentially non-contiguous data) or `span<T>` (for contiguous data)
     ///   and understand to convert the data from the given `T` type.
     ///
     /// * The image_span (in either case) includes the memory bounds and
diff --git a/src/libOpenImageIO/imageinput.cpp b/src/libOpenImageIO/imageinput.cpp
index 766701e759..efdc8548ac 100644
--- a/src/libOpenImageIO/imageinput.cpp
+++ b/src/libOpenImageIO/imageinput.cpp
@@ -285,6 +285,37 @@ ImageInput::read_scanline(int y, int z, TypeDesc format, void* data,
 
 
 
+bool
+ImageInput::read_scanlines(int subimage, int miplevel, int ybegin, int yend,
+                           int chbegin, int chend, TypeDesc format,
+                           const image_span<std::byte>& data)
+{
+    ImageSpec spec = spec_dimensions(subimage, miplevel);
+    if (chend < 0 || chend > spec.nchannels)
+        chend = spec.nchannels;
+    if (chbegin < 0 || chbegin >= chend) {
+        errorfmt("read_scanlines: invalid channel range [{},{})", chbegin,
+                 chend);
+        return false;
+    }
+    size_t isize = (format == TypeUnknown
+                        ? spec.pixel_bytes(chbegin, chend, true /*native*/)
+                        : format.size() * (chend - chbegin))
+                   * size_t(spec.width);
+    if (isize != data.size_bytes()) {
+        errorfmt(
+            "read_scanlines: Buffer size is incorrect ({} bytes vs {} needed)",
+            isize, data.size_bytes());
+        return false;
+    }
+
+    // Default implementation (for now): call the old pointer+stride
+    return read_scanlines(subimage, miplevel, ybegin, yend, 0, chbegin, chend,
+                          format, data.data(), data.xstride());
+}
+
+
+
 bool
 ImageInput::read_scanlines(int subimage, int miplevel, int ybegin, int yend,
                            int z, int chbegin, int chend, TypeDesc format,
@@ -610,6 +641,38 @@ ImageInput::read_tile(int x, int y, int z, TypeDesc format, void* data,
 
 
 
+bool
+ImageInput::read_tiles(int subimage, int miplevel, int xbegin, int xend,
+                       int ybegin, int yend, int zbegin, int zend, int chbegin,
+                       int chend, TypeDesc format,
+                       const image_span<std::byte>& data)
+{
+    ImageSpec spec = spec_dimensions(subimage, miplevel);
+    if (chend < 0 || chend > spec.nchannels)
+        chend = spec.nchannels;
+    if (chbegin < 0 || chbegin >= chend) {
+        errorfmt("read_tiles: invalid channel range [{},{})", chbegin, chend);
+        return false;
+    }
+    size_t isize = (format == TypeUnknown
+                        ? spec.pixel_bytes(chbegin, chend, true /*native*/)
+                        : format.size() * (chend - chbegin))
+                   * size_t(xend - xbegin) * size_t(yend - ybegin)
+                   * size_t(zend - zbegin);
+    if (isize != data.size_bytes()) {
+        errorfmt("read_tiles: Buffer size is incorrect ({} bytes vs {} needed)",
+                 isize, data.size_bytes());
+        return false;
+    }
+
+    // Default implementation (for now): call the old pointer+stride
+    return read_tiles(subimage, miplevel, ybegin, yend, xbegin, xend, zbegin,
+                      zend, chbegin, chend, format, data.data(),
+                      data.xstride());
+}
+
+
+
 bool
 ImageInput::read_tiles(int subimage, int miplevel, int xbegin, int xend,
                        int ybegin, int yend, int zbegin, int zend, int chbegin,
@@ -1095,6 +1158,110 @@ ImageInput::read_image(int subimage, int miplevel, int chbegin, int chend,
 
 
 
+bool
+ImageInput::read_image(int subimage, int miplevel, int chbegin, int chend,
+                       TypeDesc format, const image_span<std::byte>& data)
+{
+#if 0
+    ImageSpec spec = spec_dimensions(subimage, miplevel);
+    if (chend < 0 || chend > spec.nchannels)
+        chend = spec.nchannels;
+    size_t isize = (format == TypeUnknown
+                        ? spec.pixel_bytes(chbegin, chend, true /*native*/)
+                        : format.size() * (chend - chbegin))
+                   * spec.image_pixels();
+    if (isize != data.size_bytes()) {
+        errorfmt("read_image: Buffer size is incorrect ({} bytes vs {} needed)",
+                 sz, data.size_bytes());
+        return false;
+    }
+
+    // Default implementation (for now): call the old pointer+stride
+    return read_image(subimage, miplevel, chbegin, chend, format, data.data(),
+                      data.xstride(), data.ystride(), data.zstride());
+#else
+    pvt::LoggedTimer logtime("II::read_image");
+    ImageSpec spec;
+    int rps = 0;
+    {
+        // We need to lock briefly to retrieve rps from the spec
+        lock_guard lock(*this);
+        if (!seek_subimage(subimage, miplevel))
+            return false;
+        // Copying the dimensions of the designated subimage/miplevel to a
+        // local `spec` means that we can release the lock!  (Calls to
+        // read_native_* will internally lock again if necessary.)
+        spec.copy_dimensions(m_spec);
+        // For scanline files, we also need one piece of metadata
+        if (!spec.tile_width)
+            rps = m_spec.get_int_attribute("tiff:RowsPerStrip", 64);
+    }
+    if (spec.image_bytes() < 1) {
+        errorfmt("Invalid image size {} x {} ({} chans)", m_spec.width,
+                 m_spec.height, m_spec.nchannels);
+        return false;
+    }
+
+    if (chend < 0 || chend > spec.nchannels)
+        chend = spec.nchannels;
+    if (chbegin < 0 || chbegin >= chend) {
+        errorfmt("read_image: invalid channel range [{},{})", chbegin, chend);
+        return false;
+    }
+    int nchans         = chend - chbegin;
+    bool native        = (format == TypeUnknown);
+    size_t pixel_bytes = native ? spec.pixel_bytes(chbegin, chend, native)
+                                : (format.size() * nchans);
+    size_t isize       = pixel_bytes * spec.image_pixels();
+    if (isize != data.size_bytes()) {
+        errorfmt("read_image: Buffer size is incorrect ({} bytes vs {} needed)",
+                 isize, data.size_bytes());
+        return false;
+    }
+
+    bool ok = true;
+    if (spec.tile_width) {  // Tiled image -- rely on read_tiles
+        // Read in chunks of a whole row of tiles at once. If tiles are
+        // 64x64, a 2k image has 32 tiles across. That's fine for now (for
+        // parallelization purposes), but as typical core counts increase,
+        // we may someday want to revisit this to batch multiple rows.
+        for (int z = 0; z < spec.depth; z += spec.tile_depth) {
+            int zend = std::min(z + spec.z + spec.tile_depth,
+                                spec.z + spec.depth);
+            for (int y = 0; y < spec.height && ok; y += spec.tile_height) {
+                int yend = std::min(y + spec.y + spec.tile_height,
+                                    spec.y + spec.height);
+                ok &= read_tiles(subimage, miplevel, spec.x,
+                                 spec.x + spec.width, y + spec.y, yend,
+                                 z + spec.z, zend, chbegin, chend, format,
+                                 data.subspan(spec.x, spec.x + spec.width,
+                                              y + spec.y, yend, z + spec.z,
+                                              zend));
+            }
+        }
+    } else {  // Scanline image -- rely on read_scanlines.
+        // Split into reasonable chunks -- try to use around 64 MB or the
+        // oiio_read_chunk value, which ever is bigger, but also round up to
+        // a multiple of the TIFF rows per strip (or 64).
+        int chunk = std::max(1, (1 << 26) / int(spec.scanline_bytes(true)));
+        chunk     = std::max(chunk, int(oiio_read_chunk));
+        chunk     = round_to_multiple(chunk, rps);
+        for (int z = 0; z < spec.depth; ++z) {
+            for (int y = 0; y < spec.height && ok; y += chunk) {
+                int yend = std::min(y + spec.y + chunk, spec.y + spec.height);
+                ok &= read_scanlines(subimage, miplevel, y + spec.y, yend,
+                                     chbegin, chend, format,
+                                     data.subspan(spec.x, spec.x + spec.width,
+                                                  y + spec.y, yend));
+            }
+        }
+    }
+    return ok;
+#endif
+}
+
+
+
 bool
 ImageInput::read_native_deep_scanlines(int /*subimage*/, int /*miplevel*/,
                                        int /*ybegin*/, int /*yend*/, int /*z*/,
diff --git a/testsuite/docs-examples-cpp/ref/out-arm.txt b/testsuite/docs-examples-cpp/ref/out-arm.txt
index 72373f6058..bbb53c3af6 100644
--- a/testsuite/docs-examples-cpp/ref/out-arm.txt
+++ b/testsuite/docs-examples-cpp/ref/out-arm.txt
@@ -144,3 +144,5 @@ Comparing "simple.tif" and "ref/simple.tif"
 PASS
 Comparing "scanlines.tif" and "ref/scanlines.tif"
 PASS
+Comparing "tiles.tif" and "ref/tiles.tif"
+PASS
diff --git a/testsuite/docs-examples-cpp/ref/out.txt b/testsuite/docs-examples-cpp/ref/out.txt
index 58889fdab0..5854393bba 100644
--- a/testsuite/docs-examples-cpp/ref/out.txt
+++ b/testsuite/docs-examples-cpp/ref/out.txt
@@ -144,3 +144,5 @@ Comparing "simple.tif" and "ref/simple.tif"
 PASS
 Comparing "scanlines.tif" and "ref/scanlines.tif"
 PASS
+Comparing "tiles.tif" and "ref/tiles.tif"
+PASS
diff --git a/testsuite/docs-examples-cpp/ref/tiles.tif b/testsuite/docs-examples-cpp/ref/tiles.tif
new file mode 100644
index 0000000000..1624c6ca92
Binary files /dev/null and b/testsuite/docs-examples-cpp/ref/tiles.tif differ
diff --git a/testsuite/docs-examples-cpp/run.py b/testsuite/docs-examples-cpp/run.py
index ffae73e247..e856d190ef 100755
--- a/testsuite/docs-examples-cpp/run.py
+++ b/testsuite/docs-examples-cpp/run.py
@@ -100,7 +100,7 @@
 # and need the images checked into the ref directory.
 outputs = [
     # Outputs from the ImageOutput chapter:
-    "simple.tif", "scanlines.tif",
+    "simple.tif", "scanlines.tif", "tiles.tif",
     # Outputs from the ImageInput chapter:
 
     # Outputs from the ImageBuf chapter:
diff --git a/testsuite/docs-examples-cpp/src/docs-examples-imageinput.cpp b/testsuite/docs-examples-cpp/src/docs-examples-imageinput.cpp
index 7a062ef102..453d3d4cd2 100644
--- a/testsuite/docs-examples-cpp/src/docs-examples-imageinput.cpp
+++ b/testsuite/docs-examples-cpp/src/docs-examples-imageinput.cpp
@@ -49,13 +49,14 @@ simple_read()
     int xres              = spec.width;
     int yres              = spec.height;
     int nchannels         = spec.nchannels;
-    auto pixels           = std::unique_ptr<unsigned char[]>(
-        new unsigned char[xres * yres * nchannels]);
-    inp->read_image(0, 0, 0, nchannels, TypeDesc::UINT8, &pixels[0]);
+    std::vector<uint8_t> pixels(xres * yres * nchannels);
+    inp->read_image(0 /*subimage*/, 0 /*miplevel*/, 0 /*chbegin*/,
+                    nchannels /*chend*/, make_span(pixels));
     inp->close();
 }
 // END-imageinput-simple
 
+
 void
 scanlines_read()
 {
@@ -65,10 +66,11 @@ scanlines_read()
     auto inp              = ImageInput::open(filename);
     const ImageSpec& spec = inp->spec();
     if (spec.tile_width == 0) {
-        auto scanline = std::unique_ptr<unsigned char[]>(
-            new unsigned char[spec.width * spec.nchannels]);
-        for (int y = 0; y < spec.height; ++y) {
-            inp->read_scanline(y, 0, TypeDesc::UINT8, &scanline[0]);
+        std::vector<uint8_t> scanline(spec.width * spec.nchannels);
+        for (int y = spec.y; y < spec.y + spec.height; ++y) {
+            inp->read_scanlines(0 /*subimage*/, 0 /*miplevel*/, y, y + 1,
+                                0 /*chbegin*/, spec.nchannels /*chend*/,
+                                make_span(scanline));
             // ... process data in scanline[0..width*channels-1] ...
         }
     } else {
@@ -93,12 +95,22 @@ tiles_read()
     } else {
         // Tiles
         int tilesize = spec.tile_width * spec.tile_height;
-        auto tile    = std::unique_ptr<unsigned char[]>(
-            new unsigned char[tilesize * spec.nchannels]);
-        for (int y = 0; y < spec.height; y += spec.tile_height) {
-            for (int x = 0; x < spec.width; x += spec.tile_width) {
-                inp->read_tile(x, y, 0, TypeDesc::UINT8, &tile[0]);
+        std::vector<uint8_t> tile(tilesize * spec.nchannels);
+        for (int y = spec.y; y < spec.y + spec.height; y += spec.tile_height) {
+            for (int x = spec.x; x < spec.x + spec.width;
+                 x += spec.tile_width) {
+                inp->read_tiles(0 /*subimage*/, 0 /*miplevel*/, x,
+                                std::min(x + spec.tile_width, spec.width), y,
+                                std::min(y + spec.tile_height, spec.height), 0,
+                                1, 0 /*chbegin*/, spec.nchannels /*chend*/,
+                                make_span(tile));
                 // ... process the pixels in tile[] ..
+                // Watch out for "edge tiles" that are smaller than the full
+                // tile size.
+                // For example, if the image is 100x100 and the tile size is
+                // 32x32, the last tile in each row will be 4x32, the bottom
+                // row of tiles will be 32x4, and the very last
+                // tile of the whole images will be 4x4.
             }
         }
     }
@@ -106,6 +118,8 @@ tiles_read()
     // END-imageinput-tiles
 }
 
+
+
 void
 unassociatedalpha()
 {
diff --git a/testsuite/docs-examples-cpp/src/docs-examples-imageoutput.cpp b/testsuite/docs-examples-cpp/src/docs-examples-imageoutput.cpp
index 825954e67b..d771cb92c4 100644
--- a/testsuite/docs-examples-cpp/src/docs-examples-imageoutput.cpp
+++ b/testsuite/docs-examples-cpp/src/docs-examples-imageoutput.cpp
@@ -65,9 +65,9 @@ scanlines_write()
     std::unique_ptr<ImageOutput> out = ImageOutput::create(filename);
     if (!out)
         return;  // error
-    ImageSpec spec(xres, yres, channels, TypeDesc::UINT8);
 
     // BEGIN-imageoutput-scanlines
+    ImageSpec spec(xres, yres, channels, TypeDesc::UINT8);
     std::vector<unsigned char> scanline(xres * channels);
     out->open(filename, spec);
     for (int y = 0; y < yres; ++y) {
@@ -80,10 +80,56 @@ scanlines_write()
 
 
 
+void
+tiles_write()
+{
+    const char* filename = "tiles.tif";
+    const int xres = 320, yres = 240, channels = 3;
+    const int tilesize = 64;
+
+    // BEGIN-imageoutput-tiles-create
+    std::unique_ptr<ImageOutput> out = ImageOutput::create(filename);
+    if (!out)
+        return;  // error: could not create output at all
+    if (!out->supports("tiles")) {
+        // Tiles are not supported
+    }
+    // END-imageoutput-tiles-create
+
+    // BEGIN-imageoutput-tiles-make-spec-open
+    ImageSpec spec(xres, yres, channels, TypeDesc::UINT8);
+    spec.tile_width  = tilesize;
+    spec.tile_height = tilesize;
+    out->open(filename, spec);
+    // END-imageoutput-tiles-make-spec-open
+
+    // BEGIN-imageoutput-tiles
+    std::vector<uint8_t> tile(tilesize * tilesize * spec.nchannels);
+    for (int y = 0; y < yres; y += tilesize) {
+        for (int x = 0; x < xres; x += tilesize) {
+            out->write_tiles(x, std::min(x + spec.tile_width, spec.width), y,
+                             std::min(y + spec.tile_height, spec.height), 0, 1,
+                             make_span(tile));
+            // ... process the pixels in tile[] ..
+            // Watch out for "edge tiles" that are smaller than the full
+            // tile size.
+            // For example, if the image is 100x100 and the tile size is
+            // 32x32, the last tile in each row will be 4x32, the bottom
+            // row of tiles will be 32x4, and the very last
+            // tile of the whole images will be 4x4.
+        }
+    }
+    out->close();
+    // END-imageoutput-tiles
+}
+
+
+
 int
 main(int /*argc*/, char** /*argv*/)
 {
     simple_write();
     scanlines_write();
+    tiles_write();
     return 0;
 }
diff --git a/testsuite/docs-examples-python/ref/out-arm.txt b/testsuite/docs-examples-python/ref/out-arm.txt
index cfd8a4b071..90a3f94f2b 100644
--- a/testsuite/docs-examples-python/ref/out-arm.txt
+++ b/testsuite/docs-examples-python/ref/out-arm.txt
@@ -144,3 +144,5 @@ Comparing "simple.tif" and "../docs-examples-cpp/ref/simple.tif"
 PASS
 Comparing "scanlines.tif" and "../docs-examples-cpp/ref/scanlines.tif"
 PASS
+Comparing "tiles.tif" and "ref/tiles.tif"
+PASS
diff --git a/testsuite/docs-examples-python/ref/out.txt b/testsuite/docs-examples-python/ref/out.txt
index 502bbda4cd..6c6c628cb6 100644
--- a/testsuite/docs-examples-python/ref/out.txt
+++ b/testsuite/docs-examples-python/ref/out.txt
@@ -144,3 +144,5 @@ Comparing "simple.tif" and "../docs-examples-cpp/ref/simple.tif"
 PASS
 Comparing "scanlines.tif" and "../docs-examples-cpp/ref/scanlines.tif"
 PASS
+Comparing "tiles.tif" and "ref/tiles.tif"
+PASS
diff --git a/testsuite/docs-examples-python/ref/tiles.tif b/testsuite/docs-examples-python/ref/tiles.tif
new file mode 100644
index 0000000000..2462e7f667
Binary files /dev/null and b/testsuite/docs-examples-python/ref/tiles.tif differ
diff --git a/testsuite/docs-examples-python/run.py b/testsuite/docs-examples-python/run.py
index 3bab2e8957..01993676a5 100755
--- a/testsuite/docs-examples-python/run.py
+++ b/testsuite/docs-examples-python/run.py
@@ -95,7 +95,7 @@
 # and need the images checked into the ref directory.
 outputs = [
     # Outputs from the ImageOutput chapter:
-    "simple.tif", "scanlines.tif",
+    "simple.tif", "scanlines.tif", "tiles.tif",
     # Outputs from the ImageInput chapter:
 
     # Outputs from the ImageBuf chapter:
diff --git a/testsuite/docs-examples-python/src/docs-examples-imageoutput.py b/testsuite/docs-examples-python/src/docs-examples-imageoutput.py
index 899eafc5cb..e76d8716cb 100644
--- a/testsuite/docs-examples-python/src/docs-examples-imageoutput.py
+++ b/testsuite/docs-examples-python/src/docs-examples-imageoutput.py
@@ -56,9 +56,6 @@ def simple_write() :
 # END-imageoutput-simple
 
 
-import OpenImageIO as oiio
-import numpy as np
-
 def scanlines_write() :
     filename = "scanlines.tif"
     xres = 320
@@ -80,9 +77,35 @@ def scanlines_write() :
         # END-imageoutput-scanlines
 
 
+def tiles_write() :
+    filename = "tiles.tif"
+    xres = 320
+    yres = 240
+    channels = 3  # RGB
+    tilesize = 64
+    spec = oiio.ImageSpec(xres, yres, channels, 'uint8')
+    spec.tile_width = tilesize
+    spec.tile_height = tilesize
+
+    out = oiio.ImageOutput.create (filename)
+    if out:
+        # BEGIN-imageoutput-tiles
+        z = 0   # Always zero for 2D images
+        out.open (filename, spec)
+        for y in range(0, yres, tilesize) :
+            for x in range(0, xres, tilesize) :
+                # Generate pixel array for one tile.
+                # As an example, we are just making a zero-filled tile
+                tile = np.zeros((tilesize, tilesize, channels), dtype=np.uint8)
+                out.write_tile (x, y, z, tile)
+        out.close ()
+        # END-imageoutput-tiles
+
+
 
 if __name__ == '__main__':
     # Each example function needs to get called here, or it won't execute
     # as part of the test.
     simple_write()
     scanlines_write()
+    tiles_write()