如何使用python计算地球表面多边形的面积?

2022-01-14 00:00:00 python geometry geospatial geolocation

问题描述

标题基本上说明了一切.我需要使用 Python 计算地球表面多边形内的面积.计算地球表面任意多边形包围的面积说关于它的一些东西，但在技术细节上仍然含糊不清:

The title basically says it all. I need to calculate the area inside a polygon on the Earth's surface using Python. Calculating area enclosed by arbitrary polygon on Earth's surface says something about it, but remains vague on the technical details:

如果你想用更多GIS"风味，那么你需要选择您所在地区的计量单位和找到一个合适的投影保留区域(并非所有人都这样做).自从你正在谈论计算任意多边形，我会使用类似兰伯特方位角的东西等面积投影.设置投影的原点/中心是多边形的中心，项目多边形到新坐标系统，然后使用计算面积标准平面技术.

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.

那么，我如何在 Python 中做到这一点?

So, how do I do this in Python?

解决方案

假设您有一个以 GeoJSON 格式表示的科罗拉多州

Let's say you have a representation of the state of Colorado in GeoJSON format

{"type": "Polygon", 
 "coordinates": [[
   [-102.05, 41.0], 
   [-102.05, 37.0], 
   [-109.05, 37.0], 
   [-109.05, 41.0]
 ]]}

所有坐标都是经度、纬度.您可以使用 pyproj 来投影坐标和 Shapely 找到任何投影多边形的面积:

All coordinates are longitude, latitude. You can use pyproj to project the coordinates and Shapely to find the area of any projected polygon:

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 表示，变成 Shapely 几何对象，并取面积:

That's an equal area projection centered on and bracketing the area of interest. Now make new projected GeoJSON representation, turn into a Shapely geometric object, and take the area:

x, y = pa(lon, lat)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
from shapely.geometry import shape
shape(cop).area  # 268952044107.43506

这是一个非常接近调查区域的近似值.对于更复杂的特征，您需要沿顶点之间的边缘进行采样，以获得准确的值.以上关于日期变更线等的所有警告均适用.如果您只对区域感兴趣，您可以在投影之前将您的要素从日期线移开.

It's a very close approximation to the surveyed area. For more complex features, you'll need to sample along the edges, between the vertices, to get accurate values. All caveats above about datelines, etc, apply. If you're only interested in area, you can translate your feature away from the dateline before projecting.

相关文章