python-gdal write a GeoTIFF in colour from a data array

Question

I have a 14400 x 7200 (360 x 180 at .025 deg resolution) numpy array containing data which I want to make into a GeoTIFF. I have done the following which creates a GeoTIFF fine:

# Set file vars
output_file = "out.tif"
# Create gtif
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, WIDTH, HEIGHT, 1, gdal.GDT_Byte )
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform( [ -180, 0.025, 0, 90, 0, -0.025 ] )
# set the reference info 
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection( srs.ExportToWkt() )
# write the band
dst_ds.GetRasterBand(1).SetNoDataValue(255) #I set my nodata values in array to be 255
dst_ds.GetRasterBand(1).WriteArray(data_array)

But the resulting tiff when opened looks to be binary in colour (black & white). How can I assign a colourmap to this image? Even just in grayscale? My values have a small range (1.7 - 1.9). Should I just remap this to 0-255?

Edit:

I ended up scaling my data array from ~ 1.7-1.9 to 0-255 and that worked nicely.

import numpy as np
data_array_scaled = np.interp(data_array, (data_array.min(), data_array.max()), (0, 255))

This gives me grayscale only though -- how do I generate in colour?

Edit 2: Updated code

import numpy as np
import os
from osgeo import osr, gdal
data_array = np.load('fini.tif.npy')
data_array = data_array[::-1]
data_array_scaled = np.interp(data_array, (data_array.min(), data_array.max()), (0, 255))
data_array_scaled[data_array_scaled == 156.91450597804007] = -9999
r = data_array_scaled
RES = 0.025
WIDTH = data_array_scaled.shape[1]
HEIGHT = data_array_scaled.shape[0]
output_file = "out.tif"
# Create GeoTIFF
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, WIDTH, HEIGHT, 1, gdal.GDT_Byte)
# Upper Left x, Eeast-West px resolution, rotation, Upper Left y, rotation, North-South px resolution
dst_ds.SetGeoTransform( [ -180, 0.025, 0, 90, 0, -0.025 ] )
# Set CRS
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection( srs.ExportToWkt() )
# Write the band
dst_ds.GetRasterBand(1).SetNoDataValue(-9999) #optional if no-data transparent
dst_ds.GetRasterBand(1).WriteArray(r)
band = dst_ds.GetRasterBand(1)
colors = gdal.ColorTable()
colors.SetColorEntry(1, (112, 153, 89))
colors.SetColorEntry(2, (242, 238, 162))
colors.SetColorEntry(3, (242, 206, 133))
colors.SetColorEntry(4, (194, 140, 124))
colors.SetColorEntry(5, (214, 193, 156))
band.SetRasterColorTable(colors)
band.SetRasterColorInterpretation(gdal.GCI_PaletteIndex)
del band, dst_ds

• user30184: You can write a true color image with three bands RGB or otherwise you must attach a color palette to the single band image.
• mathwiz97: @user30184 How do I attach a color palette?
• Marcelo Villa: Might want to create a color table: gis.stackexchange.com/a/325751/86131
• mathwiz97: @MarceloVilla Thanks! that looks promising but I am getting the following error: ERROR 6: SetColorTable() not supported for multi-sample TIFF files.
• Marcelo Villa: Could'nt find too much information on the error besides the source code that triggers it. Are you sure you are applying it to a one band raster? Can you update your question with the code you are using that yields this error?

Answer 1

You can create a color ramp for the color table instead of having to specify a color for every single value in the raster. The CreateColorRamp() method takes four arguments:

• Start value
• Start color
• Stop value
• Stop color

For example, imagine you want to build a color ramp with three colors (e.g. three different shades of green). The first green will be for the value 0, the second for the value 127 and the third for the value 254. The colors for the rest of the values will be interpolated colors between two of the already mentioned colors. Feel free to toy with this example and add as many stops as you want and change the colors depending on your needs.

colors = gdal.ColorTable()
colors.CreateColorRamp(0, (229, 245, 249), 127, (153, 216, 201))
colors.CreateColorRamp(127, (153, 216, 201), 254, (44, 162, 95))
band.SetRasterColorTable(colors)
band.SetRasterColorInterpretation(gdal.GCI_PaletteIndex)

---

On another note, it doesn't make sense to set the NoData value to -9999 if you have a TIFF file with gdal.GDT_Byte as its data type. A byte can only store positive integers between 0 and 255. Consider using another value for NoData or changing the data type of the TIFF file so it can actually store -9999.

• mathwiz97: Thank you, CreateColorRamp() is the method I was looking for and didn't know existed! Also thanks for pointing out the NoData issue -- I wasn't exactly sure what GDT_byte was, but I've changed my NoData value to 0.