4


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!

2

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
  • Where does that get_colormap function come from? You need to write complete snippets, please. – j08lue Jan 23 '18 at 8:21
  • 1
    At the request of j08lue, I've added the code for method get_colormap. – Dobedani Jan 24 '18 at 9:58
0

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.

  • Yes, I've tried opening such a grid with different arguments but without success. I understand that QGIS also uses GDAL. I suppose that by using QGIS, we can study the behaviour of GDAL a bit. QGIS works with a different colormap file format - similar to CSV - with extension *.txt. Now, QGIS only seems to consider such colormap files in certain conditions, viz. when you add the Arc/Info Grid raster to a new project. It's unfortunate. I'll have to study the predefined colormaps. Maybe I can adapt things so that one of those becomes suitable for my purpose. – Dobedani Dec 7 '17 at 12:41
  • No, the colormaps of matplotlib are not useful for me. I tried to write a temporary GeoTIFF file, combining 2 examples on the – Dobedani Dec 7 '17 at 16:57
  • ... website hosting the rasterio Python User manual - at least that was what I was hoping to do. First I tried to write colormap, with the meta property from the input file as kwarg but I didn't see the colormap appear when I run gdalinfo from the command prompt. I also tried the example with the profile property as kwarg - with photometric set to "RGB" - but I didn't even manage to alter the ColorInterp of the GeoTIFF. The linked documentation seems rather complicated - sorry. Worrying ... – Dobedani Dec 7 '17 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.