GTiff/COG: add support for INTERLEAVE=TILE; COG: add support for INTE…

…RLEAVE=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

@@ -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"],
+        )

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>

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