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

Question

基本上，我想知道一个0.25°纬度x 0.25°经度的区域在世界各地各种多边形中适合多少次。后者的大小大约为3°纬度x 10°经度或2°纬度x 4°经度。

我有多边形角落的纬度/经度值，并像这样计算它们的面积：

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

我采用了这里的方法：如何使用Python计算地球表面多边形的面积？ 现在的问题是，为了让多边形的面积大小可比较，我应该选择哪种等面积投影。在这样的投影中，小区域的面积始终相同，无论它位于地球的哪个位置。

采用 Albers等面积投影 (aea) 的结果如下所示：

1. 240993868.90978813
2. 699931593.1047173
3. 212092562.5238676

采用 Lambert方位等面积投影(laea) 的结果如下所示：

1. 148709452.69292444
2. 409253749.5468254
3. 106218747.36092758

为什么两种投影的面积之间的关系不同？第一组比例为1:3 = 0.344；第二组比例为1:3 = 0.363。由于两者都是等面积投影，它们应该是相同的？！

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

- HyperCube

---

2°纬度 x 3°经度相当大，任何错误都将成比例放大。不同的投影方式之所以存在是因为没有一种投影方式是完美的，它根本无法做到这一点。 - Mark Ransom

我只是为了记录更新了补丁/多边形的大小。嗯，你是说这种方法误差变得太大了？还有什么其他的尝试方法吗？ - HyperCube

唯一完美的方法是计算球体切口的表面积，而不是投影。否则，您将不得不根据对您重要的属性来决定投影。 - Mark Ransom

“2x3补丁区域的面积始终相同” - 这是不正确的。子午线长度大约为111公里，赤道附近的一度纬度也是如此。但随着位置从赤道移向极地，后者变得越来越小。 - Igor

也许你可以从这个问题中找到一些有用的东西。 - nadya

Answer 1

如果您关心实际相对面积，仅计算瓷砖数量不能给出正确答案。如果您需要实际表面积，请使用椭球几何、球形几何或以下思路。

通用向量法：生成全球范围的0.25度网格作为多边形，并与多边形相交，计算结果。每个多边形可以是笛卡尔平方或实际椭圆形。为每个瓦片创建自定义理想投影，然后计算和存储每个瓦片的面积作为属性。您只需要这样做一次 :)

- Chris