How to calculate the area of self-intersecting polygons on Earth

Question

I am trying to find the area of self-intersecting (complex) polygons on Earth. Is there any library that implements correctly area computation for such geometries?

So far I have tried geographiclib and Polygon2 (based on gpc) but they give wrong results for complex polygons.

In alternative, is there a simple way to convert self-intersecting polygons into a set of simple polygons, so that I can use algorithms for simple geometries and then sum all the polygons? I know that a solution would be to implement the Bentley-Ottman algorithm to find intersection points and then divide the polygon into a set of simple polygons, but if there is a library I could use I would gladly avoid reinventing the wheel (and possibly introducing bugs).

Update: The coordinates I am using do not explicit include intersection points. Therefore clipping libraries seem not to handle them properly. For example, in the following code coords_1 and coords_2 define the same polygon. Using Polygon2 clipping library.

>>> coords_1 = [(0,0), (1,1), (1,0), (0,1)]
>>> coords_2 = [(0,0), (0.5, 0.5), (1,0), (1,1), (0.5,0.5), (0,1)]
>>> Polygon(coords_1).area()
0
>>> Polygon(coords_2).area()
0.5

I would like to get the second result, by using coords_1. I have also tried other clipping libraries but no luck so far.

Answer 1

An approximate solution would be to map the polygons to the plane, and use a plane geometry library to do the computational geometry. This approximation should work really well for polygons which are small relative to the entire sphere; but even for large polygons, you can be arbitrarily accurate by sufficiently refining your polygon.

I don't think there is much library functionality for intrinsically spherical operations, unfortunately.

Edit: here is a shapely solution which I think does what you want:

from shapely.geometry.polygon import  LinearRing, Polygon

coords_1 = [(0,0), (1,1), (1,0), (0,1)]
##coords_1 = [(0,0), (0,1), (1,1), (1,0)]

lr = LinearRing(coords_1)

if lr.is_valid:
    print Polygon(lr).area
else:
    print Polygon(coords_1).buffer(0).area + Polygon(coords_1[::-1]).buffer(0).area

Answer 2

Convert your lat/lon polygons to play nice on the Euclidean plane. Then use something like: http://sourceforge.net/projects/polyclipping/.