> > 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 listqgis-u...@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 > --------------------------------------------- > > 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.comwww.wheregroup.com > 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 >
_______________________________________________ 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