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 <[email protected]> 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/[email protected] > > 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 <[email protected]> *On Behalf Of *Richard > McDonnell via QGIS-User > *Sent:* 09 May 2023 09:27 > *To:* andrea antonello <[email protected]> > *Cc:* [email protected] > *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/[email protected] > > 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 <[email protected]> *On Behalf Of *andrea > antonello via QGIS-User > *Sent:* 09 May 2023 07:28 > *To:* [email protected] > *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 > > > > >
_______________________________________________ QGIS-User mailing list [email protected] List info: https://lists.osgeo.org/mailman/listinfo/qgis-user Unsubscribe: https://lists.osgeo.org/mailman/listinfo/qgis-user
