计算地球表面任意多边形包围的面积

Posted 2023-02-16

【中文标题】计算地球表面任意多边形包围的面积

【英文标题】：Calculating area enclosed by arbitrary polygon on Earth's surface

【发布时间】：2010-11-23 08:05:45

【问题描述】：

假设我有一组任意的纬度和经度对，表示一些简单的封闭曲线上的点。在笛卡尔空间中，我可以使用格林定理轻松计算出这样一条曲线所包围的面积。计算球体表面面积的类似方法是什么？我想我所追求的是（甚至是一些近似）Matlab's areaint function 背后的算法。

【参考方案1】：

这是一个 Python 3 实现，大致受到上述答案的启发：

def polygon_area(lats, lons, algorithm = 0, radius = 6378137):
    """
    Computes area of spherical polygon, assuming spherical Earth. 
    Returns result in ratio of the sphere's area if the radius is specified.
    Otherwise, in the units of provided radius.
    lats and lons are in degrees.
    """
    from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
    lats = np.deg2rad(lats)
    lons = np.deg2rad(lons)

    # Line integral based on Green's Theorem, assumes spherical Earth

    #close polygon
    if lats[0]!=lats[-1]:
        lats = append(lats, lats[0])
        lons = append(lons, lons[0])

    #colatitudes relative to (0,0)
    a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
    colat = 2*arctan2( sqrt(a), sqrt(1-a) )

    #azimuths relative to (0,0)
    az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)

    # Calculate diffs
    # daz = diff(az) % (2*pi)
    daz = diff(az)
    daz = (daz + pi) % (2 * pi) - pi

    deltas=diff(colat)/2
    colat=colat[0:-1]+deltas

    # Perform integral
    integrands = (1-cos(colat)) * daz

    # Integrate 
    area = abs(sum(integrands))/(4*pi)

    area = min(area,1-area)
    if radius is not None: #return in units of radius
        return area * 4*pi*radius**2
    else: #return in ratio of sphere total area
        return area

请找到一个更明确的版本（以及更多的参考资料和待办事项...）here。

【参考方案2】：

您还可以查看spherical_geometry 包的代码：Here 和here。它确实提供了两种不同的方法来计算球面多边形的面积。

【参考方案3】：

有几种方法可以做到这一点。

1) 整合来自纬度带的贡献。这里每条带的面积为(Rcos(A)(B1-B0))(RdA)，其中A为纬度，B1和B0为起点和终点经度，所有角度均以弧度为单位。

2) 将曲面分解为spherical triangles，并使用Girard's Theorem 计算面积，并将它们相加。

3) 正如 James Schek 所建议的，在 GIS 工作中，他们使用区域保留投影到平坦空间并计算其中的面积。

根据您的数据描述，听起来第一种方法可能是最简单的。 （当然，可能还有其他我不知道的更简单的方法。）

编辑——比较这两种方法：

第一次检查时，球面三角形方法似乎最简单，但通常情况并非如此。问题在于，不仅需要将该区域分解为三角形，还需要将其分解为球形三角形，即边为大圆弧的三角形。例如，纬度边界不符合条件，因此需要将这些边界分解为更接近大圆弧的边。对于大圆需要特定球面角组合的任意边缘，这变得更加困难。例如，考虑如何将球体周围的中间带分割成球面三角形，例如 0 到 45 度之间的所有区域。

最后，如果要正确执行此操作，并且每种方法的错误相似，方法 2 将给出更少的三角形，但它们将更难确定。方法 1 给出了更多条带，但它们很容易确定。因此，我建议方法 1 是更好的方法。

【讨论】：

我的回答是对您的 (2) 的详细说明。在计算上，向量数学将比集成便宜得多，而且很可能更容易编码。请注意，所有向量运算都可以使用球坐标向量完成，纬度/经度本质上是。

@Jefromi：我认为你的评论不正确，我已经编辑了我的答案来解决这个问题。

谢谢汤姆。我假设 Matlab 函数的功能类似于您的 (1)。我看看能不能拿到那张纸。关于您对球形三角形的反对意见，我的问题在这一点上可能并不完全清楚，但我所拥有的只是顶点——一组有序的纬度/经度对。边缘只是隐含的，因此我们不妨假设它们是用于任何计算的大圆。

保罗...这是有道理的，尤其是当您的点靠得很近时。

我设法找到了那篇论文。而且，相当令人惊讶的是，由于文章中提到的 FTP 服务器不见了，相关的代码。所以我会提高我的 Fortran 技能并检查一下。

【参考方案4】：

我用java重写了MATLAB的“areaint”函数，结果完全一样。 “areaint”计算“单位面积”，所以我将答案乘以地球表面积（5.10072e14 平方米）。

private double area(ArrayList<Double> lats,ArrayList<Double> lons)
       
    double sum=0;
    double prevcolat=0;
    double prevaz=0;
    double colat0=0;
    double az0=0;
    for (int i=0;i<lats.size();i++)
    
        double colat=2*Math.atan2(Math.sqrt(Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)+ Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)),Math.sqrt(1-  Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)- Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)));
        double az=0;
        if (lats.get(i)>=90)
        
            az=0;
        
        else if (lats.get(i)<=-90)
        
            az=Math.PI;
        
        else
        
            az=Math.atan2(Math.cos(lats.get(i)*Math.PI/180) * Math.sin(lons.get(i)*Math.PI/180),Math.sin(lats.get(i)*Math.PI/180))% (2*Math.PI);
        
        if(i==0)
        
             colat0=colat;
             az0=az;
                   
        if(i>0 && i<lats.size())
        
            sum=sum+(1-Math.cos(prevcolat  + (colat-prevcolat)/2))*Math.PI*((Math.abs(az-prevaz)/Math.PI)-2*Math.ceil(((Math.abs(az-prevaz)/Math.PI)-1)/2))* Math.signum(az-prevaz);
        
        prevcolat=colat;
        prevaz=az;
    
    sum=sum+(1-Math.cos(prevcolat  + (colat0-prevcolat)/2))*(az0-prevaz);
    return 5.10072E14* Math.min(Math.abs(sum)/4/Math.PI,1-Math.abs(sum)/4/Math.PI);

【讨论】：

我需要同样的东西，但在 php 中，但代码似乎太复杂，我无法理解。你能帮我解决这个问题吗？

【参考方案5】：

您在其中一个标签中提到了“地理”，所以我只能假设您在大地水准面表面上的多边形区域之后。通常，这是使用投影坐标系而不是地理坐标系（即经度/纬度）来完成的。如果你用 lon/lat 来做，那么我会假设返回的测量单位是球体表面的百分比。

如果您想以更“GIS”的风格来执行此操作，那么您需要为您的区域选择一个测量单位并找到一个合适的投影来保留区域（并非所有人都这样做）。由于您正在谈论计算任意多边形，因此我将使用 Lambert Azimuthal Equal Area 投影之类的东西。将投影的原点/中心设置为多边形的中心，将多边形投影到新坐标系，然后使用标准平面技术计算面积。

如果您需要在一个地理区域内制作多个多边形，则可能有其他投影可以工作（或足够接近）。例如，如果您的所有多边形都聚集在一条子午线上，UTM 就是一个极好的近似值。

我不确定这是否与 Matlab 的 areaint 函数的工作方式有关。

【讨论】：

谢谢詹姆斯。我想知道首先将多边形投影到平面上是否可行。我看到投影保留了区域，所以这可能是理想的。

+1... 对，与一位也从事大量 GIS 工作的朋友交谈，她告诉我他们就是这样做的。这种方法有原因吗？

@Paul——您可能已经知道这一点，但请注意您选择的投影。有些投影保留面积，有些则没有。大多数地图上使用的常见 Web Mercator 仅保留形状。

@tom 不知道为什么...我的猜测是使用笛卡尔/平面系统更容易。如果您需要做的不仅仅是计算多边形的面积，那么拥有平面表示会使生活更轻松。另外——USGS 等提供大多数主要投影技术的“参考”实现。

@James：从计算的角度来看：哪个等面积投影在计算面积时最便宜？我的意思是哪个投影具有最简单的转换公式？

【参考方案6】：

我对 Matlab 的功能一无所知，但我们开始吧。考虑将球面多边形分割成球面三角形，比如从顶点绘制对角线。球面三角形的表面积由下式给出

R^2 * ( A + B + C - \pi)

其中R 是球体的半径，A、B 和C 是三角形的内角（以弧度为单位）。括号中的数量称为“球形过剩”。

您的n 边多边形将被分割成n-2 三角形。对所有三角形求和，提取R^2 的公因数，并将所有\pi 放在一起，多边形的面积为

R^2 * ( S - (n-2)\pi )

其中S 是多边形的角度和。括号中的数量又是多边形的球面余量。

[edit] 无论多边形是否凸面都是如此。重要的是它可以分解成三角形。

您可以通过一些矢量数学来确定角度。假设您有三个顶点A、B、C，并且对B 处的角度感兴趣。因此，我们必须从点B沿大圆段（多边形边缘）找到两个与球体相切的向量（它们的大小无关）。让我们为BA 解决这个问题。大圆位于OA 和OB 定义的平面内，其中O 是球体的中心，因此它应该垂直于法向量OA x OB。它也应该垂直于OB，因为它在那里相切。因此，这样的向量由OB x (OA x OB) 给出。您可以使用右手定则来验证这是否在适当的方向上。另请注意，这简化为OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA)。

然后您可以使用良好的 ol' 点积来找到边之间的角度：BA'.BC' = |BA'|*|BC'|*cos(B)，其中 BA' 和 BC' 是从 B 沿边到 A 和 C 的切向量.

[编辑清楚这些是切向量，而不是点之间的文字]

【讨论】：

吉拉德定理的证明非常优雅——如果你想完全理解你在这里所做的事情，看看math.rice.edu/~pcmi/sphere/gos3.html和math.rice.edu/~pcmi/sphere/gos4.html

第二个方程（涉及S的那个）是否要求多边形是凸的？

谢谢杰弗罗米。非凸多边形也会使初始分割成球形三角形变得复杂。是否有众所周知的算法来实现这一目标？

等等，我们为什么要分解它？面积公式仍然有效！证明不依赖于凸性。多边形的面积仍然是三角形面积的总和，无论你如何切割它。

抱歉，我不是在质疑证据，而是在质疑切割本身。在某些时候，我希望能够以编程方式执行此操作，并且显然从顶点绘制对角线仅适用于凸多边形。我要问的是是否有另一种不会被非凸形状难倒的分割算法。