[Qgis-user] [pyqgis] create new raster

[andrea antonello]

Hello,
I am trying to find out the best workflow to create a new raster.
As an example I take an existing elevation model raster and loop over it to
set values between 1500 and 2000 to novalue and write the result to a new
raster.

This is the only way I found to do so:

# create new raster with novalues between 1500 and 2000
dataType = dtmLayer.dataProvider().dataType(1)
crs = QgsCoordinateReferenceSystem('EPSG:3003')
params = {
'EXTENT': dtmLayer.extent(),
'TARGET_CRS': crs,
'PIXEL_SIZE': dtmLayer.rasterUnitsPerPixelX(),
'NUMBER': -9999.0,
# 'OUTPUT_TYPE': dataType,
'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
}
newRaster = processing.run('qgis:createconstantrasterlayer', params)[
'OUTPUT']
newRasterLayer = QgsRasterLayer(newRaster, 'temp', 'gdal')
newRasterProvider = newRasterLayer.dataProvider()

block = QgsRasterBlock(dataType, cols, rows)

for row in range(rows):
for col in range(cols):
point = dtmLayer.dataProvider().transformCoordinates(QgsPoint(col, row),
transformType)
value, res = dtmLayer.dataProvider().sample(QgsPointXY(point.x(), point.y()),
1)
if res and value != -9999.0:
if value < 1000 or value > 2000:
block.setValue(row, col, value)

newRasterProvider.setEditable(True)
newRasterProvider.writeBlock(block, band=1)
newRasterProvider.setEditable(False)


This code has two main issues:
1. if I uncomment the line containing OUTPUT_TYPE, I am getting an error
about the type passed. But I can't find the right type needed there, it
should be the one taken from the original provider.
2. the resulting raster is scrambled as if there was a shift in the setting
of the values. But the QgsRasterBlock seems to be built correctly (rows,
cols) and the values set properly (col, row).
3. in the above example, the dtmLayer has an epsg 3033 crs and when loaded
manually into QGIS, it is recognized. But when I read the layer's metadata
crs with pyQGIS , it is not able to read it and tells me it is invalid.

Has anyone a hint about what I am doing wrong?

Thanks,
Andrea
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20230509/caaeedf6/attachment.htm>