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

Question

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

26

假设我有一组任意的纬度和经度对，表示某个简单封闭曲线上的点。在笛卡尔空间中，我可以使用格林定理轻松计算由这样的曲线所包围的面积。那么在球面上计算面积的类比方法是什么？我猜想我想知道的是（即使是一些近似）Matlab的areaint函数背后的算法。

- Paul A. Hoadley

6个回答

11

我用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);
}

- user2548538

我需要用PHP实现同样的功能，但是这段代码对我来说太复杂了，我无法理解。你能帮我吗？ - arsal khan

9

您在标签中提到了“地理”，因此我只能假设您想要计算地球上多边形的面积。通常，这是使用投影坐标系而不是地理坐标系（即经度/纬度）完成的。如果您使用经度/纬度，则返回的单位将是球体表面的百分比。

如果您想以更加“GIS”的方式进行此操作，则需要选择一个面积单位，并找到一个保留面积的适当投影（并非所有投影都有此功能）。由于您要计算任意多边形，因此建议使用类似于Lambert Azimuthal Equal Area投影。将投影的原点/中心设置为多边形的中心，将多边形投影到新的坐标系中，然后使用标准平面技术计算面积。

如果您需要在地理区域内计算许多多边形，则可能有其他适用的投影（或足够接近的投影）。例如，如果您的所有多边形都聚集在单个子午线周围，则UTM是一个很好的近似值。

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

- James Schek

谢谢詹姆斯。我一直在想先将多边形投影到平面上是否可行。我发现投影可以保持面积，所以这可能是最理想的方法。 - Paul A. Hoadley

+1...对，我和一个也做很多GIS工作的朋友交流后，她告诉我这就是他们的做法。这种方法有什么原因吗？ - tom10

@Paul--你可能已经知道了，但要小心选择哪种投影方式。有些投影方式可以保留面积，而有些则不行。大多数地图使用的常见Web Mercator投影方式只能保留形状。 - James Schek

@tom 不确定为什么... 我猜是因为使用笛卡尔/平面系统更容易。如果你需要做的不仅仅是计算多边形的面积，那么拥有一个平面表示会让生活更轻松。此外，美国地质调查局等机构提供了大多数主要投影技术的“参考”实现。 - James Schek

@James：从计算的角度来看，哪种等面积投影是计算面积最便宜的？我的意思是，哪种投影有最简单的转换公式？ - Igor Brejc

@Igor：很遗憾，我不了解投影算法的具体细节——对我来说，这都是一个黑盒子。我会打听一下，看看能否找到一个确定的答案。 - James Schek

5

我不了解Matlab的函数，但是我们可以来看看。考虑通过从一个顶点画对角线将你的球面多边形分成球面三角形。球面三角形的表面积由以下公式给出:

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

其中 R 是球体的半径，A、B 和 C 分别为三角形的内角（以弧度表示）。括号中的量被称为“球面过剩量”。

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

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

其中S是你的多边形的角度和。括号中的量再次是该多边形的球面超额。

[编辑]这适用于多边形是否是凸的。重要的是它可以被分解成三角形。

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

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

[编辑以明确这些是切向量，而不是字面上的两点之间]

- Cascabel

第二个方程（涉及S的那个）需要多边形是凸多边形吗？ - Justin

谢谢Jefromi。非凸多边形还会使最初的分割成为球面三角形变得复杂。是否有一个众所周知的算法可以实现这一点？ - Paul A. Hoadley

@Paul - 把三角形内部与中心顶点最远的边平行的小区域称为“远边”。然后，如果您沿着周长走，并在映射区域内时添加三角形面积，在其外时减去它，我认为这个球面三角形算法将适用于凹凸（甚至是不连通的）形状。 - tom10

证明为您提供了一个多边形面积的公式，您可以直接使用。您不必将其分成三角形。那只是我证明公式有效的方法。 - Cascabel

作为最后的定论，请看这里。公式已经简单地打印出来了。我只是认为值得看一下它为什么有效。http://mathworld.wolfram.com/SphericalPolygon.html - Cascabel

显示剩余6条评论

0

这里是一个 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

请在这里找到一个更加明确的版本（并且有更多的参考和待办事项...）这里。

- Yellows

-1

你也可以查看这个 spherical_geometry 包的代码：这里和这里。它提供了两种不同的方法来计算球面多边形的面积。

- Felix D.

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- tom10 · Accepted Answer

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

1) 整合纬度带的贡献。这里每个带的面积将为(Rcos(A)(B1-B0))(RdA)，其中A是纬度，B1和B0是起始和结束经度，所有角度都以弧度表示。

2) 将表面分成球面三角形，并使用Girard定理计算面积，并将其相加。

3) 正如James Schek在此处建议的那样，在GIS工作中，他们使用一个面积保持投影到平面空间中，并在其中计算面积。

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

编辑 - 比较这两种方法：

初看起来，球面三角形方法似乎最容易，但通常情况下并非如此。问题在于不仅需要将区域分割成三角形，还需要分割成球面三角形，即边是大圆弧的三角形。例如，纬度边界不符合要求，因此需要将这些边界分解为更好地近似大圆弧的边缘。当大圆需要特定的球面角的组合时，任意边缘变得更加难以处理。例如，考虑如何将围绕球体的中间带子进行分解，例如介于0度和45度之间的所有区域。

最终，如果要按照每种方法的相似误差正确完成此操作，则方法2将给出更少的三角形，但它们更难确定。方法1会产生更多条带，但它们很容易确定。因此，我建议使用方法1作为更好的方法。