[QGIS-Developer] Difference in area calculation using memory vs output layer

[Rudi von Staden]

Hi all,

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:

### from output layer

        habitatModelLayer = QgsVectorLayer(habitatModel,genspec,'ogr')
        iterator =
habitatModelLayer.getFeatures(QgsFeatureRequest().setFilterFid(0))
        feature = next(iterator)
        spatialiteFeature = QgsFeature(feature)     # make a copy
        spatialiteFields = spatialiteLayer.fields()
        spatialiteFeature.setFields(spatialiteFields)  # replace source
with destination fields
        for f in spatialiteFeature.fields().names():
            if f == 'm_area':
                spatialiteFeature[f] = feature.geometry().area()/1000000  #
result 2294.51 km2

### from memory layer

        habitatModelLayer = processingResult['model:zonal stats for
model_1:habitat model']
        iterator = habitatModelLayer.getFeatures()    # setFilterFid(0)
results in StopIteration
        feature = next(iterator)
        spatialiteFeature = QgsFeature(feature)     # make a copy
        spatialiteFields = spatialiteLayer.fields()
        spatialiteFeature.setFields(spatialiteFields)  # replace source
with destination fields
        for f in spatialiteFeature.fields().names():
            if f ==  'm_area':
                spatialiteFeature[f] = feature.geometry().area()/1000000 #
result 2156.43 km2

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.

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:

"+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"

If this is not expected behaviour I could try to create a simpler script to
reproduce the issue.

Thanks,
Rudi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/qgis-developer/attachments/20180808/09d71290/attachment.html>