0

I am trying to use the method writeRaster from qgis.core.writeRaster to create a singleBand raster of float and Nans but according to the documentation, I need to provide theses inputs:

writeRaster(
    self,  # OK
    pipe: QgsRasterPipe,  # Q1
    nCols: int,  # OK
    nRows: int,  # OK
    outputExtent: QgsRectangle,  # Q2
    crs: QgsCoordinateReferenceSystem,  # OK
    feedback: QgsRasterBlockFeedback = None  # OK
) → QgsRasterFileWriter.WriterError

I have 2 questions here:

  • Q1: What is a QgsRasterPipe, how to use it and what is its purpose? The documentation says: Constructor for QgsRasterPipe. Base class for processing modules. Few examples online of writeRaster just initialize this object. So what do I need to provide in the argument pipe ?

  • Q2: The argument outputExtent of type QgsRectangle seems to be the bounding area of my raster: QgsRectangle(x_min, y_min, x_max, y_max). But here is my question: Where do I declare the values of pixels?

Here is the script (not working) I have for the moment:

import os
import numpy
from qgis.core import (
    QgsMapLayer,
    QgsRasterFileWriter,
    QgsCoordinateReferenceSystem,
    QgsRasterPipe,
)


def write_to_geotiff(data: list, filename: str, epsg: str, layer: str=None) -> None:
    
    x_data = data[0]
    y_data = data[1]
    z_data = data[2]
    
    nx, ny = len(x_data), len(y_data)
    
    QgsRasterFileWriter.writeRaster(
        QgsRasterPipe(),
        nCols=nx,
        nRows=ny,
        QgsRectangle(
            min(x_data),
            min(y_data),
            max(x_data),
            max(y_data)
        ),
        crs = QgsCoordinateReferenceSystem(f"epsg:{epsg}"),
    )
    

if __name__ == "__main__":
    filename = r"C:\Users\vince\Downloads\test.gpkg"
    x_data = numpy.asarray([0, 1, 2])
    y_data = numpy.asarray([0, 1])
    z_data = numpy.asarray([
        [0.1, numpy.nan],
        [0.5, 139.5],
        [150.98, numpy.nan],
    ])
    epsg = "4326"
    write_to_geotiff(
        [x_data, y_data, z_data], 
        filename, 
        epsg
    )

I saw this answer for Q1, the data is in the pipe variable. But I don t know how to create a qgsRasterBlock from my numpy array...

Vincent Bénet
  • 1,212
  • 6
  • 20

1 Answers1

0

I get it using the method QgsRasterFileWriter.createOneBandRaster creating a provider. You can get the bloc of the provider of type QgsRasterBlock and use the method setValue to associate values.

    writer = QgsRasterFileWriter(filename)
    provider = QgsRasterFileWriter.createOneBandRaster(
        writer,
        dataType=Qgis.Float32,
        width=nx,
        height=ny,
        extent=QgsRectangle(
            min(x_data),
            min(y_data),
            max(x_data),
            max(y_data)
        ),
        crs = QgsCoordinateReferenceSystem(f"epsg:{epsg}"),
    )
    provider.setNoDataValue(1, -1)
    provider.setEditable(True)
    
    block = provider.block(
        bandNo=1,
        boundingBox=provider.extent(),
        width=provider.xSize(),
        height=provider.ySize()
    )
    
    for ix in range(nx):
        for iy in range(ny):
            value = z_data[ix][iy]
            if value == numpy.nan:
                continue
            block.setValue(iy, ix, value)
    
    provider.writeBlock(
        block=block, 
        band=1,
        xOffset=0,
        yOffset=0
    )
    provider.setEditable(False)
    

This will create a tiffile:

enter image description here

Vincent Bénet
  • 1,212
  • 6
  • 20