问题
标题基本上说明了一切。我需要在 Python 中计算地球表面多边形内的面积。计算由任意多边形包围的地球表面的面积说明了这一点,但在技术细节上仍然含糊不清:
那么,我怎样才能在 Python 中实现这一点呢?
回答
假设您有一个以 GeoJSON 格式表示的科罗拉多州
- {"type": "Polygon",
- "coordinates": [[
- [-102.05, 41.0],
- [-102.05, 37.0],
- [-109.05, 37.0],
- [-109.05, 41.0]
- ]]}
复制代码
所有坐标都是经度,纬度。您可以使用 pyproj 投影坐标,并使用 Shapely 查找任何投影多边形的面积:
- co = {"type": "Polygon", "coordinates": [
- [(-102.05, 41.0),
- (-102.05, 37.0),
- (-109.05, 37.0),
- (-109.05, 41.0)]]}
- lon, lat = zip(*co['coordinates'][0])
- from pyproj import Proj
- pa = Proj("+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55")
复制代码
这是一个以感兴趣区域为中心并用括号括起来的等面积投影。现在创建一个新的投影 GeoJSON 表示,将其转换为形状几何对象,并获取面积:
- x, y = pa(lon, lat)
- cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
- from shapely.geometry import shape
- shape(cop).area # 268952044107.43506
复制代码
这是一个非常接近的测量区域。对于更复杂的特征,您需要在顶点之间沿边进行采样,以获得准确的值。以上关于日期线等的所有警告都适用。如果您只对该区域感兴趣,则可以在投影之前将要素转换为日期线之外。
|