Rasterio: writing tiff with cmap on float values

Question

I have this very simple code:

import numpy as np
from rasterio.transform import Affine

nx = 5
maxx = 4.0
minx = -4.0
res = (maxx - minx) / nx
maxy = 3.0
miny = -3.0
ny = int((maxy - miny) / res)

x = np.linspace(minx, maxx, nx)
y = np.linspace(miny, maxy, ny)
z = numpy.array([
    [-1, 10, 15.1, 6.3, 50.4],
    [26.7, -1, 15.7, 40.7, 5],
    [5, -1, 9.0, 38, 40.3],
])
cmap = plt.get_cmap("nipy_spectral")
with rasterio.open(
    os.path.join(os.path.dirname(__file__), "test.tiff"),
    "w",
    driver='GTiff',
    height=z.shape[0],
    width=z.shape[1],
    count=1,
    dtype=z.dtype,
    crs='+proj=latlong',
    transform=Affine.translation(x[0]-res/2, y[0]-res/2) * Affine.scale(res, res),
    nodata=-1,
) as df:
    df.colorinterp = [ColorInterp.palette]
    # df.write_colormap(1, cmap)
    df.write(z, 1)

It create a basic image when drag and drop in QGIS:

I would like to drag and drop this file in Qgis and it have the cmap working from matplotlib named nipy_spectral:

The line # df.write_colormap(1, cmap) is working only for uint8 data (when cmap is a dictionary using int values as keys) according to the documentation, but there is no documentation about float data...

My question and need is simple but there is nothing in documentation: how to apply this cmap to my df rasterio object in the python code?

For the moment it is working when I force data to be uint8 but i can have only 256 values, which is not enough...

.

An other solution is to add manually in qgis a predefined cmap like this:

Then it is possible to export the style as a folder. Maybe it is possible to automatically apply this style to tiff file using qgis.core python module?

It's highly likely that the tiff format itself introduces the restriction that cmap only works with uint8 pixels. — Mark Ransom, Jan 04 '22 at 01:47
Ok, to display only 256 values is not a problem but how to do this storing float values in the tiff file? — Vincent Bénet, Jan 04 '22 at 01:50
And what about using the qgis.core to apply automatically a predefined style for cmap? (Updating question...) ? — Vincent Bénet, Jan 04 '22 at 08:51

score 0 · Accepted Answer · answered Jan 05 '22 at 08:42

I finaly have done the trick using theses lines of code:

import os
import numpy
from qgis.core import (
    QgsCoordinateReferenceSystem, QgsSingleBandPseudoColorRenderer, QgsColorRampShader, QgsStyle, QgsRasterBandStats,
    QgsRasterShader, QgsApplication, QgsProject, QgsRasterLayer)
from calcul import conversion_to_geotiff



def create_project(path_qgz, tiffs=None, epsg=2154):
    qgs = QgsApplication([], False)
    qgs.initQgis()
    project = QgsProject.instance()
    project.setTitle('test')
    project.setCrs(QgsCoordinateReferenceSystem(epsg))
    for data in tiffs:
        x, y, z, path_tif, colormap, precision = data
        conversion_to_geotiff.create_tiff(path_tif, x, y, z)
        layer = QgsRasterLayer(path_tif, os.path.splitext(os.path.basename(path_tif))[0])
        stats = layer.dataProvider().bandStatistics(1, QgsRasterBandStats.All)
        minimum = stats.minimumValue
        maximum = stats.maximumValue
        delta = maximum - minimum
        nclass = max(2, int(delta / precision))
        fractional_steps = [i / (nclass - 1) for i in range(nclass)]
        ramp = QgsStyle().defaultStyle().colorRamp(colormap)
        colors = [ramp.color(f) for f in fractional_steps]
        steps = [minimum + f * delta for f in fractional_steps]
        ramp_items = [
            QgsColorRampShader.ColorRampItem(step, color, str(step))
            for step, color in zip(steps, colors)
        ]
        shader_function = QgsColorRampShader()
        shader_function.setClassificationMode(QgsColorRampShader.EqualInterval)
        shader_function.setColorRampItemList(ramp_items)
        raster_shader = QgsRasterShader()
        raster_shader.setRasterShaderFunction(shader_function)
        renderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), 1, raster_shader)
        layer.setRenderer(renderer)
        layer.triggerRepaint()
        project.addMapLayer(layer)
    project.write(path_qgz)
    qgs.exitQgis()


if __name__ == "__main__":
    x = np.linspace(minx, maxx, nx)
    y = np.linspace(miny, maxy, ny)
    z = numpy.array([
        [-1, 10, 15.1, 6.3, 50.4],
        [26.7, -1, 15.7, 40.7, 5],
        [5, -1, 9.0, 38, 40.3],
    ])

    create_project(
        "/home/vince/test.qgz",
        tiffs=[
            [x, y, z, "/home/vince/test.tif", "Turbo", 5]
        ]
    )

The hard part is to add the qgis module to your PYTHONPATH. It is very hard in windows because of DLLs, much easier in linux. (Just make sure your python version of qgis and your python version of your code is the same.

Rasterio: writing tiff with cmap on float values

1 Answers1