and*_*s-h 30 python geometry geolocation geospatial
标题基本上都说明了一切.我需要使用Python计算地球表面多边形内的区域.计算由地球表面任意多边形包围的区域说明了一些内容,但对技术细节仍然模糊:
如果你想用更"GIS"的味道来做这件事,那么你需要为你的区域选择一个度量单位,找到一个保留区域的适当投影(不是所有的).既然你在谈论计算任意多边形,我会使用类似Lambert Azimuthal等面积投影的东西.将投影的原点/中心设置为多边形的中心,将多边形投影到新坐标系,然后使用标准平面技术计算面积.
那么,我如何在Python中执行此操作？
sgi*_*ies 36
假设您以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表示,转换为Shapely几何对象,并取以下区域:
x, y = pa(lon, lat)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
from shapely.geometry import shape
shape(cop).area # 268952044107.43506
它与调查区域非常接近.对于更复杂的特征,您需要沿顶点之间的边缘采样,以获得准确的值.以上关于日期线等的所有警告均适用.如果您只对区域感兴趣,则可以在投影前将您的功能从日期线转换出来.
Joe*_*ton 24
(在我看来)最简单的方法是将事物投射到(一个非常简单的)等面积投影中,并使用一种常用的平面技术来计算面积.
首先,如果你问这个问题,我会假设球形地球足够接近你的目的.如果没有,那么你需要使用适当的椭球重新投影你的数据,在这种情况下你将需要使用一个实际的投影库(这些天在幕后使用proj4),例如对GDAL/OGR的python绑定或者(更友好的)pyproj.
但是,如果你对球形地球没问题,那么没有任何专门的库就可以很容易地做到这一点.
要计算的最简单的等面积投影是正弦投影.基本上,您只需将纬度乘以一个纬度的长度,经度乘以纬度的长度和纬度的余弦.
def reproject(latitude, longitude):
"""Returns the x & y coordinates in meters using a sinusoidal projection"""
from math import pi, cos, radians
earth_radius = 6371009 # in meters
lat_dist = pi * earth_radius / 180.0
y = [lat * lat_dist for lat in latitude]
x = [long * lat_dist * cos(radians(lat))
for lat, long in zip(latitude, longitude)]
return x, y
好的......现在我们要做的就是计算平面中任意多边形的面积.
有很多方法可以做到这一点.我将使用这里最常见的一个.
def area_of_polygon(x, y):
"""Calculates the area of an arbitrary polygon given its verticies"""
area = 0.0
for i in range(-1, len(x)-1):
area += x[i] * (y[i+1] - y[i-1])
return abs(area) / 2.0
无论如何,希望这将指向正确的方向......
或者干脆使用一个库：https : //github.com/scisco/area
from area import area
>>> obj = {'type':'Polygon','coordinates':[[[-180,-90],[-180,90],[180,90],[180,-90],[-180,-90]]]}
>>> area(obj)
511207893395811.06
...以平方米为单位返回面积。
也许有点晚了,但这是一种不同的方法,使用吉拉德定理.它指出,大圆形多边形的面积是多边形之间的角度之和减去(N-2)*pi的R**2倍,其中N是角的数量.
我认为这是值得发布的,因为它不依赖于任何其他库而不是numpy,并且它是一种与其他库完全不同的方法.当然,这仅适用于球体,因此在将其应用于地球时会有一些不准确之处.
首先,我定义了一个函数来计算从点1沿大圆到第2点的方位角:
import numpy as np
from numpy import cos, sin, arctan2
d2r = np.pi/180
def greatCircleBearing(lon1, lat1, lon2, lat2):
dLong = lon1 - lon2
s = cos(d2r*lat2)*sin(d2r*dLong)
c = cos(d2r*lat1)*sin(d2r*lat2) - sin(lat1*d2r)*cos(d2r*lat2)*cos(d2r*dLong)
return np.arctan2(s, c)
现在我可以用它来找到角度,然后是区域(在下面,当然应该指定lons和lats,它们应该是正确的顺序.另外,应该指定球体的半径.)
N = len(lons)
angles = np.empty(N)
for i in range(N):
phiB1, phiA, phiB2 = np.roll(lats, i)[:3]
LB1, LA, LB2 = np.roll(lons, i)[:3]
# calculate angle with north (eastward)
beta1 = greatCircleBearing(LA, phiA, LB1, phiB1)
beta2 = greatCircleBearing(LA, phiA, LB2, phiB2)
# calculate angle between the polygons and add to angle array
angles[i] = np.arccos(cos(-beta1)*cos(-beta2) + sin(-beta1)*sin(-beta2))
area = (sum(angles) - (N-2)*np.pi)*R**2
在另一个回复中给出科罗拉多坐标,并且地球半径为6371 km,我得到的区域是268930758560.74808
这是一个使用basemap, 而不是pyproj和shapely进行坐标转换的解决方案。这个想法与@sgillies 的建议相同。请注意，我已经添加了第 5 个点，因此路径是一个闭环。
import numpy
from mpl_toolkits.basemap import Basemap
coordinates=numpy.array([
[-102.05, 41.0],
[-102.05, 37.0],
[-109.05, 37.0],
[-109.05, 41.0],
[-102.05, 41.0]])
lats=coordinates[:,1]
lons=coordinates[:,0]
lat1=numpy.min(lats)
lat2=numpy.max(lats)
lon1=numpy.min(lons)
lon2=numpy.max(lons)
bmap=Basemap(projection='cea',llcrnrlat=lat1,llcrnrlon=lon1,urcrnrlat=lat2,urcrnrlon=lon2)
xs,ys=bmap(lons,lats)
area=numpy.abs(0.5*numpy.sum(ys[:-1]*numpy.diff(xs)-xs[:-1]*numpy.diff(ys)))
area=area/1e6
print area
结果是 268993.609651，以 km^2 为单位。
更新：底图已被弃用，因此您可能需要首先考虑替代解决方案。
您可以直接在球体上计算面积，而不是使用等积投影。
\n\n此外，根据此讨论，似乎吉拉德定理（sulkeh 的答案）在某些情况下并没有给出准确的结果，例如“从极点到极点由 30\xc2\xba 月牙形包围的区域和以本初子午线和 30\xc2\xbaE" 为界（参见此处）。
\n\n更精确的解决方案是直接在球体上执行线积分。从下面的比较可以看出，这种方法更加精确。
\n\n与所有其他答案一样，我应该提到我们假设地球是球形的警告，但我认为对于非关键目的这已经足够了。
\n\n这是一个使用线积分和格林定理的 Python 3 实现：
\n\ndef polygon_area(lats, lons, radius = 6378137):\n """\n Computes area of spherical polygon, assuming spherical Earth. \n Returns result in ratio of the sphere\'s area if the radius is specified.\n Otherwise, in the units of provided radius.\n lats and lons are in degrees.\n """\n from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad\n lats = np.deg2rad(lats)\n lons = np.deg2rad(lons)\n\n # Line integral based on Green\'s Theorem, assumes spherical Earth\n\n #close polygon\n if lats[0]!=lats[-1]:\n lats = append(lats, lats[0])\n lons = append(lons, lons[0])\n\n #colatitudes relative to (0,0)\n a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2\n colat = 2*arctan2( sqrt(a), sqrt(1-a) )\n\n #azimuths relative to (0,0)\n az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)\n\n # Calculate diffs\n # daz = diff(az) % (2*pi)\n daz = diff(az)\n daz = (daz + pi) % (2 * pi) - pi\n\n deltas=diff(colat)/2\n colat=colat[0:-1]+deltas\n\n # Perform integral\n integrands = (1-cos(colat)) * daz\n\n # Integrate \n area = abs(sum(integrands))/(4*pi)\n\n area = min(area,1-area)\n if radius is not None: #return in units of radius\n return area * 4*pi*radius**2\n else: #return in ratio of sphere total area\n return area\n我在sphericalgeometry包中编写了一个更明确的版本（并且有更多参考文献和 TODO...）。
\n\n科罗拉多州将成为参考，因为之前的所有答案都是在其地区进行评估的。其精确总面积为 104,093.67 平方英里（来自美国人口普查局，第 89 页，另见此处），即 269601367661 平方米。我没有找到 USCB 实际方法的来源，但我认为它是基于对地面实际测量值进行求和，或使用 WGS84/EGM2008 进行精确计算。
\n\nMethod | Author | Result | Variation from ground truth\n--------------------------------------------------------------------------------\nAlbers Equal Area | sgillies | 268952044107 | -0.24%\nSinusoidal | J. Kington | 268885360163 | -0.26%\nGirard\'s theorem | sulkeh | 268930758560 | -0.25%\nEqual Area Cylindrical | Jason | 268993609651 | -0.22%\nLine integral | Yellows | 269397764066 | **-0.07%**\n结论：使用直接积分更精确。
\n\n我还没有对不同的方法进行基准测试，将纯 Python 代码与编译的 PROJ 投影进行比较是没有意义的。直观上需要更少的计算。另一方面，三角函数可能需要大量计算。
\n我知道十年后回答有一些优势，但对于今天看这个问题的人来说，提供更新的答案似乎是公平的。
pyproj直接计算面积，无需调用shapely：
# Modules:
from pyproj import Geod
import numpy as np
# Define WGS84 as CRS:
geod = Geod('+a=6378137 +f=0.0033528106647475126')
# Data for Colorado (no need to close the polygon):
coordinates = np.array([
[-102.05, 41.0],
[-102.05, 37.0],
[-109.05, 37.0],
[-109.05, 41.0]])
lats = coordinates[:,1]
lons = coordinates[:,0]
# Compute:
area, perim = geod.polygon_area_perimeter(lons, lats)
print(abs(area)) # Positive is counterclockwise, the data is clockwise.
结果为：269154.54988400977 km2，或报告正确值 (269601.367661 km2) 的 -0.17%。
| 归档时间： |
|
| 查看次数： |
24796 次 |
| 最近记录： |