[Qgis-user] [pyqgis] create new raster

[Richard McDonnell]

Hi Andrea,

If you are ok substituting the method you are using for GDAL, then you can use the following..

gdal_calc.py -A input.tif --outfile=result.tif --calc="(A>=-2000)*(A<=1500)" --NoDataValue=-9999

You can also implement this in QGIS using GDAL Calculator from the Processing Toolbox, It can also be ran as a batch in QGIS or  built into a model if you are looking to batch it or it’s something that’s regularly done.
The output type is Float32, as you are using DTM, probably with decimal places a NoDATA value of -9999

Regards,

Richard



——
Richard McDonnell MSc GIS, FME Certified Professional
FRM Data Management

——
Oifig na nOibreacha Poiblí
Office of Public Works

Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
Jonathan Swift Street, Trim, Co Meath, C15 NX36
——
M +353 87 688 5964 T +353 46 942 2409
https://gov.ie/opw

——
To send me files larger than 30MB, please use the link below https://filetransfer.opw.ie/filedrop/richard.mcdonnell@opw.ie

Email Disclaimer: https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/

——
MSc GIS, FME Certified Professional

——
Oifig na nOibreacha Poiblí
Office of Public Works

Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
Jonathan Swift Street, Trim, Co Meath, C15 NX36
——
M +353 87 688 5964 T +353 46 942 2409
https://https://gov.ie/opw<https://www.opw.ie>

——
Email Disclaimer: https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/<https://www.opw.ie/en/disclaimer/>
From: QGIS-User <qgis-user-bounces at lists.osgeo.org> On Behalf Of andrea antonello via QGIS-User
Sent: 09 May 2023 07:28
To: qgis-user at lists.osgeo.org
Subject: [Qgis-user] [pyqgis] create new raster

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/bd710c4e/attachment-0001.htm>