计算地球表面上任意多边形所包围的区域

Question

假设我有一组任意的纬度和经度对代表一些简单的闭合曲线上的点.在笛卡尔空间中,我可以使用格林定理轻松计算出这条曲线所包围的区域.计算球体表面区域的类似方法是什么？我想我所追求的是Matlab areaint函数背后的算法(甚至是一些近似).

Answer 1

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

1)整合来自纬度带的贡献.这里每个条带的面积为(Rcos(A)(B1-B0))(RdA),其中A是纬度,B1和B0是起始和结束经度,所有角度都是弧度.

2)将曲面分解为球形三角形,并使用Girard定理计算面积,然后将它们相加.

3)正如James Schek所建议的那样,在GIS工作中,他们使用一个区域保留投影到平面空间并计算那里的面积.

从您的数据描述中,听起来像第一种方法可能是最简单的.(当然,可能还有其他一些我不知道的简单方法.)

编辑 - 比较这两种方法:

在第一次检查时,似乎球形三角形方法最简单,但一般情况下并非如此.问题在于,不仅需要将区域分成三角形,而且需要分成球形三角形,即边长为大圆弧的三角形.例如,纬度边界不合格,因此需要将这些边界分解为更接近大圆弧的边缘.并且对于任意边缘来说这变得更加困难,其中大圆圈需要球面角度的特定组合.例如,考虑如何分解球体周围的中间带,将纬度0和45度之间的所有区域称为球形三角形.

最后,如果要对每种方法使用类似的错误正确地执行此操作,方法2将提供更少的三角形,但它们将更难确定.方法1给出了更多条带,但它们很容易确定.因此,我建议方法1作为更好的方法.

Answer 2

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

Answer 3

你在其中一个标签中提到"地理",所以我只能假设你是在大地水准面表面上的多边形区域之后.通常,这是使用投影坐标系而不是地理坐标系(即lon/lat)来完成的.如果您是在lon/lat中进行的,那么我会假设返回的度量单位是球体表面的百分比.

如果你想用更"GIS"的味道来做这件事,那么你需要为你的区域选择一个度量单位,找到一个保留区域的适当投影(不是所有的).既然你在谈论计算任意多边形,我会使用类似Lambert Azimuthal等面积投影的东西.将投影的原点/中心设置为多边形的中心,将多边形投影到新坐标系,然后使用标准平面技术计算面积.

如果你需要在一个地理区域做很多多边形,可能还有其他预测可以使用(或者足够接近).例如,如果所有多边形都围绕单个子午线聚集,则UTM是一个很好的近似值.

我不确定这是否与Matlab的isaint函数如何工作有任何关系.

Answer 4

我对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).

然后你可以使用好的'点积来找到边之间的角度:BA'.BC' = |BA'|*|BC'|*cos(B),在哪里BA'和BC'是从B边到边的切线向量A和C.

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