[Qgis-user] [pyqgis] create new raster

Tue May 9 04:53:43 PDT 2023

Hi and thanks for the replies.

As for my other email related to reading a raster, I wanted to investigate
the possibility of creating a raster using pure pyQGIS API.
It seems that it isn't possible without using at least the
constantmap processing algo.
But I was hoping that it would be possible to use at least the Block class
to work the raster as a matrix.

I was about to open an issue for the thing, but suspect I am using the API
in a wrong way, since the results are scrambled.

Cheers,
Andrea

On Tue, May 9, 2023 at 10:57 AM Richard McDonnell <richard.mcdonnell at opw.ie>
wrote:

> Hi Andrea,
>
> Please disregard my last email the formula is incorrect.
>
> Apologies..
>
>
>
> Richard
>
>
>
> ——
> Richard McDonnell MSc GIS, FME Certified Professional
> *FRM Data Management*
>
> ——
> Oifig na nOibreacha Poiblí
> Office of Public Works
>
> Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
> Jonathan Swift Street, Trim, Co Meath, C15 NX36
> ——
> M +353 87 688 5964 T +353 46 942 2409
> https://gov.ie/opw
>
> ——
> To send me files larger than 30MB, please use the link below
> https://filetransfer.opw.ie/filedrop/richard.mcdonnell@opw.ie
>
> Email Disclaimer:
> https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/
>
> ——
> MSc GIS, FME Certified Professional
>
> ——
> Oifig na nOibreacha Poiblí
> Office of Public Works
>
> Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36
> Jonathan Swift Street, Trim, Co Meath, C15 NX36
> ——
> M +353 87 688 5964 T +353 46 942 2409
> https://https://gov.ie/opw <https://www.opw.ie>
>
> ——
> Email Disclaimer:
> https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/
> <https://www.opw.ie/en/disclaimer/>
>
> *From:* QGIS-User <qgis-user-bounces at lists.osgeo.org> *On Behalf Of *Richard
> McDonnell via QGIS-User
> *Sent:* 09 May 2023 09:27
> *To:* andrea antonello <andrea.antonello at gmail.com>
> *Cc:* qgis-user at lists.osgeo.org
> *Subject:* Re: [Qgis-user] [pyqgis] create new raster
>
>
>
> Hi Andrea,
>
>
>
> If you are ok substituting the method you are using for GDAL, then you can
> use the following..
>
>
>
> gdal_calc.py -A input.tif --outfile=result.tif
> --calc="(A>=-2000)*(A<=1500)" --NoDataValue=-9999
>
>
>
> You can also implement this in QGIS using GDAL Calculator from the
> Processing Toolbox, It can also be ran as a batch in QGIS or  built into a
> model if you are looking to batch it or it’s something that’s regularly
> done.
>
> The output type is Float32, as you are using DTM, probably with decimal
> places a NoDATA value of -9999
>
>
>
> Regards,
>
>
>
> Richard
>
>
>
>
>
>
> * ——*
> *Richard McDonnell MSc GIS, FME Certified Professional*
> *FRM Data Management*
>
> ——
> *Oifig na nOibreacha Poiblí*
> Office of Public Works
>
> *Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36 *
> Jonathan Swift Street, Trim, Co Meath, C15 NX36
> ——
> M +353 87 688 5964 T +353 46 942 2409
> https://gov.ie/opw
>
> ——
> To send me files larger than 30MB, please use the link below
> https://filetransfer.opw.ie/filedrop/richard.mcdonnell@opw.ie
>
> Email Disclaimer:
> https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/
>
>
>
>
> * ——*
> *MSc GIS, FME Certified Professional*
>
> ——
> *Oifig na nOibreacha Poiblí*
> Office of Public Works
>
> *Sráid Jonathan Swift, Baile Átha Troim, Co na Mí, C15 NX36 *
> Jonathan Swift Street, Trim, Co Meath, C15 NX36
> ——
> M +353 87 688 5964 T +353 46 942 2409
> https://https://gov.ie/opw <https://www.opw.ie>
>
> ——
> Email Disclaimer:
> https://www.gov.ie/en/organisation-information/439daf-email-disclaimer/
> <https://www.opw.ie/en/disclaimer/>
>
> *From:* QGIS-User <qgis-user-bounces at lists.osgeo.org> *On Behalf Of *andrea
> antonello via QGIS-User
> *Sent:* 09 May 2023 07:28
> *To:* qgis-user at lists.osgeo.org
> *Subject:* [Qgis-user] [pyqgis] create new raster
>
>
>
> 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
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-user/attachments/20230509/bce757e1/attachment-0001.htm>