[Qgis-user] [pyqgis] create new raster

[Johannes Kröger (WhereGroup)]

What do you mean by "new raster"? What is supposed to be in it and how 
should its extents, CRS etc be defined?

Do you have existing data that you want to use to create it? Or why is a 
constant value raster not what you want?

If you want to base it on an existing raster, best use a raster 
calculator to create a plain copy.

Cheers, Hannes

Am 09.05.23 um 08:28 schrieb andrea antonello via QGIS-User:
> 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
>
>
>
> _______________________________________________
> QGIS-User mailing list
> QGIS-User at lists.osgeo.org
> List info:https://lists.osgeo.org/mailman/listinfo/qgis-user
> Unsubscribe:https://lists.osgeo.org/mailman/listinfo/qgis-user

-- 
Johannes Kröger / GIS-Entwickler/-Berater

---------------------------------------------
Aufwind durch Wissen!
Web-Seminare und Online-Schulungen
bei derwww.foss-academy.com
---------------------------------------------

WhereGroup GmbH
c/o KK03 GmbH
Lange Reihe 29
20099 Hamburg
Germany

Tel: +49 (0)228 / 90 90 38 - 36
Fax: +49 (0)228 / 90 90 38 - 11

johannes.kroeger at wheregroup.com
www.wheregroup.com
Geschäftsführer:
Olaf Knopp, Peter Stamm
Amtsgericht Bonn, HRB 9885
-------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20230509/8292cb5c/attachment.htm>