<div dir="ltr" gmail_original="1">Instead of manually iterating over the raster, you can use the 'Raster Calculator' Algorithm to set the pixel values which will be much easier and will resolve your issues. Test it with the processing toolbox and then all via PyQGIS with the required parameters.<div><br></div><div>I tested and the following expression works. Replace <span style="white-space: pre-wrap;">N27E086 with the name of your layer</span></div><div><span style="white-space: pre-wrap;"><br></span></div><div><p style="margin: 0px; white-space: pre-wrap;">"N27E086@1"*("N27E086@1"<1500 or "N27E086@1">2000) + -9999*("N27E086@1">=1500 and "N27E086@1"<=2000)</p><div><br></div><div><br clear="all"><div><div dir="ltr" class="gmail_signature"><div dir="ltr">---<div>Ujaval Gandhi</div><div>Spatial Thoughts</div><div><a href="https://mailtrack.io/trace/link/c359cb3ba6621429914ba43aa93bbeae2d52326f?w=cWdpcy11c2VyQGxpc3RzLm9zZ2VvLm9yZw&url=http%3A%2F%2Fwww.spatialthoughts.com&userId=8747767&signature=aae92c96b65f2608" target="_blank">www.spatialthoughts.com</a></div><div><br></div></div></div></div><br></div></div></div><br><img width="0" height="0" class="mailtrack-img" alt="" style="display:flex" src="https://mailtrack.io/trace/mail/w/cWdpcy11c2VyQGxpc3RzLm9zZ2VvLm9yZw/b0840d5e4df992bdd6f13e84c4b17e1db97e8888.png?u=8747767" ><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, May 9, 2023 at 11:58 AM andrea antonello via QGIS-User <qgis-user@lists.osgeo.org> wrote:<br></div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><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-wrap;"><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);">dataTy<wbr>pe</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);">'<wbr>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>(<wbr>),</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:<wbr>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);">transf<wbr>ormCoordinates</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><wbr>(<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);"><wbr>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);">b<wbr>lock</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);"><wbr>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>
______________________________<wbr>_________________<br>
QGIS-User mailing list<br>
<a href="mailto:QGIS-User@lists.osgeo.org" target="_blank">QGIS-User@lists.osgeo.org</a><br>
List info: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-user" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailma<wbr>n/listinfo/qgis-user</a><br>
Unsubscribe: <a href="https://lists.osgeo.org/mailman/listinfo/qgis-user" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailma<wbr>n/listinfo/qgis-user</a><br>
</blockquote></div>