<div dir="ltr">Hello,<br><div>I am trying to find out the best workflow to create a new raster.</div><div>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.</div><div><br></div><div>This is the only way I found to do so:</div><div><br></div><div><div style="color:rgb(0,0,0);font-family:"Droid Sans Mono","monospace",monospace;font-size:14px;line-height:19px;white-space:pre"><div><span style="color:rgb(0,128,0)"># create new raster with novalues between 1500 and 2000</span></div><div><span style="color:rgb(0,16,128)">dataType</span> = <span style="color:rgb(0,16,128)">dtmLayer</span>.<span style="color:rgb(121,94,38)">dataProvider</span>().<span style="color:rgb(121,94,38)">dataType</span>(<span style="color:rgb(9,134,88)">1</span>)</div><div><span style="color:rgb(0,16,128)">crs</span> = <span style="color:rgb(38,127,153)">QgsCoordinateReferenceSystem</span>(<span style="color:rgb(163,21,21)">'EPSG:3003'</span>)</div><div><span style="color:rgb(0,16,128)">params</span> = {</div><div>    <span style="color:rgb(163,21,21)">'EXTENT'</span>: <span style="color:rgb(0,16,128)">dtmLayer</span>.<span style="color:rgb(121,94,38)">extent</span>(),</div><div>    <span style="color:rgb(163,21,21)">'TARGET_CRS'</span>: <span style="color:rgb(0,16,128)">crs</span>,</div><div>    <span style="color:rgb(163,21,21)">'PIXEL_SIZE'</span>: <span style="color:rgb(0,16,128)">dtmLayer</span>.<span style="color:rgb(121,94,38)">rasterUnitsPerPixelX</span>(),</div><div>    <span style="color:rgb(163,21,21)">'NUMBER'</span>: -<span style="color:rgb(9,134,88)">9999.0</span>,</div><div>    <span style="color:rgb(0,128,0)"># 'OUTPUT_TYPE': dataType,</span></div><div>    <span style="color:rgb(163,21,21)">'OUTPUT'</span>: <span style="color:rgb(38,127,153)">QgsProcessing</span>.<span style="color:rgb(0,16,128)">TEMPORARY_OUTPUT</span></div><div>}</div><div><span style="color:rgb(0,16,128)">newRaster</span> = <span style="color:rgb(38,127,153)">processing</span>.run(<span style="color:rgb(163,21,21)">'qgis:createconstantrasterlayer'</span>, <span style="color:rgb(0,16,128)">params</span>)[<span style="color:rgb(163,21,21)">'OUTPUT'</span>]</div><div><span style="color:rgb(0,16,128)">newRasterLayer</span> = <span style="color:rgb(38,127,153)">QgsRasterLayer</span>(<span style="color:rgb(0,16,128)">newRaster</span>, <span style="color:rgb(163,21,21)">'temp'</span>, <span style="color:rgb(163,21,21)">'gdal'</span>)</div><div><span style="color:rgb(0,16,128)">newRasterProvider</span> = <span style="color:rgb(0,16,128)">newRasterLayer</span>.<span style="color:rgb(121,94,38)">dataProvider</span>()</div><br><div><span style="color:rgb(0,16,128)">block</span> = <span style="color:rgb(38,127,153)">QgsRasterBlock</span>(<span style="color:rgb(0,16,128)">dataType</span>, <span style="color:rgb(0,16,128)">cols</span>, <span style="color:rgb(0,16,128)">rows</span>)</div><div><br></div><div><span style="color:rgb(175,0,219)">for</span> <span style="color:rgb(0,16,128)">row</span> <span style="color:rgb(0,0,255)">in</span> <span style="color:rgb(38,127,153)">range</span>(<span style="color:rgb(0,16,128)">rows</span>):</div><div>    <span style="color:rgb(175,0,219)">for</span> <span style="color:rgb(0,16,128)">col</span> <span style="color:rgb(0,0,255)">in</span> <span style="color:rgb(38,127,153)">range</span>(<span style="color:rgb(0,16,128)">cols</span>):</div><div>        <span style="color:rgb(0,16,128)">point</span> = <span style="color:rgb(0,16,128)">dtmLayer</span>.<span style="color:rgb(121,94,38)">dataProvider</span>().<span style="color:rgb(121,94,38)">transformCoordinates</span>(<span style="color:rgb(38,127,153)">QgsPoint</span>(<span style="color:rgb(0,16,128)">col</span>, <span style="color:rgb(0,16,128)">row</span>), <span style="color:rgb(0,16,128)">transformType</span>)</div><div>        <span style="color:rgb(0,16,128)">value</span>, <span style="color:rgb(0,16,128)">res</span> = <span style="color:rgb(0,16,128)">dtmLayer</span>.<span style="color:rgb(121,94,38)">dataProvider</span>().<span style="color:rgb(121,94,38)">sample</span>(<span style="color:rgb(38,127,153)">QgsPointXY</span>(<span style="color:rgb(0,16,128)">point</span>.<span style="color:rgb(121,94,38)">x</span>(), <span style="color:rgb(0,16,128)">point</span>.<span style="color:rgb(121,94,38)">y</span>()), <span style="color:rgb(9,134,88)">1</span>)</div><div>        </div><div>        <span style="color:rgb(175,0,219)">if</span> <span style="color:rgb(0,16,128)">res</span> <span style="color:rgb(0,0,255)">and</span> <span style="color:rgb(0,16,128)">value</span> != -<span style="color:rgb(9,134,88)">9999.0</span>:</div><div>            <span style="color:rgb(175,0,219)">if</span> <span style="color:rgb(0,16,128)">value</span> < <span style="color:rgb(9,134,88)">1000</span> <span style="color:rgb(0,0,255)">or</span> <span style="color:rgb(0,16,128)">value</span> > <span style="color:rgb(9,134,88)">2000</span>:</div><div>                <span style="color:rgb(0,16,128)">block</span>.<span style="color:rgb(121,94,38)">setValue</span>(<span style="color:rgb(0,16,128)">row</span>, <span style="color:rgb(0,16,128)">col</span>, <span style="color:rgb(0,16,128)">value</span>)</div><div><br></div><div><span style="color:rgb(0,16,128)">newRasterProvider</span>.<span style="color:rgb(121,94,38)">setEditable</span>(<span style="color:rgb(0,0,255)">True</span>)</div><div><span style="color:rgb(0,16,128)">newRasterProvider</span>.<span style="color:rgb(121,94,38)">writeBlock</span>(<span style="color:rgb(0,16,128)">block</span>, <span style="color:rgb(0,16,128)">band</span>=<span style="color:rgb(9,134,88)">1</span>)</div><div><span style="color:rgb(0,16,128)">newRasterProvider</span>.<span style="color:rgb(121,94,38)">setEditable</span>(<span style="color:rgb(0,0,255)">False</span>)</div></div></div><div><br></div><div><br></div><div>This code has two main issues:</div><div>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.</div><div>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).<br>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. </div><div><br></div><div>Has anyone a hint about what I am doing wrong? </div><div><br></div><div>Thanks,</div><div>Andrea</div><div><br></div><div><br></div></div>