Changing color of raster images based on their data values - gdal [duplicate]

Question

Let's say I have 1 band raster images with Gray Color Interpretation and I want to change their color based on their data values. Data values -1, will be Blue, -0.5 will be Light Blue, 0 will be Yellow, 0.5 will be Orange and 1 will be Red. If the data values are <-3000, I want to have that colored Black and if the data values are >3.402, that should be colored White . Will that be possible in gdal? And process that in batch?

Answer 1

c1 to c4 in the vrt file are the RGBA values you see in a 4-channel RGBA tif file. Every row corresponds to an integer palette value in the source file.

---

EDIT

gdaldem seems to be the right tool for your task:

1. Take a raster of your choice (like SRTM as an example)
2. Create the following text file named col.txt:

0 black

100 blue

200 yellow

300 orange

400 red

500 white

3. On command line, run

gdaldem color-relief N51E007.hgt col.txt out.tif

and the result will look like this:

I have added white dashed contour lines to show that the colours get interpolated.

Instead of colour names, you can use RGB values as well, see http://www.gdal.org/gdaldem.html#gdaldem_color_relief for details:

The supported list is : white, black, red, green, blue, yellow, magenta, cyan, aqua, grey/gray, orange, brown, purple/violet and indigo

You can create a vrt file as output too and look into that on how the interpolating is done with <LUT> keys:

gdaldem color-relief -of VRT N51E007.hgt col.txt out.vrt

  ...
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Red</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">N51E007.hgt</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1201" RasterYSize="1201" DataType="Int16" BlockXSize="1201" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="1201" ySize="1201"/>
      <DstRect xOff="0" yOff="0" xSize="1201" ySize="1201"/>
      <LUT>0:0,100:0,200:0,300:255,400:255,500:255</LUT>
    </ComplexSource>
  </VRTRasterBand>

---

update

To apply different transparencies, I got around with this batch:

gdaldem color-relief -of GTiff N51E006.hgt color.txt output6.tif -alpha 
gdaldem color-relief -of GTiff N51E007.hgt color.txt output7.tif -alpha 

and this color.txt:

0. 255 255 255 0
50. 255 0 0 50
100. 0 255 0 100
150. 0 0 255 150

with elevation in meters, and RGBA values between 0 and 255.

Answer 2

Here is a auto color script based on AndreJ's code:

================================================================================
Usage:
-------------
Custom color:
python gdaldem.py input_tif.tif color.txt output_color.tif
Auto color:
python gdaldem.py input_tif.tif auto output_color.tif

"""
import subprocess
import sys
import os
import tempfile
import numpy as np


def main(input_file, color_file, output_file):
    cmd = "gdaldem color-relief " + input_file \
          + ' ' + color_file + ' ' + output_file
    subprocess.check_call(cmd, shell=True)


def gen_color_file(input_file):
    fp, temp_file = tempfile.mkstemp(suffix='.txt')

    dem = DEM(input_file)
    dem.open()
    phase_data = dem.height_band.ReadAsArray()

    max_ph = np.nanmax(phase_data)
    min_ph = np.nanmin(phase_data)
    range_ph = max_ph-min_ph
    colors = ['black', 'blue', 'yellow', 'orange', 'red', 'white']
    with open(temp_file, 'w') as f:
        for i, c in enumerate(colors[:-1]):
            f.write(str(int(min_ph + (i + 1)*range_ph/len(colors))) +
                    ' ' + c + '\n')
        f.write(str(int(max_ph - range_ph/len(colors))) +
                ' ' + colors[-1] + '\n')
    os.close(fp)
    return temp_file


if __name__ == '__main__':
    input_file = sys.argv[1]
    color_file = sys.argv[2]
    output_file = sys.argv[3]
    if color_file == 'auto':
        print '\nauto generating color file'
        color_file = gen_color_file(input_file)
    with open(color_file, 'r') as f:
        print '\ncolor file contents'
        print '='*50
        for l in f.readlines():
            print l
    print '='*50
    main(input_file, color_file, output_file)