<div dir="ltr">Hi all,<div><br></div><div>I'm building a script to iterate over a relatively complex processing model. I started out using a parameterAsOutputLayer to get the result of the algorithm, and then loaded that as a QgsVectorLayer() using the 'ogr' provider. I then refactored that to use a 'memory:' layer instead. The strange thing is that one of the steps in the script is to calculate the area of the single feature produced by the algorithm, and the area result is different when calculated from the shapefile loaded in vs calculated directly from the memory layer. My code below:</div><div><br></div><div>### from output layer</div><div><br></div><div><div>        habitatModelLayer = QgsVectorLayer(habitatModel,genspec,'ogr')</div><div>        iterator = habitatModelLayer.getFeatures(QgsFeatureRequest().setFilterFid(0))</div><div>        feature = next(iterator)</div><div>        spatialiteFeature = QgsFeature(feature)     # make a copy        </div><div>        spatialiteFields = spatialiteLayer.fields()</div><div>        spatialiteFeature.setFields(spatialiteFields)  # replace source with destination fields</div><div>        for f in spatialiteFeature.fields().names():</div><div>            if f == 'm_area':<br></div><div>                spatialiteFeature[f] = feature.geometry().area()/1000000  # result 2294.51 km2</div></div><div><br></div><div>### from memory layer</div><div><br></div><div><div>        habitatModelLayer = processingResult['model:zonal stats for model_1:habitat model']</div><div>        iterator = habitatModelLayer.getFeatures()    # setFilterFid(0) results in StopIteration</div><div>        feature = next(iterator)</div><div>        spatialiteFeature = QgsFeature(feature)     # make a copy        </div><div>        spatialiteFields = spatialiteLayer.fields()</div><div>        spatialiteFeature.setFields(spatialiteFields)  # replace source with destination fields</div><div>        for f in spatialiteFeature.fields().names():</div><div>            if f ==  'm_area':</div><div>                spatialiteFeature[f] = feature.geometry().area()/1000000 # result 2156.43 km2</div><div><br></div><div>I did also calculate a range of zonal stats on the layer (which has only one multipolygon feature), but those results were consistent between output / memory layers.</div><div><br></div><div>Is the area difference expected because of how shapefiles store geometry vs how they are handled by memory layers, or could there be a bug somewhere? Maybe this isn't the recommended way to calculate area? In case it matters, I'm using a custom CRS:</div><div><br></div><div>"+proj=aea +lat_1=-24 +lat_2=-32 +lat_0=0 +lon_0=24 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"</div><div><br></div></div><div>If this is not expected behaviour I could try to create a simpler script to reproduce the issue.</div><div><br></div><div>Thanks,</div><div>Rudi</div><div><br></div></div>