在非凸多边形上通过地理坐标计算面积

Question

我想计算从GPS轨迹中获取的多边形的面积。基本上，我在一段时间后（比如5秒）存储设备/用户的位置。

对于这个轨迹多边形，我想计算它所在的区域的面积。对于凸多边形来说，这应该不是问题，因为我猜我只需要计算三角形的面积（当每个三角形的起点在第一个点时）。基本上就像左图所示。（黄色多边形是由GPS位置组成的多边形，黑色线条显示用于计算面积的三角形，浅黄色是所需面积）

但昨晚我发现了这个想法的一个缺陷，即当多边形不是凸多边形时。不仅会计算在多边形外部的一部分面积（左上侧），还会测量多边形的某些区域超过一次（查看左下方的重叠三角形）。

有人有想法吗？我的意思是，如果我的多边形是S形的，甚至知道哪个区域应该被计算仍然很困难...（但我可以接受那样...只要它对几乎封闭的多边形得到足够公正的结果即可）。

我计算多边形的凸包，然后在那上面进行面积计算的另一个想法也行不通，如果多边形是非凸多边形。我不会重复计算某些区域，但是如右图所示，我将计算比实际更大的区域。

如果有人能帮我解决这个问题就太好了！谢谢！

- florianbaethge

---

我建议将多边形分割成许多凸多边形，这样面积就可以轻松计算。然而，我想不出任何好的分割方法。 - Mikulas Dite

如果您选择分割路线，您可能需要考虑使用Delaunay三角剖分。 - fmark

注意：你的方法几乎与Howard的方法相同。只是你没有使用负面区域来表示“回溯”的线段。 - Tom Sirgedas

有一个带有公式的Excel电子表格，可以使用笛卡尔坐标计算正多边形和不规则多边形。它可用于计算土地面积等... http://maruzar.blogspot.com/2011/12/irregular-and-regular-polygon-area.html - user1115212

实际上，赛道占用零面积，但通过将最后一个点连接回第一个点，您可以将其变成多边形。假设赛道从所有细线相交的点开始，如果您展示的任何一条赛道下一个是向右急转弯，以至于它进入图中所示的多边形内部，则您的多边形面积将减少。如果您遵循MathWorld页面中的顺时针-逆时针约定，则您展示的两个多边形的面积都为负数。这些结果可能是您想要的 - 我们无法确定是否需要更多上下文。 - David K

Answer 1

你可以查看多边形面积的通用公式：http://mathworld.wolfram.com/PolygonArea.html。这也适用于非凸多边形（只要它们不自交）。

- Howard

---

非常感谢，这看起来太棒了，正是我需要的！而且它完美地运行了...计算2x2矩阵的行列式很容易，所以这个方法很适用...谢谢！ - florianbaethge

我们在这里忽略了地球的曲率吗？ - Manu Manjunath

Answer 2

我经常这样做，但我不知道算法，因为我使用类似于C / C ++中的GEOS，Java中的JTS或Python中的Shapely库。

如果您可以承担额外的依赖项，我强烈建议使用它，因为计算是稳健的，经过测试的，采用开放标准输入格式（Well-Known Text），并且可以处理奇怪和不寻常的几何图形（例如具有孔的多边形等）。一旦您完成此操作，就可以执行各种奇怪而美妙的几何操作。

- fmark

Answer 3

有点晚了，但是我刚刚在Java中实现了这个GPS坐标。您必须按照此处描述的方式规范化GPS坐标：使用从笛卡尔空间和世界文件生成的纬度和经度计算多边形面积。效果很好，但由于笛卡尔坐标是基于理想地球半径的近似值，所以存在约0.005％的误差。请注意，下面的代码假定geojson样式[longitude，latitude]对成对出现，而不是相反的情况。

public static double area(double[][] polygon) {
    Validate.isTrue(polygon.length > 3,"polygon should have at least three elements");

    double total=0;
    double[] previous=polygon[0];

    double[] center = polygonCenter(polygon);
    double xRef=center[0];
    double yRef=center[1];


    for(int i=1; i< polygon.length;i++) {
        double[] current = polygon[i];
        // convert to cartesian coordinates in meters, note this not very exact
        double x1 = ((previous[0]-xRef)*( 6378137*PI/180 ))*Math.cos( yRef*PI/180 );
        double y1 = (previous[1]-yRef)*( Math.toRadians( 6378137 ) );
        double x2 = ((current[0]-xRef)*( 6378137*PI/180 ))*Math.cos( yRef*PI/180 );
        double y2 = (current[1]-yRef)*( Math.toRadians( 6378137 ) );

        // calculate crossproduct
        total += x1*y2 - x2*y1;
        previous=current;
    }

    return 0.5 * Math.abs(total);
}

- Jilles van Gurp

Answer 4

你想查看一种空间填充曲线，例如z曲线或希尔伯特曲线。SFC曲线将表面细分为许多瓷砖，并将二维问题简化为一维问题。你想搜索Nick关于空间索引希尔伯特曲线和四叉树的博客。

- Micromega