Rendering single-band raster in Python with set color ramp

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()

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.

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

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

Look this:

ndvi = (nir - red) / (nir + red)

ndvi[ndvi > 1] = 1
ndvi[ndvi < -1] = -1
ndvi