通过任意多边形地球表面封闭面积计算(Calculating area enclosed by arb

说我有代表对一些简单的，封闭的曲线点的纬度和经度对任意一组。在笛卡尔空间，我可以很容易地计算出使用格林定理这样的曲线围成的区域。什么是类似的方法，以一个球的表面上计算的面积？我猜我后（甚至有些近似）背后的算法Matlab的areaint功能。

Answer 1:

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

1）从集成纬度带的贡献。在此，每个条带的区域将是（RCOS（A）（B1-B0））（RDA），其中A是纬度，B1和B0是起始和结束经度，并且所有的角度以弧度。

2）表面打入球面三角形，并使用计算面积Girard的定理，并添加这些起来。

3）由于这里由詹姆斯Schek建议，在GIS的工作，他们使用的区域保持投影到平整的空间，并在那里计算面积。

从你的数据像第一种方法的描述，声音可能是最容易的。（当然，可能还有其他更容易的方法，我不知道。）

编辑-比较这两种方法：

乍一看，它可能看起来球面三角的方法是最简单的，但是，在一般情况下，这种情况并非如此。问题是，一个不仅需要打破区域成三角形，但成球形三角形 ，也就是三角形，其边大圆弧。例如， 纬度边界不符合条件 ，所以这些边界需要被分解成边缘，更好地近似大圆弧。这变得更加困难，其中大圆需要球角度的特定组合的任意边做。举个例子，怎么一会分手中间带围绕一个球，说所有的区域纬度0和45deg之间成球形三角形。

最后，如果一个人有每个方法类似的错误，正确地做到这一点，方法2将给予较少的三角形，但他们将很难确定。方法1提供了更多的条，但他们是微不足道的决定。因此，我建议的方法1的更好的方法。

Answer 2:

我重写了Java中的MATLAB的“areaint”的功能，它具有完全相同的结果。 “areaint”计算“单位suface”，所以我再乘以地球表面积（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 /纬度，那么我将承担返回的结果是球体表面的百分之单位的小节。

如果你想与更多的“GIS”的味道要做到这一点，那么你需要选择一个单位的，测量你的区域，发现保存面积适当的投影（不是所有的事）。既然你是在谈论计算任意多边形，我会用像一兰勃特方位等积投影。将投影原点/中心是你的多边形的中心，项目的多边形到新的坐标系，然后计算出使用标准平面技术领域。

如果你需要做的多边形在一个地理区域，也有可能其他的，将工作（或将是足够接近）的预测。 UTM，例如，如果是你所有的多边形都围绕一个经络集群极好的近似。

我不知道是否有任何这有什么用Matlab的areaint功能是如何工作的。

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是你的多边形的角度总和。括号中的数量又是球形多余的多边形。

[编辑]这是真实的多边形是否为凸的。所有重要的是，它可以被划分成三角形。

您可以从一个位矢量数学的判断的角度。假设你有三个顶点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) OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA)