From 624533986beb123b59087c2eae90d2569965ae78 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 23 Dec 2024 02:35:11 +0100
Subject: [PATCH 1/7] GTiff/COG: add support for INTERLEAVE=TILE; COG: add
 support for INTERLEAVE=BAND

Fixes #10859

Up to now the COG driver only produced INTERLEAVE=PIXEL /
PlanarConfiguration=Contiguous files where values for all bands are in the
same tile. This PR add supports for INTERLEAVE=BAND where the tiles of
band 1 are placed first, followed by the ones of band 2, etc. It also
introduces a INTERLEAVE=TILE mode, which is similar to the BIL
(Band Interleave by Line), but generalize to tiles. That is you put first
tile (x,y)=(0,0) of band 1, then tile (0, 0) of band 2, ... tile (0, 0) of
band N, tile (1, 0) of band 1, ... tile (1, 0) of band N, etc. Both modes
can be useful for hyper-spectral datasets for example.
---
 autotest/gcore/cog.py                         | 298 ++++++++++
 autotest/gcore/test_driver_metadata.py        |   1 +
 autotest/gcore/tiff_write.py                  |  52 ++
 doc/source/drivers/raster/cog.rst             | 126 ++++-
 doc/source/drivers/raster/gtiff.rst           |  94 +++-
 frmts/gtiff/cogdriver.cpp                     |  31 +-
 frmts/gtiff/geotiff.cpp                       |   5 +-
 frmts/gtiff/gtiffdataset.cpp                  |  23 +-
 frmts/gtiff/gtiffdataset.h                    |   8 +-
 frmts/gtiff/gtiffdataset_read.cpp             | 525 +++++++++++++++++-
 frmts/gtiff/gtiffdataset_write.cpp            | 508 ++++++++++-------
 frmts/gtiff/gtiffrasterband.cpp               |  19 +-
 frmts/gtiff/gtiffrasterband.h                 |   4 -
 frmts/gtiff/gtiffrasterband_read.cpp          | 482 ----------------
 port/cpl_known_config_options.h               |   2 +-
 .../validate_cloud_optimized_geotiff.py       | 139 ++++-
 16 files changed, 1586 insertions(+), 731 deletions(-)

diff --git a/autotest/gcore/cog.py b/autotest/gcore/cog.py
index 8056e244caed..4b6e4b6a6aba 100755
--- a/autotest/gcore/cog.py
+++ b/autotest/gcore/cog.py
@@ -18,6 +18,7 @@
 
 import gdaltest
 import pytest
+import webserver
 from test_py_scripts import samples_path
 
 from osgeo import gdal, osr
@@ -1990,3 +1991,300 @@ def test_cog_preserve_ALPHA_PREMULTIPLIED_on_copy(tmp_vsimem):
             ds.GetRasterBand(4).GetMetadataItem("ALPHA", "IMAGE_STRUCTURE")
             == "PREMULTIPLIED"
         )
+
+
+###############################################################################
+#
+
+
+@gdaltest.enable_exceptions()
+def test_cog_write_interleave_tile(tmp_vsimem):
+    out_filename = str(tmp_vsimem / "out.tif")
+
+    with gdal.quiet_errors():
+        gdal.GetDriverByName("COG").CreateCopy(
+            out_filename,
+            gdal.Open("data/rgbsmall.tif"),
+            options=["INTERLEAVE=TILE", "BLOCKSIZE=32"],
+        )
+
+    ds = gdal.Open(out_filename)
+    assert ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE") == "TILE"
+    assert ds.GetMetadataItem("LAYOUT", "IMAGE_STRUCTURE") == "COG"
+
+    assert [ds.GetRasterBand(band + 1).Checksum() for band in range(3)] == [
+        21212,
+        21053,
+        21349,
+    ]
+
+    _check_cog(out_filename)
+
+
+###############################################################################
+#
+
+
+@gdaltest.enable_exceptions()
+def test_cog_write_interleave_band(tmp_vsimem):
+    out_filename = str(tmp_vsimem / "out.tif")
+
+    with gdal.quiet_errors():
+        gdal.GetDriverByName("COG").CreateCopy(
+            out_filename,
+            gdal.Open("data/rgbsmall.tif"),
+            options=["INTERLEAVE=BAND", "BLOCKSIZE=32"],
+        )
+
+    ds = gdal.Open(out_filename)
+    assert ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE") == "BAND"
+    assert ds.GetMetadataItem("LAYOUT", "IMAGE_STRUCTURE") == "COG"
+
+    assert [ds.GetRasterBand(band + 1).Checksum() for band in range(3)] == [
+        21212,
+        21053,
+        21349,
+    ]
+
+    _check_cog(out_filename)
+
+
+###############################################################################
+#
+
+
+@gdaltest.enable_exceptions()
+def test_cog_write_interleave_tile_with_mask(tmp_vsimem):
+    out_filename = str(tmp_vsimem / "out.tif")
+
+    with gdal.quiet_errors():
+        gdal.GetDriverByName("COG").CreateCopy(
+            out_filename,
+            gdal.Translate(
+                "", "data/stefan_full_rgba.tif", options="-f MEM -b 1 -b 2 -b 3 -mask 4"
+            ),
+            options=["INTERLEAVE=TILE", "BLOCKSIZE=32"],
+        )
+
+    ds = gdal.Open(out_filename)
+    assert ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE") == "TILE"
+    assert ds.GetMetadataItem("LAYOUT", "IMAGE_STRUCTURE") == "COG"
+
+    assert [ds.GetRasterBand(band + 1).Checksum() for band in range(3)] == [
+        12603,
+        58561,
+        36064,
+    ]
+    assert ds.GetRasterBand(1).GetMaskBand().Checksum() == 22499
+
+    _check_cog(out_filename)
+
+    # Check that the tiles are in the expected order in the file
+    last_offset = 0
+    for y in range(2):
+        for x in range(2):
+            for band in range(3):
+                offset = int(
+                    ds.GetRasterBand(band + 1).GetMetadataItem(
+                        f"BLOCK_OFFSET_{x}_{y}", "TIFF"
+                    )
+                )
+                assert offset > last_offset
+                last_offset = offset
+            offset = int(
+                ds.GetRasterBand(1)
+                .GetMaskBand()
+                .GetMetadataItem(f"BLOCK_OFFSET_{x}_{y}", "TIFF")
+            )
+            assert offset > last_offset
+            last_offset = offset
+
+
+###############################################################################
+#
+
+
+@gdaltest.enable_exceptions()
+def test_cog_write_interleave_tile_with_mask_and_ovr(tmp_vsimem):
+    out_filename = str(tmp_vsimem / "out.tif")
+    out2_filename = str(tmp_vsimem / "out2.tif")
+
+    ds = gdal.Translate(
+        out_filename,
+        "data/stefan_full_rgba.tif",
+        options="-b 1 -b 2 -b 3 -mask 4 -outsize 1024 0",
+    )
+    ds.BuildOverviews("NEAR", [2])
+    ds.Close()
+
+    ds = gdal.Open(out_filename)
+    expected_md = [ds.GetRasterBand(band + 1).Checksum() for band in range(3)]
+    expected_md += [ds.GetRasterBand(1).GetMaskBand().Checksum()]
+    expected_ovr_md = [
+        ds.GetRasterBand(band + 1).GetOverview(0).Checksum() for band in range(3)
+    ]
+    expected_ovr_md += [ds.GetRasterBand(1).GetMaskBand().GetOverview(0).Checksum()]
+
+    gdal.GetDriverByName("COG").CreateCopy(
+        out2_filename,
+        ds,
+        options=["INTERLEAVE=TILE", "OVERVIEW_RESAMPLING=NEAREST"],
+    )
+
+    ds = gdal.Open(out2_filename)
+    assert ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE") == "TILE"
+    assert ds.GetMetadataItem("LAYOUT", "IMAGE_STRUCTURE") == "COG"
+
+    _check_cog(out2_filename)
+
+    got_md = [ds.GetRasterBand(band + 1).Checksum() for band in range(3)]
+    got_md += [ds.GetRasterBand(1).GetMaskBand().Checksum()]
+    assert got_md == expected_md
+    got_ovr_md = [
+        ds.GetRasterBand(band + 1).GetOverview(0).Checksum() for band in range(3)
+    ]
+    got_ovr_md += [ds.GetRasterBand(1).GetMaskBand().GetOverview(0).Checksum()]
+    assert got_ovr_md == expected_ovr_md
+
+
+###############################################################################
+# Check that our reading of a COG with /vsicurl is efficient
+
+
+@pytest.mark.require_curl()
+@pytest.mark.skipif(
+    not check_libtiff_internal_or_at_least(4, 0, 11),
+    reason="libtiff >= 4.0.11 required",
+)
+@pytest.mark.parametrize("INTERLEAVE", ["BAND", "TILE"])
+def test_cog_interleave_tile_or_band_vsicurl(tmp_vsimem, INTERLEAVE):
+
+    gdal.VSICurlClearCache()
+
+    webserver_process = None
+    webserver_port = 0
+
+    (webserver_process, webserver_port) = webserver.launch(
+        handler=webserver.DispatcherHttpHandler
+    )
+    if webserver_port == 0:
+        pytest.skip()
+
+    in_filename = str(tmp_vsimem / "in.tif")
+    cog_filename = str(tmp_vsimem / "cog.tif")
+
+    ds = gdal.Translate(
+        in_filename,
+        "data/stefan_full_rgba.tif",
+        options="-b 1 -b 2 -b 3 -mask 4 -outsize 1024 0",
+    )
+    ds.BuildOverviews("NEAR", [2])
+    ds.Close()
+
+    src_ds = gdal.Open(in_filename)
+    gdal.GetDriverByName("COG").CreateCopy(
+        cog_filename,
+        src_ds,
+        options=["INTERLEAVE=" + INTERLEAVE, "OVERVIEW_RESAMPLING=NEAREST"],
+    )
+
+    def extract(offset, size):
+        f = gdal.VSIFOpenL(cog_filename, "rb")
+        gdal.VSIFSeekL(f, offset, 0)
+        data = gdal.VSIFReadL(size, 1, f)
+        gdal.VSIFCloseL(f)
+        return data
+
+    try:
+        filesize = gdal.VSIStatL(cog_filename).size
+
+        handler = webserver.SequentialHandler()
+        handler.add("HEAD", "/cog.tif", 200, {"Content-Length": "%d" % filesize})
+        handler.add(
+            "GET",
+            "/cog.tif",
+            206,
+            {"Content-Length": "16384"},
+            extract(0, 16384),
+            expected_headers={"Range": "bytes=0-16383"},
+        )
+        with webserver.install_http_handler(handler):
+            ds = gdal.Open("/vsicurl/http://localhost:%d/cog.tif" % webserver_port)
+
+        handler = webserver.SequentialHandler()
+
+        def method(request):
+            # sys.stderr.write('%s\n' % str(request.headers))
+
+            if request.headers["Range"].startswith("bytes="):
+                rng = request.headers["Range"][len("bytes=") :]
+                assert len(rng.split("-")) == 2
+                start = int(rng.split("-")[0])
+                end = int(rng.split("-")[1])
+
+                request.protocol_version = "HTTP/1.1"
+                request.send_response(206)
+                request.send_header("Content-type", "application/octet-stream")
+                request.send_header(
+                    "Content-Range", "bytes %d-%d/%d" % (start, end, filesize)
+                )
+                request.send_header("Content-Length", end - start + 1)
+                request.send_header("Connection", "close")
+                request.end_headers()
+
+                request.wfile.write(extract(start, end - start + 1))
+
+        handler.add("GET", "/cog.tif", custom_method=method)
+        with webserver.install_http_handler(handler):
+            ret = ds.ReadRaster()
+        assert ret == src_ds.ReadRaster()
+
+    finally:
+        webserver.server_stop(webserver_process, webserver_port)
+
+        gdal.VSICurlClearCache()
+
+
+###############################################################################
+
+
+@pytest.mark.require_creation_option("COG", "JPEG")
+@gdaltest.enable_exceptions()
+def test_cog_write_interleave_tile_jpeg(tmp_vsimem):
+    out_filename = str(tmp_vsimem / "out.tif")
+
+    gdal.GetDriverByName("GTiff").CreateCopy(
+        out_filename,
+        gdal.Open("data/rgbsmall.tif"),
+        options=["INTERLEAVE=BAND", "COMPRESS=JPEG"],
+    )
+    with gdal.Open(out_filename) as ds:
+        expected_md = [ds.GetRasterBand(band + 1).Checksum() for band in range(3)]
+
+    gdal.GetDriverByName("COG").CreateCopy(
+        out_filename,
+        gdal.Open("data/rgbsmall.tif"),
+        options=["INTERLEAVE=TILE", "COMPRESS=JPEG"],
+    )
+    with gdal.Open(out_filename) as ds:
+        got_md = [ds.GetRasterBand(band + 1).Checksum() for band in range(3)]
+
+    assert got_md == expected_md
+
+
+###############################################################################
+
+
+@pytest.mark.require_creation_option("COG", "WEBP")
+@gdaltest.enable_exceptions()
+def test_cog_write_interleave_tile_webp_error(tmp_vsimem):
+    out_filename = str(tmp_vsimem / "out.tif")
+
+    with pytest.raises(
+        Exception, match="COMPRESS=WEBP only supported for INTERLEAVE=PIXEL"
+    ):
+        gdal.GetDriverByName("COG").CreateCopy(
+            out_filename,
+            gdal.Open("data/rgbsmall.tif"),
+            options=["INTERLEAVE=TILE", "COMPRESS=WEBP"],
+        )
diff --git a/autotest/gcore/test_driver_metadata.py b/autotest/gcore/test_driver_metadata.py
index d36628be4055..b8907c47e938 100644
--- a/autotest/gcore/test_driver_metadata.py
+++ b/autotest/gcore/test_driver_metadata.py
@@ -130,6 +130,7 @@
         <xs:extension base="xs:string">
           <xs:attribute type="xs:string" name="alias" use="optional"/>
           <xs:attribute type="xs:string" name="aliasOf" use="optional"/>
+          <xs:attribute type="xs:string" name="remark" use="optional"/>
         </xs:extension>
       </xs:simpleContent>
     </xs:complexType>
diff --git a/autotest/gcore/tiff_write.py b/autotest/gcore/tiff_write.py
index 25f58a827f9c..c589fdf6b200 100755
--- a/autotest/gcore/tiff_write.py
+++ b/autotest/gcore/tiff_write.py
@@ -11988,3 +11988,55 @@ def test_tiff_write_warn_ignore_predictor_option(tmp_vsimem):
             out_filename, 1, 1, options=["PREDICTOR=2"]
         )
     assert "PREDICTOR option is ignored" in gdal.GetLastErrorMsg()
+
+
+###############################################################################
+#
+
+
+@gdaltest.enable_exceptions()
+@pytest.mark.parametrize("COPY_SRC_OVERVIEWS", ["YES", "NO"])
+def test_tiff_write_interleave_tile(tmp_vsimem, COPY_SRC_OVERVIEWS):
+    out_filename = str(tmp_vsimem / "out.tif")
+    with pytest.raises(
+        Exception, match="INTERLEAVE=TILE is only supported for CreateCopy"
+    ):
+        gdal.GetDriverByName("GTiff").Create(
+            out_filename, 1, 1, 2, options=["INTERLEAVE=TILE"]
+        )
+
+    ds = gdal.GetDriverByName("GTiff").CreateCopy(
+        out_filename,
+        gdal.Open("data/rgbsmall.tif"),
+        options=[
+            "INTERLEAVE=TILE",
+            "TILED=YES",
+            "BLOCKXSIZE=32",
+            "BLOCKYSIZE=32",
+            "COPY_SRC_OVERVIEWS=" + COPY_SRC_OVERVIEWS,
+        ],
+    )
+    assert ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE") == "TILE"
+    ds.Close()
+
+    ds = gdal.Open(out_filename)
+    assert ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE") == "TILE"
+
+    assert [ds.GetRasterBand(band + 1).Checksum() for band in range(3)] == [
+        21212,
+        21053,
+        21349,
+    ]
+
+    # Check that the tiles are in the expected order in the file
+    last_offset = 0
+    for y in range(2):
+        for x in range(2):
+            for band in range(3):
+                offset = int(
+                    ds.GetRasterBand(band + 1).GetMetadataItem(
+                        f"BLOCK_OFFSET_{x}_{y}", "TIFF"
+                    )
+                )
+                assert offset > last_offset
+                last_offset = offset
diff --git a/doc/source/drivers/raster/cog.rst b/doc/source/drivers/raster/cog.rst
index c46bac1f5703..a84a592725d6 100644
--- a/doc/source/drivers/raster/cog.rst
+++ b/doc/source/drivers/raster/cog.rst
@@ -43,6 +43,105 @@ General creation options
 
       Sets the tile width and height in pixels. Must be divisible by 16.
 
+-  .. co:: INTERLEAVE
+      :choices: BAND, PIXEL, TILE
+      :since: 3.11
+
+      Set the interleaving to use
+
+      * ``PIXEL``: for each spatial block, one TIFF tile/strip gathering values for
+        all bands is used . This matches the ``contiguous`` planar configuration in
+        TIFF terminology.
+        This is also known as a ``BIP (Band Interleaved Per Pixel)`` organization.
+        Such method is the default, and may be slightly less
+        efficient than BAND interleaving for some purposes, but some applications
+        only support pixel interleaved TIFF files. On the other hand, image-based
+        compression methods may perform better using PIXEL interleaving. For JPEG
+        PHOTOMETRIC=YCbCr, pixel interleaving is required. It is also required for
+        WebP compression.
+
+        Prior to GDAL 3.11, INTERLEAVE=PIXEL was the only possible interleaving
+        output by the COG driver.
+
+        Assuming a pixel[band][y][x] indexed array, when using CreateCopy(),
+        the pseudo code for the file disposition is:
+
+        ::
+
+          for y in 0 ... numberOfBlocksInHeight - 1:
+              for x in 0 ... numberOfTilesInWidth - 1:
+                  for j in 0 ... blockHeight - 1:
+                      for i in 0 ... blockWidth -1:
+                          start_new_strip_or_tile()
+                          for band in 0 ... numberBands -1:
+                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
+                          end_new_strip_or_tile()
+                      end_for
+                  end_for
+              end_for
+          end_for
+
+
+      * ``BAND``: for each spatial block, one TIFF tile/strip is used for each band.
+        This matches the contiguous ``separate`` configuration in TIFF terminology.
+        This is also known as a ``BSQ (Band SeQuential)`` organization.
+
+        In addition to that, when using CreateCopy(), data for the first band is
+        written first, followed by data for the second band, etc.
+        The pseudo code for the file disposition is:
+
+        ::
+
+          for y in 0 ... numberOfBlocksInHeight - 1:
+              for x in 0 ... numberOfTilesInWidth - 1:
+                  start_new_strip_or_tile()
+                  for band in 0 ... numberBands -1:
+                      for j in 0 ... blockHeight - 1:
+                          for i in 0 ... blockWidth -1:
+                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
+                          end_for
+                      end_for
+                  end_new_strip_or_tile()
+              end_for
+          end_for
+
+
+      * ``TILE`` (added in 3.11): this is a sort of
+        compromise between PIXEL and BAND, using the ``separate`` configuration,
+        but where data for a same spatial block is written for all bands, before
+        the data of the next spatial block.
+        When the block height is 1, this is also known as a
+        ``BIL (Band Interleaved per Line)`` organization.
+
+        Such a layout may be useful for writing hyperspectral datasets (several
+        hundred of bands), to get a compromise between efficient spatial query,
+        and partial band selection.
+
+        Assuming a pixel[band][y][x] indexed array, when using CreateCopy(),
+        the pseudo code for the file disposition is:
+
+        ::
+
+          for y in 0 ... numberOfBlocksInHeight - 1:
+              for x in 0 ... numberOfTilesInWidth - 1:
+                  for band in 0 ... numberBands -1:
+                      start_new_strip_or_tile()
+                      for j in 0 ... blockHeight - 1:
+                          for i in 0 ... blockWidth -1:
+                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
+                          end_for
+                      end_for
+                      end_new_strip_or_tile()
+                  end_for
+              end_for
+          end_for
+
+
+      Starting with GDAL 3.5, when copying from a source dataset with multiple bands
+      which advertises a INTERLEAVE metadata item, if the INTERLEAVE creation option
+      is not specified, the source dataset INTERLEAVE will be automatically taken
+      into account, unless the COMPRESS creation option is specified.
+
 -  .. co:: COMPRESS
       :choices: NONE, LZW, JPEG, DEFLATE, ZSTD, WEBP, LERC, LERC_DEFLATE, LERC_ZSTD, LZMA
       :default: LZW
@@ -56,7 +155,7 @@ General creation options
         files (seen as UInt16 bands with NBITS=12).
         For the COG driver, JPEG compression for 3 or 4-band images automatically
         selects the PHOTOMETRIC=YCBCR colorspace with a 4:2:2 subsampling of the Y,Cb,Cr
-        components.
+        components with the default INTERLEAVE=PIXEL.
         For a input dataset (single-band or 3-band), plus an alpha band,
         the alpha band will be converted as a 1-bit DEFLATE compressed mask.
 
@@ -65,6 +164,10 @@ General creation options
       * ``ZSTD`` is available when using internal libtiff and if GDAL built against
         libzstd >=1.0, or if built against external libtiff with zstd support.
 
+      * ``WEBP`` is available when using internal libtiff and if GDAL built against
+        libwebp, or if built against external libtiff with WebP support.
+        It can only be used with the default INTERLEAVE=PIXEL.
+
       * ``LERC`` is available when using internal libtiff.
 
       * ``LERC_ZSTD`` is available when ``LERC`` and ``ZSTD`` are available.
@@ -559,11 +662,18 @@ line).
      mask.TileByteCounts[i], but none of them actually need to be read)
    * trailer of mask data (4 bytes)
 
+   This is only written if INTERLEAVE=PIXEL.
+
+- ``INTERLEAVE=BAND`` or ``INTERLEAVE=TILE``: (GDAL >= 3.11)
+  Reflects the value of the INTERLEAVE creation option.
+  Omission implies INTERLEAVE=PIXEL.
+
+
 .. note::
 
     The content of the header ghost area can be retrieved by getting the
     ``GDAL_STRUCTURAL_METADATA`` metadata item of the ``TIFF`` metadata domain
-    on the datasett object (with GetMetadataItem())
+    on the dataset object (with GetMetadataItem())
 
 .. _cog.tile_data_leader_trailer:
 
@@ -576,7 +686,8 @@ that follows it. This leader is *ghost* in the sense that the
 TileOffsets[] array does not point to it, but points to the real payload. Hence
 the offset of the leader is TileOffsets[i]-4.
 
-An optimized reader seeing the ``BLOCK_LEADER=SIZE_AS_UINT4`` metadata item will thus look for TileOffset[i]
+For INTERLEAVE=PIXEL or INTERLEAVE=TILE, an optimized reader seeing the
+``BLOCK_LEADER=SIZE_AS_UINT4`` metadata item will thus look for TileOffset[i]
 and TileOffset[i+1] to deduce it must fetch the data starting at
 offset=TileOffset[i] - 4 and of size=TileOffset[i+1]-TileOffset[i]+4. It then
 checks the 4 first bytes to see if the size in this leader marker is
@@ -587,6 +698,9 @@ MASK_INTERLEAVED_WITH_IMAGERY=YES, then the tile size indicated in the leader
 will be < TileOffset[i+1]-TileOffset[i] since the data for the mask will
 follow the imagery data (see MASK_INTERLEAVED_WITH_IMAGERY=YES)
 
+For INTERLEAVE=BAND, the above paragraph applies but the successor of tile i
+is not tile i+1, but tile i+nTilesPerBand.
+
 Each tile data is immediately followed by a trailer, consisting of the repetition
 of the last 4 bytes of the payload of the tile data. The size of this trailer is
 *not* included in the TileByteCounts[] array. The purpose of this trailer is forces
@@ -613,3 +727,9 @@ See Also
 - If your source dataset is an internally tiled geotiff with the desired georeferencing and compression,
   using `cogger <https://github.com/airbusgeo/cogger>`__ (possibly along with gdaladdo to create overviews) will
   be much faster than the COG driver.
+
+
+.. below is an allow-list for spelling checker.
+
+.. spelling:word-list::
+      nTilesPerBand
diff --git a/doc/source/drivers/raster/gtiff.rst b/doc/source/drivers/raster/gtiff.rst
index 31f2d4a46ece..c10a625852cc 100644
--- a/doc/source/drivers/raster/gtiff.rst
+++ b/doc/source/drivers/raster/gtiff.rst
@@ -423,13 +423,95 @@ This driver supports the following creation options:
       RPC information is available.
 
 -  .. co:: INTERLEAVE
-      :choices: BAND, PIXEL
+      :choices: BAND, PIXEL, TILE
+
+      Set the interleaving to use
+
+      * ``PIXEL``: for each spatial block, one TIFF tile/strip gathering values for
+        all bands is used . This matches the ``contiguous`` planar configuration in
+        TIFF terminology.
+        This is also known as a ``BIP (Band Interleaved Per Pixel)`` organization.
+        Such method is the default, and may be slightly less
+        efficient than BAND interleaving for some purposes, but some applications
+        only support pixel interleaved TIFF files. On the other hand, image-based
+        compression methods may perform better using PIXEL interleaving. For JPEG
+        PHOTOMETRIC=YCbCr, pixel interleaving is required. It is also required for
+        WebP compression.
+
+        Assuming a pixel[band][y][x] indexed array, when using CreateCopy(),
+        the pseudo code for the file disposition is:
+
+        ::
+
+          for y in 0 ... numberOfBlocksInHeight - 1:
+              for x in 0 ... numberOfTilesInWidth - 1:
+                  for j in 0 ... blockHeight - 1:
+                      for i in 0 ... blockWidth -1:
+                          start_new_strip_or_tile()
+                          for band in 0 ... numberBands -1:
+                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
+                          end_new_strip_or_tile()
+                      end_for
+                  end_for
+              end_for
+          end_for
+
+
+      * ``BAND``: for each spatial block, one TIFF tile/strip is used for each band.
+        This matches the contiguous ``separate`` configuration in TIFF terminology.
+        This is also known as a ``BSQ (Band SeQuential)`` organization.
+
+        In addition to that, when using CreateCopy(), data for the first band is
+        written first, followed by data for the second band, etc.
+        The pseudo code for the file disposition is:
+
+        ::
+
+          for y in 0 ... numberOfBlocksInHeight - 1:
+              for x in 0 ... numberOfTilesInWidth - 1:
+                  start_new_strip_or_tile()
+                  for band in 0 ... numberBands -1:
+                      for j in 0 ... blockHeight - 1:
+                          for i in 0 ... blockWidth -1:
+                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
+                          end_for
+                      end_for
+                  end_new_strip_or_tile()
+              end_for
+          end_for
+
+
+      * ``TILE`` (added in 3.11, only for CreateCopy()): this is a sort of
+        compromise between PIXEL and BAND, using the ``separate`` configuration,
+        but where data for a same spatial block is written for all bands, before
+        the data of the next spatial block.
+        When the block height is 1, this is also known as a
+        ``BIL (Band Interleaved per Line)`` organization.
+
+        Such a layout may be useful for writing hyperspectral datasets (several
+        hundred of bands), to get a compromise between efficient spatial query,
+        and partial band selection.
+
+        Assuming a pixel[band][y][x] indexed array, when using CreateCopy(),
+        the pseudo code for the file disposition is:
+
+        ::
+
+          for y in 0 ... numberOfBlocksInHeight - 1:
+              for x in 0 ... numberOfTilesInWidth - 1:
+                  for band in 0 ... numberBands -1:
+                      start_new_strip_or_tile()
+                      for j in 0 ... blockHeight - 1:
+                          for i in 0 ... blockWidth -1:
+                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
+                          end_for
+                      end_for
+                      end_new_strip_or_tile()
+                  end_for
+              end_for
+          end_for
+
 
-      By default TIFF files with pixel
-      interleaving (PLANARCONFIG_CONTIG in TIFF terminology) are created.
-      These are slightly less efficient than BAND interleaving for some
-      purposes, but some applications only support pixel interleaved TIFF
-      files.
       Starting with GDAL 3.5, when copying from a source dataset with multiple bands
       which advertises a INTERLEAVE metadata item, if the INTERLEAVE creation option
       is not specified, the source dataset INTERLEAVE will be automatically taken
diff --git a/frmts/gtiff/cogdriver.cpp b/frmts/gtiff/cogdriver.cpp
index 310204a64127..73ced0419555 100644
--- a/frmts/gtiff/cogdriver.cpp
+++ b/frmts/gtiff/cogdriver.cpp
@@ -763,6 +763,21 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
         return nullptr;
     }
 
+    const CPLString osCompress = CSLFetchNameValueDef(
+        papszOptions, "COMPRESS", gbHasLZW ? "LZW" : "NONE");
+
+    const char *pszInterleave =
+        CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "PIXEL");
+    if (EQUAL(osCompress, "WEBP"))
+    {
+        if (!EQUAL(pszInterleave, "PIXEL"))
+        {
+            CPLError(CE_Failure, CPLE_NotSupported,
+                     "COMPRESS=WEBP only supported for INTERLEAVE=PIXEL");
+            return nullptr;
+        }
+    }
+
     CPLConfigOptionSetter oSetterReportDirtyBlockFlushing(
         "GDAL_REPORT_DIRTY_BLOCK_FLUSHING", "NO", true);
 
@@ -877,9 +892,7 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
         }
     }
 
-    CPLString osCompress = CSLFetchNameValueDef(papszOptions, "COMPRESS",
-                                                gbHasLZW ? "LZW" : "NONE");
-    if (EQUAL(osCompress, "JPEG") &&
+    if (EQUAL(osCompress, "JPEG") && EQUAL(pszInterleave, "PIXEL") &&
         (poCurDS->GetRasterCount() == 2 || poCurDS->GetRasterCount() == 4) &&
         poCurDS->GetRasterBand(poCurDS->GetRasterCount())
                 ->GetColorInterpretation() == GCI_AlphaBand)
@@ -1198,7 +1211,7 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
     if (EQUAL(osCompress, "JPEG"))
     {
         aosOptions.SetNameValue("JPEG_QUALITY", pszQuality);
-        if (nBands == 3)
+        if (nBands == 3 && EQUAL(pszInterleave, "PIXEL"))
             aosOptions.SetNameValue("PHOTOMETRIC", "YCBCR");
     }
     else if (EQUAL(osCompress, "WEBP"))
@@ -1316,7 +1329,8 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
     }
 
     std::unique_ptr<CPLConfigOptionSetter> poPhotometricSetter;
-    if (nBands == 3 && EQUAL(pszOverviewCompress, "JPEG"))
+    if (nBands == 3 && EQUAL(pszOverviewCompress, "JPEG") &&
+        EQUAL(pszInterleave, "PIXEL"))
     {
         poPhotometricSetter.reset(
             new CPLConfigOptionSetter("PHOTOMETRIC_OVERVIEW", "YCBCR", true));
@@ -1347,6 +1361,8 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
         aosOptions.AddNameValue("SRC_MDD", *papszSrcMDDIter);
     CSLDestroy(papszSrcMDD);
 
+    aosOptions.SetNameValue("INTERLEAVE", pszInterleave);
+
     CPLDebug("COG", "Generating final product: start");
     auto poRet =
         poGTiffDrv->CreateCopy(pszFilename, poCurDS, false, aosOptions.List(),
@@ -1525,6 +1541,11 @@ void GDALCOGDriver::InitializeCreationOptionList()
         "(16)'/>"
         "   <Option name='BLOCKSIZE' type='int' "
         "description='Tile size in pixels' min='128' default='512'/>"
+        "   <Option name='INTERLEAVE' type='string-select' default='PIXEL'>"
+        "       <Value>BAND</Value>"
+        "       <Value>PIXEL</Value>"
+        "       <Value>TILE</Value>"
+        "   </Option>"
         "   <Option name='BIGTIFF' type='string-select' description='"
         "Force creation of BigTIFF file'>"
         "     <Value>YES</Value>"
diff --git a/frmts/gtiff/geotiff.cpp b/frmts/gtiff/geotiff.cpp
index bd27e8bfe681..0e2f6d3b63f8 100644
--- a/frmts/gtiff/geotiff.cpp
+++ b/frmts/gtiff/geotiff.cpp
@@ -646,11 +646,13 @@ void GTiffWriteJPEGTables(TIFF *hTIFF, const char *pszPhotometric,
     papszLocalParameters =
         CSLSetNameValue(papszLocalParameters, "WRITE_JPEGTABLE_TAG", "NO");
 
+    bool bTileInterleaving;
     TIFF *hTIFFTmp =
         GTiffDataset::CreateLL(osTmpFilenameIn, nInMemImageWidth,
                                nInMemImageHeight, (nBands <= 4) ? nBands : 1,
                                (l_nBitsPerSample <= 8) ? GDT_Byte : GDT_UInt16,
-                               0.0, 0, papszLocalParameters, &fpTmp, osTmp);
+                               0.0, 0, papszLocalParameters, &fpTmp, osTmp,
+                               /* bCreateCopy=*/false, bTileInterleaving);
     CSLDestroy(papszLocalParameters);
     if (hTIFFTmp)
     {
@@ -1338,6 +1340,7 @@ void GDALRegister_GTiff()
         "   <Option name='INTERLEAVE' type='string-select' default='PIXEL'>"
         "       <Value>BAND</Value>"
         "       <Value>PIXEL</Value>"
+        "       <Value remark='Only for CreateCopy()'>TILE</Value>"
         "   </Option>"
         "   <Option name='TILED' type='boolean' description='Switch to tiled "
         "format'/>"
diff --git a/frmts/gtiff/gtiffdataset.cpp b/frmts/gtiff/gtiffdataset.cpp
index 26d1f917d28c..75fe9a07cad3 100644
--- a/frmts/gtiff/gtiffdataset.cpp
+++ b/frmts/gtiff/gtiffdataset.cpp
@@ -99,7 +99,7 @@ GTiffDataset::GTiffDataset()
       m_bLeaderSizeAsUInt4(false), m_bTrailerRepeatedLast4BytesRepeated(false),
       m_bMaskInterleavedWithImagery(false), m_bKnownIncompatibleEdition(false),
       m_bWriteKnownIncompatibleEdition(false), m_bHasUsedReadEncodedAPI(false),
-      m_bWriteCOGLayout(false)
+      m_bWriteCOGLayout(false), m_bTileInterleave(false)
 {
     // CPLDebug("GDAL", "sizeof(GTiffDataset) = %d bytes", static_cast<int>(
     //     sizeof(GTiffDataset)));
@@ -483,8 +483,8 @@ CPLErr GTiffDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
         const int nBlockY2 = (nYOff + nYSize - 1) / m_nBlockYSize;
         const int nXBlocks = nBlockX2 - nBlockX1 + 1;
         const int nYBlocks = nBlockY2 - nBlockY1 + 1;
-        const int nBlocks =
-            nXBlocks * nYBlocks *
+        const size_t nBlocks =
+            static_cast<size_t>(nXBlocks) * nYBlocks *
             (m_nPlanarConfig == PLANARCONFIG_CONTIG ? 1 : nBandCount);
         if (nBlocks > 1)
         {
@@ -497,13 +497,23 @@ CPLErr GTiffDataset::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
     const auto eDataType = poFirstBand->GetRasterDataType();
 
     if (eAccess == GA_ReadOnly && eRWFlag == GF_Read &&
-        (nBands == 1 || m_nPlanarConfig == PLANARCONFIG_CONTIG) &&
         HasOptimizedReadMultiRange() &&
         !(bCanUseMultiThreadedRead &&
           VSI_TIFFGetVSILFile(TIFFClientdata(m_hTIFF))->HasPRead()))
     {
-        pBufferedData = poFirstBand->CacheMultiRange(
-            nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize, psExtraArg);
+        if (nBands == 1 || m_nPlanarConfig == PLANARCONFIG_CONTIG)
+        {
+            const int nBandOne = 1;
+            pBufferedData =
+                CacheMultiRange(nXOff, nYOff, nXSize, nYSize, nBufXSize,
+                                nBufYSize, &nBandOne, 1, psExtraArg);
+        }
+        else
+        {
+            pBufferedData =
+                CacheMultiRange(nXOff, nYOff, nXSize, nYSize, nBufXSize,
+                                nBufYSize, panBandMap, nBandCount, psExtraArg);
+        }
     }
     else if (bCanUseMultiThreadedRead)
     {
@@ -1004,6 +1014,7 @@ void GTiffDataset::SetStructuralMDFromParent(GTiffDataset *poParentDS)
         poParentDS->m_bTrailerRepeatedLast4BytesRepeated;
     m_bMaskInterleavedWithImagery = poParentDS->m_bMaskInterleavedWithImagery;
     m_bWriteEmptyTiles = poParentDS->m_bWriteEmptyTiles;
+    m_bTileInterleave = poParentDS->m_bTileInterleave;
 }
 
 /************************************************************************/
diff --git a/frmts/gtiff/gtiffdataset.h b/frmts/gtiff/gtiffdataset.h
index 8d8a543279ec..e2babec77259 100644
--- a/frmts/gtiff/gtiffdataset.h
+++ b/frmts/gtiff/gtiffdataset.h
@@ -305,6 +305,7 @@ class GTiffDataset final : public GDALPamDataset
     bool m_bWriteKnownIncompatibleEdition : 1;
     bool m_bHasUsedReadEncodedAPI : 1;  // for debugging
     bool m_bWriteCOGLayout : 1;
+    bool m_bTileInterleave : 1;
 
     void ScanDirectories();
     bool ReadStrile(int nBlockId, void *pOutputBuffer,
@@ -439,6 +440,10 @@ class GTiffDataset final : public GDALPamDataset
 
     bool IsWholeBlock(int nXOff, int nYOff, int nXSize, int nYSize) const;
 
+    void *CacheMultiRange(int nXOff, int nYOff, int nXSize, int nYSize,
+                          int nBufXSize, int nBufYSize, const int *panBandMap,
+                          int nBandCount, GDALRasterIOExtraArg *psExtraArg);
+
     static void ThreadDecompressionFunc(void *pData);
 
     static GTIF *GTIFNew(TIFF *hTIFF);
@@ -541,7 +546,8 @@ class GTiffDataset final : public GDALPamDataset
                           int nBands, GDALDataType eType,
                           double dfExtraSpaceForOverviews,
                           int nColorTableMultiplier, char **papszParamList,
-                          VSILFILE **pfpL, CPLString &osTmpFilename);
+                          VSILFILE **pfpL, CPLString &osTmpFilename,
+                          bool bCreateCopy, bool &bTileInterleavingOut);
 
     CPLErr WriteEncodedTileOrStrip(uint32_t tile_or_strip, void *data,
                                    int bPreserveDataBuffer);
diff --git a/frmts/gtiff/gtiffdataset_read.cpp b/frmts/gtiff/gtiffdataset_read.cpp
index 1ba913fa4e23..2393629eac00 100644
--- a/frmts/gtiff/gtiffdataset_read.cpp
+++ b/frmts/gtiff/gtiffdataset_read.cpp
@@ -3908,13 +3908,20 @@ GDALDataset *GTiffDataset::Open(GDALOpenInfo *poOpenInfo)
         poDS->m_bBlockOrderRowMajor =
             strstr(pszStructuralMD, "BLOCK_ORDER=ROW_MAJOR") != nullptr;
         poDS->m_bLeaderSizeAsUInt4 =
-            strstr(pszStructuralMD, "BLOCK_LEADER=SIZE_AS_UINT4") != nullptr;
+            strstr(pszStructuralMD, "BLOCK_LEADER=SIZE_AS_UINT4") != nullptr &&
+            (strstr(pszStructuralMD, "INTERLEAVE=") == nullptr ||
+             strstr(pszStructuralMD, "INTERLEAVE=BAND") != nullptr ||
+             strstr(pszStructuralMD, "INTERLEAVE=TILE") != nullptr);
         poDS->m_bTrailerRepeatedLast4BytesRepeated =
             strstr(pszStructuralMD, "BLOCK_TRAILER=LAST_4_BYTES_REPEATED") !=
-            nullptr;
+                nullptr &&
+            (strstr(pszStructuralMD, "INTERLEAVE=") == nullptr ||
+             strstr(pszStructuralMD, "INTERLEAVE=BAND") != nullptr ||
+             strstr(pszStructuralMD, "INTERLEAVE=TILE") != nullptr);
         poDS->m_bMaskInterleavedWithImagery =
             strstr(pszStructuralMD, "MASK_INTERLEAVED_WITH_IMAGERY=YES") !=
-            nullptr;
+                nullptr &&
+            strstr(pszStructuralMD, "INTERLEAVE=") == nullptr;
         poDS->m_bKnownIncompatibleEdition =
             strstr(pszStructuralMD, "KNOWN_INCOMPATIBLE_EDITION=YES") !=
             nullptr;
@@ -5469,8 +5476,24 @@ CPLErr GTiffDataset::OpenOffset(TIFF *hTIFFIn, toff_t nDirOffsetIn,
                 continue;
             if (EQUAL(pszDomain, "IMAGE_STRUCTURE"))
             {
-                if (m_nCompression == COMPRESSION_WEBP &&
-                    EQUAL(pszKey, "COMPRESSION_REVERSIBILITY"))
+                if (EQUAL(pszKey, "INTERLEAVE"))
+                {
+                    if (EQUAL(pszValue, "TILE"))
+                    {
+                        m_bTileInterleave = true;
+                        m_oGTiffMDMD.SetMetadataItem("INTERLEAVE", "TILE",
+                                                     "IMAGE_STRUCTURE");
+                    }
+                    else
+                    {
+                        CPLDebug("GTiff",
+                                 "Unhandled INTERLEAVE=%s found in "
+                                 "GDAL_METADATA tag",
+                                 pszValue);
+                    }
+                }
+                else if (m_nCompression == COMPRESSION_WEBP &&
+                         EQUAL(pszKey, "COMPRESSION_REVERSIBILITY"))
                 {
                     if (EQUAL(pszValue, "LOSSLESS"))
                         m_bWebPLossless = true;
@@ -6644,3 +6667,495 @@ bool GTiffDataset::HasOptimizedReadMultiRange()
                "GTIFF_HAS_OPTIMIZED_READ_MULTI_RANGE", "NO")));
     return m_nHasOptimizedReadMultiRange != 0;
 }
+
+/************************************************************************/
+/*                         CacheMultiRange()                            */
+/************************************************************************/
+
+static bool CheckTrailer(const GByte *strileData, vsi_l_offset nStrileSize)
+{
+    GByte abyTrailer[4];
+    memcpy(abyTrailer, strileData + nStrileSize, 4);
+    GByte abyLastBytes[4] = {};
+    if (nStrileSize >= 4)
+        memcpy(abyLastBytes, strileData + nStrileSize - 4, 4);
+    else
+    {
+        // The last bytes will be zero due to the above {} initialization,
+        // and that's what should be in abyTrailer too when the trailer is
+        // correct.
+        memcpy(abyLastBytes, strileData, static_cast<size_t>(nStrileSize));
+    }
+    return memcmp(abyTrailer, abyLastBytes, 4) == 0;
+}
+
+void *GTiffDataset::CacheMultiRange(int nXOff, int nYOff, int nXSize,
+                                    int nYSize, int nBufXSize, int nBufYSize,
+                                    const int *panBandMap, int nBandCount,
+                                    GDALRasterIOExtraArg *psExtraArg)
+{
+    void *pBufferedData = nullptr;
+    // Same logic as in GDALRasterBand::IRasterIO()
+    double dfXOff = nXOff;
+    double dfYOff = nYOff;
+    double dfXSize = nXSize;
+    double dfYSize = nYSize;
+    if (psExtraArg->bFloatingPointWindowValidity)
+    {
+        dfXOff = psExtraArg->dfXOff;
+        dfYOff = psExtraArg->dfYOff;
+        dfXSize = psExtraArg->dfXSize;
+        dfYSize = psExtraArg->dfYSize;
+    }
+    const double dfSrcXInc = dfXSize / static_cast<double>(nBufXSize);
+    const double dfSrcYInc = dfYSize / static_cast<double>(nBufYSize);
+    const double EPS = 1e-10;
+    const int nBlockX1 =
+        static_cast<int>(std::max(0.0, (0 + 0.5) * dfSrcXInc + dfXOff + EPS)) /
+        m_nBlockXSize;
+    const int nBlockY1 =
+        static_cast<int>(std::max(0.0, (0 + 0.5) * dfSrcYInc + dfYOff + EPS)) /
+        m_nBlockYSize;
+    const int nBlockX2 =
+        static_cast<int>(
+            std::min(static_cast<double>(nRasterXSize - 1),
+                     (nBufXSize - 1 + 0.5) * dfSrcXInc + dfXOff + EPS)) /
+        m_nBlockXSize;
+    const int nBlockY2 =
+        static_cast<int>(
+            std::min(static_cast<double>(nRasterYSize - 1),
+                     (nBufYSize - 1 + 0.5) * dfSrcYInc + dfYOff + EPS)) /
+        m_nBlockYSize;
+
+    struct StrileData
+    {
+        vsi_l_offset nOffset;
+        vsi_l_offset nByteCount;
+        bool bTryMask;
+    };
+
+    std::map<int, StrileData> oMapStrileToOffsetByteCount;
+
+    // Dedicated method to retrieved the offset and size in an efficient way
+    // when m_bBlockOrderRowMajor and m_bLeaderSizeAsUInt4 conditions are
+    // met.
+    // Except for the last block, we just read the offset from the TIFF offset
+    // array, and retrieve the size in the leader 4 bytes that come before the
+    // payload.
+    auto OptimizedRetrievalOfOffsetSize =
+        [&](int nBand, int nBlockId, vsi_l_offset &nOffset, vsi_l_offset &nSize,
+            size_t nTotalSize, size_t nMaxRawBlockCacheSize)
+    {
+        bool bTryMask = m_bMaskInterleavedWithImagery;
+        nOffset = TIFFGetStrileOffset(m_hTIFF, nBlockId);
+        if (nOffset >= 4)
+        {
+            if ((m_nPlanarConfig == PLANARCONFIG_CONTIG &&
+                 nBlockId == m_nBlocksPerBand - 1) ||
+                (m_nPlanarConfig == PLANARCONFIG_SEPARATE &&
+                 nBlockId == m_nBlocksPerBand * nBands - 1))
+            {
+                // Special case for the last block. As there is no next block
+                // from which to retrieve an offset, use the good old method
+                // that consists in reading the ByteCount array.
+                if (m_nPlanarConfig == PLANARCONFIG_CONTIG && bTryMask &&
+                    GetRasterBand(1)->GetMaskBand() && m_poMaskDS)
+                {
+                    auto nMaskOffset =
+                        TIFFGetStrileOffset(m_poMaskDS->m_hTIFF, nBlockId);
+                    if (nMaskOffset)
+                    {
+                        nSize = nMaskOffset +
+                                TIFFGetStrileByteCount(m_poMaskDS->m_hTIFF,
+                                                       nBlockId) -
+                                nOffset;
+                    }
+                    else
+                    {
+                        bTryMask = false;
+                    }
+                }
+                if (nSize == 0)
+                {
+                    nSize = TIFFGetStrileByteCount(m_hTIFF, nBlockId);
+                }
+                if (nSize && m_bTrailerRepeatedLast4BytesRepeated)
+                {
+                    nSize += 4;
+                }
+            }
+            else
+            {
+                const auto nNextBlockId =
+                    (m_bTileInterleave && nBand < nBands)
+                        ? nBlockId + m_nBlocksPerBand
+                    : (m_bTileInterleave && nBand == nBands)
+                        ? nBlockId - (nBand - 1) * m_nBlocksPerBand + 1
+                        : nBlockId + 1;
+                auto nOffsetNext = TIFFGetStrileOffset(m_hTIFF, nNextBlockId);
+                if (nOffsetNext > nOffset)
+                {
+                    nSize = nOffsetNext - nOffset;
+                }
+                else
+                {
+                    // Shouldn't happen for a compliant file
+                    if (nOffsetNext != 0)
+                    {
+                        CPLDebug("GTiff", "Tile %d is not located after %d",
+                                 nNextBlockId, nBlockId);
+                    }
+                    bTryMask = false;
+                    nSize = TIFFGetStrileByteCount(m_hTIFF, nBlockId);
+                    if (m_bTrailerRepeatedLast4BytesRepeated)
+                        nSize += 4;
+                }
+            }
+            if (nSize)
+            {
+                nOffset -= 4;
+                nSize += 4;
+                if (nTotalSize + nSize < nMaxRawBlockCacheSize)
+                {
+                    StrileData data;
+                    data.nOffset = nOffset;
+                    data.nByteCount = nSize;
+                    data.bTryMask = bTryMask;
+                    oMapStrileToOffsetByteCount[nBlockId] = data;
+                }
+            }
+        }
+        else
+        {
+            // Sparse tile
+            StrileData data;
+            data.nOffset = 0;
+            data.nByteCount = 0;
+            data.bTryMask = false;
+            oMapStrileToOffsetByteCount[nBlockId] = data;
+        }
+    };
+
+    // This lambda fills m_poDS->m_oCacheStrileToOffsetByteCount (and
+    // m_poDS->m_poMaskDS->m_oCacheStrileToOffsetByteCount, when there is a
+    // mask) from the temporary oMapStrileToOffsetByteCount.
+    auto FillCacheStrileToOffsetByteCount =
+        [&](const std::vector<vsi_l_offset> &anOffsets,
+            const std::vector<size_t> &anSizes,
+            const std::vector<void *> &apData)
+    {
+        CPLAssert(m_bLeaderSizeAsUInt4);
+        size_t i = 0;
+        vsi_l_offset nLastOffset = 0;
+        for (const auto &entry : oMapStrileToOffsetByteCount)
+        {
+            const auto nBlockId = entry.first;
+            const auto nOffset = entry.second.nOffset;
+            const auto nSize = entry.second.nByteCount;
+            if (nOffset == 0)
+            {
+                // Sparse tile
+                m_oCacheStrileToOffsetByteCount.insert(nBlockId,
+                                                       std::pair(0, 0));
+                continue;
+            }
+
+            if (nOffset < nLastOffset)
+            {
+                // shouldn't happen normally if tiles are sorted
+                i = 0;
+            }
+            nLastOffset = nOffset;
+            while (i < anOffsets.size() &&
+                   !(nOffset >= anOffsets[i] &&
+                     nOffset + nSize <= anOffsets[i] + anSizes[i]))
+            {
+                i++;
+            }
+            CPLAssert(i < anOffsets.size());
+            CPLAssert(nOffset >= anOffsets[i]);
+            CPLAssert(nOffset + nSize <= anOffsets[i] + anSizes[i]);
+            GUInt32 nSizeFromLeader;
+            memcpy(&nSizeFromLeader,
+                   // cppcheck-suppress containerOutOfBounds
+                   static_cast<GByte *>(apData[i]) + nOffset - anOffsets[i],
+                   sizeof(nSizeFromLeader));
+            CPL_LSBPTR32(&nSizeFromLeader);
+            bool bOK = true;
+            constexpr int nLeaderSize = 4;
+            const int nTrailerSize =
+                (m_bTrailerRepeatedLast4BytesRepeated ? 4 : 0);
+            if (nSizeFromLeader > nSize - nLeaderSize - nTrailerSize)
+            {
+                CPLDebug("GTiff",
+                         "Inconsistent block size from in leader of block %d",
+                         nBlockId);
+                bOK = false;
+            }
+            else if (m_bTrailerRepeatedLast4BytesRepeated)
+            {
+                // Check trailer consistency
+                const GByte *strileData = static_cast<GByte *>(apData[i]) +
+                                          nOffset - anOffsets[i] + nLeaderSize;
+                if (!CheckTrailer(strileData, nSizeFromLeader))
+                {
+                    CPLDebug("GTiff", "Inconsistent trailer of block %d",
+                             nBlockId);
+                    bOK = false;
+                }
+            }
+            if (!bOK)
+            {
+                return false;
+            }
+
+            {
+                const vsi_l_offset nRealOffset = nOffset + nLeaderSize;
+                const vsi_l_offset nRealSize = nSizeFromLeader;
+#ifdef DEBUG_VERBOSE
+                CPLDebug("GTiff",
+                         "Block %d found at offset " CPL_FRMT_GUIB
+                         " with size " CPL_FRMT_GUIB,
+                         nBlockId, nRealOffset, nRealSize);
+#endif
+                m_oCacheStrileToOffsetByteCount.insert(
+                    nBlockId, std::pair(nRealOffset, nRealSize));
+            }
+
+            // Processing of mask
+            if (!(entry.second.bTryMask && m_bMaskInterleavedWithImagery &&
+                  GetRasterBand(1)->GetMaskBand() && m_poMaskDS))
+            {
+                continue;
+            }
+
+            bOK = false;
+            const vsi_l_offset nMaskOffsetWithLeader =
+                nOffset + nLeaderSize + nSizeFromLeader + nTrailerSize;
+            if (nMaskOffsetWithLeader + nLeaderSize <=
+                anOffsets[i] + anSizes[i])
+            {
+                GUInt32 nMaskSizeFromLeader;
+                memcpy(&nMaskSizeFromLeader,
+                       static_cast<GByte *>(apData[i]) + nMaskOffsetWithLeader -
+                           anOffsets[i],
+                       sizeof(nMaskSizeFromLeader));
+                CPL_LSBPTR32(&nMaskSizeFromLeader);
+                if (nMaskOffsetWithLeader + nLeaderSize + nMaskSizeFromLeader +
+                        nTrailerSize <=
+                    anOffsets[i] + anSizes[i])
+                {
+                    bOK = true;
+                    if (m_bTrailerRepeatedLast4BytesRepeated)
+                    {
+                        // Check trailer consistency
+                        const GByte *strileMaskData =
+                            static_cast<GByte *>(apData[i]) + nOffset -
+                            anOffsets[i] + nLeaderSize + nSizeFromLeader +
+                            nTrailerSize + nLeaderSize;
+                        if (!CheckTrailer(strileMaskData, nMaskSizeFromLeader))
+                        {
+                            CPLDebug("GTiff",
+                                     "Inconsistent trailer of mask of block %d",
+                                     nBlockId);
+                            bOK = false;
+                        }
+                    }
+                }
+                if (bOK)
+                {
+                    const vsi_l_offset nRealOffset = nOffset + nLeaderSize +
+                                                     nSizeFromLeader +
+                                                     nTrailerSize + nLeaderSize;
+                    const vsi_l_offset nRealSize = nMaskSizeFromLeader;
+#ifdef DEBUG_VERBOSE
+                    CPLDebug("GTiff",
+                             "Mask of block %d found at offset " CPL_FRMT_GUIB
+                             " with size " CPL_FRMT_GUIB,
+                             nBlockId, nRealOffset, nRealSize);
+#endif
+
+                    m_poMaskDS->m_oCacheStrileToOffsetByteCount.insert(
+                        nBlockId, std::pair(nRealOffset, nRealSize));
+                }
+            }
+            if (!bOK)
+            {
+                CPLDebug("GTiff",
+                         "Mask for block %d is not properly interleaved with "
+                         "imagery block",
+                         nBlockId);
+            }
+        }
+        return true;
+    };
+
+    thandle_t th = TIFFClientdata(m_hTIFF);
+    if (!VSI_TIFFHasCachedRanges(th))
+    {
+        std::vector<std::pair<vsi_l_offset, size_t>> aOffsetSize;
+        size_t nTotalSize = 0;
+        const unsigned int nMaxRawBlockCacheSize = atoi(
+            CPLGetConfigOption("GDAL_MAX_RAW_BLOCK_CACHE_SIZE", "10485760"));
+        bool bGoOn = true;
+        for (int iBand = 0; iBand < nBandCount; ++iBand)
+        {
+            const int nBand = panBandMap[iBand];
+            GTiffRasterBand *poBand =
+                cpl::down_cast<GTiffRasterBand *>(papoBands[nBand - 1]);
+            for (int iY = nBlockY1; bGoOn && iY <= nBlockY2; iY++)
+            {
+                for (int iX = nBlockX1; bGoOn && iX <= nBlockX2; iX++)
+                {
+                    GDALRasterBlock *poBlock =
+                        poBand->TryGetLockedBlockRef(iX, iY);
+                    if (poBlock != nullptr)
+                    {
+                        poBlock->DropLock();
+                        continue;
+                    }
+                    int nBlockId = iX + iY * m_nBlocksPerRow;
+                    if (m_nPlanarConfig == PLANARCONFIG_SEPARATE)
+                        nBlockId += (nBand - 1) * m_nBlocksPerBand;
+                    vsi_l_offset nOffset = 0;
+                    vsi_l_offset nSize = 0;
+
+                    if (!m_bStreamingIn && m_bBlockOrderRowMajor &&
+                        m_bLeaderSizeAsUInt4)
+                    {
+                        OptimizedRetrievalOfOffsetSize(nBand, nBlockId, nOffset,
+                                                       nSize, nTotalSize,
+                                                       nMaxRawBlockCacheSize);
+                    }
+                    else
+                    {
+                        CPL_IGNORE_RET_VAL(IsBlockAvailable(nBlockId, &nOffset,
+                                                            &nSize, nullptr));
+                    }
+                    if (nSize)
+                    {
+                        if (nTotalSize + nSize < nMaxRawBlockCacheSize)
+                        {
+#ifdef DEBUG_VERBOSE
+                            CPLDebug(
+                                "GTiff",
+                                "Precaching for block (%d, %d), " CPL_FRMT_GUIB
+                                "-" CPL_FRMT_GUIB,
+                                iX, iY, nOffset,
+                                nOffset + static_cast<size_t>(nSize) - 1);
+#endif
+                            aOffsetSize.push_back(
+                                std::pair(nOffset, static_cast<size_t>(nSize)));
+                            nTotalSize += static_cast<size_t>(nSize);
+                        }
+                        else
+                        {
+                            bGoOn = false;
+                        }
+                    }
+                }
+            }
+        }
+
+        std::sort(aOffsetSize.begin(), aOffsetSize.end());
+
+        if (nTotalSize > 0)
+        {
+            pBufferedData = VSI_MALLOC_VERBOSE(nTotalSize);
+            if (pBufferedData)
+            {
+                std::vector<vsi_l_offset> anOffsets;
+                std::vector<size_t> anSizes;
+                std::vector<void *> apData;
+                anOffsets.push_back(aOffsetSize[0].first);
+                apData.push_back(static_cast<GByte *>(pBufferedData));
+                size_t nChunkSize = aOffsetSize[0].second;
+                size_t nAccOffset = 0;
+                // Try to merge contiguous or slightly overlapping ranges
+                for (size_t i = 0; i < aOffsetSize.size() - 1; i++)
+                {
+                    if (aOffsetSize[i].first < aOffsetSize[i + 1].first &&
+                        aOffsetSize[i].first + aOffsetSize[i].second >=
+                            aOffsetSize[i + 1].first)
+                    {
+                        const auto overlap = aOffsetSize[i].first +
+                                             aOffsetSize[i].second -
+                                             aOffsetSize[i + 1].first;
+                        // That should always be the case for well behaved
+                        // TIFF files.
+                        if (aOffsetSize[i + 1].second > overlap)
+                        {
+                            nChunkSize += static_cast<size_t>(
+                                aOffsetSize[i + 1].second - overlap);
+                        }
+                    }
+                    else
+                    {
+                        // terminate current block
+                        anSizes.push_back(nChunkSize);
+#ifdef DEBUG_VERBOSE
+                        CPLDebug("GTiff",
+                                 "Requesting range [" CPL_FRMT_GUIB
+                                 "-" CPL_FRMT_GUIB "]",
+                                 anOffsets.back(),
+                                 anOffsets.back() + anSizes.back() - 1);
+#endif
+                        nAccOffset += nChunkSize;
+                        // start a new range
+                        anOffsets.push_back(aOffsetSize[i + 1].first);
+                        apData.push_back(static_cast<GByte *>(pBufferedData) +
+                                         nAccOffset);
+                        nChunkSize = aOffsetSize[i + 1].second;
+                    }
+                }
+                // terminate last block
+                anSizes.push_back(nChunkSize);
+#ifdef DEBUG_VERBOSE
+                CPLDebug(
+                    "GTiff",
+                    "Requesting range [" CPL_FRMT_GUIB "-" CPL_FRMT_GUIB "]",
+                    anOffsets.back(), anOffsets.back() + anSizes.back() - 1);
+#endif
+
+                VSILFILE *fp = VSI_TIFFGetVSILFile(th);
+
+                // An error in VSIFReadMultiRangeL() will not be critical,
+                // as this method is an optimization, and if it fails,
+                // tile-by-tile data acquisition will be done, so we can
+                // temporary turn failures into warnings.
+                CPLTurnFailureIntoWarning(true);
+                const bool ok =
+                    VSIFReadMultiRangeL(static_cast<int>(anSizes.size()),
+                                        &apData[0], &anOffsets[0], &anSizes[0],
+                                        fp) == 0;
+                CPLTurnFailureIntoWarning(false);
+                if (ok)
+                {
+                    if (!oMapStrileToOffsetByteCount.empty() &&
+                        !FillCacheStrileToOffsetByteCount(anOffsets, anSizes,
+                                                          apData))
+                    {
+                        // Retry without optimization
+                        CPLFree(pBufferedData);
+                        m_bLeaderSizeAsUInt4 = false;
+                        void *pRet = CacheMultiRange(
+                            nXOff, nYOff, nXSize, nYSize, nBufXSize, nBufYSize,
+                            panBandMap, nBandCount, psExtraArg);
+                        m_bLeaderSizeAsUInt4 = true;
+                        return pRet;
+                    }
+
+                    VSI_TIFFSetCachedRanges(
+                        th, static_cast<int>(anSizes.size()), &apData[0],
+                        &anOffsets[0], &anSizes[0]);
+                }
+                else
+                {
+                    CPLFree(pBufferedData);
+                    pBufferedData = nullptr;
+                }
+            }
+        }
+    }
+    return pBufferedData;
+}
diff --git a/frmts/gtiff/gtiffdataset_write.cpp b/frmts/gtiff/gtiffdataset_write.cpp
index aa45e94280fc..5d0dd63ed98f 100644
--- a/frmts/gtiff/gtiffdataset_write.cpp
+++ b/frmts/gtiff/gtiffdataset_write.cpp
@@ -598,7 +598,6 @@ void GTiffDataset::WriteDealWithLercAndNan(T *pBuffer, int nActualBlockWidth,
 bool GTiffDataset::WriteEncodedTile(uint32_t tile, GByte *pabyData,
                                     int bPreserveDataBuffer)
 {
-
     const int iColumn = (tile % m_nBlocksPerBand) % m_nBlocksPerRow;
     const int iRow = (tile % m_nBlocksPerBand) / m_nBlocksPerRow;
 
@@ -4285,6 +4284,14 @@ bool GTiffDataset::WriteMetadata(GDALDataset *poSrcDS, TIFF *l_hTIFF,
     if (CPLTestBool(
             CPLGetConfigOption("GTIFF_WRITE_IMAGE_STRUCTURE_METADATA", "YES")))
     {
+        const char *pszInterleave =
+            CSLFetchNameValue(papszCreationOptions, "INTERLEAVE");
+        if (pszInterleave && EQUAL(pszInterleave, "TILE"))
+        {
+            AppendMetadataItem(&psRoot, &psTail, "INTERLEAVE", "TILE", 0,
+                               nullptr, "IMAGE_STRUCTURE");
+        }
+
         const char *pszCompress =
             CSLFetchNameValue(papszCreationOptions, "COMPRESS");
         if (pszCompress && EQUAL(pszCompress, "WEBP"))
@@ -5038,9 +5045,12 @@ TIFF *GTiffDataset::CreateLL(const char *pszFilename, int nXSize, int nYSize,
                              int l_nBands, GDALDataType eType,
                              double dfExtraSpaceForOverviews,
                              int nColorTableMultiplier, char **papszParamList,
-                             VSILFILE **pfpL, CPLString &l_osTmpFilename)
+                             VSILFILE **pfpL, CPLString &l_osTmpFilename,
+                             bool bCreateCopy, bool &bTileInterleavingOut)
 
 {
+    bTileInterleavingOut = false;
+
     GTiffOneTimeInit();
 
     /* -------------------------------------------------------------------- */
@@ -5120,11 +5130,30 @@ TIFF *GTiffDataset::CreateLL(const char *pszFilename, int nXSize, int nYSize,
         {
             nPlanar = PLANARCONFIG_SEPARATE;
         }
+        else if (EQUAL(pszValue, "BAND"))
+        {
+            nPlanar = PLANARCONFIG_SEPARATE;
+        }
+        else if (EQUAL(pszValue, "TILE"))
+        {
+            bTileInterleavingOut = true;
+            nPlanar = PLANARCONFIG_SEPARATE;
+            if (!bCreateCopy)
+            {
+                CPLError(CE_Failure, CPLE_NotSupported,
+                         "INTERLEAVE=TILE is only supported for CreateCopy(), "
+                         "not Create()");
+                return nullptr;
+            }
+        }
         else
         {
             ReportError(
                 pszFilename, CE_Failure, CPLE_IllegalArg,
-                "INTERLEAVE=%s unsupported, value must be PIXEL or BAND.",
+                bCreateCopy
+                    ? "INTERLEAVE=%s unsupported, value must be PIXEL, BAND or "
+                      "TILE."
+                    : "INTERLEAVE=%s unsupported, value must be PIXEL or BAND.",
                 pszValue);
             return nullptr;
         }
@@ -6176,10 +6205,11 @@ int GTiffDataset::GuessJPEGQuality(bool &bOutHasQuantizationTable,
 
         CPLPushErrorHandler(CPLQuietErrorHandler);
         CPLString osTmp;
-        TIFF *hTIFFTmp =
-            CreateLL(osTmpFilenameIn, 16, 16, (nBands <= 4) ? nBands : 1,
-                     GetRasterBand(1)->GetRasterDataType(), 0.0, 0,
-                     papszLocalParameters, &fpTmp, osTmp);
+        bool bTileInterleaving;
+        TIFF *hTIFFTmp = CreateLL(
+            osTmpFilenameIn, 16, 16, (nBands <= 4) ? nBands : 1,
+            GetRasterBand(1)->GetRasterDataType(), 0.0, 0, papszLocalParameters,
+            &fpTmp, osTmp, /* bCreateCopy=*/false, bTileInterleaving);
         CPLPopErrorHandler();
         if (!hTIFFTmp)
         {
@@ -6338,9 +6368,11 @@ GDALDataset *GTiffDataset::Create(const char *pszFilename, int nXSize,
     /* -------------------------------------------------------------------- */
     /*      Create the underlying TIFF file.                                */
     /* -------------------------------------------------------------------- */
-    TIFF *l_hTIFF = CreateLL(pszFilename, nXSize, nYSize, l_nBands, eType, 0,
-                             nColorTableMultiplier, papszParamList, &l_fpL,
-                             l_osTmpFilename);
+    bool bTileInterleaving;
+    TIFF *l_hTIFF =
+        CreateLL(pszFilename, nXSize, nYSize, l_nBands, eType, 0,
+                 nColorTableMultiplier, papszParamList, &l_fpL, l_osTmpFilename,
+                 /* bCreateCopy=*/false, bTileInterleaving);
     const bool bStreaming = !l_osTmpFilename.empty();
 
     if (l_hTIFF == nullptr)
@@ -6572,19 +6604,15 @@ CPLErr GTiffDataset::CopyImageryAndMask(GTiffDataset *poDstDS,
     const auto eType = poDstDS->GetRasterBand(1)->GetRasterDataType();
     const int nDataTypeSize = GDALGetDataTypeSizeBytes(eType);
     const int l_nBands = poDstDS->GetRasterCount();
-    void *pBlockBuffer =
+    GByte *pBlockBuffer = static_cast<GByte *>(
         VSI_MALLOC3_VERBOSE(poDstDS->m_nBlockXSize, poDstDS->m_nBlockYSize,
-                            cpl::fits_on<int>(l_nBands * nDataTypeSize));
+                            cpl::fits_on<int>(l_nBands * nDataTypeSize)));
     if (pBlockBuffer == nullptr)
     {
         eErr = CE_Failure;
     }
     const int nYSize = poDstDS->nRasterYSize;
     const int nXSize = poDstDS->nRasterXSize;
-    const int nBlocks = poDstDS->m_nBlocksPerBand;
-
-    CPLAssert(l_nBands == 1 || poDstDS->m_nPlanarConfig == PLANARCONFIG_CONTIG);
-
     const bool bIsOddBand =
         dynamic_cast<GTiffOddBitsBand *>(poDstDS->GetRasterBand(1)) != nullptr;
 
@@ -6594,134 +6622,276 @@ CPLErr GTiffDataset::CopyImageryAndMask(GTiffDataset *poDstDS,
         CPLAssert(poDstDS->m_poMaskDS->m_nBlockYSize == poDstDS->m_nBlockYSize);
     }
 
-    int iBlock = 0;
-    for (int iY = 0, nYBlock = 0; iY < nYSize && eErr == CE_None;
-         iY = ((nYSize - iY < poDstDS->m_nBlockYSize)
-                   ? nYSize
-                   : iY + poDstDS->m_nBlockYSize),
-             nYBlock++)
+    if (poDstDS->m_nPlanarConfig == PLANARCONFIG_SEPARATE &&
+        !poDstDS->m_bTileInterleave)
     {
-        const int nReqYSize = std::min(nYSize - iY, poDstDS->m_nBlockYSize);
-        for (int iX = 0, nXBlock = 0; iX < nXSize && eErr == CE_None;
-             iX = ((nXSize - iX < poDstDS->m_nBlockXSize)
-                       ? nXSize
-                       : iX + poDstDS->m_nBlockXSize),
-                 nXBlock++)
+        int iBlock = 0;
+        const int nBlocks = poDstDS->m_nBlocksPerBand *
+                            (l_nBands + (poDstDS->m_poMaskDS ? 1 : 0));
+        for (int i = 0; eErr == CE_None && i < l_nBands; i++)
         {
-            const int nReqXSize = std::min(nXSize - iX, poDstDS->m_nBlockXSize);
-            if (nReqXSize < poDstDS->m_nBlockXSize ||
-                nReqYSize < poDstDS->m_nBlockYSize)
+            for (int iY = 0, nYBlock = 0; iY < nYSize && eErr == CE_None;
+                 iY = ((nYSize - iY < poDstDS->m_nBlockYSize)
+                           ? nYSize
+                           : iY + poDstDS->m_nBlockYSize),
+                     nYBlock++)
             {
-                memset(pBlockBuffer, 0,
-                       static_cast<size_t>(poDstDS->m_nBlockXSize) *
-                           poDstDS->m_nBlockYSize * l_nBands * nDataTypeSize);
-            }
-
-            if (!bIsOddBand)
-            {
-                eErr = poSrcDS->RasterIO(
-                    GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
-                    nReqXSize, nReqYSize, eType, l_nBands, nullptr,
-                    static_cast<GSpacing>(nDataTypeSize) * l_nBands,
-                    static_cast<GSpacing>(nDataTypeSize) * l_nBands *
-                        poDstDS->m_nBlockXSize,
-                    nDataTypeSize, nullptr);
-                if (eErr == CE_None)
+                const int nReqYSize =
+                    std::min(nYSize - iY, poDstDS->m_nBlockYSize);
+                for (int iX = 0, nXBlock = 0; iX < nXSize && eErr == CE_None;
+                     iX = ((nXSize - iX < poDstDS->m_nBlockXSize)
+                               ? nXSize
+                               : iX + poDstDS->m_nBlockXSize),
+                         nXBlock++)
                 {
-                    eErr = poDstDS->WriteEncodedTileOrStrip(
-                        iBlock, pBlockBuffer, false);
+                    const int nReqXSize =
+                        std::min(nXSize - iX, poDstDS->m_nBlockXSize);
+                    if (nReqXSize < poDstDS->m_nBlockXSize ||
+                        nReqYSize < poDstDS->m_nBlockYSize)
+                    {
+                        memset(pBlockBuffer, 0,
+                               static_cast<size_t>(poDstDS->m_nBlockXSize) *
+                                   poDstDS->m_nBlockYSize * nDataTypeSize);
+                    }
+                    eErr = poSrcDS->GetRasterBand(i + 1)->RasterIO(
+                        GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
+                        nReqXSize, nReqYSize, eType, nDataTypeSize,
+                        static_cast<GSpacing>(nDataTypeSize) *
+                            poDstDS->m_nBlockXSize,
+                        nullptr);
+                    if (eErr == CE_None)
+                    {
+                        eErr = poDstDS->WriteEncodedTileOrStrip(
+                            iBlock, pBlockBuffer, false);
+                    }
+
+                    iBlock++;
+                    if (pfnProgress &&
+                        !pfnProgress(static_cast<double>(iBlock) / nBlocks,
+                                     nullptr, pProgressData))
+                    {
+                        eErr = CE_Failure;
+                    }
+
+                    if (poDstDS->m_bWriteError)
+                        eErr = CE_Failure;
                 }
             }
-            else
+        }
+        if (poDstDS->m_poMaskDS && eErr == CE_None)
+        {
+            int iBlockMask = 0;
+            for (int iY = 0, nYBlock = 0; iY < nYSize && eErr == CE_None;
+                 iY = ((nYSize - iY < poDstDS->m_nBlockYSize)
+                           ? nYSize
+                           : iY + poDstDS->m_nBlockYSize),
+                     nYBlock++)
             {
-                // In the odd bit case, this is a bit messy to ensure
-                // the strile gets written synchronously.
-                // We load the content of the n-1 bands in the cache,
-                // and for the last band we invoke WriteBlock() directly
-                // We also force FlushBlockBuf()
-                std::vector<GDALRasterBlock *> apoLockedBlocks;
-                for (int i = 0; eErr == CE_None && i < l_nBands - 1; i++)
+                const int nReqYSize =
+                    std::min(nYSize - iY, poDstDS->m_nBlockYSize);
+                for (int iX = 0, nXBlock = 0; iX < nXSize && eErr == CE_None;
+                     iX = ((nXSize - iX < poDstDS->m_nBlockXSize)
+                               ? nXSize
+                               : iX + poDstDS->m_nBlockXSize),
+                         nXBlock++)
                 {
-                    auto poBlock =
-                        poDstDS->GetRasterBand(i + 1)->GetLockedBlockRef(
-                            nXBlock, nYBlock, TRUE);
-                    if (poBlock)
+                    const int nReqXSize =
+                        std::min(nXSize - iX, poDstDS->m_nBlockXSize);
+                    if (nReqXSize < poDstDS->m_nBlockXSize ||
+                        nReqYSize < poDstDS->m_nBlockYSize)
                     {
-                        eErr = poSrcDS->GetRasterBand(i + 1)->RasterIO(
-                            GF_Read, iX, iY, nReqXSize, nReqYSize,
-                            poBlock->GetDataRef(), nReqXSize, nReqYSize, eType,
-                            nDataTypeSize,
-                            static_cast<GSpacing>(nDataTypeSize) *
-                                poDstDS->m_nBlockXSize,
-                            nullptr);
-                        poBlock->MarkDirty();
-                        apoLockedBlocks.emplace_back(poBlock);
+                        memset(pBlockBuffer, 0,
+                               static_cast<size_t>(poDstDS->m_nBlockXSize) *
+                                   poDstDS->m_nBlockYSize * nDataTypeSize);
                     }
-                    else
+                    eErr = poSrcMaskBand->RasterIO(
+                        GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
+                        nReqXSize, nReqYSize, eType, nDataTypeSize,
+                        static_cast<GSpacing>(nDataTypeSize) *
+                            poDstDS->m_nBlockXSize,
+                        nullptr);
+                    if (eErr == CE_None)
+                    {
+                        eErr = poDstDS->m_poMaskDS->WriteEncodedTileOrStrip(
+                            iBlockMask, pBlockBuffer, false);
+                    }
+
+                    iBlockMask++;
+                    if (pfnProgress &&
+                        !pfnProgress(static_cast<double>(iBlock + iBlockMask) /
+                                         nBlocks,
+                                     nullptr, pProgressData))
                     {
                         eErr = CE_Failure;
                     }
+
+                    if (poDstDS->m_poMaskDS->m_bWriteError)
+                        eErr = CE_Failure;
+                }
+            }
+        }
+    }
+    else
+    {
+        int iBlock = 0;
+        const int nBlocks = poDstDS->m_nBlocksPerBand;
+        for (int iY = 0, nYBlock = 0; iY < nYSize && eErr == CE_None;
+             iY = ((nYSize - iY < poDstDS->m_nBlockYSize)
+                       ? nYSize
+                       : iY + poDstDS->m_nBlockYSize),
+                 nYBlock++)
+        {
+            const int nReqYSize = std::min(nYSize - iY, poDstDS->m_nBlockYSize);
+            for (int iX = 0, nXBlock = 0; iX < nXSize && eErr == CE_None;
+                 iX = ((nXSize - iX < poDstDS->m_nBlockXSize)
+                           ? nXSize
+                           : iX + poDstDS->m_nBlockXSize),
+                     nXBlock++)
+            {
+                const int nReqXSize =
+                    std::min(nXSize - iX, poDstDS->m_nBlockXSize);
+                if (nReqXSize < poDstDS->m_nBlockXSize ||
+                    nReqYSize < poDstDS->m_nBlockYSize)
+                {
+                    memset(pBlockBuffer, 0,
+                           static_cast<size_t>(poDstDS->m_nBlockXSize) *
+                               poDstDS->m_nBlockYSize * l_nBands *
+                               nDataTypeSize);
                 }
-                if (eErr == CE_None)
+
+                if (poDstDS->m_bTileInterleave)
                 {
-                    eErr = poSrcDS->GetRasterBand(l_nBands)->RasterIO(
+                    eErr = poSrcDS->RasterIO(
                         GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
-                        nReqXSize, nReqYSize, eType, nDataTypeSize,
+                        nReqXSize, nReqYSize, eType, l_nBands, nullptr,
+                        nDataTypeSize,
                         static_cast<GSpacing>(nDataTypeSize) *
                             poDstDS->m_nBlockXSize,
+                        static_cast<GSpacing>(nDataTypeSize) *
+                            poDstDS->m_nBlockXSize * poDstDS->m_nBlockYSize,
                         nullptr);
+                    if (eErr == CE_None)
+                    {
+                        for (int i = 0; eErr == CE_None && i < l_nBands; i++)
+                        {
+                            eErr = poDstDS->WriteEncodedTileOrStrip(
+                                iBlock + i * poDstDS->m_nBlocksPerBand,
+                                pBlockBuffer + static_cast<size_t>(i) *
+                                                   poDstDS->m_nBlockXSize *
+                                                   poDstDS->m_nBlockYSize *
+                                                   nDataTypeSize,
+                                false);
+                        }
+                    }
                 }
-                if (eErr == CE_None)
+                else if (!bIsOddBand)
                 {
-                    // Avoid any attempt to load from disk
-                    poDstDS->m_nLoadedBlock = iBlock;
-                    eErr = poDstDS->GetRasterBand(l_nBands)->WriteBlock(
-                        nXBlock, nYBlock, pBlockBuffer);
+                    eErr = poSrcDS->RasterIO(
+                        GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
+                        nReqXSize, nReqYSize, eType, l_nBands, nullptr,
+                        static_cast<GSpacing>(nDataTypeSize) * l_nBands,
+                        static_cast<GSpacing>(nDataTypeSize) * l_nBands *
+                            poDstDS->m_nBlockXSize,
+                        nDataTypeSize, nullptr);
                     if (eErr == CE_None)
-                        eErr = poDstDS->FlushBlockBuf();
+                    {
+                        eErr = poDstDS->WriteEncodedTileOrStrip(
+                            iBlock, pBlockBuffer, false);
+                    }
                 }
-                for (auto poBlock : apoLockedBlocks)
+                else
                 {
-                    poBlock->MarkClean();
-                    poBlock->DropLock();
+                    // In the odd bit case, this is a bit messy to ensure
+                    // the strile gets written synchronously.
+                    // We load the content of the n-1 bands in the cache,
+                    // and for the last band we invoke WriteBlock() directly
+                    // We also force FlushBlockBuf()
+                    std::vector<GDALRasterBlock *> apoLockedBlocks;
+                    for (int i = 0; eErr == CE_None && i < l_nBands - 1; i++)
+                    {
+                        auto poBlock =
+                            poDstDS->GetRasterBand(i + 1)->GetLockedBlockRef(
+                                nXBlock, nYBlock, TRUE);
+                        if (poBlock)
+                        {
+                            eErr = poSrcDS->GetRasterBand(i + 1)->RasterIO(
+                                GF_Read, iX, iY, nReqXSize, nReqYSize,
+                                poBlock->GetDataRef(), nReqXSize, nReqYSize,
+                                eType, nDataTypeSize,
+                                static_cast<GSpacing>(nDataTypeSize) *
+                                    poDstDS->m_nBlockXSize,
+                                nullptr);
+                            poBlock->MarkDirty();
+                            apoLockedBlocks.emplace_back(poBlock);
+                        }
+                        else
+                        {
+                            eErr = CE_Failure;
+                        }
+                    }
+                    if (eErr == CE_None)
+                    {
+                        eErr = poSrcDS->GetRasterBand(l_nBands)->RasterIO(
+                            GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
+                            nReqXSize, nReqYSize, eType, nDataTypeSize,
+                            static_cast<GSpacing>(nDataTypeSize) *
+                                poDstDS->m_nBlockXSize,
+                            nullptr);
+                    }
+                    if (eErr == CE_None)
+                    {
+                        // Avoid any attempt to load from disk
+                        poDstDS->m_nLoadedBlock = iBlock;
+                        eErr = poDstDS->GetRasterBand(l_nBands)->WriteBlock(
+                            nXBlock, nYBlock, pBlockBuffer);
+                        if (eErr == CE_None)
+                            eErr = poDstDS->FlushBlockBuf();
+                    }
+                    for (auto poBlock : apoLockedBlocks)
+                    {
+                        poBlock->MarkClean();
+                        poBlock->DropLock();
+                    }
                 }
-            }
 
-            if (eErr == CE_None && poDstDS->m_poMaskDS)
-            {
-                if (nReqXSize < poDstDS->m_nBlockXSize ||
-                    nReqYSize < poDstDS->m_nBlockYSize)
+                if (eErr == CE_None && poDstDS->m_poMaskDS)
                 {
-                    memset(pBlockBuffer, 0,
-                           static_cast<size_t>(poDstDS->m_nBlockXSize) *
-                               poDstDS->m_nBlockYSize);
-                }
-                eErr = poSrcMaskBand->RasterIO(
-                    GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
-                    nReqXSize, nReqYSize, GDT_Byte, 1, poDstDS->m_nBlockXSize,
-                    nullptr);
-                if (eErr == CE_None)
-                {
-                    // Avoid any attempt to load from disk
-                    poDstDS->m_poMaskDS->m_nLoadedBlock = iBlock;
-                    eErr = poDstDS->m_poMaskDS->GetRasterBand(1)->WriteBlock(
-                        nXBlock, nYBlock, pBlockBuffer);
+                    if (nReqXSize < poDstDS->m_nBlockXSize ||
+                        nReqYSize < poDstDS->m_nBlockYSize)
+                    {
+                        memset(pBlockBuffer, 0,
+                               static_cast<size_t>(poDstDS->m_nBlockXSize) *
+                                   poDstDS->m_nBlockYSize);
+                    }
+                    eErr = poSrcMaskBand->RasterIO(
+                        GF_Read, iX, iY, nReqXSize, nReqYSize, pBlockBuffer,
+                        nReqXSize, nReqYSize, GDT_Byte, 1,
+                        poDstDS->m_nBlockXSize, nullptr);
                     if (eErr == CE_None)
-                        eErr = poDstDS->m_poMaskDS->FlushBlockBuf();
+                    {
+                        // Avoid any attempt to load from disk
+                        poDstDS->m_poMaskDS->m_nLoadedBlock = iBlock;
+                        eErr =
+                            poDstDS->m_poMaskDS->GetRasterBand(1)->WriteBlock(
+                                nXBlock, nYBlock, pBlockBuffer);
+                        if (eErr == CE_None)
+                            eErr = poDstDS->m_poMaskDS->FlushBlockBuf();
+                    }
                 }
-            }
-            if (poDstDS->m_bWriteError)
-                eErr = CE_Failure;
+                if (poDstDS->m_bWriteError)
+                    eErr = CE_Failure;
 
-            iBlock++;
-            if (pfnProgress &&
-                !pfnProgress(static_cast<double>(iBlock) / nBlocks, nullptr,
-                             pProgressData))
-            {
-                eErr = CE_Failure;
+                iBlock++;
+                if (pfnProgress &&
+                    !pfnProgress(static_cast<double>(iBlock) / nBlocks, nullptr,
+                                 pProgressData))
+                {
+                    eErr = CE_Failure;
+                }
             }
         }
     }
+
     poDstDS->FlushCache(false);  // mostly to wait for thread completion
     VSIFree(pBlockBuffer);
 
@@ -7044,9 +7214,11 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
                      papszOptions, "COLOR_TABLE_MULTIPLIER",
                      CPLSPrintf("%d", DEFAULT_COLOR_TABLE_MULTIPLIER_257)))));
 
+    bool bTileInterleaving = false;
     TIFF *l_hTIFF = CreateLL(pszFilename, nXSize, nYSize, l_nBands, eType,
                              dfExtraSpaceForOverviews, nColorTableMultiplier,
-                             papszCreateOptions, &l_fpL, l_osTmpFilename);
+                             papszCreateOptions, &l_fpL, l_osTmpFilename,
+                             /* bCreateCopy = */ true, bTileInterleaving);
     const bool bStreaming = !l_osTmpFilename.empty();
 
     CSLDestroy(papszCreateOptions);
@@ -7258,13 +7430,20 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
     const int nMaskFlags = poSrcDS->GetRasterBand(1)->GetMaskFlags();
     bool bCreateMask = false;
     CPLString osHiddenStructuralMD;
-    if ((l_nBands == 1 || l_nPlanarConfig == PLANARCONFIG_CONTIG) &&
-        bCopySrcOverviews)
+    const char *pszInterleave =
+        CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "PIXEL");
+    if (bCopySrcOverviews)
     {
         osHiddenStructuralMD += "LAYOUT=IFDS_BEFORE_DATA\n";
         osHiddenStructuralMD += "BLOCK_ORDER=ROW_MAJOR\n";
         osHiddenStructuralMD += "BLOCK_LEADER=SIZE_AS_UINT4\n";
         osHiddenStructuralMD += "BLOCK_TRAILER=LAST_4_BYTES_REPEATED\n";
+        if (l_nBands > 1 && !EQUAL(pszInterleave, "PIXEL"))
+        {
+            osHiddenStructuralMD += "INTERLEAVE=";
+            osHiddenStructuralMD += CPLString(pszInterleave).toupper();
+            osHiddenStructuralMD += "\n";
+        }
         osHiddenStructuralMD +=
             "KNOWN_INCOMPATIBLE_EDITION=NO\n ";  // Final space intended, so
                                                  // this can be replaced by YES
@@ -7274,7 +7453,7 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
     {
         bCreateMask = true;
         if (GTiffDataset::MustCreateInternalMask() &&
-            !osHiddenStructuralMD.empty())
+            !osHiddenStructuralMD.empty() && EQUAL(pszInterleave, "PIXEL"))
         {
             osHiddenStructuralMD += "MASK_INTERLEAVED_WITH_IMAGERY=YES\n";
         }
@@ -7627,6 +7806,13 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
     poDS->m_fpL = l_fpL;
     poDS->m_bIMDRPCMetadataLoaded = true;
     poDS->m_nColorTableMultiplier = nColorTableMultiplier;
+    poDS->m_bTileInterleave = bTileInterleaving;
+
+    if (bTileInterleaving)
+    {
+        poDS->m_oGTiffMDMD.SetMetadataItem("INTERLEAVE", "TILE",
+                                           "IMAGE_STRUCTURE");
+    }
 
     const bool bAppend = CPLFetchBool(papszOptions, "APPEND_SUBDATASET", false);
     if (poDS->OpenOffset(l_hTIFF,
@@ -8060,10 +8246,6 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
                 }
             }
 
-            char *papszCopyWholeRasterOptions[2] = {nullptr, nullptr};
-            if (l_nCompression != COMPRESSION_NONE)
-                papszCopyWholeRasterOptions[0] =
-                    const_cast<char *>("COMPRESSED=YES");
             // Now copy the imagery.
             // Begin with the smallest overview.
             for (int iOvrLevel = nSrcOverviews - 1;
@@ -8099,6 +8281,7 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
                 poDstDS->m_bFillEmptyTilesAtClosing =
                     poDS->m_bFillEmptyTilesAtClosing;
                 poDstDS->m_bWriteEmptyTiles = poDS->m_bWriteEmptyTiles;
+                poDstDS->m_bTileInterleave = poDS->m_bTileInterleave;
                 GDALRasterBand *poSrcMaskBand = nullptr;
                 if (poDstDS->m_poMaskDS)
                 {
@@ -8120,68 +8303,22 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
                             : poSrcOvrBand->GetMaskBand();
                 }
 
-                if (l_nBands == 1 ||
-                    poDstDS->m_nPlanarConfig == PLANARCONFIG_CONTIG)
+                if (poDstDS->m_poMaskDS)
                 {
-                    if (poDstDS->m_poMaskDS)
-                    {
-                        dfNextCurPixels +=
-                            static_cast<double>(poSrcOvrBand->GetXSize()) *
-                            poSrcOvrBand->GetYSize();
-                    }
-                    void *pScaledData = GDALCreateScaledProgress(
-                        dfCurPixels / dfTotalPixels,
-                        dfNextCurPixels / dfTotalPixels, pfnProgress,
-                        pProgressData);
-
-                    eErr =
-                        CopyImageryAndMask(poDstDS, poSrcOvrDS, poSrcMaskBand,
-                                           GDALScaledProgress, pScaledData);
-
-                    dfCurPixels = dfNextCurPixels;
-                    GDALDestroyScaledProgress(pScaledData);
+                    dfNextCurPixels +=
+                        static_cast<double>(poSrcOvrBand->GetXSize()) *
+                        poSrcOvrBand->GetYSize();
                 }
-                else
-                {
-                    void *pScaledData = GDALCreateScaledProgress(
-                        dfCurPixels / dfTotalPixels,
-                        dfNextCurPixels / dfTotalPixels, pfnProgress,
-                        pProgressData);
-
-                    eErr = GDALDatasetCopyWholeRaster(
-                        GDALDataset::ToHandle(poSrcOvrDS),
-                        GDALDataset::ToHandle(poDstDS),
-                        papszCopyWholeRasterOptions, GDALScaledProgress,
-                        pScaledData);
-
-                    dfCurPixels = dfNextCurPixels;
-                    GDALDestroyScaledProgress(pScaledData);
+                void *pScaledData =
+                    GDALCreateScaledProgress(dfCurPixels / dfTotalPixels,
+                                             dfNextCurPixels / dfTotalPixels,
+                                             pfnProgress, pProgressData);
 
-                    poDstDS->FlushCache(false);
+                eErr = CopyImageryAndMask(poDstDS, poSrcOvrDS, poSrcMaskBand,
+                                          GDALScaledProgress, pScaledData);
 
-                    // Copy mask of the overview.
-                    if (eErr == CE_None &&
-                        (poMaskOvrDS ||
-                         poSrcOvrBand->GetMaskFlags() == GMF_PER_DATASET) &&
-                        poDstDS->m_poMaskDS != nullptr)
-                    {
-                        dfNextCurPixels +=
-                            static_cast<double>(poSrcOvrBand->GetXSize()) *
-                            poSrcOvrBand->GetYSize();
-                        pScaledData = GDALCreateScaledProgress(
-                            dfCurPixels / dfTotalPixels,
-                            dfNextCurPixels / dfTotalPixels, pfnProgress,
-                            pProgressData);
-                        eErr = GDALRasterBandCopyWholeRaster(
-                            poSrcMaskBand,
-                            poDstDS->m_poMaskDS->GetRasterBand(1),
-                            papszCopyWholeRasterOptions, GDALScaledProgress,
-                            pScaledData);
-                        dfCurPixels = dfNextCurPixels;
-                        GDALDestroyScaledProgress(pScaledData);
-                        poDstDS->m_poMaskDS->FlushCache(false);
-                    }
-                }
+                dfCurPixels = dfNextCurPixels;
+                GDALDestroyScaledProgress(pScaledData);
 
                 if (poSrcOvrDS != poOvrDS.get())
                     delete poSrcOvrDS;
@@ -8336,6 +8473,7 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
         {
             papszCopyWholeRasterOptions[iNextOption++] = "COMPRESSED=YES";
         }
+
         // For streaming with separate, we really want that bands are written
         // after each other, even if the source is pixel interleaved.
         else if (bStreaming && poDS->m_nPlanarConfig == PLANARCONFIG_SEPARATE)
@@ -8343,21 +8481,17 @@ GDALDataset *GTiffDataset::CreateCopy(const char *pszFilename,
             papszCopyWholeRasterOptions[iNextOption++] = "INTERLEAVE=BAND";
         }
 
-        if (bCopySrcOverviews &&
-            (l_nBands == 1 || poDS->m_nPlanarConfig == PLANARCONFIG_CONTIG))
+        if (bCopySrcOverviews || bTileInterleaving)
         {
             poDS->m_bBlockOrderRowMajor = true;
-            poDS->m_bLeaderSizeAsUInt4 = true;
-            poDS->m_bTrailerRepeatedLast4BytesRepeated = true;
+            poDS->m_bLeaderSizeAsUInt4 = bCopySrcOverviews;
+            poDS->m_bTrailerRepeatedLast4BytesRepeated = bCopySrcOverviews;
             if (poDS->m_poMaskDS)
             {
                 poDS->m_poMaskDS->m_bBlockOrderRowMajor = true;
-                poDS->m_poMaskDS->m_bLeaderSizeAsUInt4 = true;
-                poDS->m_poMaskDS->m_bTrailerRepeatedLast4BytesRepeated = true;
-            }
-
-            if (poDS->m_poMaskDS)
-            {
+                poDS->m_poMaskDS->m_bLeaderSizeAsUInt4 = bCopySrcOverviews;
+                poDS->m_poMaskDS->m_bTrailerRepeatedLast4BytesRepeated =
+                    bCopySrcOverviews;
                 GDALDestroyScaledProgress(pScaledData);
                 pScaledData =
                     GDALCreateScaledProgress(dfCurPixels / dfTotalPixels, 1.0,
diff --git a/frmts/gtiff/gtiffrasterband.cpp b/frmts/gtiff/gtiffrasterband.cpp
index cbb2a451ccda..ab77f9cb910d 100644
--- a/frmts/gtiff/gtiffrasterband.cpp
+++ b/frmts/gtiff/gtiffrasterband.cpp
@@ -315,19 +315,21 @@ CPLErr GTiffRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
         else
         {
             bCanUseMultiThreadedRead = false;
-            GTiffRasterBand *poBandForCache = this;
+            GTiffDataset *poDSForCache = m_poGDS;
+            int nBandForCache = nBand;
             if (!m_poGDS->m_bStreamingIn && m_poGDS->m_bBlockOrderRowMajor &&
                 m_poGDS->m_bLeaderSizeAsUInt4 &&
                 m_poGDS->m_bMaskInterleavedWithImagery &&
                 m_poGDS->m_poImageryDS)
             {
-                poBandForCache = cpl::down_cast<GTiffRasterBand *>(
-                    m_poGDS->m_poImageryDS->GetRasterBand(1));
+                poDSForCache = m_poGDS->m_poImageryDS;
+                nBandForCache = 1;
             }
-            bufferedDataFreer.Init(poBandForCache->CacheMultiRange(
-                                       nXOff, nYOff, nXSize, nYSize, nBufXSize,
-                                       nBufYSize, psExtraArg),
-                                   poBandForCache->m_poGDS->m_hTIFF);
+            bufferedDataFreer.Init(
+                poDSForCache->CacheMultiRange(nXOff, nYOff, nXSize, nYSize,
+                                              nBufXSize, nBufYSize,
+                                              &nBandForCache, 1, psExtraArg),
+                poDSForCache->m_hTIFF);
         }
     }
 
@@ -379,7 +381,8 @@ CPLErr GTiffRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
              !m_poGDS->m_bLoadedBlockDirty &&
              (m_poGDS->nBands == 1 ||
               m_poGDS->m_nPlanarConfig == PLANARCONFIG_SEPARATE) &&
-             (nXOff % nBlockXSize) == 0 && (nYOff % nBlockYSize) == 0 &&
+             !m_poGDS->m_bLeaderSizeAsUInt4 && (nXOff % nBlockXSize) == 0 &&
+             (nYOff % nBlockYSize) == 0 &&
              (nXOff + nXSize == nRasterXSize || (nXSize % nBlockXSize) == 0) &&
              (nYOff + nYSize == nRasterYSize || (nYSize % nBlockYSize) == 0))
     {
diff --git a/frmts/gtiff/gtiffrasterband.h b/frmts/gtiff/gtiffrasterband.h
index 8aab99209a79..ff250d51f5c8 100644
--- a/frmts/gtiff/gtiffrasterband.h
+++ b/frmts/gtiff/gtiffrasterband.h
@@ -55,10 +55,6 @@ class GTiffRasterBand CPL_NON_FINAL : public GDALPamRasterBand
                                              GIntBig *pnLineSpace,
                                              char **papszOptions);
 
-    void *CacheMultiRange(int nXOff, int nYOff, int nXSize, int nYSize,
-                          int nBufXSize, int nBufYSize,
-                          GDALRasterIOExtraArg *psExtraArg);
-
   protected:
     GTiffDataset *m_poGDS = nullptr;
     GDALMultiDomainMetadata m_oGTiffMDMD{};
diff --git a/frmts/gtiff/gtiffrasterband_read.cpp b/frmts/gtiff/gtiffrasterband_read.cpp
index 44097ff8fc31..6d36652121a6 100644
--- a/frmts/gtiff/gtiffrasterband_read.cpp
+++ b/frmts/gtiff/gtiffrasterband_read.cpp
@@ -630,488 +630,6 @@ CPLVirtualMem *GTiffRasterBand::GetVirtualMemAutoInternal(GDALRWFlag eRWFlag,
     return pVMem;
 }
 
-/************************************************************************/
-/*                         CacheMultiRange()                            */
-/************************************************************************/
-
-static bool CheckTrailer(const GByte *strileData, vsi_l_offset nStrileSize)
-{
-    GByte abyTrailer[4];
-    memcpy(abyTrailer, strileData + nStrileSize, 4);
-    GByte abyLastBytes[4] = {};
-    if (nStrileSize >= 4)
-        memcpy(abyLastBytes, strileData + nStrileSize - 4, 4);
-    else
-    {
-        // The last bytes will be zero due to the above {} initialization,
-        // and that's what should be in abyTrailer too when the trailer is
-        // correct.
-        memcpy(abyLastBytes, strileData, static_cast<size_t>(nStrileSize));
-    }
-    return memcmp(abyTrailer, abyLastBytes, 4) == 0;
-}
-
-void *GTiffRasterBand::CacheMultiRange(int nXOff, int nYOff, int nXSize,
-                                       int nYSize, int nBufXSize, int nBufYSize,
-                                       GDALRasterIOExtraArg *psExtraArg)
-{
-    void *pBufferedData = nullptr;
-    // Same logic as in GDALRasterBand::IRasterIO()
-    double dfXOff = nXOff;
-    double dfYOff = nYOff;
-    double dfXSize = nXSize;
-    double dfYSize = nYSize;
-    if (psExtraArg->bFloatingPointWindowValidity)
-    {
-        dfXOff = psExtraArg->dfXOff;
-        dfYOff = psExtraArg->dfYOff;
-        dfXSize = psExtraArg->dfXSize;
-        dfYSize = psExtraArg->dfYSize;
-    }
-    const double dfSrcXInc = dfXSize / static_cast<double>(nBufXSize);
-    const double dfSrcYInc = dfYSize / static_cast<double>(nBufYSize);
-    const double EPS = 1e-10;
-    const int nBlockX1 =
-        static_cast<int>(std::max(0.0, (0 + 0.5) * dfSrcXInc + dfXOff + EPS)) /
-        nBlockXSize;
-    const int nBlockY1 =
-        static_cast<int>(std::max(0.0, (0 + 0.5) * dfSrcYInc + dfYOff + EPS)) /
-        nBlockYSize;
-    const int nBlockX2 =
-        static_cast<int>(
-            std::min(static_cast<double>(nRasterXSize - 1),
-                     (nBufXSize - 1 + 0.5) * dfSrcXInc + dfXOff + EPS)) /
-        nBlockXSize;
-    const int nBlockY2 =
-        static_cast<int>(
-            std::min(static_cast<double>(nRasterYSize - 1),
-                     (nBufYSize - 1 + 0.5) * dfSrcYInc + dfYOff + EPS)) /
-        nBlockYSize;
-
-    const int nBlockCount = nBlocksPerRow * nBlocksPerColumn;
-
-    struct StrileData
-    {
-        vsi_l_offset nOffset;
-        vsi_l_offset nByteCount;
-        bool bTryMask;
-    };
-
-    std::map<int, StrileData> oMapStrileToOffsetByteCount;
-
-    // Dedicated method to retrieved the offset and size in an efficient way
-    // when m_bBlockOrderRowMajor and m_bLeaderSizeAsUInt4 conditions are
-    // met.
-    // Except for the last block, we just read the offset from the TIFF offset
-    // array, and retrieve the size in the leader 4 bytes that come before the
-    // payload.
-    auto OptimizedRetrievalOfOffsetSize =
-        [&](int nBlockId, vsi_l_offset &nOffset, vsi_l_offset &nSize,
-            size_t nTotalSize, size_t nMaxRawBlockCacheSize)
-    {
-        bool bTryMask = m_poGDS->m_bMaskInterleavedWithImagery;
-        nOffset = TIFFGetStrileOffset(m_poGDS->m_hTIFF, nBlockId);
-        if (nOffset >= 4)
-        {
-            if (nBlockId == nBlockCount - 1)
-            {
-                // Special case for the last block. As there is no next block
-                // from which to retrieve an offset, use the good old method
-                // that consists in reading the ByteCount array.
-                if (bTryMask && m_poGDS->GetRasterBand(1)->GetMaskBand() &&
-                    m_poGDS->m_poMaskDS)
-                {
-                    auto nMaskOffset = TIFFGetStrileOffset(
-                        m_poGDS->m_poMaskDS->m_hTIFF, nBlockId);
-                    if (nMaskOffset)
-                    {
-                        nSize = nMaskOffset +
-                                TIFFGetStrileByteCount(
-                                    m_poGDS->m_poMaskDS->m_hTIFF, nBlockId) -
-                                nOffset;
-                    }
-                    else
-                    {
-                        bTryMask = false;
-                    }
-                }
-                if (nSize == 0)
-                {
-                    nSize = TIFFGetStrileByteCount(m_poGDS->m_hTIFF, nBlockId);
-                }
-                if (nSize && m_poGDS->m_bTrailerRepeatedLast4BytesRepeated)
-                {
-                    nSize += 4;
-                }
-            }
-            else
-            {
-                auto nOffsetNext =
-                    TIFFGetStrileOffset(m_poGDS->m_hTIFF, nBlockId + 1);
-                if (nOffsetNext > nOffset)
-                {
-                    nSize = nOffsetNext - nOffset;
-                }
-                else
-                {
-                    // Shouldn't happen for a compliant file
-                    if (nOffsetNext != 0)
-                    {
-                        CPLDebug("GTiff", "Tile %d is not located after %d",
-                                 nBlockId + 1, nBlockId);
-                    }
-                    bTryMask = false;
-                    nSize = TIFFGetStrileByteCount(m_poGDS->m_hTIFF, nBlockId);
-                    if (m_poGDS->m_bTrailerRepeatedLast4BytesRepeated)
-                        nSize += 4;
-                }
-            }
-            if (nSize)
-            {
-                nOffset -= 4;
-                nSize += 4;
-                if (nTotalSize + nSize < nMaxRawBlockCacheSize)
-                {
-                    StrileData data;
-                    data.nOffset = nOffset;
-                    data.nByteCount = nSize;
-                    data.bTryMask = bTryMask;
-                    oMapStrileToOffsetByteCount[nBlockId] = data;
-                }
-            }
-        }
-        else
-        {
-            // Sparse tile
-            StrileData data;
-            data.nOffset = 0;
-            data.nByteCount = 0;
-            data.bTryMask = false;
-            oMapStrileToOffsetByteCount[nBlockId] = data;
-        }
-    };
-
-    // This lambda fills m_poDS->m_oCacheStrileToOffsetByteCount (and
-    // m_poDS->m_poMaskDS->m_oCacheStrileToOffsetByteCount, when there is a
-    // mask) from the temporary oMapStrileToOffsetByteCount.
-    auto FillCacheStrileToOffsetByteCount =
-        [&](const std::vector<vsi_l_offset> &anOffsets,
-            const std::vector<size_t> &anSizes,
-            const std::vector<void *> &apData)
-    {
-        CPLAssert(m_poGDS->m_bLeaderSizeAsUInt4);
-        size_t i = 0;
-        vsi_l_offset nLastOffset = 0;
-        for (const auto &entry : oMapStrileToOffsetByteCount)
-        {
-            const auto nBlockId = entry.first;
-            const auto nOffset = entry.second.nOffset;
-            const auto nSize = entry.second.nByteCount;
-            if (nOffset == 0)
-            {
-                // Sparse tile
-                m_poGDS->m_oCacheStrileToOffsetByteCount.insert(
-                    nBlockId, std::pair(0, 0));
-                continue;
-            }
-
-            if (nOffset < nLastOffset)
-            {
-                // shouldn't happen normally if tiles are sorted
-                i = 0;
-            }
-            nLastOffset = nOffset;
-            while (i < anOffsets.size() &&
-                   !(nOffset >= anOffsets[i] &&
-                     nOffset + nSize <= anOffsets[i] + anSizes[i]))
-            {
-                i++;
-            }
-            CPLAssert(i < anOffsets.size());
-            CPLAssert(nOffset >= anOffsets[i]);
-            CPLAssert(nOffset + nSize <= anOffsets[i] + anSizes[i]);
-            GUInt32 nSizeFromLeader;
-            memcpy(&nSizeFromLeader,
-                   // cppcheck-suppress containerOutOfBounds
-                   static_cast<GByte *>(apData[i]) + nOffset - anOffsets[i],
-                   sizeof(nSizeFromLeader));
-            CPL_LSBPTR32(&nSizeFromLeader);
-            bool bOK = true;
-            constexpr int nLeaderSize = 4;
-            const int nTrailerSize =
-                (m_poGDS->m_bTrailerRepeatedLast4BytesRepeated ? 4 : 0);
-            if (nSizeFromLeader > nSize - nLeaderSize - nTrailerSize)
-            {
-                CPLDebug("GTiff",
-                         "Inconsistent block size from in leader of block %d",
-                         nBlockId);
-                bOK = false;
-            }
-            else if (m_poGDS->m_bTrailerRepeatedLast4BytesRepeated)
-            {
-                // Check trailer consistency
-                const GByte *strileData = static_cast<GByte *>(apData[i]) +
-                                          nOffset - anOffsets[i] + nLeaderSize;
-                if (!CheckTrailer(strileData, nSizeFromLeader))
-                {
-                    CPLDebug("GTiff", "Inconsistent trailer of block %d",
-                             nBlockId);
-                    bOK = false;
-                }
-            }
-            if (!bOK)
-            {
-                return false;
-            }
-
-            {
-                const vsi_l_offset nRealOffset = nOffset + nLeaderSize;
-                const vsi_l_offset nRealSize = nSizeFromLeader;
-#ifdef DEBUG_VERBOSE
-                CPLDebug("GTiff",
-                         "Block %d found at offset " CPL_FRMT_GUIB
-                         " with size " CPL_FRMT_GUIB,
-                         nBlockId, nRealOffset, nRealSize);
-#endif
-                m_poGDS->m_oCacheStrileToOffsetByteCount.insert(
-                    nBlockId, std::pair(nRealOffset, nRealSize));
-            }
-
-            // Processing of mask
-            if (!(entry.second.bTryMask &&
-                  m_poGDS->m_bMaskInterleavedWithImagery &&
-                  m_poGDS->GetRasterBand(1)->GetMaskBand() &&
-                  m_poGDS->m_poMaskDS))
-            {
-                continue;
-            }
-
-            bOK = false;
-            const vsi_l_offset nMaskOffsetWithLeader =
-                nOffset + nLeaderSize + nSizeFromLeader + nTrailerSize;
-            if (nMaskOffsetWithLeader + nLeaderSize <=
-                anOffsets[i] + anSizes[i])
-            {
-                GUInt32 nMaskSizeFromLeader;
-                memcpy(&nMaskSizeFromLeader,
-                       static_cast<GByte *>(apData[i]) + nMaskOffsetWithLeader -
-                           anOffsets[i],
-                       sizeof(nMaskSizeFromLeader));
-                CPL_LSBPTR32(&nMaskSizeFromLeader);
-                if (nMaskOffsetWithLeader + nLeaderSize + nMaskSizeFromLeader +
-                        nTrailerSize <=
-                    anOffsets[i] + anSizes[i])
-                {
-                    bOK = true;
-                    if (m_poGDS->m_bTrailerRepeatedLast4BytesRepeated)
-                    {
-                        // Check trailer consistency
-                        const GByte *strileMaskData =
-                            static_cast<GByte *>(apData[i]) + nOffset -
-                            anOffsets[i] + nLeaderSize + nSizeFromLeader +
-                            nTrailerSize + nLeaderSize;
-                        if (!CheckTrailer(strileMaskData, nMaskSizeFromLeader))
-                        {
-                            CPLDebug("GTiff",
-                                     "Inconsistent trailer of mask of block %d",
-                                     nBlockId);
-                            bOK = false;
-                        }
-                    }
-                }
-                if (bOK)
-                {
-                    const vsi_l_offset nRealOffset = nOffset + nLeaderSize +
-                                                     nSizeFromLeader +
-                                                     nTrailerSize + nLeaderSize;
-                    const vsi_l_offset nRealSize = nMaskSizeFromLeader;
-#ifdef DEBUG_VERBOSE
-                    CPLDebug("GTiff",
-                             "Mask of block %d found at offset " CPL_FRMT_GUIB
-                             " with size " CPL_FRMT_GUIB,
-                             nBlockId, nRealOffset, nRealSize);
-#endif
-
-                    m_poGDS->m_poMaskDS->m_oCacheStrileToOffsetByteCount.insert(
-                        nBlockId, std::pair(nRealOffset, nRealSize));
-                }
-            }
-            if (!bOK)
-            {
-                CPLDebug("GTiff",
-                         "Mask for block %d is not properly interleaved with "
-                         "imagery block",
-                         nBlockId);
-            }
-        }
-        return true;
-    };
-
-    thandle_t th = TIFFClientdata(m_poGDS->m_hTIFF);
-    if (!VSI_TIFFHasCachedRanges(th))
-    {
-        std::vector<std::pair<vsi_l_offset, size_t>> aOffsetSize;
-        size_t nTotalSize = 0;
-        const unsigned int nMaxRawBlockCacheSize = atoi(
-            CPLGetConfigOption("GDAL_MAX_RAW_BLOCK_CACHE_SIZE", "10485760"));
-        bool bGoOn = true;
-        for (int iY = nBlockY1; bGoOn && iY <= nBlockY2; iY++)
-        {
-            for (int iX = nBlockX1; bGoOn && iX <= nBlockX2; iX++)
-            {
-                GDALRasterBlock *poBlock = TryGetLockedBlockRef(iX, iY);
-                if (poBlock != nullptr)
-                {
-                    poBlock->DropLock();
-                    continue;
-                }
-                int nBlockId = iX + iY * nBlocksPerRow;
-                if (m_poGDS->m_nPlanarConfig == PLANARCONFIG_SEPARATE)
-                    nBlockId += (nBand - 1) * m_poGDS->m_nBlocksPerBand;
-                vsi_l_offset nOffset = 0;
-                vsi_l_offset nSize = 0;
-
-                if ((m_poGDS->m_nPlanarConfig == PLANARCONFIG_CONTIG ||
-                     m_poGDS->nBands == 1) &&
-                    !m_poGDS->m_bStreamingIn &&
-                    m_poGDS->m_bBlockOrderRowMajor &&
-                    m_poGDS->m_bLeaderSizeAsUInt4)
-                {
-                    OptimizedRetrievalOfOffsetSize(nBlockId, nOffset, nSize,
-                                                   nTotalSize,
-                                                   nMaxRawBlockCacheSize);
-                }
-                else
-                {
-                    CPL_IGNORE_RET_VAL(m_poGDS->IsBlockAvailable(
-                        nBlockId, &nOffset, &nSize, nullptr));
-                }
-                if (nSize)
-                {
-                    if (nTotalSize + nSize < nMaxRawBlockCacheSize)
-                    {
-#ifdef DEBUG_VERBOSE
-                        CPLDebug("GTiff",
-                                 "Precaching for block (%d, %d), " CPL_FRMT_GUIB
-                                 "-" CPL_FRMT_GUIB,
-                                 iX, iY, nOffset,
-                                 nOffset + static_cast<size_t>(nSize) - 1);
-#endif
-                        aOffsetSize.push_back(
-                            std::pair(nOffset, static_cast<size_t>(nSize)));
-                        nTotalSize += static_cast<size_t>(nSize);
-                    }
-                    else
-                    {
-                        bGoOn = false;
-                    }
-                }
-            }
-        }
-
-        std::sort(aOffsetSize.begin(), aOffsetSize.end());
-
-        if (nTotalSize > 0)
-        {
-            pBufferedData = VSI_MALLOC_VERBOSE(nTotalSize);
-            if (pBufferedData)
-            {
-                std::vector<vsi_l_offset> anOffsets;
-                std::vector<size_t> anSizes;
-                std::vector<void *> apData;
-                anOffsets.push_back(aOffsetSize[0].first);
-                apData.push_back(static_cast<GByte *>(pBufferedData));
-                size_t nChunkSize = aOffsetSize[0].second;
-                size_t nAccOffset = 0;
-                // Try to merge contiguous or slightly overlapping ranges
-                for (size_t i = 0; i < aOffsetSize.size() - 1; i++)
-                {
-                    if (aOffsetSize[i].first < aOffsetSize[i + 1].first &&
-                        aOffsetSize[i].first + aOffsetSize[i].second >=
-                            aOffsetSize[i + 1].first)
-                    {
-                        const auto overlap = aOffsetSize[i].first +
-                                             aOffsetSize[i].second -
-                                             aOffsetSize[i + 1].first;
-                        // That should always be the case for well behaved
-                        // TIFF files.
-                        if (aOffsetSize[i + 1].second > overlap)
-                        {
-                            nChunkSize += static_cast<size_t>(
-                                aOffsetSize[i + 1].second - overlap);
-                        }
-                    }
-                    else
-                    {
-                        // terminate current block
-                        anSizes.push_back(nChunkSize);
-#ifdef DEBUG_VERBOSE
-                        CPLDebug("GTiff",
-                                 "Requesting range [" CPL_FRMT_GUIB
-                                 "-" CPL_FRMT_GUIB "]",
-                                 anOffsets.back(),
-                                 anOffsets.back() + anSizes.back() - 1);
-#endif
-                        nAccOffset += nChunkSize;
-                        // start a new range
-                        anOffsets.push_back(aOffsetSize[i + 1].first);
-                        apData.push_back(static_cast<GByte *>(pBufferedData) +
-                                         nAccOffset);
-                        nChunkSize = aOffsetSize[i + 1].second;
-                    }
-                }
-                // terminate last block
-                anSizes.push_back(nChunkSize);
-#ifdef DEBUG_VERBOSE
-                CPLDebug(
-                    "GTiff",
-                    "Requesting range [" CPL_FRMT_GUIB "-" CPL_FRMT_GUIB "]",
-                    anOffsets.back(), anOffsets.back() + anSizes.back() - 1);
-#endif
-
-                VSILFILE *fp = VSI_TIFFGetVSILFile(th);
-
-                // An error in VSIFReadMultiRangeL() will not be critical,
-                // as this method is an optimization, and if it fails,
-                // tile-by-tile data acquisition will be done, so we can
-                // temporary turn failures into warnings.
-                CPLTurnFailureIntoWarning(true);
-                const bool ok =
-                    VSIFReadMultiRangeL(static_cast<int>(anSizes.size()),
-                                        &apData[0], &anOffsets[0], &anSizes[0],
-                                        fp) == 0;
-                CPLTurnFailureIntoWarning(false);
-                if (ok)
-                {
-                    if (!oMapStrileToOffsetByteCount.empty() &&
-                        !FillCacheStrileToOffsetByteCount(anOffsets, anSizes,
-                                                          apData))
-                    {
-                        // Retry without optimization
-                        CPLFree(pBufferedData);
-                        m_poGDS->m_bLeaderSizeAsUInt4 = false;
-                        void *pRet =
-                            CacheMultiRange(nXOff, nYOff, nXSize, nYSize,
-                                            nBufXSize, nBufYSize, psExtraArg);
-                        m_poGDS->m_bLeaderSizeAsUInt4 = true;
-                        return pRet;
-                    }
-
-                    VSI_TIFFSetCachedRanges(
-                        th, static_cast<int>(anSizes.size()), &apData[0],
-                        &anOffsets[0], &anSizes[0]);
-                }
-                else
-                {
-                    CPLFree(pBufferedData);
-                    pBufferedData = nullptr;
-                }
-            }
-        }
-    }
-    return pBufferedData;
-}
-
 /************************************************************************/
 /*                       IGetDataCoverageStatus()                       */
 /************************************************************************/
diff --git a/port/cpl_known_config_options.h b/port/cpl_known_config_options.h
index c866bad5b3b5..6826761b480b 100644
--- a/port/cpl_known_config_options.h
+++ b/port/cpl_known_config_options.h
@@ -324,7 +324,7 @@ constexpr static const char* const apszKnownConfigOptions[] =
    "GDAL_MAX_CONNECTIONS", // from gdalogcapidataset.cpp, gdalwmsdataset.cpp
    "GDAL_MAX_DATASET_POOL_RAM_USAGE", // from gdalproxypool.cpp
    "GDAL_MAX_DATASET_POOL_SIZE", // from gdal_translate_bin.cpp, gdalproxypool.cpp, gdalwarp_bin.cpp
-   "GDAL_MAX_RAW_BLOCK_CACHE_SIZE", // from gtiffrasterband_read.cpp
+   "GDAL_MAX_RAW_BLOCK_CACHE_SIZE", // from gtiffdataset_read.cpp
    "GDAL_MEM_ENABLE_OPEN", // from memdataset.cpp
    "GDAL_NETCDF_ASSUME_LONGLAT", // from netcdfdataset.cpp
    "GDAL_NETCDF_BOTTOMUP", // from netcdfdataset.cpp
diff --git a/swig/python/gdal-utils/osgeo_utils/samples/validate_cloud_optimized_geotiff.py b/swig/python/gdal-utils/osgeo_utils/samples/validate_cloud_optimized_geotiff.py
index 22ccaa6d2179..e2b346b04881 100644
--- a/swig/python/gdal-utils/osgeo_utils/samples/validate_cloud_optimized_geotiff.py
+++ b/swig/python/gdal-utils/osgeo_utils/samples/validate_cloud_optimized_geotiff.py
@@ -38,6 +38,7 @@ class ValidateCloudOptimizedGeoTIFFException(Exception):
 
 def full_check_band(
     f,
+    interleave,
     band_name,
     band,
     errors,
@@ -45,6 +46,7 @@ def full_check_band(
     block_leader_size_as_uint4,
     block_trailer_last_4_bytes_repeated,
     mask_interleaved_with_imagery,
+    last_offset,
 ):
 
     block_size = band.GetBlockSize()
@@ -60,7 +62,6 @@ def full_check_band(
 
     yblocks = (band.YSize + block_size[1] - 1) // block_size[1]
     xblocks = (band.XSize + block_size[0] - 1) // block_size[0]
-    last_offset = 0
     for y in range(yblocks):
         for x in range(xblocks):
 
@@ -98,7 +99,13 @@ def full_check_band(
                                 % (x, y)
                             ]
 
-            if mask_band:
+            if mask_band and (
+                interleave == "PIXEL"
+                or (
+                    interleave == "TILE"
+                    and band.GetBand() == band.GetDataset().RasterCount
+                )
+            ):
                 offset_mask = mask_band.GetMetadataItem(
                     "BLOCK_OFFSET_%d_%d" % (x, y), "TIFF"
                 )
@@ -131,6 +138,33 @@ def full_check_band(
 
             last_offset = offset
 
+    return last_offset
+
+
+def check_tile_interleave(ds, ds_name, block_order_row_major, errors):
+
+    block_size = ds.GetRasterBand(1).GetBlockSize()
+    yblocks = (ds.RasterYSize + block_size[1] - 1) // block_size[1]
+    xblocks = (ds.RasterXSize + block_size[0] - 1) // block_size[0]
+    last_offset = 0
+    for y in range(yblocks):
+        for x in range(xblocks):
+            for band_idx in range(ds.RasterCount):
+                offset = ds.GetRasterBand(band_idx + 1).GetMetadataItem(
+                    "BLOCK_OFFSET_%d_%d" % (x, y), "TIFF"
+                )
+                offset = int(offset) if offset is not None else 0
+
+                if offset > 0:
+                    if block_order_row_major and offset < last_offset:
+                        errors += [
+                            ds_name
+                            + ": offset of block (%d, %d) is smaller than previous block"
+                            % (x, y)
+                        ]
+
+                    last_offset = offset
+
 
 def validate(ds, check_tiled=True, full_check=False):
     """Check if a file is a (Geo)TIFF with cloud optimized compatible structure.
@@ -331,6 +365,8 @@ def get_block_offset(band):
             "should be after the one of the overview of index %d" % (ovr_count - 1)
         ]
 
+    interleave = ds.GetMetadataItem("INTERLEAVE", "IMAGE_STRUCTURE")
+
     if full_check and (
         block_order_row_major
         or block_leader_size_as_uint4
@@ -341,22 +377,49 @@ def get_block_offset(band):
         if not f:
             raise ValidateCloudOptimizedGeoTIFFException("Cannot open file")
 
-        full_check_band(
-            f,
-            "Main resolution image",
-            main_band,
-            errors,
-            block_order_row_major,
-            block_leader_size_as_uint4,
-            block_trailer_last_4_bytes_repeated,
-            mask_interleaved_with_imagery,
-        )
+        if interleave == "PIXEL":
+            full_check_band(
+                f,
+                interleave,
+                "Main resolution image",
+                main_band,
+                errors,
+                block_order_row_major,
+                block_leader_size_as_uint4,
+                block_trailer_last_4_bytes_repeated,
+                mask_interleaved_with_imagery,
+                0,
+            )
+        else:
+            last_offset = 0
+            for band_idx in range(ds.RasterCount):
+                if interleave != "BAND":
+                    last_offset = 0
+                last_offset = full_check_band(
+                    f,
+                    interleave,
+                    "Band %d of main resolution image" % (band_idx + 1),
+                    ds.GetRasterBand(band_idx + 1),
+                    errors,
+                    block_order_row_major,
+                    block_leader_size_as_uint4,
+                    block_trailer_last_4_bytes_repeated,
+                    mask_interleaved_with_imagery,
+                    last_offset,
+                )
+
+        if interleave == "TILE":
+            check_tile_interleave(
+                ds, "Main resolution image", block_order_row_major, errors
+            )
+
         if (
             main_band.GetMaskFlags() == gdal.GMF_PER_DATASET
             and (filename + ".msk") not in ds.GetFileList()
         ):
             full_check_band(
                 f,
+                interleave,
                 "Mask band of main resolution image",
                 main_band.GetMaskBand(),
                 errors,
@@ -364,25 +427,56 @@ def get_block_offset(band):
                 block_leader_size_as_uint4,
                 block_trailer_last_4_bytes_repeated,
                 False,
+                0,
             )
         for i in range(ovr_count):
             ovr_band = ds.GetRasterBand(1).GetOverview(i)
-            full_check_band(
-                f,
-                "Overview %d" % i,
-                ovr_band,
-                errors,
-                block_order_row_major,
-                block_leader_size_as_uint4,
-                block_trailer_last_4_bytes_repeated,
-                mask_interleaved_with_imagery,
-            )
+            if interleave == "PIXEL":
+                full_check_band(
+                    f,
+                    interleave,
+                    "Overview %d" % i,
+                    ovr_band,
+                    errors,
+                    block_order_row_major,
+                    block_leader_size_as_uint4,
+                    block_trailer_last_4_bytes_repeated,
+                    mask_interleaved_with_imagery,
+                    0,
+                )
+            else:
+                last_offset = 0
+                for band_idx in range(ds.RasterCount):
+                    if interleave != "BAND":
+                        last_offset = 0
+                    last_offset = full_check_band(
+                        f,
+                        interleave,
+                        "Band %d of overview %d" % (band_idx + 1, i),
+                        ds.GetRasterBand(band_idx + 1).GetOverview(i),
+                        errors,
+                        block_order_row_major,
+                        block_leader_size_as_uint4,
+                        block_trailer_last_4_bytes_repeated,
+                        mask_interleaved_with_imagery,
+                        last_offset,
+                    )
+
+            if interleave == "TILE":
+                check_tile_interleave(
+                    ds.GetRasterBand(1).GetDataset(),
+                    "Overview %d" % i,
+                    block_order_row_major,
+                    errors,
+                )
+
             if (
                 ovr_band.GetMaskFlags() == gdal.GMF_PER_DATASET
                 and (filename + ".msk") not in ds.GetFileList()
             ):
                 full_check_band(
                     f,
+                    interleave,
                     "Mask band of overview %d" % i,
                     ovr_band.GetMaskBand(),
                     errors,
@@ -390,6 +484,7 @@ def get_block_offset(band):
                     block_leader_size_as_uint4,
                     block_trailer_last_4_bytes_repeated,
                     False,
+                    0,
                 )
         gdal.VSIFCloseL(f)
 

From a3d090ef300c0870ca38bc80be2bb84e98776d94 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 10 Jan 2025 15:50:11 +0100
Subject: [PATCH 2/7] GTiff: no longer allow INTERLEAVE=TILE (only in a secret
 way use by the COG driver)

---
 autotest/gcore/tiff_write.py        |  8 +---
 doc/source/drivers/raster/gtiff.rst | 38 +++--------------
 frmts/gtiff/cogdriver.cpp           | 10 ++++-
 frmts/gtiff/geotiff.cpp             |  1 -
 frmts/gtiff/gtiffdataset_write.cpp  | 65 ++++++++++++++---------------
 5 files changed, 47 insertions(+), 75 deletions(-)

diff --git a/autotest/gcore/tiff_write.py b/autotest/gcore/tiff_write.py
index c589fdf6b200..30cdc78e6e24 100755
--- a/autotest/gcore/tiff_write.py
+++ b/autotest/gcore/tiff_write.py
@@ -11998,18 +11998,12 @@ def test_tiff_write_warn_ignore_predictor_option(tmp_vsimem):
 @pytest.mark.parametrize("COPY_SRC_OVERVIEWS", ["YES", "NO"])
 def test_tiff_write_interleave_tile(tmp_vsimem, COPY_SRC_OVERVIEWS):
     out_filename = str(tmp_vsimem / "out.tif")
-    with pytest.raises(
-        Exception, match="INTERLEAVE=TILE is only supported for CreateCopy"
-    ):
-        gdal.GetDriverByName("GTiff").Create(
-            out_filename, 1, 1, 2, options=["INTERLEAVE=TILE"]
-        )
 
     ds = gdal.GetDriverByName("GTiff").CreateCopy(
         out_filename,
         gdal.Open("data/rgbsmall.tif"),
         options=[
-            "INTERLEAVE=TILE",
+            "@TILE_INTERLEAVE=YES",
             "TILED=YES",
             "BLOCKXSIZE=32",
             "BLOCKYSIZE=32",
diff --git a/doc/source/drivers/raster/gtiff.rst b/doc/source/drivers/raster/gtiff.rst
index c10a625852cc..1d59fd384e46 100644
--- a/doc/source/drivers/raster/gtiff.rst
+++ b/doc/source/drivers/raster/gtiff.rst
@@ -423,7 +423,7 @@ This driver supports the following creation options:
       RPC information is available.
 
 -  .. co:: INTERLEAVE
-      :choices: BAND, PIXEL, TILE
+      :choices: BAND, PIXEL
 
       Set the interleaving to use
 
@@ -481,42 +481,16 @@ This driver supports the following creation options:
           end_for
 
 
-      * ``TILE`` (added in 3.11, only for CreateCopy()): this is a sort of
-        compromise between PIXEL and BAND, using the ``separate`` configuration,
-        but where data for a same spatial block is written for all bands, before
-        the data of the next spatial block.
-        When the block height is 1, this is also known as a
-        ``BIL (Band Interleaved per Line)`` organization.
-
-        Such a layout may be useful for writing hyperspectral datasets (several
-        hundred of bands), to get a compromise between efficient spatial query,
-        and partial band selection.
-
-        Assuming a pixel[band][y][x] indexed array, when using CreateCopy(),
-        the pseudo code for the file disposition is:
-
-        ::
-
-          for y in 0 ... numberOfBlocksInHeight - 1:
-              for x in 0 ... numberOfTilesInWidth - 1:
-                  for band in 0 ... numberBands -1:
-                      start_new_strip_or_tile()
-                      for j in 0 ... blockHeight - 1:
-                          for i in 0 ... blockWidth -1:
-                              write(pixel[band][y*blockHeight+j][x*blockWidth+i])
-                          end_for
-                      end_for
-                      end_new_strip_or_tile()
-                  end_for
-              end_for
-          end_for
-
-
       Starting with GDAL 3.5, when copying from a source dataset with multiple bands
       which advertises a INTERLEAVE metadata item, if the INTERLEAVE creation option
       is not specified, the source dataset INTERLEAVE will be automatically taken
       into account, unless the COMPRESS creation option is specified.
 
+      .. note::
+
+          Starting with GDAL 3.11, a ``TILE`` interleaving is available,
+          but only when using the :ref:`raster.cog` driver.
+
 -  .. co:: TILED
       :choices: YES, NO
       :default: NO
diff --git a/frmts/gtiff/cogdriver.cpp b/frmts/gtiff/cogdriver.cpp
index 73ced0419555..1786798a21b8 100644
--- a/frmts/gtiff/cogdriver.cpp
+++ b/frmts/gtiff/cogdriver.cpp
@@ -1361,7 +1361,15 @@ GDALDataset *GDALCOGCreator::Create(const char *pszFilename,
         aosOptions.AddNameValue("SRC_MDD", *papszSrcMDDIter);
     CSLDestroy(papszSrcMDD);
 
-    aosOptions.SetNameValue("INTERLEAVE", pszInterleave);
+    if (EQUAL(pszInterleave, "TILE"))
+    {
+        aosOptions.SetNameValue("INTERLEAVE", "BAND");
+        aosOptions.SetNameValue("@TILE_INTERLEAVE", "YES");
+    }
+    else
+    {
+        aosOptions.SetNameValue("INTERLEAVE", pszInterleave);
+    }
 
     CPLDebug("COG", "Generating final product: start");
     auto poRet =
diff --git a/frmts/gtiff/geotiff.cpp b/frmts/gtiff/geotiff.cpp
index 0e2f6d3b63f8..225d00dfa304 100644
--- a/frmts/gtiff/geotiff.cpp
+++ b/frmts/gtiff/geotiff.cpp
@@ -1340,7 +1340,6 @@ void GDALRegister_GTiff()
         "   <Option name='INTERLEAVE' type='string-select' default='PIXEL'>"
         "       <Value>BAND</Value>"
         "       <Value>PIXEL</Value>"
-        "       <Value remark='Only for CreateCopy()'>TILE</Value>"
         "   </Option>"
         "   <Option name='TILED' type='boolean' description='Switch to tiled "
         "format'/>"
diff --git a/frmts/gtiff/gtiffdataset_write.cpp b/frmts/gtiff/gtiffdataset_write.cpp
index 5d0dd63ed98f..cf7124482bb3 100644
--- a/frmts/gtiff/gtiffdataset_write.cpp
+++ b/frmts/gtiff/gtiffdataset_write.cpp
@@ -4284,9 +4284,9 @@ bool GTiffDataset::WriteMetadata(GDALDataset *poSrcDS, TIFF *l_hTIFF,
     if (CPLTestBool(
             CPLGetConfigOption("GTIFF_WRITE_IMAGE_STRUCTURE_METADATA", "YES")))
     {
-        const char *pszInterleave =
-            CSLFetchNameValue(papszCreationOptions, "INTERLEAVE");
-        if (pszInterleave && EQUAL(pszInterleave, "TILE"))
+        const char *pszTileInterleave =
+            CSLFetchNameValue(papszCreationOptions, "@TILE_INTERLEAVE");
+        if (pszTileInterleave && CPLTestBool(pszTileInterleave))
         {
             AppendMetadataItem(&psRoot, &psTail, "INTERLEAVE", "TILE", 0,
                                nullptr, "IMAGE_STRUCTURE");
@@ -5121,47 +5121,44 @@ TIFF *GTiffDataset::CreateLL(const char *pszFilename, int nXSize, int nYSize,
     }
 
     int nPlanar = 0;
-    pszValue = CSLFetchNameValue(papszParamList, "INTERLEAVE");
-    if (pszValue != nullptr)
+
+    // Hidden @TILE_INTERLEAVE=YES parameter used by the COG driver
+    if (bCreateCopy &&
+        (pszValue = CSLFetchNameValue(papszParamList, "@TILE_INTERLEAVE")) &&
+        CPLTestBool(pszValue))
     {
-        if (EQUAL(pszValue, "PIXEL"))
-            nPlanar = PLANARCONFIG_CONTIG;
-        else if (EQUAL(pszValue, "BAND"))
-        {
-            nPlanar = PLANARCONFIG_SEPARATE;
-        }
-        else if (EQUAL(pszValue, "BAND"))
-        {
-            nPlanar = PLANARCONFIG_SEPARATE;
-        }
-        else if (EQUAL(pszValue, "TILE"))
+        bTileInterleavingOut = true;
+        nPlanar = PLANARCONFIG_SEPARATE;
+    }
+    else
+    {
+        pszValue = CSLFetchNameValue(papszParamList, "INTERLEAVE");
+        if (pszValue != nullptr)
         {
-            bTileInterleavingOut = true;
-            nPlanar = PLANARCONFIG_SEPARATE;
-            if (!bCreateCopy)
+            if (EQUAL(pszValue, "PIXEL"))
+                nPlanar = PLANARCONFIG_CONTIG;
+            else if (EQUAL(pszValue, "BAND"))
+            {
+                nPlanar = PLANARCONFIG_SEPARATE;
+            }
+            else if (EQUAL(pszValue, "BAND"))
+            {
+                nPlanar = PLANARCONFIG_SEPARATE;
+            }
+            else
             {
-                CPLError(CE_Failure, CPLE_NotSupported,
-                         "INTERLEAVE=TILE is only supported for CreateCopy(), "
-                         "not Create()");
+                ReportError(
+                    pszFilename, CE_Failure, CPLE_IllegalArg,
+                    "INTERLEAVE=%s unsupported, value must be PIXEL or BAND.",
+                    pszValue);
                 return nullptr;
             }
         }
         else
         {
-            ReportError(
-                pszFilename, CE_Failure, CPLE_IllegalArg,
-                bCreateCopy
-                    ? "INTERLEAVE=%s unsupported, value must be PIXEL, BAND or "
-                      "TILE."
-                    : "INTERLEAVE=%s unsupported, value must be PIXEL or BAND.",
-                pszValue);
-            return nullptr;
+            nPlanar = PLANARCONFIG_CONTIG;
         }
     }
-    else
-    {
-        nPlanar = PLANARCONFIG_CONTIG;
-    }
 
     int l_nCompression = COMPRESSION_NONE;
     pszValue = CSLFetchNameValue(papszParamList, "COMPRESS");

From 8221e1287712265e5630a24707fb9aeb1866b414 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 17 Jan 2025 16:58:47 +0100
Subject: [PATCH 3/7] DXF: do not error out on INSERT blocks with row count or
 column count equal to 0

Fixes #11591
---
 .../data/dxf/insert_only_col_count_zero.dxf   | 200 ++++++++++++++++++
 autotest/ogr/ogr_dxf.py                       |  21 ++
 ogr/ogrsf_frmts/dxf/ogrdxfdatasource.cpp      |   2 +
 ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp           |  10 +-
 4 files changed, 231 insertions(+), 2 deletions(-)
 create mode 100644 autotest/ogr/data/dxf/insert_only_col_count_zero.dxf

diff --git a/autotest/ogr/data/dxf/insert_only_col_count_zero.dxf b/autotest/ogr/data/dxf/insert_only_col_count_zero.dxf
new file mode 100644
index 000000000000..04b8de6fda85
--- /dev/null
+++ b/autotest/ogr/data/dxf/insert_only_col_count_zero.dxf
@@ -0,0 +1,200 @@
+  0
+SECTION
+  2
+TABLES
+  0
+ENDSEC
+  0
+SECTION
+  2
+BLOCKS
+  0
+BLOCK
+  5
+44
+100
+AcDbEntity
+  8
+0
+100
+AcDbBlockBegin
+  2
+STAR
+ 70
+0
+ 10
+0.0
+ 20
+0.0
+ 30
+0.0
+  3
+STAR
+  1
+
+  0
+LINE
+  5
+45
+100
+AcDbEntity
+  8
+0
+ 62
+256
+100
+AcDbLine
+ 10
+-0.0281474976710656
+ 20
+1.0414574138294284
+ 30
+0.0
+ 11
+0.6192449487634439
+ 21
+-1.0696049115004942
+ 31
+0.0
+  0
+LINE
+  5
+46
+100
+AcDbEntity
+  8
+0
+ 62
+256
+100
+AcDbLine
+ 10
+0.6192449487634439
+ 20
+-1.0696049115004942
+ 30
+0.0
+ 11
+-0.9570149208162315
+ 21
+0.4785074604081158
+ 31
+0.0
+  0
+LINE
+  5
+47
+100
+AcDbEntity
+  8
+0
+ 62
+256
+100
+AcDbLine
+ 10
+-0.9570149208162315
+ 20
+0.4785074604081158
+ 30
+0.0
+ 11
+1.0414574138294284
+ 21
+0.3659174697238533
+ 31
+0.0
+  0
+LINE
+  5
+48
+100
+AcDbEntity
+  8
+0
+ 62
+256
+100
+AcDbLine
+ 10
+1.0414574138294284
+ 20
+0.3659174697238533
+ 30
+0.0
+ 11
+-0.4785074604081158
+ 21
+-1.0414574138294284
+ 31
+0.0
+  0
+LINE
+  5
+49
+100
+AcDbEntity
+  8
+0
+ 62
+256
+100
+AcDbLine
+ 10
+-0.4785074604081158
+ 20
+-1.0414574138294284
+ 30
+0.0
+ 11
+-0.0562949953421313
+ 21
+1.0133099161583627
+ 31
+0.0
+  0
+ENDBLK
+  5
+4A
+100
+AcDbEntity
+  8
+0
+100
+AcDbBlockEnd
+  0
+ENDSEC
+  0
+SECTION
+  2
+ENTITIES
+  0
+INSERT
+  5
+55
+100
+AcDbEntity
+  8
+0
+ 62
+256
+100
+AcDbBlockReference
+  2
+STAR
+ 10
+79.0976537766561876
+ 20
+119.9621950624433424
+ 30
+0.0
+ 44
+1.0
+ 45
+1.0
+ 70
+0
+  0
+ENDSEC
+  0
+EOF
diff --git a/autotest/ogr/ogr_dxf.py b/autotest/ogr/ogr_dxf.py
index 16cd74c362f8..9932e5288454 100644
--- a/autotest/ogr/ogr_dxf.py
+++ b/autotest/ogr/ogr_dxf.py
@@ -4124,3 +4124,24 @@ def test_ogr_dxf_write_MEASUREMENT(tmp_vsimem):
             pass
     with ogr.Open(filename) as ds:
         assert ds.GetMetadataItem("$MEASUREMENT", "DXF_HEADER_VARIABLES") == "     0"
+
+
+###############################################################################
+# Use case of https://github.com/OSGeo/gdal/issues/11591
+# Test reading a INSERT block whose column count is zero.
+# Not totally sure about the exact behavior we should do. Currently
+# the zero count is strictly honored and no geometry from the block will be
+# inserted into the regular geometries. LibreCAD 2.1.3 does the same thing
+
+
+@gdaltest.enable_exceptions()
+def test_ogr_dxf_insert_col_count_zero():
+
+    with ogr.Open("data/dxf/insert_only_col_count_zero.dxf") as ds:
+        lyr = ds.GetLayer(0)
+        assert lyr.GetFeatureCount() == 0
+
+    with gdal.config_option("DXF_INLINE_BLOCKS", "NO"):
+        with ogr.Open("data/dxf/insert_only_col_count_zero.dxf") as ds:
+            lyr = ds.GetLayerByName("blocks")
+            assert lyr.GetFeatureCount() == 1
diff --git a/ogr/ogrsf_frmts/dxf/ogrdxfdatasource.cpp b/ogr/ogrsf_frmts/dxf/ogrdxfdatasource.cpp
index 5c3d26b5e0df..cbdb9dd6173d 100644
--- a/ogr/ogrsf_frmts/dxf/ogrdxfdatasource.cpp
+++ b/ogr/ogrsf_frmts/dxf/ogrdxfdatasource.cpp
@@ -88,6 +88,8 @@ int OGRDXFDataSource::Open(const char *pszFilename, bool bHeaderOnly,
                            CSLConstList papszOptionsIn)
 
 {
+    SetDescription(pszFilename);
+
     osEncoding = CPL_ENC_ISO8859_1;
 
     bInlineBlocks = CPLTestBool(
diff --git a/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp b/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp
index 719d3bb85324..64aa386f81bc 100644
--- a/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp
+++ b/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp
@@ -3167,7 +3167,7 @@ bool OGRDXFLayer::TranslateINSERT()
 
             case 70:
                 m_oInsertState.m_nColumnCount = atoi(szLineBuf);
-                if (m_oInsertState.m_nColumnCount <= 0)
+                if (m_oInsertState.m_nColumnCount < 0)
                 {
                     DXF_LAYER_READER_ERROR();
                     m_oInsertState.m_nRowCount = 0;
@@ -3178,7 +3178,7 @@ bool OGRDXFLayer::TranslateINSERT()
 
             case 71:
                 m_oInsertState.m_nRowCount = atoi(szLineBuf);
-                if (m_oInsertState.m_nRowCount <= 0)
+                if (m_oInsertState.m_nRowCount < 0)
                 {
                     DXF_LAYER_READER_ERROR();
                     m_oInsertState.m_nRowCount = 0;
@@ -3205,6 +3205,12 @@ bool OGRDXFLayer::TranslateINSERT()
         return false;
     }
 
+    if (m_oInsertState.m_nRowCount == 0 || m_oInsertState.m_nColumnCount == 0)
+    {
+        m_oInsertState.m_nRowCount = 0;
+        m_oInsertState.m_nColumnCount = 0;
+    }
+
     /* -------------------------------------------------------------------- */
     /*      Process any attribute entities.                                 */
     /* -------------------------------------------------------------------- */

From 1326594e6d097d54c79be16bd611b7f2f7e7f75b Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 17 Jan 2025 01:07:45 +0100
Subject: [PATCH 4/7] Add OGRPolygon::OGRPolygon(double x1, double y1, double
 x2, double y2)

---
 autotest/cpp/test_ogr.cpp |  9 +++++++++
 ogr/ogr_geometry.h        |  2 ++
 ogr/ogrpolygon.cpp        | 21 +++++++++++++++++++++
 3 files changed, 32 insertions(+)

diff --git a/autotest/cpp/test_ogr.cpp b/autotest/cpp/test_ogr.cpp
index 3faf8f0791b9..0ce901592594 100644
--- a/autotest/cpp/test_ogr.cpp
+++ b/autotest/cpp/test_ogr.cpp
@@ -4628,4 +4628,13 @@ TEST_F(test_ogr, GetDefaultArcStepSize)
     }
 }
 
+TEST_F(test_ogr, OGRPolygon_two_vertex_constructor)
+{
+    OGRPolygon p(1, 2, 3, 4);
+    char *outWKT = nullptr;
+    p.exportToWkt(&outWKT, wkbVariantIso);
+    EXPECT_STREQ(outWKT, "POLYGON ((1 2,1 4,3 4,3 2,1 2))");
+    CPLFree(outWKT);
+}
+
 }  // namespace
diff --git a/ogr/ogr_geometry.h b/ogr/ogr_geometry.h
index e2227f3443f0..ddfdec51bd2a 100644
--- a/ogr/ogr_geometry.h
+++ b/ogr/ogr_geometry.h
@@ -2716,6 +2716,8 @@ class CPL_DLL OGRPolygon : public OGRCurvePolygon
     /** Create an empty polygon. */
     OGRPolygon() = default;
 
+    OGRPolygon(double x1, double y1, double x2, double y2);
+
     OGRPolygon(const OGRPolygon &other);
     /** Move constructor */
     OGRPolygon(OGRPolygon &&other) = default;
diff --git a/ogr/ogrpolygon.cpp b/ogr/ogrpolygon.cpp
index 1daceec6e6c4..084ec0617859 100644
--- a/ogr/ogrpolygon.cpp
+++ b/ogr/ogrpolygon.cpp
@@ -27,6 +27,27 @@
 #include "ogr_sfcgal.h"
 #include "ogr_p.h"
 
+/************************************************************************/
+/*        OGRPolygon(double x1, double y1, double x2, double y2)        */
+/************************************************************************/
+
+/**
+ * \brief Construct a rectangular polygon from opposite corner coordinates.
+ *
+ * @since GDAL 3.11
+ */
+
+OGRPolygon::OGRPolygon(double x1, double y1, double x2, double y2)
+{
+    auto poLR = std::make_unique<OGRLinearRing>();
+    poLR->addPoint(x1, y1);
+    poLR->addPoint(x1, y2);
+    poLR->addPoint(x2, y2);
+    poLR->addPoint(x2, y1);
+    poLR->addPoint(x1, y1);
+    addRingDirectly(poLR.release());
+}
+
 /************************************************************************/
 /*                     OGRPolygon( const OGRPolygon& )                  */
 /************************************************************************/

From fa9f8038ae6ae8de4358a8426705ecf74a22add5 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Mon, 20 Jan 2025 22:16:07 +0100
Subject: [PATCH 5/7] gdal_raster_clip.rst: improve description of clip

---
 doc/source/programs/gdal_raster_clip.rst | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/doc/source/programs/gdal_raster_clip.rst b/doc/source/programs/gdal_raster_clip.rst
index 8a79c0c67418..f7eddd87c1ce 100644
--- a/doc/source/programs/gdal_raster_clip.rst
+++ b/doc/source/programs/gdal_raster_clip.rst
@@ -89,7 +89,10 @@ Standard options
 
 .. option:: --like <DATASET>
 
-    Raster dataset to use as a template for bounds. This option is mutually
+    Raster dataset to use as a template for bounds, forming a rectangular shape
+    following the geotransformation matrix (and thus potentially including
+    nodata collar).
+    This option is mutually
     exclusive with :option:`--bbox` and :option:`--bbox-crs`.
 
 Advanced options

From 2b585ffc3aae90532a81efa8ad92c70a91ac91e5 Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Fri, 17 Jan 2025 01:44:16 +0100
Subject: [PATCH 6/7] Add 'gdal vector clip' and 'gdal vector pipeline read ...
 ! clip ... ! write ...'

---
 apps/CMakeLists.txt                           |   1 +
 apps/gdalalg_vector.cpp                       |   2 +
 apps/gdalalg_vector_clip.cpp                  | 482 +++++++++
 apps/gdalalg_vector_clip.h                    |  67 ++
 apps/gdalalg_vector_pipeline.cpp              |   2 +
 .../utilities/test_gdalalg_vector_clip.py     | 924 ++++++++++++++++++
 doc/source/conf.py                            |   7 +
 .../programs/gdal_options/co_vector.rst       |  18 +
 doc/source/programs/gdal_vector.rst           |   4 +-
 doc/source/programs/gdal_vector_clip.rst      | 170 ++++
 doc/source/programs/gdal_vector_convert.rst   |   2 +-
 doc/source/programs/gdal_vector_pipeline.rst  |  22 +
 doc/source/programs/index.rst                 |   2 +
 13 files changed, 1701 insertions(+), 2 deletions(-)
 create mode 100644 apps/gdalalg_vector_clip.cpp
 create mode 100644 apps/gdalalg_vector_clip.h
 create mode 100755 autotest/utilities/test_gdalalg_vector_clip.py
 create mode 100644 doc/source/programs/gdal_options/co_vector.rst
 create mode 100644 doc/source/programs/gdal_vector_clip.rst

diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index a080179aa656..a89f8c456f15 100644
--- a/apps/CMakeLists.txt
+++ b/apps/CMakeLists.txt
@@ -24,6 +24,7 @@ add_library(
   gdalalg_raster_write.cpp
   gdalalg_vector.cpp
   gdalalg_vector_info.cpp
+  gdalalg_vector_clip.cpp
   gdalalg_vector_convert.cpp
   gdalalg_vector_pipeline.cpp
   gdalalg_vector_read.cpp
diff --git a/apps/gdalalg_vector.cpp b/apps/gdalalg_vector.cpp
index 41a15da845ef..75191f2c01ba 100644
--- a/apps/gdalalg_vector.cpp
+++ b/apps/gdalalg_vector.cpp
@@ -13,6 +13,7 @@
 #include "gdalalgorithm.h"
 
 #include "gdalalg_vector_info.h"
+#include "gdalalg_vector_clip.h"
 #include "gdalalg_vector_convert.h"
 #include "gdalalg_vector_pipeline.h"
 #include "gdalalg_vector_filter.h"
@@ -37,6 +38,7 @@ class GDALVectorAlgorithm final : public GDALAlgorithm
     GDALVectorAlgorithm() : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
     {
         RegisterSubAlgorithm<GDALVectorInfoAlgorithm>();
+        RegisterSubAlgorithm<GDALVectorClipAlgorithmStandalone>();
         RegisterSubAlgorithm<GDALVectorConvertAlgorithm>();
         RegisterSubAlgorithm<GDALVectorPipelineAlgorithm>();
         RegisterSubAlgorithm<GDALVectorFilterAlgorithmStandalone>();
diff --git a/apps/gdalalg_vector_clip.cpp b/apps/gdalalg_vector_clip.cpp
new file mode 100644
index 000000000000..8e529f1fda0d
--- /dev/null
+++ b/apps/gdalalg_vector_clip.cpp
@@ -0,0 +1,482 @@
+/******************************************************************************
+ *
+ * Project:  GDAL
+ * Purpose:  "clip" step of "vector pipeline", or "gdal vector clip" standalone
+ * Author:   Even Rouault <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#include "gdalalg_vector_clip.h"
+
+#include "gdal_priv.h"
+#include "gdal_utils.h"
+#include "ogrsf_frmts.h"
+
+#include <algorithm>
+
+//! @cond Doxygen_Suppress
+
+#ifndef _
+#define _(x) (x)
+#endif
+
+/************************************************************************/
+/*          GDALVectorClipAlgorithm::GDALVectorClipAlgorithm()          */
+/************************************************************************/
+
+GDALVectorClipAlgorithm::GDALVectorClipAlgorithm(bool standaloneStep)
+    : GDALVectorPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
+                                      standaloneStep)
+{
+    AddBBOXArg(&m_bbox, _("Clipping bounding box as xmin,ymin,xmax,ymax"))
+        .SetMutualExclusionGroup("bbox-geometry-like");
+    AddArg("bbox-crs", 0, _("CRS of clipping bounding box"), &m_bboxCrs)
+        .SetIsCRSArg()
+        .AddHiddenAlias("bbox_srs");
+    AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
+        .SetMutualExclusionGroup("bbox-geometry-like");
+    AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
+        .SetIsCRSArg()
+        .AddHiddenAlias("geometry_srs");
+    AddArg("like", 0, _("Dataset to use as a template for bounds"),
+           &m_likeDataset, GDAL_OF_RASTER | GDAL_OF_VECTOR)
+        .SetMetaVar("DATASET")
+        .SetMutualExclusionGroup("bbox-geometry-like");
+    AddArg("like-sql", 0, ("SELECT statement to run on the 'like' dataset"),
+           &m_likeSQL)
+        .SetMetaVar("SELECT-STATEMENT")
+        .SetMutualExclusionGroup("sql-where");
+    AddArg("like-layer", 0, ("Name of the layer of the 'like' dataset"),
+           &m_likeLayer)
+        .SetMetaVar("LAYER-NAME");
+    AddArg("like-where", 0, ("WHERE SQL clause to run on the 'like' dataset"),
+           &m_likeWhere)
+        .SetMetaVar("WHERE-EXPRESSION")
+        .SetMutualExclusionGroup("sql-where");
+}
+
+/************************************************************************/
+/*                   GDALVectorClipAlgorithmDataset                     */
+/************************************************************************/
+
+namespace
+{
+class GDALVectorClipAlgorithmDataset final : public GDALDataset
+{
+    std::vector<std::unique_ptr<OGRLayer>> m_layers{};
+
+  public:
+    GDALVectorClipAlgorithmDataset() = default;
+
+    void AddLayer(std::unique_ptr<OGRLayer> poLayer)
+    {
+        m_layers.push_back(std::move(poLayer));
+    }
+
+    int GetLayerCount() override
+    {
+        return static_cast<int>(m_layers.size());
+    }
+
+    OGRLayer *GetLayer(int idx) override
+    {
+        return idx >= 0 && idx < GetLayerCount() ? m_layers[idx].get()
+                                                 : nullptr;
+    }
+};
+
+class GDALVectorClipAlgorithmLayer final : public OGRLayer
+{
+  public:
+    GDALVectorClipAlgorithmLayer(OGRLayer *poSrcLayer,
+                                 std::unique_ptr<OGRGeometry> poClipGeom)
+        : m_poSrcLayer(poSrcLayer), m_poClipGeom(std::move(poClipGeom)),
+          m_eSrcLayerGeomType(m_poSrcLayer->GetGeomType()),
+          m_eFlattenSrcLayerGeomType(wkbFlatten(m_eSrcLayerGeomType)),
+          m_bSrcLayerGeomTypeIsCollection(OGR_GT_IsSubClassOf(
+              m_eFlattenSrcLayerGeomType, wkbGeometryCollection))
+    {
+        SetDescription(poSrcLayer->GetDescription());
+        poSrcLayer->SetSpatialFilter(m_poClipGeom.get());
+    }
+
+    OGRFeatureDefn *GetLayerDefn() override
+    {
+        return m_poSrcLayer->GetLayerDefn();
+    }
+
+    void ResetReading() override
+    {
+        m_poSrcLayer->ResetReading();
+        m_poSrcFeature.reset();
+        m_poCurGeomColl.reset();
+        m_idxInCurGeomColl = 0;
+    }
+
+    OGRFeature *GetNextFeature() override
+    {
+        if (m_poSrcFeature && m_poCurGeomColl)
+        {
+            while (m_idxInCurGeomColl < m_poCurGeomColl->getNumGeometries())
+            {
+                const auto poGeom =
+                    m_poCurGeomColl->getGeometryRef(m_idxInCurGeomColl);
+                ++m_idxInCurGeomColl;
+                if (m_eFlattenSrcLayerGeomType == wkbUnknown ||
+                    m_eFlattenSrcLayerGeomType ==
+                        wkbFlatten(poGeom->getGeometryType()))
+                {
+                    auto poDstFeature =
+                        std::unique_ptr<OGRFeature>(m_poSrcFeature->Clone());
+                    poDstFeature->SetGeometry(poGeom);
+                    return poDstFeature.release();
+                }
+            }
+            m_poSrcFeature.reset();
+            m_poCurGeomColl.reset();
+            m_idxInCurGeomColl = 0;
+        }
+
+        while (auto poFeature =
+                   std::unique_ptr<OGRFeature>(m_poSrcLayer->GetNextFeature()))
+        {
+            auto poGeom = poFeature->GetGeometryRef();
+            if (!poGeom)
+                continue;
+
+            auto poIntersection = std::unique_ptr<OGRGeometry>(
+                poGeom->Intersection(m_poClipGeom.get()));
+            if (!poIntersection)
+                continue;
+
+            const auto eFeatGeomType =
+                wkbFlatten(poIntersection->getGeometryType());
+            if (m_eFlattenSrcLayerGeomType != wkbUnknown &&
+                m_eFlattenSrcLayerGeomType != eFeatGeomType)
+            {
+                // If the intersection is a collection of geometry and the
+                // layer geometry type is of non-collection type, create
+                // one feature per element of the collection.
+                if (!m_bSrcLayerGeomTypeIsCollection &&
+                    OGR_GT_IsSubClassOf(eFeatGeomType, wkbGeometryCollection))
+                {
+                    m_poSrcFeature = std::move(poFeature);
+                    m_poCurGeomColl.reset(
+                        poIntersection.release()->toGeometryCollection());
+                    m_idxInCurGeomColl = 0;
+                    return GetNextFeature();
+                }
+                else if (OGR_GT_GetCollection(eFeatGeomType) ==
+                         m_eFlattenSrcLayerGeomType)
+                {
+                    poIntersection.reset(OGRGeometryFactory::forceTo(
+                        poIntersection.release(), m_eSrcLayerGeomType));
+                    poFeature->SetGeometryDirectly(poIntersection.release());
+                    return poFeature.release();
+                }
+                // else discard geometries of incompatible type with the
+                // layer geometry type
+            }
+            else
+            {
+                poFeature->SetGeometryDirectly(poIntersection.release());
+                return poFeature.release();
+            }
+        }
+        return nullptr;
+    }
+
+    int TestCapability(const char *pszCap) override
+    {
+        if (EQUAL(pszCap, OLCStringsAsUTF8) ||
+            EQUAL(pszCap, OLCCurveGeometries) || EQUAL(pszCap, OLCZGeometries))
+            return m_poSrcLayer->TestCapability(pszCap);
+        return false;
+    }
+
+  private:
+    OGRLayer *m_poSrcLayer = nullptr;
+    std::unique_ptr<OGRGeometry> m_poClipGeom{};
+    const OGRwkbGeometryType m_eSrcLayerGeomType;
+    const OGRwkbGeometryType m_eFlattenSrcLayerGeomType;
+    const bool m_bSrcLayerGeomTypeIsCollection;
+    std::unique_ptr<OGRFeature> m_poSrcFeature{};
+    std::unique_ptr<OGRGeometryCollection> m_poCurGeomColl{};
+    int m_idxInCurGeomColl = 0;
+
+    CPL_DISALLOW_COPY_ASSIGN(GDALVectorClipAlgorithmLayer)
+};
+
+}  // namespace
+
+/************************************************************************/
+/*                           LoadGeometry()                             */
+/************************************************************************/
+
+static std::unique_ptr<OGRGeometry> LoadGeometry(GDALDataset *poDS,
+                                                 const std::string &osSQL,
+                                                 const std::string &osLyr,
+                                                 const std::string &osWhere)
+{
+    OGRLayer *poLyr = nullptr;
+    if (!osSQL.empty())
+        poLyr = poDS->ExecuteSQL(osSQL.c_str(), nullptr, nullptr);
+    else if (!osLyr.empty())
+        poLyr = poDS->GetLayerByName(osLyr.c_str());
+    else
+        poLyr = poDS->GetLayer(0);
+
+    if (poLyr == nullptr)
+    {
+        CPLError(CE_Failure, CPLE_AppDefined,
+                 "Failed to identify source layer from clipping dataset.");
+        return nullptr;
+    }
+
+    if (!osWhere.empty())
+        poLyr->SetAttributeFilter(osWhere.c_str());
+
+    OGRGeometryCollection oGC;
+
+    const auto poSRSSrc = poLyr->GetSpatialRef();
+    if (poSRSSrc)
+    {
+        auto poSRSClone = poSRSSrc->Clone();
+        oGC.assignSpatialReference(poSRSClone);
+        poSRSClone->Release();
+    }
+
+    for (auto &poFeat : poLyr)
+    {
+        auto poSrcGeom = std::unique_ptr<OGRGeometry>(poFeat->StealGeometry());
+        if (poSrcGeom)
+        {
+            // Only take into account areal geometries.
+            if (poSrcGeom->getDimension() == 2)
+            {
+                if (!poSrcGeom->IsValid())
+                {
+                    CPLError(CE_Failure, CPLE_AppDefined,
+                             "Geometry of feature " CPL_FRMT_GIB " of %s "
+                             "is invalid.",
+                             poFeat->GetFID(), poDS->GetDescription());
+                    return nullptr;
+                }
+                else
+                {
+                    oGC.addGeometry(std::move(poSrcGeom));
+                }
+            }
+        }
+    }
+
+    if (!osSQL.empty())
+        poDS->ReleaseResultSet(poLyr);
+
+    if (oGC.IsEmpty())
+    {
+        CPLError(CE_Failure, CPLE_AppDefined, "No clipping geometry found");
+        return nullptr;
+    }
+
+    return std::unique_ptr<OGRGeometry>(oGC.UnaryUnion());
+}
+
+/************************************************************************/
+/*                 GDALVectorClipAlgorithm::RunStep()                   */
+/************************************************************************/
+
+bool GDALVectorClipAlgorithm::RunStep(GDALProgressFunc, void *)
+{
+    CPLAssert(m_inputDataset.GetDatasetRef());
+    CPLAssert(m_outputDataset.GetName().empty());
+    CPLAssert(!m_outputDataset.GetDatasetRef());
+
+    auto poSrcDS = m_inputDataset.GetDatasetRef();
+
+    std::unique_ptr<OGRGeometry> poClipGeom;
+
+    const int nLayerCount = poSrcDS->GetLayerCount();
+    bool bSrcLayerHasSRS = false;
+    for (int i = 0; i < nLayerCount; ++i)
+    {
+        auto poSrcLayer = poSrcDS->GetLayer(i);
+        if (poSrcLayer && poSrcLayer->GetSpatialRef())
+        {
+            bSrcLayerHasSRS = true;
+            break;
+        }
+    }
+
+    if (!m_bbox.empty())
+    {
+        poClipGeom = std::make_unique<OGRPolygon>(m_bbox[0], m_bbox[1],
+                                                  m_bbox[2], m_bbox[3]);
+
+        if (!m_bboxCrs.empty())
+        {
+            auto poSRS = new OGRSpatialReference();
+            poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+            CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_bboxCrs.c_str()));
+            poClipGeom->assignSpatialReference(poSRS);
+            poSRS->Release();
+        }
+    }
+    else if (!m_geometry.empty())
+    {
+        {
+            CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
+            auto [poGeom, eErr] =
+                OGRGeometryFactory::createFromWkt(m_geometry.c_str());
+            if (eErr == OGRERR_NONE)
+            {
+                poClipGeom = std::move(poGeom);
+            }
+            else
+            {
+                poClipGeom.reset(
+                    OGRGeometryFactory::createFromGeoJson(m_geometry.c_str()));
+                if (poClipGeom)
+                {
+                    auto poSRS = new OGRSpatialReference();
+                    poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+                    CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput("WGS84"));
+                    poClipGeom->assignSpatialReference(poSRS);
+                    poSRS->Release();
+                }
+            }
+        }
+        if (!poClipGeom)
+        {
+            ReportError(
+                CE_Failure, CPLE_AppDefined,
+                "Clipping geometry is neither a valid WKT or GeoJSON geometry");
+            return false;
+        }
+
+        if (!m_geometryCrs.empty())
+        {
+            auto poSRS = new OGRSpatialReference();
+            poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+            CPL_IGNORE_RET_VAL(poSRS->SetFromUserInput(m_geometryCrs.c_str()));
+            poClipGeom->assignSpatialReference(poSRS);
+            poSRS->Release();
+        }
+    }
+    else if (auto poLikeDS = m_likeDataset.GetDatasetRef())
+    {
+        if (poLikeDS->GetLayerCount() > 1 && m_likeLayer.empty() &&
+            m_likeSQL.empty())
+        {
+            ReportError(
+                CE_Failure, CPLE_AppDefined,
+                "Only single layer dataset can be specified with --like when "
+                "neither --like-layer or --like-sql have been specified");
+            return false;
+        }
+        else if (poLikeDS->GetLayerCount() > 0)
+        {
+            poClipGeom =
+                LoadGeometry(poLikeDS, m_likeSQL, m_likeLayer, m_likeWhere);
+            if (!poClipGeom)
+                return false;
+        }
+        else if (poLikeDS->GetRasterCount() > 0)
+        {
+            double adfGT[6];
+            if (poLikeDS->GetGeoTransform(adfGT) != CE_None)
+            {
+                ReportError(
+                    CE_Failure, CPLE_AppDefined,
+                    "Dataset '%s' has no geotransform matrix. Its bounds "
+                    "cannot be established.",
+                    poLikeDS->GetDescription());
+                return false;
+            }
+            auto poLikeSRS = poLikeDS->GetSpatialRef();
+            if (bSrcLayerHasSRS && !poLikeSRS)
+            {
+                ReportError(CE_Warning, CPLE_AppDefined,
+                            "Dataset '%s' has no SRS. Assuming its SRS is the "
+                            "same as the input vector.",
+                            poLikeDS->GetDescription());
+            }
+            const double dfTLX = adfGT[0];
+            const double dfTLY = adfGT[3];
+            const double dfTRX =
+                adfGT[0] + poLikeDS->GetRasterXSize() * adfGT[1];
+            const double dfTRY =
+                adfGT[3] + poLikeDS->GetRasterXSize() * adfGT[4];
+            const double dfBLX =
+                adfGT[0] + poLikeDS->GetRasterYSize() * adfGT[2];
+            const double dfBLY =
+                adfGT[3] + poLikeDS->GetRasterYSize() * adfGT[5];
+            const double dfBRX = adfGT[0] +
+                                 poLikeDS->GetRasterXSize() * adfGT[1] +
+                                 poLikeDS->GetRasterYSize() * adfGT[2];
+            const double dfBRY = adfGT[3] +
+                                 poLikeDS->GetRasterXSize() * adfGT[4] +
+                                 poLikeDS->GetRasterYSize() * adfGT[5];
+
+            auto poPoly = std::make_unique<OGRPolygon>();
+            auto poLR = std::make_unique<OGRLinearRing>();
+            poLR->addPoint(dfTLX, dfTLY);
+            poLR->addPoint(dfTRX, dfTRY);
+            poLR->addPoint(dfBRX, dfBRY);
+            poLR->addPoint(dfBLX, dfBLY);
+            poLR->addPoint(dfTLX, dfTLY);
+            poPoly->addRingDirectly(poLR.release());
+            poPoly->assignSpatialReference(poLikeSRS);
+            poClipGeom = std::move(poPoly);
+        }
+        else
+        {
+            ReportError(CE_Failure, CPLE_AppDefined,
+                        "Cannot get extent from clip dataset");
+            return false;
+        }
+    }
+    else
+    {
+        ReportError(CE_Failure, CPLE_AppDefined,
+                    "--bbox, --geometry or --like must be specified");
+        return false;
+    }
+
+    auto outDS = std::make_unique<GDALVectorClipAlgorithmDataset>();
+    outDS->SetDescription(poSrcDS->GetDescription());
+
+    bool ret = true;
+    for (int i = 0; ret && i < nLayerCount; ++i)
+    {
+        auto poSrcLayer = poSrcDS->GetLayer(i);
+        ret = (poSrcLayer != nullptr);
+        if (ret)
+        {
+            auto poClipGeomForLayer =
+                std::unique_ptr<OGRGeometry>(poClipGeom->clone());
+            if (poClipGeomForLayer->getSpatialReference() &&
+                poSrcLayer->GetSpatialRef())
+            {
+                ret = poClipGeomForLayer->transformTo(
+                          poSrcLayer->GetSpatialRef()) == OGRERR_NONE;
+            }
+            if (ret)
+            {
+                outDS->AddLayer(std::make_unique<GDALVectorClipAlgorithmLayer>(
+                    poSrcLayer, std::move(poClipGeomForLayer)));
+            }
+        }
+    }
+
+    if (ret)
+        m_outputDataset.Set(std::move(outDS));
+
+    return ret;
+}
+
+//! @endcond
diff --git a/apps/gdalalg_vector_clip.h b/apps/gdalalg_vector_clip.h
new file mode 100644
index 000000000000..edc30c0e5fd1
--- /dev/null
+++ b/apps/gdalalg_vector_clip.h
@@ -0,0 +1,67 @@
+/******************************************************************************
+ *
+ * Project:  GDAL
+ * Purpose:  "clip" step of "vector pipeline", or "gdal vector clip" standalone
+ * Author:   Even Rouault <even dot rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
+ *
+ * SPDX-License-Identifier: MIT
+ ****************************************************************************/
+
+#ifndef GDALALG_VECTOR_CLIP_INCLUDED
+#define GDALALG_VECTOR_CLIP_INCLUDED
+
+#include "gdalalg_vector_pipeline.h"
+
+//! @cond Doxygen_Suppress
+
+/************************************************************************/
+/*                       GDALVectorClipAlgorithm                        */
+/************************************************************************/
+
+class GDALVectorClipAlgorithm /* non final */
+    : public GDALVectorPipelineStepAlgorithm
+{
+  public:
+    static constexpr const char *NAME = "clip";
+    static constexpr const char *DESCRIPTION = "Clip a vector dataset.";
+    static constexpr const char *HELP_URL = "/programs/gdal_vector_clip.html";
+
+    static std::vector<std::string> GetAliases()
+    {
+        return {};
+    }
+
+    explicit GDALVectorClipAlgorithm(bool standaloneStep = false);
+
+  private:
+    bool RunStep(GDALProgressFunc pfnProgress, void *pProgressData) override;
+
+    std::vector<double> m_bbox{};
+    std::string m_bboxCrs{};
+    std::string m_geometry{};
+    std::string m_geometryCrs{};
+    GDALArgDatasetValue m_likeDataset{};
+    std::string m_likeLayer{};
+    std::string m_likeSQL{};
+    std::string m_likeWhere{};
+};
+
+/************************************************************************/
+/*                   GDALVectorClipAlgorithmStandalone                  */
+/************************************************************************/
+
+class GDALVectorClipAlgorithmStandalone final : public GDALVectorClipAlgorithm
+{
+  public:
+    GDALVectorClipAlgorithmStandalone()
+        : GDALVectorClipAlgorithm(/* standaloneStep = */ true)
+    {
+    }
+};
+
+//! @endcond
+
+#endif /* GDALALG_VECTOR_CLIP_INCLUDED */
diff --git a/apps/gdalalg_vector_pipeline.cpp b/apps/gdalalg_vector_pipeline.cpp
index 7bee6af9a19a..eba4c0ae10ab 100644
--- a/apps/gdalalg_vector_pipeline.cpp
+++ b/apps/gdalalg_vector_pipeline.cpp
@@ -12,6 +12,7 @@
 
 #include "gdalalg_vector_pipeline.h"
 #include "gdalalg_vector_read.h"
+#include "gdalalg_vector_clip.h"
 #include "gdalalg_vector_filter.h"
 #include "gdalalg_vector_reproject.h"
 #include "gdalalg_vector_write.h"
@@ -174,6 +175,7 @@ GDALVectorPipelineAlgorithm::GDALVectorPipelineAlgorithm()
 
     m_stepRegistry.Register<GDALVectorReadAlgorithm>();
     m_stepRegistry.Register<GDALVectorWriteAlgorithm>();
+    m_stepRegistry.Register<GDALVectorClipAlgorithm>();
     m_stepRegistry.Register<GDALVectorReprojectAlgorithm>();
     m_stepRegistry.Register<GDALVectorFilterAlgorithm>();
 }
diff --git a/autotest/utilities/test_gdalalg_vector_clip.py b/autotest/utilities/test_gdalalg_vector_clip.py
new file mode 100755
index 000000000000..9c5afbe0379a
--- /dev/null
+++ b/autotest/utilities/test_gdalalg_vector_clip.py
@@ -0,0 +1,924 @@
+#!/usr/bin/env pytest
+# -*- coding: utf-8 -*-
+###############################################################################
+# Project:  GDAL/OGR Test Suite
+# Purpose:  'gdal vector clip' testing
+# Author:   Even Rouault <even dot rouault @ spatialys.com>
+#
+###############################################################################
+# Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
+#
+# SPDX-License-Identifier: MIT
+###############################################################################
+
+import ogrtest
+import pytest
+
+from osgeo import gdal, ogr, osr
+
+pytestmark = pytest.mark.require_geos
+
+
+def get_clip_alg():
+    reg = gdal.GetGlobalAlgorithmRegistry()
+    vector = reg.InstantiateAlg("vector")
+    return vector.InstantiateSubAlgorithm("clip")
+
+
+@pytest.mark.require_driver("GPKG")
+def test_gdalalg_vector_clip_general_behavior(tmp_vsimem):
+
+    out_filename = str(tmp_vsimem / "out.gpkg")
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        ["--bbox", "-1e10,-1e10,1e10,1e10", "../ogr/data/poly.shp", out_filename]
+    )
+
+    with gdal.OpenEx(out_filename, gdal.OF_UPDATE) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 10
+        for i in range(10):
+            ds.GetLayer(0).DeleteFeature(i + 1)
+
+    clip = get_clip_alg()
+    with pytest.raises(
+        Exception, match="already exists. Specify the --overwrite option"
+    ):
+        clip.ParseRunAndFinalize(
+            ["--bbox", "-1e10,-1e10,1e10,1e10", "../ogr/data/poly.shp", out_filename]
+        )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 0
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--overwrite",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 10
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--append",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 20
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--update",
+            "--nln",
+            "layer2",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayerByName("poly").GetFeatureCount() == 20
+        assert ds.GetLayerByName("layer2").GetFeatureCount() == 10
+
+    clip = get_clip_alg()
+    assert clip.ParseRunAndFinalize(
+        [
+            "--bbox",
+            "-1e10,-1e10,1e10,1e10",
+            "--of",
+            "GPKG",
+            "--overwrite-layer",
+            "--nln",
+            "poly",
+            "../ogr/data/poly.shp",
+            out_filename,
+        ]
+    )
+
+    with gdal.OpenEx(out_filename) as ds:
+        assert ds.GetLayerByName("poly").GetFeatureCount() == 10
+        assert ds.GetLayerByName("layer2").GetFeatureCount() == 10
+
+
+def test_gdalalg_vector_clip_bbox():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 -1,-1 -1,-1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--bbox", "0.2,0.3,0.7,0.8", "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_bbox_srs():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        [
+            "--bbox",
+            "22263.8981586547,33395.9998333802,77923.6435552915,89058.4864167096",
+            "--bbox-crs",
+            "EPSG:3857",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    assert out_lyr.GetSpatialRef().GetAuthorityCode(None) == "4326"
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_split_multipart():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test", geom_type=ogr.wkbPolygon)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt(
+            "POLYGON ((0 0,0 1,1 1,1 0,0.8 0,0.8 0.8,0.2 0.8,0.2 0,0 0))"
+        )
+    )
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "baz"
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.4 0.4,0.4 0.6,0.6 0.6,0.6 0.4,0.4 0.4))")
+    )
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--bbox", "-0.1,0.3,1.1,0.7", "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+
+    ref1_geom = ogr.CreateGeometryFromWkt(
+        "POLYGON ((0.0 0.7,0.2 0.7,0.2 0.3,0.0 0.3,0.0 0.7))"
+    )
+    ref2_geom = ogr.CreateGeometryFromWkt(
+        "POLYGON ((1.0 0.3,0.8 0.3,0.8 0.7,1.0 0.7,1.0 0.3))"
+    )
+
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    g = out_f.GetGeometryRef()
+    was_one = g.Within(ref1_geom) and ref1_geom.Within(g)
+    assert was_one or (g.Within(ref2_geom) and ref2_geom.Within(g))
+
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    g = out_f.GetGeometryRef()
+    if was_one:
+        assert g.Within(ref2_geom) and ref2_geom.Within(g)
+    else:
+        assert g.Within(ref1_geom) and ref1_geom.Within(g)
+
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "baz"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.4 0.4,0.4 0.6,0.6 0.6,0.6 0.4,0.4 0.4))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+@pytest.mark.parametrize(
+    "clip_geom",
+    [
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+        '{"type":"Polygon","coordinates":[[[0.2,0.8],[0.7,0.8],[0.7,0.3],[0.2,0.3],[0.2,0.8]]]}',
+    ],
+)
+def test_gdalalg_vector_clip_geom(clip_geom):
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 -1,-1 -1,-1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--geometry", clip_geom, "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+@pytest.mark.parametrize(
+    "clip_geom",
+    [
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+        '{"type":"Polygon","coordinates":[[[0.2,0.8],[0.7,0.8],[0.7,0.3],[0.2,0.3],[0.2,0.8]]]}',
+    ],
+)
+def test_gdalalg_vector_clip_geom_srs(clip_geom):
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip_geom = clip_geom.replace("0.2", "22263.8981586547")
+    clip_geom = clip_geom.replace("0.3", "33395.9998333802")
+    clip_geom = clip_geom.replace("0.7", "77923.6435552915")
+    clip_geom = clip_geom.replace("0.8", "89058.4864167096")
+    assert clip.ParseCommandLineArguments(
+        [
+            "--geometry",
+            clip_geom,
+            "--geometry-crs",
+            "EPSG:3857",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_geom_not_rectangle():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    # Bounding box of this polygon intersects the bounding box of the clipping
+    # geometry, but their full intersection is empty
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.4 0.4,0.4 0.5,0.5 0.5,0.4 0.4))")
+    )
+    src_lyr.CreateFeature(f)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((-1 -1,-1 2,2 2,2 -1,-1 -1))"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        [
+            "--geometry",
+            "POLYGON ((0 0,0 1,1 1,1 0,0.8 0,0.8 0.8,0.2 0.8,0.2 0,0 0))",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0 0,0 1,1 1,1 0,0.8 0,0.8 0.8,0.2 0.8,0.2 0,0 0))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_intersection_incompatible_geometry_type():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test", geom_type=ogr.wkbLineString)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f.SetGeometry(ogr.CreateGeometryFromWkt("LINESTRING(1 1,2 2)"))
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        [
+            "--geometry",
+            "POLYGON ((0 0,0 1,1 1,1 0,0 0))",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_intersection_promote_simple_type_to_multi():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test", geom_type=ogr.wkbMultiPolygon)
+
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt(
+            "MULTIPOLYGON(((-0.5 -0.5,-0.5 1.5,1.5 1.5,1.5 -0.5,-0.5 -0.5)))"
+        )
+    )
+    src_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        [
+            "--geometry",
+            "POLYGON ((0 0,0 1,1 1,1 0,0 0))",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "MULTIPOLYGON (((0 1,1 1,1 0,0 0,0 1)))",
+    )
+
+
+def test_gdalalg_vector_clip_like_vector():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_lyr = like_ds.CreateLayer("test")
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))")
+    )
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_vector_invalid_geom():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create(
+        "clip_ds", 0, 0, 0, gdal.GDT_Unknown
+    )
+    like_lyr = like_ds.CreateLayer("test")
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,1 1,0 1,1 0,0 0))"))
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    with gdal.quiet_errors(), pytest.raises(
+        Exception, match="Geometry of feature 0 of clip_ds is invalid"
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_vector_like_layer():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_ds.CreateLayer("unused")
+    like_lyr = like_ds.CreateLayer("my_layer")
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))")
+    )
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--like-layer", "my_layer", "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_vector_like_layer_invalid():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_ds.CreateLayer("test")
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--like-layer", "invalid", "--of", "Memory", "--output", "memory_ds"]
+    )
+    with pytest.raises(
+        Exception, match="Failed to identify source layer from clipping dataset."
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_vector_like_sql():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_ds.CreateLayer("unused")
+    like_lyr = like_ds.CreateLayer("my_layer")
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))")
+    )
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(
+        [
+            "--like-sql",
+            "SELECT * FROM my_layer",
+            "--of",
+            "Memory",
+            "--output",
+            "memory_ds",
+        ]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_vector_like_where():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_lyr = like_ds.CreateLayer("my_layer")
+    like_lyr.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f["id"] = 0
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    like_lyr.CreateFeature(f)
+
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    f["id"] = 1
+    f.SetGeometry(
+        ogr.CreateGeometryFromWkt("POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))")
+    )
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--like-where", "id=1", "--of", "Memory", "--output", "memory_ds"]
+    )
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_vector_like_where_empty():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_lyr = like_ds.CreateLayer("my_layer")
+    like_lyr.CreateField(ogr.FieldDefn("id", ogr.OFTInteger))
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--like-where", "id=1", "--of", "Memory", "--output", "memory_ds"]
+    )
+    with pytest.raises(Exception, match="No clipping geometry found"):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_vector_srs():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_3857 = osr.SpatialReference()
+    srs_3857.SetFromUserInput("EPSG:3857")
+    srs_3857.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    like_lyr = like_ds.CreateLayer("test", srs=srs_3857)
+    f = ogr.Feature(like_lyr.GetLayerDefn())
+    wkt = "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    wkt = wkt.replace("0.2", "22263.8981586547")
+    wkt = wkt.replace("0.3", "33395.9998333802")
+    wkt = wkt.replace("0.7", "77923.6435552915")
+    wkt = wkt.replace("0.8", "89058.4864167096")
+    f.SetGeometry(ogr.CreateGeometryFromWkt(wkt))
+    like_lyr.CreateFeature(f)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    ogrtest.check_feature_geometry(
+        ogr.CreateGeometryFromWkt(out_f.GetGeometryRef().ExportToWkt()),
+        "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))",
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_raster():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_lyr = src_ds.CreateLayer("test")
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1)
+    like_ds.SetGeoTransform([0.2, 0.5, 0, 0.3, 0, 0.5])
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_like_raster_srs():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    srs_wgs84 = osr.SpatialReference()
+    srs_wgs84.SetFromUserInput("WGS84")
+    srs_wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    src_lyr = src_ds.CreateLayer("test", srs=srs_wgs84)
+    src_lyr.CreateField(ogr.FieldDefn("foo"))
+    f = ogr.Feature(src_lyr.GetLayerDefn())
+    f["foo"] = "bar"
+    f.SetGeometry(ogr.CreateGeometryFromWkt("POLYGON ((0 0,0 1,1 1,1 0,0 0))"))
+    src_lyr.CreateFeature(f)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1)
+    like_ds.SetGeoTransform(
+        [22263.8981586547, 55660.4518654215, 0, 33395.9998333802, 0, 55659.7453966368]
+    )
+    srs_3857 = osr.SpatialReference()
+    srs_3857.SetFromUserInput("EPSG:3857")
+    srs_3857.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
+    like_ds.SetSpatialRef(srs_3857)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    assert clip.Run()
+
+    out_ds = clip.GetArg("output").Get().GetDataset()
+    out_lyr = out_ds.GetLayer(0)
+    out_f = out_lyr.GetNextFeature()
+    assert out_f["foo"] == "bar"
+    ogrtest.check_feature_geometry(
+        out_f, "POLYGON ((0.2 0.8,0.7 0.8,0.7 0.3,0.2 0.3,0.2 0.8))"
+    )
+
+    assert out_lyr.GetNextFeature() is None
+
+
+def test_gdalalg_vector_clip_missing_arg(tmp_vsimem):
+
+    out_filename = str(tmp_vsimem / "out.shp")
+
+    clip = get_clip_alg()
+    with pytest.raises(
+        Exception, match="clip: --bbox, --geometry or --like must be specified"
+    ):
+        clip.ParseRunAndFinalize(["../ogr/data/poly.shp", out_filename])
+
+
+def test_gdalalg_vector_clip_geometry_invalid():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_ds.CreateLayer("test")
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+
+    assert clip.ParseCommandLineArguments(
+        ["--geometry", "invalid", "--of", "Memory", "--output", "memory_ds"]
+    )
+    with pytest.raises(
+        Exception,
+        match="clip: Clipping geometry is neither a valid WKT or GeoJSON geometry",
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_vector_too_many_layers():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_ds.CreateLayer("test")
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("MEM").Create("", 1, 1, 1)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    with pytest.raises(
+        Exception,
+        match="clip: Dataset '' has no geotransform matrix. Its bounds cannot be established",
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_raster_no_geotransform():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    src_ds.CreateLayer("test")
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+    like_ds.CreateLayer("test1")
+    like_ds.CreateLayer("test2")
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    with pytest.raises(
+        Exception, match="clip: Only single layer dataset can be specified with --like"
+    ):
+        clip.Run()
+
+
+def test_gdalalg_vector_clip_like_neither_raster_no_vector():
+
+    src_ds = gdal.GetDriverByName("Memory").Create("", 0, 0, 0, gdal.GDT_Unknown)
+
+    # Create "like" dataset
+    like_ds = gdal.GetDriverByName("Memory").Create("clip", 0, 0, 0, gdal.GDT_Unknown)
+
+    clip = get_clip_alg()
+    clip.GetArg("input").Get().SetDataset(src_ds)
+    clip.GetArg("like").Get().SetDataset(like_ds)
+
+    assert clip.ParseCommandLineArguments(["--of", "Memory", "--output", "memory_ds"])
+    with pytest.raises(
+        Exception,
+        match="clip: Cannot get extent from clip dataset",
+    ):
+        clip.Run()
diff --git a/doc/source/conf.py b/doc/source/conf.py
index cb9096e2d997..51be7d65f45c 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -292,6 +292,13 @@
         [author_evenr],
         1,
     ),
+    (
+        "programs/gdal_vector_clip",
+        "gdal-vector-clip",
+        "Clip a vector dataset",
+        [author_evenr],
+        1,
+    ),
     (
         "programs/gdal_vector_convert",
         "gdal-vector-convert",
diff --git a/doc/source/programs/gdal_options/co_vector.rst b/doc/source/programs/gdal_options/co_vector.rst
new file mode 100644
index 000000000000..329bb104623a
--- /dev/null
+++ b/doc/source/programs/gdal_options/co_vector.rst
@@ -0,0 +1,18 @@
+.. option:: --co <NAME>=<VALUE>
+
+    Many formats have one or more optional dataset creation options that can be
+    used to control particulars about the file created. For instance,
+    the GeoPackage driver supports creation options to control the version.
+
+    May be repeated.
+
+    The dataset creation options available vary by format driver, and some
+    simple formats have no creation options at all. A list of options
+    supported for a format can be listed with the
+    :ref:`--formats <vector_common_options>`
+    command line option but the documentation for the format is the
+    definitive source of information on driver creation options.
+    See :ref:`vector_drivers` format
+    specific documentation for legal creation options for each format.
+
+    Note that dataset creation options are different from layer creation options.
diff --git a/doc/source/programs/gdal_vector.rst b/doc/source/programs/gdal_vector.rst
index 49d1a8ae48a7..9b503db16b7a 100644
--- a/doc/source/programs/gdal_vector.rst
+++ b/doc/source/programs/gdal_vector.rst
@@ -19,6 +19,7 @@ Synopsis
 
     Usage: gdal vector <SUBCOMMAND>
     where <SUBCOMMAND> is one of:
+      - clip:      Clip a vector dataset.
       - convert:   Convert a vector dataset.
       - filter:    Filter a vector dataset.
       - info:      Return information on a vector dataset.
@@ -28,8 +29,9 @@ Synopsis
 Available sub-commands
 ----------------------
 
-- :ref:`gdal_vector_info_subcommand`
+- :ref:`gdal_vector_clip_subcommand`
 - :ref:`gdal_vector_convert_subcommand`
+- :ref:`gdal_vector_info_subcommand`
 - :ref:`gdal_vector_pipeline_subcommand`
 
 Examples
diff --git a/doc/source/programs/gdal_vector_clip.rst b/doc/source/programs/gdal_vector_clip.rst
new file mode 100644
index 000000000000..b8d51e1da013
--- /dev/null
+++ b/doc/source/programs/gdal_vector_clip.rst
@@ -0,0 +1,170 @@
+.. _gdal_vector_clip_subcommand:
+
+================================================================================
+"gdal vector clip" sub-command
+================================================================================
+
+.. versionadded:: 3.11
+
+.. only:: html
+
+    Clip a vector dataset.
+
+.. Index:: gdal vector clip
+
+Synopsis
+--------
+
+.. code-block::
+
+    Usage: gdal vector clip [OPTIONS] <INPUT> <OUTPUT>
+
+    Clip a vector dataset.
+
+    Positional arguments:
+      -i, --input <INPUT>                                  Input vector dataset [required]
+      -o, --output <OUTPUT>                                Output vector dataset [required]
+
+    Common Options:
+      -h, --help                                           Display help message and exit
+      --version                                            Display GDAL version and exit
+      --json-usage                                         Display usage as JSON document and exit
+      --drivers                                            Display driver list as JSON document and exit
+      --config <KEY>=<VALUE>                               Configuration option [may be repeated]
+      --progress                                           Display progress bar
+
+    Options:
+      -l, --layer, --input-layer <INPUT-LAYER>             Input layer name(s) [may be repeated]
+      -f, --of, --format, --output-format <OUTPUT-FORMAT>  Output format
+      --co, --creation-option <KEY>=<VALUE>                Creation option [may be repeated]
+      --lco, --layer-creation-option <KEY>=<VALUE>         Layer creation option [may be repeated]
+      --overwrite                                          Whether overwriting existing output is allowed
+      --update                                             Whether to open existing dataset in update mode
+      --overwrite-layer                                    Whether overwriting existing layer is allowed
+      --append                                             Whether appending to existing layer is allowed
+      --output-layer <OUTPUT-LAYER>                        Output layer name
+      --bbox <BBOX>                                        Clipping bounding box as xmin,ymin,xmax,ymax
+                                                           Mutually exclusive with --geometry, --like
+      --bbox-crs <BBOX-CRS>                                CRS of clipping bounding box
+      --geometry <GEOMETRY>                                Clipping geometry (WKT or GeoJSON)
+                                                           Mutually exclusive with --bbox, --like
+      --geometry-crs <GEOMETRY-CRS>                        CRS of clipping geometry
+      --like <DATASET>                                     Dataset to use as a template for bounds
+                                                           Mutually exclusive with --bbox, --geometry
+      --like-sql <SELECT-STATEMENT>                        SELECT statement to run on the 'like' dataset
+      --like-layer <LAYER-NAME>                            Name of the layer of the 'like' dataset
+      --like-where <WHERE-EXPRESSION>                      Expression for a WHERE SQL clause to run on the 'like' dataset
+
+
+    Advanced Options:
+      --if, --input-format <INPUT-FORMAT>                  Input formats [may be repeated]
+      --oo, --open-option <KEY=VALUE>                      Open options [may be repeated]
+
+
+
+Description
+-----------
+
+:program:`gdal vector clip` can be used to clip a vector dataset using
+georeferenced coordinates.
+
+Either :option:`--bbox`, :option:`--geometry` or :option:`--like` must be specified.
+
+``clip`` can also be used as a step of :ref:`gdal_vector_pipeline_subcommand`.
+
+Clipping sometimes results in geometries that are of a type not compatible
+with the geometry type. This program splits multi-geometries or geometry
+collections into their parts for layers that are of a single geometry type
+(such as point, linestring, polygon) and only keeps the parts that are of that
+type. It also promotes single geometry types (e.g. polygons) to multi
+geometry types (e.g. multi-polygons). If the user needs to preserve any type of
+intersection, the layer must use the ``wkbUnknown`` (any geometry) type.
+
+Standard options
+++++++++++++++++
+
+.. include:: gdal_options/of_vector.rst
+
+.. include:: gdal_options/co_vector.rst
+
+.. include:: gdal_options/overwrite.rst
+
+.. option:: --bbox <xmin>,<ymin>,<xmax>,<ymax>
+
+    Bounds to which to clip the dataset. They are assumed to be in the CRS of
+    the input dataset, unless :option:`--bbox-crs` is specified.
+    The X and Y axis are the "GIS friendly ones", that is X is longitude or easting,
+    and Y is latitude or northing.
+    Mutually exclusive with :option:`--like` and :option:`--geometry`.
+
+.. option:: --bbox-crs <CRS>
+
+    CRS in which the <xmin>,<ymin>,<xmax>,<ymax> values of :option:`--bbox`
+    are expressed. If not specified, it is assumed to be the CRS of the input
+    dataset.
+    Not that specifying --bbox-crs does not involve doing vector reprojection.
+    Instead, the bounds are reprojected from the bbox-crs to the CRS of the
+    input dataset.
+
+.. option:: --geometry <WKT_or_GeoJSON>
+
+    Geometry as a WKT or GeoJSON string to which to clip the dataset.
+    It is assumed to be in the CRS of
+    the input dataset, unless :option:`--geometry-crs` is specified.
+    The X and Y axis are the "GIS friendly ones", that is X is longitude or easting,
+    and Y is latitude or northing.
+    Mutually exclusive with :option:`--bbox` and :option:`--like`.
+
+.. option:: --geometry-crs <CRS>
+
+    CRS in which the coordinates values of :option:`--geometry`
+    are expressed. If not specified, it is assumed to be the CRS of the input
+    dataset.
+    Not that specifying --geometry-crs does not involve doing vector reprojection.
+    Instead, the bounds are reprojected from the geometry-crs to the CRS of the
+    input dataset.
+
+.. option:: --like <DATASET>
+
+    Vector or raster dataset to use as a template for bounds.
+    If the specified dataset is a raster, its rectangular bounds are use as
+    the clipping geometry.
+    If the specified dataset is a vector dataset, its polygonal geometries
+    are unioned together to form the clipping geometry (beware that the result
+    union might not be perfect if there are gaps between individual polygon
+    features). If several layers are present,
+    :option:`--like-sql` or :option:`--like-layer` must be specified.
+    Mutually exclusive with :option:`--bbox` and :option:`--geometry`.
+
+.. option:: --like-sql <SELECT-STATEMENT>
+
+    Select desired geometries from the vector clip dataset using an SQL query.
+    e.g ``SELECT geom FROM my_layer WHERE country = 'France'``.
+    Mutually exclusive with :option:`--like-layer` and :option:`--like-where`
+
+.. option:: --like-layer <LAYER-NAME>
+
+    Select the named layer from the vector clip dataset.
+    Mutually exclusive with :option:`--like-sql`
+
+.. option:: --like-where <WHERE-EXPRESSION>
+
+    Restrict desired geometries from vector clip dataset layer based on an attribute query.
+    e.g ``country = 'France'``.
+
+Advanced options
+++++++++++++++++
+
+.. include:: gdal_options/oo.rst
+
+.. include:: gdal_options/if.rst
+
+Examples
+--------
+
+.. example::
+   :title: Clip a GeoPackage file to the bounding box from longitude 2, latitude 49, to longitude 3, latitude 50 in WGS 84
+
+   .. code-block:: bash
+
+        $ gdal vector clip --bbox=2,49,3,50 --bbox-crs=EPSG:4326 in.gpkg out.gpkg --overwrite
diff --git a/doc/source/programs/gdal_vector_convert.rst b/doc/source/programs/gdal_vector_convert.rst
index fdd736101b28..f475ecccc818 100644
--- a/doc/source/programs/gdal_vector_convert.rst
+++ b/doc/source/programs/gdal_vector_convert.rst
@@ -59,7 +59,7 @@ Standard options
 
 .. include:: gdal_options/of_vector.rst
 
-.. include:: gdal_options/co.rst
+.. include:: gdal_options/co_vector.rst
 
 .. include:: gdal_options/overwrite.rst
 
diff --git a/doc/source/programs/gdal_vector_pipeline.rst b/doc/source/programs/gdal_vector_pipeline.rst
index 394ff0d092ae..20d6b756fccd 100644
--- a/doc/source/programs/gdal_vector_pipeline.rst
+++ b/doc/source/programs/gdal_vector_pipeline.rst
@@ -52,6 +52,28 @@ Potential steps are:
       --if, --input-format <INPUT-FORMAT>                  Input formats [may be repeated]
       --oo, --open-option <KEY=VALUE>                      Open options [may be repeated]
 
+* clip [OPTIONS]
+
+.. code-block::
+
+    Clip a vector dataset.
+
+    Options:
+      --bbox <BBOX>                                        Clipping bounding box as xmin,ymin,xmax,ymax
+                                                           Mutually exclusive with --geometry, --like
+      --bbox-crs <BBOX-CRS>                                CRS of clipping bounding box
+      --geometry <GEOMETRY>                                Clipping geometry (WKT or GeoJSON)
+                                                           Mutually exclusive with --bbox, --like
+      --geometry-crs <GEOMETRY-CRS>                        CRS of clipping geometry
+      --like <DATASET>                                     Dataset to use as a template for bounds
+                                                           Mutually exclusive with --bbox, --geometry
+      --like-sql <SELECT-STATEMENT>                        SELECT statement to run on the 'like' dataset
+      --like-layer <LAYER-NAME>                            Name of the layer of the 'like' dataset
+      --like-where <WHERE-EXPRESSION>                      Expression for a WHERE SQL clause to run on the 'like' dataset
+
+
+Details for options can be found in :ref:`gdal_vector_clip_subcommand`.
+
 * filter [OPTIONS]
 
 .. code-block::
diff --git a/doc/source/programs/index.rst b/doc/source/programs/index.rst
index 8fd8d3ab2945..aebb8a2c34b6 100644
--- a/doc/source/programs/index.rst
+++ b/doc/source/programs/index.rst
@@ -41,6 +41,7 @@ single :program:`gdal` program that accepts commands and subcommands.
    gdal_raster_reproject
    gdal_vector
    gdal_vector_info
+   gdal_vector_clip
    gdal_vector_convert
    gdal_vector_pipeline
 
@@ -62,6 +63,7 @@ single :program:`gdal` program that accepts commands and subcommands.
     - :ref:`gdal_raster_reproject_subcommand`: Reproject a raster dataset
     - :ref:`gdal_vector_command`: Entry point for vector commands
     - :ref:`gdal_vector_info_subcommand`: Get information on a vector dataset
+    - :ref:`gdal_vector_clip_subcommand`: Clip a vector dataset
     - :ref:`gdal_vector_convert_subcommand`: Convert a vector dataset
     - :ref:`gdal_vector_pipeline_subcommand`: Process a vector dataset
 

From 8056d2f0fc2bf1fa553ff11da4de8d4e16c2bd6b Mon Sep 17 00:00:00 2001
From: Even Rouault <even.rouault@spatialys.com>
Date: Tue, 21 Jan 2025 12:49:52 +0100
Subject: [PATCH 7/7] DXF: interpret INSERT blocks with row count or column
 count equal to 0 as 1

---
 autotest/ogr/ogr_dxf.py             | 6 ++----
 ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp | 6 ++++--
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/autotest/ogr/ogr_dxf.py b/autotest/ogr/ogr_dxf.py
index 9932e5288454..90c496d999b7 100644
--- a/autotest/ogr/ogr_dxf.py
+++ b/autotest/ogr/ogr_dxf.py
@@ -4129,9 +4129,7 @@ def test_ogr_dxf_write_MEASUREMENT(tmp_vsimem):
 ###############################################################################
 # Use case of https://github.com/OSGeo/gdal/issues/11591
 # Test reading a INSERT block whose column count is zero.
-# Not totally sure about the exact behavior we should do. Currently
-# the zero count is strictly honored and no geometry from the block will be
-# inserted into the regular geometries. LibreCAD 2.1.3 does the same thing
+# Interpretating it as 1, as AutoCAD does
 
 
 @gdaltest.enable_exceptions()
@@ -4139,7 +4137,7 @@ def test_ogr_dxf_insert_col_count_zero():
 
     with ogr.Open("data/dxf/insert_only_col_count_zero.dxf") as ds:
         lyr = ds.GetLayer(0)
-        assert lyr.GetFeatureCount() == 0
+        assert lyr.GetFeatureCount() == 1
 
     with gdal.config_option("DXF_INLINE_BLOCKS", "NO"):
         with ogr.Open("data/dxf/insert_only_col_count_zero.dxf") as ds:
diff --git a/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp b/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp
index 64aa386f81bc..bc3e519499c2 100644
--- a/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp
+++ b/ogr/ogrsf_frmts/dxf/ogrdxflayer.cpp
@@ -3207,8 +3207,10 @@ bool OGRDXFLayer::TranslateINSERT()
 
     if (m_oInsertState.m_nRowCount == 0 || m_oInsertState.m_nColumnCount == 0)
     {
-        m_oInsertState.m_nRowCount = 0;
-        m_oInsertState.m_nColumnCount = 0;
+        // AutoCad doesn't allow setting to 0 in its UI, but interprets 0
+        // as 1 (but other software such as LibreCAD interpret 0 as 0)
+        m_oInsertState.m_nRowCount = 1;
+        m_oInsertState.m_nColumnCount = 1;
     }
 
     /* -------------------------------------------------------------------- */