找回密码
 立即注册
问题
标题基本上说明了一切。我需要在 Python 中计算地球表面多边形内的面积。计算由任意多边形包围的地球表面的面积说明了这一点,但在技术细节上仍然含糊不清:

那么,我怎样才能在 Python 中实现这一点呢?

回答
假设您有一个以 GeoJSON 格式表示的科罗拉多州
  1. {"type": "Polygon",
  2. "coordinates": [[
  3.    [-102.05, 41.0],
  4.    [-102.05, 37.0],
  5.    [-109.05, 37.0],
  6.    [-109.05, 41.0]
  7. ]]}
复制代码

所有坐标都是经度,纬度。您可以使用 pyproj 投影坐标,并使用 Shapely 查找任何投影多边形的面积:
  1. co = {"type": "Polygon", "coordinates": [
  2.     [(-102.05, 41.0),
  3.      (-102.05, 37.0),
  4.      (-109.05, 37.0),
  5.      (-109.05, 41.0)]]}
  6. lon, lat = zip(*co['coordinates'][0])
  7. from pyproj import Proj
  8. pa = Proj("+proj=aea +lat_1=37.0 +lat_2=41.0 +lat_0=39.0 +lon_0=-106.55")
复制代码

这是一个以感兴趣区域为中心并用括号括起来的等面积投影。现在创建一个新的投影 GeoJSON 表示,将其转换为形状几何对象,并获取面积:
  1. x, y = pa(lon, lat)
  2. cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
  3. from shapely.geometry import shape
  4. shape(cop).area  # 268952044107.43506
复制代码

这是一个非常接近的测量区域。对于更复杂的特征,您需要在顶点之间沿边进行采样,以获得准确的值。以上关于日期线等的所有警告都适用。如果您只对该区域感兴趣,则可以在投影之前将要素转换为日期线之外。





上一篇:如何判断我的 R 脚本中没有使用哪些包?
下一篇:悬停时缩放图像但不超过父 div 边框