当前位置： 文江博客 话题详情

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

发布于 2024-08-03 06:45:33 字数 292 浏览 9 评论 0 原文

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

原文

Say I have an arbitrary set of latitude and longitude pairs representing points on some simple, closed curve. In Cartesian space I could easily calculate the area enclosed by such a curve using Green's Theorem. What is the analogous approach to calculating the area on the surface of a sphere? I guess what I am after is (even some approximation of) the algorithm behind Matlab's areaint function.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（8）

蓝戈者 2024-08-10 06:45:33

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

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

2) 将曲面分成球形三角形，并使用吉拉德定理，并将它们相加。

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

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

编辑 - 比较这两种方法：

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

最后，如果要正确地做到这一点，并且每种方法都有类似的错误，方法 2 将给出更少的三角形，但它们将更难确定。方法 1 提供了更多条带，但确定它们很简单。因此，我建议方法1作为更好的方法。

There several ways to do this.

1) Integrate the contributions from latitudinal strips. Here the area of each strip will be (Rcos(A)(B1-B0))(RdA), where A is the latitude, B1 and B0 are the starting and ending longitudes, and all angles are in radians.

2) Break the surface into spherical triangles, and calculate the area using Girard's Theorem, and add these up.

3) As suggested here by James Schek, in GIS work they use an area preserving projection onto a flat space and calculate the area in there.

From the description of your data, in sounds like the first method might be the easiest. (Of course, there may be other easier methods I don't know of.)

Edit – comparing these two methods:

On first inspection, it may seem that the spherical triangle approach is easiest, but, in general, this is not the case. The problem is that one not only needs to break the region up into triangles, but into spherical triangles, that is, triangles whose sides are great circle arcs. For example, latitudinal boundaries don't qualify, so these boundaries need to be broken up into edges that better approximate great circle arcs. And this becomes more difficult to do for arbitrary edges where the great circles require specific combinations of spherical angles. Consider, for example, how one would break up a middle band around a sphere, say all the area between lat 0 and 45deg into spherical triangles.

In the end, if one is to do this properly with similar errors for each method, method 2 will give fewer triangles, but they will be harder to determine. Method 1 gives more strips, but they are trivial to determine. Therefore, I suggest method 1 as the better approach.

回复 原文
我爱人 2024-08-10 06:45:33

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

I rewrote the MATLAB's "areaint" function in java, which has exactly the same result.
"areaint" calculates the "suface per unit", so I multiplied the answer by Earth's Surface Area (5.10072e14 sq m).

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);
}
回复 原文
你げ笑在眉眼 2024-08-10 06:45:33

您在标签之一中提到了“地理”，因此我只能假设您正在寻找大地水准面表面上的多边形区域。通常，这是使用投影坐标系而不是地理坐标系（即经度/纬度）来完成的。如果您要以经度/纬度进行计算，那么我会假设返回的测量单位将是球体表面的百分比。

如果您想以更具“GIS”风格的方式执行此操作，那么您需要为您的区域选择一个测量单位，并找到一个保留面积的适当投影（并非全部如此）。由于您正在谈论计算任意多边形，所以我会使用类似 兰伯特方位等积投影。将投影的原点/中心设置为多边形的中心，将多边形投影到新的坐标系，然后使用标准平面技术计算面积。

如果您需要在一个地理区域中绘制许多多边形，则可能有其他投影可以使用（或足够接近）。例如，如果所有多边形都聚集在一条子午线周围，那么 UTM 就是一个很好的近似值。

我不确定这是否与Matlab的areint函数的工作原理有关。

You mention "geography" in one of your tags so I can only assume you are after the area of a polygon on the surface of a geoid. Normally, this is done using a projected coordinate system rather than a geographic coordinate system (i.e. lon/lat). If you were to do it in lon/lat, then I would assume the unit-of-measure returned would be percent of sphere surface.

If you want to do this with a more "GIS" flavor, then you need to select an unit-of-measure for your area and find an appropriate projection that preserves area (not all do). Since you are talking about calculating an arbitrary polygon, I would use something like a Lambert Azimuthal Equal Area projection. Set the origin/center of the projection to be the center of your polygon, project the polygon to the new coordinate system, then calculate the area using standard planar techniques.

If you needed to do many polygons in a geographic area, there are likely other projections that will work (or will be close enough). UTM, for example, is an excellent approximation if all of your polygons are clustered around a single meridian.

I am not sure if any of this has anything to do with how Matlab's areaint function works.

回复 原文
绳情 2024-08-10 06:45:33

我对Matlab的功能一无所知，但我们开始吧。考虑将球形多边形分割成球形三角形，例如从顶点绘制对角线。球形三角形的表面积由下式给出，

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

其中 R 是球体的半径，A、B 和 C 是三角形的内角（以弧度为单位）。括号中的数量称为“球形过剩”。

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

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

S< /code> 是多边形的角度和。括号中的数量又是多边形的球面余量。

[编辑] 无论多边形是否为凸多边形，这都是正确的。重要的是它可以被分割成三角形。

您可以通过一些向量数学来确定角度。假设您有三个顶点 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）。

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

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

I don't know anything about Matlab's function, but here we go. Consider splitting your spherical polygon into spherical triangles, say by drawing diagonals from a vertex. The surface area of a spherical triangle is given by

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

where R is the radius of the sphere, and A, B, and C are the interior angles of the triangle (in radians). The quantity in the parentheses is known as the "spherical excess".

Your n-sided polygon will be split into n-2 triangles. Summing over all the triangles, extracting the common factor of R^2, and bringing all of the \pi together, the area of your polygon is

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

where S is the angle sum of your polygon. The quantity in parentheses is again the spherical excess of the polygon.

[edit] This is true whether or not the polygon is convex. All that matters is that it can be dissected into triangles.

You can determine the angles from a bit of vector math. Suppose you have three vertices A,B,C and are interested in the angle at B. We must therefore find two tangent vectors (their magnitudes are irrelevant) to the sphere from point B along the great circle segments (the polygon edges). Let's work it out for BA. The great circle lies in the plane defined by OA and OB, where O is the center of the sphere, so it should be perpendicular to the normal vector OA x OB. It should also be perpendicular to OB since it's tangent there. Such a vector is therefore given by OB x (OA x OB). You can use the right-hand rule to verify that this is in the appropriate direction. Note also that this simplifies to OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA).

You can then use the good ol' dot product to find the angle between sides: BA'.BC' = |BA'|*|BC'|*cos(B), where BA' and BC' are the tangent vectors from B along sides to A and C.

[edited to be clear that these are tangent vectors, not literal between the points]

回复 原文
心的憧憬 2024-08-10 06:45:33

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

请找到一个更明确的版本（以及更多参考文献和 TODO...） 此处。

Here is a Python 3 implementation, loosely inspired by the above answers:

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

Please find a somewhat more explicit version (and with many more references and TODOs...) here.

回复 原文
离笑几人歌 2024-08-10 06:45:33

如果您的多边形适合某个半球，那么您可以将多边形的每条边从其第一个顶点扫过的三角形面积相加。以下公式获取单位向量顶点 A、B 和 C 的球面三角形的面积：

area = 2 * atan( (AxB.C) / (1 + A.B + B.C + C.A) )

上式使用点积，(AxB.C) 是三重积。计算出的面积是有符号的，因此该方法将自动处理凹多边形。

如果您的多边形不适合半球，则多边形中的某些顶点可能相距 > 180 度，并且上述公式将使用大圆的错误部分来测量三角形。所以可能会失败。

注意：早期的答案建议使用吉拉德定理来获取球面三角形的面积。这对于简单的凸多边形来说是可以的，但是很难将吉拉德定理用于凹多边形，因为将多边形切割成位于多边形内部的三角形很困难。此外，如果存在重复的顶点，吉拉德定理往往会崩溃，并且/或者当顶点非常接近时会得到不准确的答案。此答案中的面积公式对于重复和接近的顶点表现更好。

If your polygon fits within some hemisphere then you can sum up the triangular areas swept out by each of the polygon's edges from its first vertex. The following formula gets the area of a spherical triangle with unit vector vertices A, B and C:

area = 2 * atan( (AxB.C) / (1 + A.B + B.C + C.A) )

The above formula uses dot products, and (AxB.C) is the triple product. The calculated area is signed, so this method will handle concave polygons automatically.

If your polygon does not fit onto a hemisphere then some vertices could be >180 degrees apart in the polygon, and the above formula will measure the triangle with the wrong part of the great circle. So it may fail.

Note: An earlier answer suggested using Girard's Theorem to get the area of a spherical triangle. This will be OK with simple convex polygons, but it is difficult to use Girard's Theorem with concave polygons as cutting the polygon into triangles that lie inside the polygon is difficult. Also, Girard's Theorem tends to blow up if there are repeated vertices, and/or get inaccurate answers when vertices are very close together. The area formula in this answer behaves better with repeated and close vertices.

回复 原文
锦爱 2024-08-10 06:45:33

您还可以查看 spherical_geometry 包的代码： 此处 和 这里。它确实提供了两种不同的方法来计算球形多边形的面积。

You could also have a look at this code of the spherical_geometry package: Here and here. It does provide two different methods for calculating the area of a spherical polygon.

回复 原文
赠意 2024-08-10 06:45:33

截至目前，谷歌地图 utils 库提供了一种计算面积的方法。

您必须添加以下依赖项

implementation("com.google.maps.android:android-maps-utils:3.8.2")

，然后调用下面的方法并为其提供您的 LatLng 列表。它将返回以平方米为单位的面积。

fun calculateSurfaceArea(latLngList: List<LatLng>): Double {
    if (latLngList.size < 3) {
        return 0.0
    }
    return SphericalUtil.computeArea(latLngList)
}

As of now google maps utils library provide a method to calculate the Area.

You have to add below dependency

implementation("com.google.maps.android:android-maps-utils:3.8.2")

And then call below method and give it yours LatLng list. It will return area in square meters.

fun calculateSurfaceArea(latLngList: List<LatLng>): Double {
    if (latLngList.size < 3) {
        return 0.0
    }
    return SphericalUtil.computeArea(latLngList)
}
回复 原文
~没有更多了~

关于作者

口干舌燥

暂无简介

文章
评论
26 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

忆悲凉

文章 0 评论 0

hgfg1645

文章 0 评论 0

qq_qLPLYi

文章 0 评论 0

戏舞

文章 0 评论 0

殊姿

文章 0 评论 0

﹂绝世的画

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文