Re: [Qgis-user] [pyqgis] create new raster

WhereGroup Tue, 09 May 2023 09:36:29 -0700

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

What does the QGIS GUI say about the original and the created layersdatatypes 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@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.kroe...@wheregroup.com
www.wheregroup.com [2]
Geschäftsführer:
Olaf Knopp, Peter Stamm
Amtsgericht Bonn, HRB 9885
-------------------------------

_______________________________________________
QGIS-User mailing list
QGIS-User@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@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

_______________________________________________
QGIS-User mailing list
QGIS-User@lists.osgeo.org
List info: https://lists.osgeo.org/mailman/listinfo/qgis-user
Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user

Re: [Qgis-user] [pyqgis] create new raster

Reply via email to