Rendering single-band raster in Python with set color ramp

Question

I have a set of single-band gray-scale rasters that I want to save as a rendered GeoTIFF/image that is colored according to a specific color ramp.

For example, if I had some NDVIs that range from a minimum and maximum value of -1 to 1 and I want them colored according to a ramp that would be similar to this: 0 (red) -> 0.5 (yellow) -> 1 (green).

I know how to do this within QGIS and ArcMap but I am looking to do it straight up in open source Python. I have looked into gdal and rasterio but I don't think I've found exactly what I need. Especially when it comes to saving them with the applied color ramp and forcing a min/max value of the ramp.

Are there any resources I can be pointed to?

Answer 1

import rasterio
import matplotlib
from matplotlib import pyplot as plt
import matplotlib.colors
rasterfile = r"/home/bera/Desktop/gistest/nvdi.tif"
dataset = rasterio.open(rasterfile, mode="r")
array = dataset.read(1)
norm=plt.Normalize(0, 0.6) #My nvdi have values from ~ 0-0.6
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", ["red","yellow","limegreen"])
plt.imshow(array, cmap=cmap, norm=norm)
plt.colorbar()

enter image description here

Answer 2

Follow the very very recent development in here https://github.com/OSGeo/gdal/pull/3133. The pull request has been submitted two days ago but is has not been accepted yet. From the description of the PR

What does this PR do?

This PR adds option to attach color table from a supported file and to more extent options (union, intersection, custom). The related functions were implemented in #3131.

Because the pull request contains new complete python scripts you may be able to just capture them and try right away. However, they may require some components from the previous PR from the same author.

Answer 3

I think using gdaldem solved my issue for now using the color-relief option

https://gdal.org/programs/gdaldem.html

Answer 4

Look this:

ndvi = (nir - red) / (nir + red)
ndvi[ndvi > 1] = 1
ndvi[ndvi < -1] = -1
ndvi