在 Python 中从纬度/经度多边形计算面积

Posted 2023-02-23

【中文标题】在 Python 中从纬度/经度多边形计算面积

【英文标题】：Calculating area from lat/lon polygons in Python

【发布时间】：2012-10-30 21:45:41

【问题描述】：

基本上，我想知道 0.25°lat x 0.25°lon 补丁在世界各地的各种多边形中拟合了多少次。后者的尺寸约为 3°lat x 10°lon 或 2°lat x 4°lon。

我有多边形角的纬度/经度值，我正在计算它们的面积，如下所示：

from pyproj import Proj
from shapely.geometry import shape

def getArea(coords):
    c = "type": "Polygon",
    "coordinates": [[ (coords[0], coords[2]), (coords[1], coords[2]),
                      (coords[0], coords[3]), (coords[1], coords[3]) ]]
    lon, lat = zip(*c['coordinates'][0])
    pro = Proj("+proj=aea")
    x, y = pro(lon, lat)
    poly = "type": "Polygon", "coordinates": [zip(x, y)]
    return shape(cop).area

我从这里采取了方法：How to calculate the area of a polygon on the earth's surface using python?

现在的问题是，我应该选择哪个等面积投影才能使多边形具有可比较的面积大小。在这种投影中，无论它位于地球上的哪个位置，小块的面积始终相同。

采用 Albers 等面积投影 (aea) 得到三个多边形的这些区域：

240993868.90978813 699931593.1047173 212092562.5238676

采用 Lambert 方位角等面积投影 (laea) 得到相同多边形的这些区域：

148709452.69292444 409253749.5468254 106218747.36092758

为什么两个投影中的面积关系不同？前 1:3 = 0.344；第二个 1:3 = 0.363；它们应该是相同的，因为它们都是相等的面积投影？！

这让我想知道将小块与任一投影中的多边形区域进行比较是否合法。你有什么建议吗？

【问题讨论】：

2°lat x 3°lon 相当大，任何错误都会成比例地变大。有不同的投影的原因是没有投影是完美的，它根本不可能。

我更新了补丁/多边形的大小只是为了记录。嗯，所以你说这种方法的错误越来越大？我还能尝试什么？

完美地做到这一点的唯一方法是在不投影的情况下计算球体切口的表面积。否则，您将不得不根据对您很重要的属性来决定投影。

“2x3 补丁的面积始终相同” - 不正确。经线度数约为 111 公里，赤道附近的纬度为 1 度。但随着位置从赤道移动到极地，后者变得更小。

也许你可以从this question找到有用的东西

【参考方案1】：

如果您关心实际的相对面积，计算瓷砖不会给您正确的答案。 如果您需要实际表面积，请使用椭球或球面几何，或以下想法。

蛮力矢量方法：将全球 0.25 度网格生成为多边形，与多边形相交，计算结果。每个多边形可以是笛卡尔正方形或实际的椭圆正方形。为每个图块创建自定义理想投影，然后计算面积并将其存储为每个图块的属性。你只需要这样做一次:)