[GIS] python-gdal write a GeoTIFF in colour from a data array

colormapgdalgeotiff-tiffpython

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

Best Answer

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.

Related Question