如何计算地球上自相交多边形的面积

Question

我试图找到地球上自相交（复杂）多边形的面积。是否有任何库可以为此类几何图形正确实现面积计算？

到目前为止，我已经尝试过 geographiclib 和 Polygon2（基于 gpc），但它们对复杂多边形给出了错误的结果。

或者，有没有一种简单的方法可以将自相交的多边形转换为一组简单的多边形，这样我就可以对简单的几何图形使用算法，然后对所有多边形求和？我知道一个解决方案是实现Bentley-Ottman algorithm 来查找交点，然后将多边形划分为一组简单的多边形，但是如果有一个我可以使用的库，我很乐意避免重新发明轮子（并可能引入错误)。

更新：我使用的坐标不明确包括交点。因此剪辑库似乎无法正确处理它们。例如，在下面的代码中coords_1 和coords_2 定义了同一个多边形。使用 Polygon2 裁剪库。

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

我想通过使用 coords_1 获得第二个结果。我也尝试过其他剪辑库，但到目前为止没有运气。

Answer 1

一个近似的解决方案是将多边形映射到平面，并使用平面几何库来进行计算几何。这种近似应该适用于相对于整个球体较小的多边形；但即使对于大型多边形，您也可以通过充分细化多边形来任意精确。

不幸的是，我认为本质上球形操作的库功能并不多。

编辑：这是一个我认为可以满足您的需求的匀称的解决方案：

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

转换您的纬度/经度多边形以在欧几里得平面上表现出色。然后使用类似：http://sourceforge.net/projects/polyclipping/.