Area calculation seems really wrong

[seav]

Steps to reproduce the error

1. Add a rectangular polygon location in a high-latitude location, like Helsinki.
2. Export the project data as XLS then inspect the area provided for the location.

Actual behavior

For a sample building I traced in Helsinki, the calculated area is ~12987 m².

Expected behavior

When I calculate the area using two other methods, I get an area of ~3194 m² and ~3198 m². The first method is the simple length×width calculation using JOSM and the second method is by using the Ellipsoidal Area script in QGIS applied on the location when exported as a shape file from the platform.

Notes

The two methods I used to calculate the area match closely, but both are way off when compared to the area computed by the platform. I have to conclude that the method used to calculate the area in PR #1534 is very wrong. Looking at the actual code, a polygon location is first transformed to the Spherical Mercator projection (EPSG:3857) and then the area in this projection is computed. But this is bound to be wrong because locations far from the equator will have a larger projected area. In the latitude of Helsinki (~60°N), we expect lengths to be twice as long as lengths on the equator (because cos 60° = 0.5). Therefore areas will be 4 times as large. This matches what we see in the test location that I created (12987 ≈ 4 × 3194).

[amplifi]

This was already identified in https://github.com/Cadasta/cadasta-platform/pull/1534#issuecomment-314882132 and https://github.com/Cadasta/cadasta-platform/pull/1534#issuecomment-314995326

[seav]

@amplifi: Not quite. The two comments only say that there's a discrepancy between the tooltip area and the database-stored area. It doesn't actually say which of the two is correct. Here, I am asserting that the database-stored area is incorrect regardless of whether the tooltip area is correct or not.

[amplifi]

@seav The error in area calculation due to projection was specifically discussed. Please read the comments linked. When I raised this issue in the PR phase, the code was decided to be acceptable as-is.

[oliverroick]

Eugene is right on this. The area is off because the projection used to calculate the area is not equal-area.

There is a way to solve this using pyproj and shapely.

The resulting code may look like this:

def calculate_area(sender, instance, **kwargs):
    geom = instance.geometry
    from django.contrib.gis.geos.polygon import Polygon

    if geom and isinstance(geom, Polygon) and geom.valid:
        import pyproj
        from shapely import ops
        from shapely.wkt import loads
        from functools import partial

        geom = loads(geom.wkt)
        geom_area = ops.transform(
            partial(
                pyproj.transform,
                pyproj.Proj(init='EPSG:4326'),
                pyproj.Proj(
                    proj='aea',
                    lat1=geom.bounds[1],
                    lat2=geom.bounds[3])),
            geom)
        instance.area = geom_area.area

For the migration, we can either use the same solution, which will be slow because we would have to update each SpatialUnit instance individually or we update our data model and the migration:

1. We need to convert geometries to geography.

2. In the migration we can do area calculations on the database without re-projecting the data:

def calculate_area_field(apps, schema_editor):
    SpatialUnit = apps.get_model('spatial', 'SpatialUnit')
    SpatialUnit.objects.exclude(geometry=None).extra(
        where=["geometrytype(geometry) LIKE 'POLYGON'"]).update(
        area=Area('geometry'))

So we'd have to add another migration, which converts to geography and re-calculates the area.

We might be able to apply a similar approach when creating new locations (instead of the pyproj and shapely approach) but I haven't looked into that.

[amplifi]

@oliverroick My point is that I did raise this issue during the PR review phase, and I was the only one to actually evaluate the area calculation as part of that review. When I raised it, it was dismissed as not being important enough to investigate and resolve before merging. There is a pattern of this situation recurring that needs to be addressed.

[seav]

Re-reading the PR comments, it seems the issue was dismissed because the conclusion was that the tooltip calculation was the one in error, and this implied that the database calculation is correct:

Basically, it appears that the Tooltip feature doesn't take into account the latitude when calculating area, meaning that the closer a polygon is to the poles, the less accurate it is.

And since the tooltip calculation is not included in SMAP, as confirmed by @bjohare, the area discrepancy was dismissed:

If we're omitting that in SMAP, I'm not fussed.

My assertion now is that the database calculation is in fact erroneous. Whether the tooltip calculation is correct or not is now irrelevant.

[alukach]

@seav Thanks for the well written bug report. Can you provide any pointers on how you manually tested area for the building footprint via QGIS or JOSM?

And yes, my earlier interpretation was that it was only the tooltip that was incorrect, not the DB values. Reading through your and @oliverroick 's comments, you do seem to be correct. I think storing the data as a geography rather than geometry should resolve this. I was under the impression that storing the data in 4326 was storing it in lat/lng, however apparently this is not completely correct...

[seav]

@alukach: storing the data in 4326 does store coordinates as lat/lng degrees. But computing the area from that without re-projecting will give you a value in square degrees (and assuming the earth is rectangularly flat), not square meters.

For the JOSM method, JOSM helpfully shows the length of a line segment in the status bar when you are drawing a way. So I used that feature to measure the length and width of the building and then multiply the two together to get the area.

For the QGIS method, I used the 3rd answer here: https://gis.stackexchange.com/questions/23355/how-to-calculate-polygon-areas-in-qgis