Repository: rasterio/rasterio (language: python)

Issue:
GeoTIFF written with v1.3.10 differs from one written with v1.3.9 in the color interpretation tag
<!--

WELCOME ABOARD!

Hi and welcome to the Rasterio project. We appreciate bug reports, questions
about documentation, and suggestions for new features. This issue template
isn't intended to ward you off; only to intercept and redirect some particular
categories of reports, and to collect a few important facts that issue reporters
often omit.

The primary forum for questions about installation and usage of Rasterio is 
https://rasterio.groups.io/g/main. The authors and other users will answer 
questions when they have expertise to share and time to explain. Please take the
time to craft a clear question and be patient about responses. Please do not
bring these questions to Rasterio's issue tracker, which we want to reserve for
bug reports and other actionable issues.

Questions about development of Rasterio, brainstorming, requests for comment,
and not-yet-actionable proposals are welcome in the project's developers 
discussion group https://rasterio.groups.io/g/dev. Issues opened in Rasterio's
GitHub repo which haven't been socialized there may be perfunctorily closed.

Please note: Rasterio contains extension modules and is thus susceptible to
C library compatibility issues. If you are reporting an installation or module
import issue, please note that this project only accepts reports about problems
with packages downloaded from the Python Package Index. Conda users should take
issues to one of the following trackers:

- https://github.com/ContinuumIO/anaconda-issues/issues
- https://github.com/conda-forge/rasterio-feedstock

You think you've found a bug? We believe you!
-->

## Expected behavior and actual behavior.

Geotiffs containing a colormap written with version `1.3.9` differ form those written with `1.3.10`. While `gdalinfo` (and therefore QGIS) see no difference, `tiffinfo` reports that the photometric interpretation of files written with `1.3.9` is `palette color (RGB from colormap)` while files written with `1.3.10` report `min-is-black`. This breaks display of such Geotiffs e.g. in https://geotiffjs.github.io/ . 

## Steps to reproduce the problem.

Minimum reproducible example writing a Geotiff of size 2x2 with a colormap of 3 values.

```
from importlib.metadata import version

import numpy as np
import rasterio
from affine import Affine
from rasterio import CRS
from rasterio.profiles import DefaultGTiffProfile

data_array = np.array([[0, 1, 2], [200, 2, 1]], dtype=np.uint8)

file_path = f"{version('rasterio')}-{version('numpy')}.tiff"

profile = DefaultGTiffProfile(
    **{
        'height': data_array.shape[0],
        'width': data_array.shape[1],
        'count': 1,
        'dtype': data_array.dtype,
        'crs': CRS.from_epsg(4326),
        'nodata': 200,
        'transform': Affine(
            0.0001374386909998465, 0.0, 8.573052230839624, 0.0, -9.041674461565189e-05, 49.45977893086753
        ),
    }
)

with rasterio.open(file_path, mode='w', **profile) as out_map_file:
    out_map_file.write(data_array, indexes=1)

    out_map_file.write_colormap(1, {0: (0, 0, 0), 1: (255, 255, 255), 2: (255, 0, 0)})

```

When examining the above created files with version 1.3.9 and 1.3.10 they differ as described above. This is independent of the numpy version used. The error still exists in version 1.4a3. 

### Notes from investigation

I am not deep into GDAL or rasterio but I will still share what I found out so far:

The relevant function seems to be https://github.com/rasterio/rasterio/blob/835c10b2289852376cb3ef3b40614502c30fa104/rasterio/_io.pyx#L1916-L1956

which seems to hot have been touched between the two versions. Because GDAL was updated between the two releases from `3.6.4` to `3.8.4` this may also cause the problem. I think 

https://github.com/rasterio/rasterio/blob/835c10b2289852376cb3ef3b40614502c30fa104/rasterio/_io.pyx#L1952

should in fact be `            GDALSetRasterColorInterpretation(hBand, <GDALColorInterp>2)` to set the correct color interpretation but I'm far outside my competences here. Maybe GDAL had some smartness (or negligence) in an older version that was removed wherefore the `min-is-black` interpretation is now set.

Here is the script that I used: [repex.zip](https://github.com/user-attachments/files/16566111/repex.zip)


#### Environment Information

```
rasterio info:
  rasterio: 1.3.9
      GDAL: 3.6.4
      PROJ: 9.0.1
      GEOS: 3.11.1
 PROJ DATA: ~/rasterio_test/v1.3.9/lib/python3.10/site-packages/rasterio/proj_data
 GDAL DATA: ~/rasterio_test/v1.3.9/lib/python3.10/site-packages/rasterio/gdal_data

System:
    python: 3.10.12 (main, Jul 29 2024, 16:56:48) [GCC 11.4.0]
executable: ~/rasterio_test/v1.3.9/bin/python3
   machine: Linux-6.5.0-45-generic-x86_64-with-glibc2.35

Python deps:
    affine: 2.4.0
     attrs: 24.2.0
   certifi: 2024.07.04
     click: 8.1.7
     cligj: 0.7.2
    cython: None
     numpy: 1.26.4
    snuggs: 1.4.7
click-plugins: None
setuptools: 59.6.0
```

and

```
rasterio info:
  rasterio: 1.3.10
      GDAL: 3.8.4
      PROJ: 9.3.1
      GEOS: 3.11.1
 PROJ DATA: ~/rasterio_test/v1.3.10-numpy1/lib/python3.10/site-packages/rasterio/proj_data
 GDAL DATA: ~/rasterio_test/v1.3.10-numpy1/lib/python3.10/site-packages/rasterio/gdal_data

System:
    python: 3.10.12 (main, Jul 29 2024, 16:56:48) [GCC 11.4.0]
executable: ~/rasterio_test/v1.3.10-numpy1/bin/python3
   machine: Linux-6.5.0-45-generic-x86_64-with-glibc2.35

Python deps:
    affine: 2.4.0
     attrs: 24.2.0
   certifi: 2024.07.04
     click: 8.1.7
     cligj: 0.7.2
    cython: None
     numpy: 1.26.4
    snuggs: 1.4.7
click-plugins: None
setuptools: 59.6.0
```



## Installation Method

The `rasterio-1.3.9-cp310-cp310-manylinux2014_x86_64.whl` and `rasterio-1.3.10-cp310-cp310-manylinux2014_x86_64.whl` installed from PyPI using pip 22.0.2.



Discussion / hints:
Thanks for the report @SlowMo24 !
@SlowMo24 With what version of tiffinfo do you see the problem? Is it older that the one linked by gdalinfo and QGIS?

Fix patch:
diff --git a/CHANGES.txt b/CHANGES.txt
index c77090ba8..9446939cb 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,14 @@
 Changes
 =======
 
+Next (TBD)
+----------
+
+Bug fixes:
+
+- Color interpretation is correctly set to "palette" after a colormap is
+  written (#3133).
+
 1.4b1 (2024-08-10)
 ------------------
 
diff --git a/rasterio/_io.pyx b/rasterio/_io.pyx
index 6b29e3c5f..6b040bbf2 100644
--- a/rasterio/_io.pyx
+++ b/rasterio/_io.pyx
@@ -2036,7 +2036,7 @@ cdef class DatasetWriterBase(DatasetReaderBase):
                 GDALSetColorEntry(hTable, i, &color)
 
             # TODO: other color interpretations?
-            GDALSetRasterColorInterpretation(hBand, <GDALColorInterp>1)
+            GDALSetRasterColorInterpretation(hBand, GCI_PaletteIndex)
             GDALSetRasterColorTable(hBand, hTable)
 
         finally:


Test patch:
diff --git a/tests/test_colormap.py b/tests/test_colormap.py
index ead496fc3..d37c8d5e3 100644
--- a/tests/test_colormap.py
+++ b/tests/test_colormap.py
@@ -1,27 +1,26 @@
+"""Colormap tests."""
+
 import rasterio
+from rasterio.enums import ColorInterp
 
 
-def test_write_colormap_warn(tmpdir, recwarn):
-    with rasterio.open('tests/data/shade.tif') as src:
-        profile = src.meta
-    tiffname = str(tmpdir.join('foo.tif'))
-    with rasterio.open(tiffname, 'w', **profile) as dst:
-        dst.write_colormap(1, {0: (255, 0, 0, 255), 255: (0, 0, 0, 0)})
+def test_write_colormap(tmp_path):
+    with rasterio.open("tests/data/shade.tif") as src:
+        shade = src.read(1)
+        profile = src.profile
 
+    profile["driver"] = "PNG"
 
-def test_write_colormap(tmpdir):
-    with rasterio.open('tests/data/shade.tif') as src:
-        shade = src.read(1)
-        meta = src.meta
-    tiffname = str(tmpdir.join('foo.png'))
-    meta['driver'] = 'PNG'
-    with rasterio.open(tiffname, 'w', **meta) as dst:
+    with rasterio.open(tmp_path / "test.tif", "w", **profile) as dst:
         dst.write(shade, indexes=1)
         dst.write_colormap(1, {0: (255, 0, 0, 255), 255: (0, 0, 0, 0)})
+        assert dst.colorinterp == (ColorInterp.palette,)
         cmap = dst.colormap(1)
         assert cmap[0] == (255, 0, 0, 255)
         assert cmap[255] == (0, 0, 0, 0)
-    with rasterio.open(tiffname) as src:
+
+    with rasterio.open(tmp_path / "test.tif") as src:
+        assert src.colorinterp == (ColorInterp.palette,)
         cmap = src.colormap(1)
         assert cmap[0] == (255, 0, 0, 255)
         assert cmap[255] == (0, 0, 0, 0)