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

Question

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

Answer 1

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

如果你想用一个更“地理信息系统”的风格来做这件事，那么你需要为你的区域选择一个度量单位，并找到一个适当的投影来保存区域（不是所有的都要做）。既然你说的是计算任意多边形，我会使用类似于Lambert Azimuthal Equal Area投影的东西。将投影的原点/中心设置为多边形的中心，将多边形投影到新的坐标系，然后使用标准平面技术计算面积。

如果你需要在一个地理区域做很多多边形，可能还有其他的投影可以工作（或者足够接近）。例如，如果所有多边形都聚集在一个子午线周围，则UTM是一个极好的近似值。

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

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

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

1）综合纬度带的贡献。这里每条带的面积是（Rcos（A）（B1-B0））（RdA），其中A是纬度，B1和B0是起始和结束经度，所有角度都是弧度。

2）将曲面分解成spherical triangles，并使用Girard's Theorem计算面积，然后将其相加。

3）正如James Schek在这里所建议的，在GIS工作中，他们使用平面空间上保留面积的投影，并计算其中的面积。

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

编辑–比较这两种方法：

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

最后，如果要在每个方法都有相似的误差的情况下正确地完成这项工作，方法2将给出更少的三角形，但它们将更难确定。方法1给出了更多的条带，但它们很容易确定。因此，我建议方法1是更好的方法。