[GIS] How to let rasterio use a colormap file when showing a TIFF

colormapgeotiff-tiffpythonrasterio

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!

Best Answer

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
Related Question