feat: Create RGB WebP example with photometric interpretation set

rasterio_generated/fixtures/uint8_rgb_webp_block64_cog.py

@@ -0,0 +1,32 @@
+"""Generate an RGB GeoTIFF with DEFLATE compression."""
+
+from pathlib import Path
+
+import numpy as np
+
+from rasterio_generated.write_utils import write_tiff
+
+
+def generate(output_path: Path) -> None:
+    """Generate a 256x256 RGB GeoTIFF with DEFLATE compression."""
+    # Create RGB gradient pattern
+    r = np.linspace(0, 127, 128, dtype=np.uint8)
+    g = np.linspace(127, 0, 128, dtype=np.uint8)
+    b = np.full(128, 128, dtype=np.uint8)
+
+    data = np.stack(
+        [
+            np.tile(r, (128, 1)),
+            np.tile(g.reshape(-1, 1), (1, 128)),
+            np.tile(b, (128, 1)),
+        ]
+    )
+
+    write_tiff(
+        output_path,
+        data,
+        driver="COG",
+        blocksize=64,
+        photometric="RGB",
+        compress="WEBP",
+    )

rasterio_generated/fixtures/uint8_rgb_webp_block64_cog_info.md

@@ -0,0 +1,55 @@
+```
+Driver: GTiff
+File: rasterio_generated/fixtures/uint8_rgb_webp_block64_cog.tif
+COG: True
+Compression: WEBP
+ColorSpace: None
+
+Profile
+    Width:            128
+    Height:           128
+    Bands:            3
+    Tiled:            True
+    Dtype:            uint8
+    NoData:           None
+    Alpha Band:       False
+    Internal Mask:    False
+    Interleave:       PIXEL
+    ColorMap:         False
+    ColorInterp:      ('red', 'green', 'blue')
+    Scales:           (1.0, 1.0, 1.0)
+    Offsets:          (0.0, 0.0, 0.0)
+
+Geo
+    Crs:              EPSG:4326
+    Origin:           (0.0, 0.0)
+    Resolution:       (0.01, -0.01)
+    BoundingBox:      (0.0, -1.28, 1.28, 0.0)
+    MinZoom:          7
+    MaxZoom:          7
+
+Image Metadata
+    AREA_OR_POINT: Area
+
+Image Structure
+    LAYOUT: COG
+    COMPRESSION: WEBP
+    INTERLEAVE: PIXEL
+    OVERVIEW_RESAMPLING: CUBIC
+    COMPRESSION_REVERSIBILITY: LOSSY
+    WEBP_LEVEL: 75
+
+Band 1
+    ColorInterp: red
+
+Band 2
+    ColorInterp: green
+
+Band 3
+    ColorInterp: blue
+
+IFD
+    Id      Size           BlockSize     Decimation
+    0       128x128        64x64         0
+    1       64x64          64x64         2
+```

rasterio_generated/write_utils.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import tempfile
 from pathlib import Path
 from typing import TYPE_CHECKING, Literal
 
@@ -33,11 +34,16 @@ def write_tiff(
     transform: Affine = from_origin(0, 0, 0.01, 0.01),
     predictor: Literal[2, 3] | None = None,
     nodata: int | float | None = None,
+    photometric: Literal["RGB", "MINISBLACK"] | None = None,
 ):
     """
     data:
       - shape (H, W)           -> 1 band
       - shape (B, H, W)        -> multi-band (e.g. RGB)
+
+    Note: The GDAL COG driver does not support the PHOTOMETRIC creation option.
+    When driver="COG" and photometric is set, a two-step process is used:
+    first write a GTiff with photometric, then convert to COG via gdal.Translate.
     """
     if data.ndim == 2:
         count = 1
@@ -47,6 +53,23 @@ def write_tiff(
     else:
         raise ValueError("data must be 2D or 3D (bands, height, width)")
 
+    # The COG driver ignores PHOTOMETRIC. Work around this by writing a
+    # temporary GTiff first, then converting to COG with gdal.Translate.
+    if driver == "COG" and photometric is not None:
+        _write_cog_with_photometric(
+            path,
+            data,
+            count=count,
+            blocksize=blocksize,
+            compress=compress,
+            crs=crs,
+            transform=transform,
+            predictor=predictor,
+            nodata=nodata,
+            photometric=photometric,
+        )
+        return
+
     profile = {
         "driver": driver,
         "width": width,
@@ -56,7 +79,7 @@ def write_tiff(
         "crs": crs,
         "transform": transform,
         "tiled": blocksize is not None,
-        "interleave": "pixel",  # explicit chunky
+        "interleave": "pixel",
     }
 
     if blocksize is not None:
@@ -71,8 +94,77 @@ def write_tiff(
     if nodata is not None:
         profile["nodata"] = nodata
 
+    if photometric is not None:
+        profile["photometric"] = photometric
+
     with rasterio.open(path, "w", **profile) as dst:
         if count == 1:
             dst.write(data, 1)
         else:
             dst.write(data)
+
+
+def _write_cog_with_photometric(
+    path: Path,
+    data: np.ndarray,
+    *,
+    count: int,
+    blocksize: int | None,
+    compress: str | None,
+    crs: str,
+    transform: Affine,
+    predictor: int | None,
+    nodata: int | float | None,
+    photometric: str,
+):
+    """Write a COG via a temporary GTiff to preserve photometric interpretation."""
+    from osgeo import gdal
+
+    height = data.shape[-2]
+    width = data.shape[-1]
+
+    with tempfile.NamedTemporaryFile(suffix=".tif", delete=True) as tmp:
+        # Step 1: Write GTiff with photometric
+        gtiff_profile = {
+            "driver": "GTiff",
+            "width": width,
+            "height": height,
+            "count": count,
+            "dtype": data.dtype,
+            "crs": crs,
+            "transform": transform,
+            "tiled": True,
+            "blockxsize": blocksize or 256,
+            "blockysize": blocksize or 256,
+            "interleave": "pixel",
+            "photometric": photometric,
+        }
+
+        if compress is not None:
+            gtiff_profile["compress"] = compress
+
+        if predictor is not None:
+            gtiff_profile["predictor"] = predictor
+
+        if nodata is not None:
+            gtiff_profile["nodata"] = nodata
+
+        with rasterio.open(tmp.name, "w", **gtiff_profile) as dst:
+            if count == 1:
+                dst.write(data, 1)
+            else:
+                dst.write(data)
+
+        # Step 2: Convert to COG via gdal.Translate
+        gdal.UseExceptions()
+        creation_options = []
+        if compress is not None:
+            creation_options.append(f"COMPRESS={compress}")
+        if blocksize is not None:
+            creation_options.append(f"BLOCKSIZE={blocksize}")
+        if predictor is not None:
+            creation_options.append(f"PREDICTOR={predictor}")
+
+        ds = gdal.Open(tmp.name)
+        gdal.Translate(str(path), ds, format="COG", creationOptions=creation_options)
+        ds = None