BUG: Fix handling of polygon holes when calculating area in Geod

[jacob-indigo]

Because geod.geometry_area_perimeter is signed, the handedness of the geometries coming in matters. There is an existing comment that says that counter-clockwise windings result in positive areas. The docs also point to the shapely.geometry.polygon.orient function which winds a polygon according to either the right-hand or left-hand rule. From the Shapely docs:

A sign of 1.0 means that the coordinates of the product’s exterior ring will be oriented counter-clockwise and the interior rings (holes) will be oriented clockwise.

That means that if you are finding the area of a polygon wound using orient, the exterior would generate a positive area and the interiors (holes) would generate negative areas. As is, Geod subtracts the holes' areas which results in area_exterior + area_interiors instead of area_exterior - area_interiors. Simply changing to an addition resolves the issue.

For reference, geos computes the area of a polygon using unsigned areas of the rings, i.e. abs(area_exterior) - sum(abs(area_interior_x)). See: https://github.com/libgeos/geos/blob/master/src/geom/Polygon.cpp#L374. It may be worthwhile considering either moving to an unsigned area in Geod or offering an unsigned function as an option.

[snowman2]

As is, Geod subtracts the holes' areas

That is the intended behavior of geometry_area_perimeter as well. The function assumes that the interiors and exteriors have the same sign. If the signs of the interiors and exteriors differ, then incorrect results will occur.

That means that if you are finding the area of a polygon wound using orient, the exterior would generate a positive area and the interiors (holes) would generate negative areas.

In shapely.geometry.polygon.orient() the interiors and exteriors are handled such that he signs change as expected. This method is used in shapely.ops.orient(). So, if the polygon's sign is changed using orient, then it should be consistent across the exterior and the interiors.

[snowman2]

For reference, geos computes the area of a polygon using unsigned areas of the rings, i.e. abs(area_exterior) - sum(abs(area_interior_x)). See: https://github.com/libgeos/geos/blob/master/src/geom/Polygon.cpp#L374. It may be worthwhile considering either moving to an unsigned area in Geod or offering an unsigned function as an option.

The sign may or may not be useful to some people. So, I think adding a signed=True kwarg to the existing function would make sense.

[snowman2]

Shapely docs:

A sign of 1.0 means that the coordinates of the product’s exterior ring will be oriented counter-clockwise and the interior rings (holes) will be oriented clockwise.

Looking at this closer, it seems like I missed the switch between using the >= and <= with the interiors and exteriors and the docstring is correct. So, the suggestion to use this function was a bad one as it always produces a polygon that is incompatible with this method 🤦‍♂️.

Your suggested solution in this PR might be the way to go. I will need to look into this and think on it 🤔.

[snowman2]

Related: pygeos/pygeos#71

[snowman2]

@jacob-indigo I think your fix in this PR is the way to go as it is consistent with PostGIS and shapely's behavior.

[jacob-indigo]

@snowman2 Thanks. Are you interested in having an unsigned variant? If so, I'm happy to add it. If not, 👍

[snowman2]

Are you interested in having an unsigned variant? If so, I'm happy to add it.

I think that would be useful. Go for it 🚀

[brendan-ward]

Having recently tripped over the pitfalls of winding order and signed areas when implementing this against pygeos geometries, an unsigned version seems like it would be useful and help avoid needing to normalize the winding order before calling. (note: in my specific case I'm using the underlying geod.polygon_area_perimeter())

[snowman2]

@jacob-indigo this looks good. Do you want to sqash your commits or do you mind if I do Squash and merge?

[jacob-indigo]

@snowman2 Thanks. Feel free to squash and merge