How to let rasterio use a colormap file when showing a TIFF

Question

I'm new to rasterio. I have tried this example. My raster file is a GeoTIFF file and I have a colormap file available in ESRI *.clr format, with 58 lines. The values in the GeoTIFF file range from 10471 thru 10534. For those who don't know the format, the first 2 lines looks like this:
10471 16 24 79
10472 104 223 247
I have read a bit about colormaps, e.g. here. I read this: "You cannot set the color interpretation on existing data". All the same, I'm hoping that I can have the colormap applied to the plotted raster in a direct way, but I have no idea how to go about it. In order to start somewhere, I wrote this code:

import csv
import os.path
import sys
cf = None
try:
    curdir = os.path.dirname(sys.argv[0])
    datapath = os.path.normpath(os.path.join(curdir, "../geodata/EUR"))
    cf = open(os.path.join(datapath, "countries.clr"), "r")
    rows = {}
    for row in (csv.reader(cf, delimiter=' ')):
        color = tuple([int(row[1]), int(row[2]), int(row[3]), 255])
        rows[int(row[0])] = color
    # Do we almost have a colormap now?
finally:
    if cf != None: cf.close()

Later on, I'd also like to indicate that rasterio should plot classes - defined somehow - in different colours. I'm afraid that it's only possible to get this working in an indirect way - with some time-consuming tricks - but if you can, please surprise me!

Answer 1

I've had partial success. At least I've been able to write a GeoTIFF - with a suitable colormap built in - to file. It was less easy than I had hoped. Anyway, these are the most important code lines:

# Now open the input raster - with wanted size - but Type=Int32, ColorInterp=Gray
datapath = os.path.normpath(os.path.join(curdir, "../geodata/EUR"))
with rasterio.open(os.path.join(datapath, "input.tif")) as rf:
    # Retrieve some particulars of the input file
    data = rf.read(1)
    meta = rf.meta

    # Prepare the new data structure - up to you which values < 255 to add
    newdata = np.empty(data.shape, dtype=np.uint8)
    newdata.fill(255)

    # Fill the new data structure cell by cell in the way you want
    for i in ...
        for k in ...
            newdata[i, k] = int(x)

    # Get the color map
    cmap = get_colormap(os.path.join(datapath, "palette.clr"))
    cmap[255] = (255, 255, 255)
    meta["dtype"] = 'uint8'
    meta["nodata"] = 255

    # Write the GeoTIFF to file
    with rasterio.open(os.path.join(datapath, "output.tif"), "w", **meta) as tf: 
        tf.nodata = 255
        tf.write(newdata, indexes=1)
        tf.write_colormap(1, cmap)

This way I obtained a GeoTIFF which is shown in the right colors in QGIS as well as in the graphics viewers commonly available on the Windows platform. Next challenge is to actually have rasterio also show the output file in the indicated colors.

P.S. Method get_colormap is as follows:

def get_colormap(fn):
    if not os.path.exists(fn):
        raise IOError("File not found!")
    with open(fn, "rb") as cf:
        colormgr = csv.reader(cf, delimiter=' ')
        result = {}
        for row in colormgr:
            color = tuple([int(row[1]), int(row[2]), int(row[3])])
            result[int(row[0])] = color
    return result

Answer 2

My understanding of Rasterio is that when you read a GeoTIFF it reads it as an array. Because Rasterio/Matplotlib display rasters as arrays it doesn't automatically use the colourmap associated with the tiff file (but it does automatically show up in ArcGIS).

Here is a slightly simpler way of setting the colourmap for a raster that I've used many times:

import rasterio

infl = r"D:/folder/raster.tif"  # Input file path  

# Read the file metadata and raster data from the tiff file                              
with rasterio.open(infl) as src:     
    data = src.read(1)
    meta = src.meta 

infl_coloured = r"D:/folder/raster_coloured.tif"     # Output file path

with rasterio.open(infl_coloured, 'w', **meta) as dst:
    dst.write(data, indexes=1)

    # Create a colourmap dictionary using RGBA values (i.e., {value1: (R,G,B,A), value2: (R,G,B,A), .....})
    dst.write_colormap(                             
        1, {0: (111, 112, 112, 1),
             1: (54, 58, 64, 1)})   

    # Add the colourmap to the destination file you're writing through rasterio
    cmap = dst.colormap(1)        

This website is very helpful for creating custom colour tables to use for the colourmap dictionary: https://www.rampgenerator.com/

Hopefully this helps!

Answer 3

As far as I can tell the GDAL library (which rasterio uses) doesn't consider .clr files when opening a GeoTIFF. It does consider a .clr file that accompanies an AIG file (see http://www.gdal.org/frmt_various.html and search for "AIG"), so you may have success if you convert the GeoTIFF to the old Arc/Info Binary Grid format.

You can in fact update the colormap of a GeoTIFF file. The documentation is out of date and misleading and will be corrected. Here's an example: https://github.com/mapbox/rasterio/blob/master/tests/test_rio_merge.py#L82-L83. I think that converting the file is the thing to try first.