[Qgis-user] [pyqgis] create new raster

Tue May 9 09:36:08 PDT 2023

What kind of error do you get if you uncomment your datatype line? What 
value does that datatype have? Does it work if you pass the number 6 as 
datatype?

What does the QGIS GUI say about the original and the created layers 
datatypes in their properties dialogs?

Am 2023-05-09 16:23, schrieb andrea antonello via QGIS-User:
>> 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 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 der www.foss-academy.com [1]
>> ---------------------------------------------
>> 
>> 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 [2]
>> 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
> 
> 
> Links:
> ------
> [1] http://www.foss-academy.com
> [2] http://www.wheregroup.com
> _______________________________________________
> 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