I have a large number of single band rasters that are colour indexed TIF files. I want to mosaic them using GDAL (or similar non-proprietary software with a Python API I can call as part of a bigger script).

The images all represent the same sort of thing but the information is embedded in the colour tables for each image and not the raster pixel value as such (i.e. all the colour tables are independent). So, sky-blue (say) always means the same thing but does not always have the same position in the different colour tables.

Mosaicing them is easy as a greyscale image, but then I lose the colour information which is vital. I can't use the '-pct' switch in gdal_merge.py because that only grabs the first colour table and causes some or all of the other tiles in the mosaic to display incorrectly.

My gut feeling is that there is nothing for it but to convert all the images to three-band RGB (and consequently balooning the data volume) and mosaic the images with a possible post process step to create a new unified colour table and swap the mode back to single band indexed colour.

However, I'd value opinion and work-flow suggestions on this (the process must be automatable). Is there a better way of doing it?

To be clear, I have no way of knowing whether any particular tile contains the full range of colours represented by all the tiles together.


If you don't mind getting your hands dirty with some Python, you could write a script that uses the GDAL Python bindings to do most of the heavy lifting.

The broad steps I'd use are:

  • Load the master raster and extract its colour table as a simple Python list.
  • Load the source raster, and for each entry in its colour table search for a corresponding colour in the source image's table.
  • Create a Python dictionary mapping source to master index values.
  • Use RasterIO to read in the data from the source (one block at a time if it's quite large), and write the data out as a new image using map() and a lambda function to do the index mapping.
  • Rinse and repeat either within the Python code, or externally as a batch file.

The biggest issue I can see is what to do if a colour in the source image doesn't exist in the master palette. I'd be inclined to use 0 or 255 and throw up a warning on the command line.

  • This is an interesting idea. 0 to 255 values should be ample. Can you define what you mean by 'master raster' though. If you mean the first raster in my mosaic, then your procedure surely won't help me much more than using the -pct switch in gdal_merge because, if that first tile happens to be sparsely populated, it's colour table will lack colours in subsequent tables. I think I would have to refine your procedure to create the dictionary and then extend the dictionary with colours I don't find in the 'master raster' but do find in other tiles. (note my edit to my question). Jun 7 '12 at 13:12
  • 2
    Ah yes, I see, I assumed there might be a file with all colours, but not necessarily on the right order. You're right, start off with an empty master palette, and as a file is opened any colours that don't exist, append to the bottom of the palette list - and hope you don't run out of entries! Alternatively, do a two-pass scan, where the first pass just gathers the colours. That'd be ideal for bounds checking, and if you want to ensure the colours are in some specific order. Jun 7 '12 at 13:21
  • 1
    Yes, I think that would work! Bounds checking is probably a good idea, though I'm pretty sure there will only be 0-255 colours though... but you never know. Jun 7 '12 at 13:43
  • +1 The two-pass approach sounds like a great idea. It could even be refined to create an optimal mapping from the input color tables to a master table for the mosaic, in case some slight changes in colors could be tolerated. This approach cleanly separates the (analytical, statistical) work of creating an output color table (and maps from the input to output tables) from the work of making the mosaic itself. This leads to a flexible engineering design and workflow. E.g., the color tables could be computed with suitable software (such as R) while the mosaicing could be done in GIS.
    – whuber
    Jun 7 '12 at 20:17
  • Yes. If necessary, I can also add a routine to handle dithering from, say, lossy compression (e.g. if the images came as JPEG) where I know that values must fall within certain ranges. This is easy for discrete data but not so easy for continuous data. However, where white equates to no-data, I can at least set limits on the colour bleeding. Jun 8 '12 at 9:54

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.