[Qgis-user] [pyqgis] create new raster

[andrea antonello]

>
> What do you mean by "new raster"? What is supposed to be in it and how
> should its extents, CRS etc be defined?
>
Well, I might want to create the extents + crs based on a different raster,
but I also might want to create it based on some user input.
For example if I write a module that starts from a vector layer and
generates a raster layer processing it.

> Do you have existing data that you want to use to create it? Or why is a
> constant value raster not what you want?
>
I agree that 99% of the times I will start from an existing raster layer
from which I take extent, resolution and crs.
And I am ok with using the processing toolbox to create a new raster
initialized with novalues.

The problem I posed in the first place (with the code snippet) is that if I
modify the raster using a QgsRasterBlock, the result is not correct.
In the snippet I am taking values from one raster and copy it into the new
raster in certain conditions. But the result is not the expected one.
So my first question was if I am doing something wrong. Because I expect a
portion following isolines to be set to novalue, but instead I get just a
noisy raster layer.

Thanks,
Andrea


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 listQGIS-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 der www.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.comwww.wheregroup.com
> Geschäftsführer:
> Olaf Knopp, Peter Stamm
> Amtsgericht Bonn, HRB 9885
> -------------------------------
>
> _______________________________________________
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20230509/101f493c/attachment.htm>