假设我有一组任意的纬度和经度对，表示一些简单的闭合曲线上的点。在笛卡尔空间中，我可以使用格林定理轻松计算出这样一条曲线所包围的面积。计算球体表面面积的类似方法是什么？我想我所追求的是（甚至是一些近似）Matlabareaint函数背后的算法。
问问题
28091 次
6 回答
13
有几种方法可以做到这一点。
1）整合来自纬度带的贡献。这里每个条带的面积将是 (Rcos(A)(B1-B0))(RdA)，其中 A 是纬度，B1 和 B0 是开始和结束经度，所有角度都以弧度为单位。
2) 将曲面分成球面三角形，并使用吉拉德定理计算面积，并将它们相加。
3) 正如 James Schek 所建议的，在 GIS 工作中，他们使用区域保留投影到平坦空间并计算其中的面积。
从您的数据描述来看，听起来第一种方法可能是最简单的。（当然，可能还有其他我不知道的更简单的方法。）
编辑——比较这两种方法：
乍一看，球面三角形方法似乎最简单，但通常情况并非如此。问题在于，不仅需要将该区域分解为三角形，还需要将其分解为球面三角形，即边为大圆弧的三角形。例如，纬度边界不符合条件，因此需要将这些边界分解为更接近大圆弧的边。对于大圆需要特定球面角组合的任意边缘，这变得更加困难。例如，考虑如何将一个球体周围的中间带分开，比如将 0 到 45 度之间的所有区域分成球形三角形。
最后，如果要正确执行此操作，并且每种方法的错误相似，方法 2 将给出更少的三角形，但它们将更难确定。方法 1 给出了更多条带，但它们很容易确定。因此，我建议方法 1 作为更好的方法。
于 2009-08-27T17:51:49.293 回答
10
我用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);
}
于 2015-02-10T05:21:17.777 回答
9
您在其中一个标签中提到了“地理”，所以我只能假设您在大地水准面表面上的多边形区域之后。通常，这是使用投影坐标系而不是地理坐标系（即经度/纬度）来完成的。如果您要在 lon/lat 中执行此操作，那么我会假设返回的测量单位将是球体表面的百分比。
如果您想以更“GIS”的风格来执行此操作，那么您需要为您的区域选择一个测量单位并找到一个合适的投影来保留区域（并非所有人都这样做）。由于您正在谈论计算任意多边形，因此我会使用Lambert Azimuthal Equal Area projection 之类的东西。将投影的原点/中心设置为多边形的中心，将多边形投影到新的坐标系，然后使用标准平面技术计算面积。
如果您需要在一个地理区域中制作许多多边形，则可能有其他投影可以工作（或足够接近）。例如，如果您的所有多边形都围绕单个子午线聚集，那么 UTM 是一个极好的近似值。
我不确定这是否与 Matlab 的 areaint 函数的工作方式有关。
于 2009-08-28T15:42:35.570 回答
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在哪里，所以它应该垂直于法向量。它也应该垂直于，因为它在那里相切。因此，这样的向量由 给出。您可以使用右手定则来验证这是否在适当的方向上。另请注意，这简化为.OOA x OBOBOB 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。
[编辑清楚这些是切向量，而不是点之间的文字]
于 2009-08-27T17:51:25.050 回答
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
请在此处找到更明确的版本（以及更多参考资料和待办事项...）。
于 2020-04-12T19:11:41.043 回答