Как рассчитать площадь многоугольника на поверхности земли с помощью питона?
Название в основном говорит обо всем. Мне нужно рассчитать площадь внутри многоугольника на поверхности Земли, используя Python. Вычисление площади, окруженной произвольным многоугольником на поверхности Земли, говорит о чем-то, но остается неясным в технических деталях:
Если вы хотите сделать это с более "ароматом ГИС", то вам нужно выбрать единицу измерения для вашей области и найти соответствующую проекцию, которая сохраняет область (не все это делают). Поскольку вы говорите о вычислении произвольного многоугольника, я бы использовал что-то вроде проекции азимутальной равной площади Ламберта. Установите начало / центр проекции как центр вашего многоугольника, спроецируйте многоугольник на новую систему координат, затем вычислите площадь, используя стандартные плоские методы.
Итак, как мне это сделать в Python?
10 ответов
Допустим, у вас есть представление о штате Колорадо в формате 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
Это очень близкое приближение к исследуемой области. Для более сложных функций вам нужно будет выполнить выборку по краям, между вершинами, чтобы получить точные значения. Применяются все оговорки, приведенные выше относительно сроков и т. Д. Если вас интересует только область, вы можете перевести свою функцию с линии времени перед проецированием.
Самый простой способ сделать это (на мой взгляд) - это проецировать вещи в (очень простую) проекцию равных площадей и использовать один из обычных плоских методов для вычисления площади.
Прежде всего, я собираюсь предположить, что сферическая земля достаточно близка для ваших целей, если вы задаете этот вопрос. Если нет, то вам нужно перепроектировать ваши данные с помощью соответствующего эллипсоида, и в этом случае вы захотите использовать реальную библиотеку проекций (в наши дни все использует proj4 за кулисами), например привязки python к GDAL/OGR или (гораздо более дружелюбный) 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
... возвращает площадь в квадратных метрах.
Возможно, немного поздно, но здесь другой метод, использующий теорему Жирара. В нем говорится, что площадь многоугольника больших кругов в R** в 2 раза больше суммы углов между многоугольниками минус (N-2)*pi, где 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 км, я получаю, что площадь составляет 268930758560.74808
Вы можете вычислить площадь непосредственно на сфере, вместо того, чтобы использовать проекцию равной площади.
Более того, согласно этому обсуждению, кажется, что теорема Жирара (ответ Сульке) не дает точных результатов в некоторых случаях, например, "область, окруженная луной 30º от полюса к полюсу и ограниченная нулевым меридианом и 30º в.д." (см. здесь).
Более точным решением было бы выполнить линейный интеграл непосредственно на сфере. Приведенное ниже сравнение показывает, что этот метод более точен.
Как и все другие ответы, я должен упомянуть оговорку, что мы предполагаем сферическую Землю, но я предполагаю, что для некритических целей этого достаточно.
Реализация Python
Вот реализация Python 3, которая использует линейный интеграл и теорему Грина:
def polygon_area(lats, lons, radius = 6378137):
"""
Computes area of spherical polygon, assuming spherical Earth.
Returns result in ratio of the sphere's area if the radius is specified.
Otherwise, in the units of provided radius.
lats and lons are in degrees.
"""
from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
lats = np.deg2rad(lats)
lons = np.deg2rad(lons)
# Line integral based on Green's Theorem, assumes spherical Earth
#close polygon
if lats[0]!=lats[-1]:
lats = append(lats, lats[0])
lons = append(lons, lons[0])
#colatitudes relative to (0,0)
a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
colat = 2*arctan2( sqrt(a), sqrt(1-a) )
#azimuths relative to (0,0)
az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)
# Calculate diffs
# daz = diff(az) % (2*pi)
daz = diff(az)
daz = (daz + pi) % (2 * pi) - pi
deltas=diff(colat)/2
colat=colat[0:-1]+deltas
# Perform integral
integrands = (1-cos(colat)) * daz
# Integrate
area = abs(sum(integrands))/(4*pi)
area = min(area,1-area)
if radius is not None: #return in units of radius
return area * 4*pi*radius**2
else: #return in ratio of sphere total area
return area
Я написал несколько более явную версию (и со многими ссылками более и несделанными...) в sphericalgeometry пакете там.
Численное сравнение
Колорадо будет эталоном, так как все предыдущие ответы были оценены на его территории. Его точная общая площадь составляет 104093,67 квадратных миль (по данным Бюро переписи населения США, стр. 89, см. Также здесь), или 269601367661 квадратный метр. Я не нашел источника фактической методологии USCB, но предполагаю, что она основана на суммировании фактических измерений на земле или точных расчетах с использованием WGS84/EGM2008.
Method | Author | Result | Variation from ground truth
--------------------------------------------------------------------------------
Albers Equal Area | sgillies | 268952044107 | -0.24%
Sinusoidal | J. Kington | 268885360163 | -0.26%
Girard's theorem | sulkeh | 268930758560 | -0.25%
Equal Area Cylindrical | Jason | 268993609651 | -0.22%
Line integral | Yellows | 269397764066 | **-0.07%**
Вывод: точнее использовать прямой интеграл.
Спектакль
Я не тестировал различные методы, и сравнение чистого кода Python с скомпилированными проекциями PROJ не имело бы смысла. Требуется интуитивно меньше вычислений. С другой стороны, тригонометрические функции могут требовать больших вычислительных ресурсов.
Вот решение, которое использует 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 в км ^2.
Я знаю, что ответ через 10 лет имеет некоторые преимущества, но тому, кто сегодня смотрит на этот вопрос, кажется справедливым дать обновленный ответ.
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 км2, или -0,17% от сообщенного правильного значения (269601,367661 км2).
Поскольку земля является замкнутой поверхностью, замкнутый многоугольник, нарисованный на ее поверхности, создает две многоугольные области. Вы также должны определить, кто из них внутри, а кто снаружи!
В большинстве случаев люди будут иметь дело с небольшими полигонами, и это "очевидно", но как только у вас появятся вещи размером с океан или континент, вам лучше убедиться, что вы все правильно поняли.
Также помните, что строки могут идти от (-179,0) до (+179,0) двумя различными способами. Один намного длиннее другого. Опять же, в основном вы будете исходить из предположения, что это линия, которая идет от (-179,0) до (-180,0), то есть (+180,0), а затем до (+179,0), но одна день... не будет.
Обращение с широтой как к простой (x,y) системе координат или даже пренебрежение тем фактом, что любая проекция координат будет иметь искажения и разрывы, может привести к серьезным ошибкам в сферах.
Вот реализация Python 3, в которой функция будет принимать список пар кортежей lats и long и возвращать область, заключенную в проецируемом многоугольнике. Он использует pyproj для проецирования координат, а затем Shapely, чтобы найти площадь любого проецируемого многоугольника.
def calc_area(lis_lats_lons):
import numpy as np
from pyproj import Proj
from shapely.geometry import shape
lons, lats = zip(*lis_lats_lons)
ll = list(set(lats))[::-1]
var = []
for i in range(len(ll)):
var.append('lat_' + str(i+1))
st = ""
for v, l in zip(var,ll):
st = st + str(v) + "=" + str(l) +" "+ "+"
st = st +"lat_0="+ str(np.mean(ll)) + " "+ "+" + "lon_0" +"=" + str(np.mean(lons))
tx = "+proj=aea +" + st
pa = Proj(tx)
x, y = pa(lons, lats)
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
return shape(cop).area
Для выборочного набора латов / лонгов он дает значение площади, близкое к исследуемому приближенному значению.
calc_area(lis_lats_lons = [(-102.05, 41.0),
(-102.05, 37.0),
(-109.05, 37.0),
(-109.05, 41.0)])
На выходе получается площадь 268952044107,4342 кв. Mts.
Согласно утверждению Желтого, прямой интеграл более точен.
Но желтые используют радиус Земли = 6378 137 м, который является эллипсоидом WGS-84, большая полуось, в то время как Сульке использует 6371 000 м.
Используя радиус = 6378 137 м в методе Сульке, получаем 269533625893 квадратных метра.
Если предположить, что истинное значение площади Колорадо (по данным Бюро переписи населения США) составляет 269601367661 квадратный метр, то отклонение от наземной истинности метода Сульке составляет: -0,025%, что лучше, чем -0,07 при использовании метода линейного интеграла.
Так что предложение Сульке пока кажется более точным.
Чтобы иметь возможность провести численное сравнение решений с предположением о сферической Земле, все расчеты должны использовать один и тот же земной радиус.