Как рассчитать площадь многоугольника на поверхности земли с помощью питона?

Question

Название в основном говорит обо всем. Мне нужно рассчитать площадь внутри многоугольника на поверхности Земли, используя Python. Расчет площади, окруженной произвольным многоугольником на поверхности Земли что-то говорит об этом, но остается неясным в технических деталях:

Если вы хотите сделать это с более «ГИС» ароматизатор, тогда нужно выбрать единица измерения для вашего района и найти подходящий прогноз, который заповедники области (не все так делают). С вами говорим о расчете произвольный многоугольник, я бы использовал что-то вроде азимутального Ламберта Равная площадь проекции. Установить начало / центр проекции, которая будет центр вашего многоугольника, проект полигон с новой координатой система, а затем рассчитать площадь, используя стандартные планарные техники.

Итак, как мне это сделать в Python?

Answer 1

Допустим, у вас есть представление о штате Колорадо в формате 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

Это очень близкое приближение к исследуемой области.Для более сложных функций вам нужно будет выполнить выборку по краям, между вершинами, чтобы получить точные значения.Применяются все оговорки, приведенные выше о сроках и т. Д.Если вас интересует только область, вы можете перевести свою функцию с даты до проецирования.

Answer 2

Самый простой способ сделать это (на мой взгляд) - это проецировать вещи в (очень простую) проекцию равной площади и использовать один из обычных плоских методов для расчета площади.

Прежде всего,Я собираюсь предположить, что сферическая земля достаточно близка для ваших целей, если вы задаете этот вопрос.Если нет, то вам необходимо перепроецировать ваши данные, используя соответствующий эллипсоид, и в этом случае вы захотите использовать реальную библиотеку проекций (в наши дни все использует 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

Надеюсь, это все равно укажет вам правильное направление ...

Answer 3

Возможно, немного поздно, но здесь другой метод, использующий теорему Жирара.В нем говорится, что площадь многоугольника больших кругов 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

Answer 4

Поскольку Земля представляет собой замкнутую поверхность, замкнутый многоугольник, нарисованный на ее поверхности, создает ДВА полигональных области.Вам также нужно определить, какая из них находится внутри, а какая снаружи!

В большинстве случаев люди будут иметь дело с небольшими полигонами, и поэтому это «очевидно», но если у вас есть вещи размером с океан или континент, вам лучшеубедитесь, что вы все правильно поняли.

Также помните, что линии могут идти от (-179,0) до (+179,0) двумя различными способами.Один намного длиннее другого.Опять же, в основном вы будете исходить из предположения, что это линия, которая идет от (-179,0) до (-180,0), то есть (+180,0), а затем до (+179,0), но однадень ... этого не будет.

Если рассматривать широту-долготу как простую (x, y) систему координат или даже игнорировать тот факт, что любая проекция координат будет иметь искажения и разрывы, это может заставить васпровалиться на шарах.

Answer 5

Или просто используйте библиотеку: 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

... возвращает площадь в квадратных метрах.

Answer 6

Вот решение, которое использует 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.