[Qgis-user] [pyqgis] create new raster

Tue May 9 23:42:56 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?
>

Yes, it does and sets the right type Float64. But I am not sure how to do
that if I want the same datatype of the original raster layer.

Interestingly enough:

dataType = dtmLayer.dataProvider().dataType(1)
print(int(dataType))

prints 7, which then results in an error:

File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line
173, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Unable to execute algorithm
Incorrect parameter value for OUTPUT_TYPE

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

Once the 6 is used, then for both properly Float64.

I have some findings on the reaster writing, but that will make more sense
in a different dedicated thread.

Cheers,
Andrea

>
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20230510/d66d2d77/attachment.htm>