Rasterio：在浮点值上用 cmap 写 tiff

Question

我有这个非常简单的代码：

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)

在 QGIS 中拖放时创建基本图像：

我想将这个文件拖放到 Qgis 中，它有来自 matplotlib 的 cmap，名为 nipy_spectral:

行 # df.write_colormap(1, cmap) 仅适用于 uint8 数据（当 cmap 是使用 int 值作为键的字典时）according to the documentation，但没有关于浮动数据...

我的问题和需求很简单，但文档中没有任何内容：如何在 python 代码中将此 cmap 应用于我的 df rasterio 对象？

目前当我强制数据为 uint8 时它正在工作，但我只能有 256 个值，这还不够...

.

另一个解决方案是在 qgis 中手动添加一个预定义的 cmap，如下所示：

然后可以将样式导出为文件夹。也许可以使用 qgis.core python 模块自动将此样式应用于 tiff 文件？

Answer 1

我终于用这些代码行完成了这个技巧：

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

困难的部分是将 qgis 模块添加到您的 PYTHONPATH。由于 DLL，在 windows 中很难，在 linux 中容易得多。（只需确保您的 python 版本的 qgis 和您的 python 版本的代码相同。

Rasterio：在浮点值上用 cmap 写 tiff

Rasterio: writing tiff with cmap on float values

python

matplotlib

tiff

qgis

rasterio