Нужно рассчитать площадь полигона, загруженного из OSM. Как лучше (правильнее) это сделать на python?
Полигон грузится из ОСМ ХМЛ? Очевидно, перевести координаты в метрическую систему, ну а далее геометрия на плоскости
а можно поподробнее (желательно со ссылками на формулы или примеры) ?
Да, спасибо, примерно это мне и требуется
Но я все еще не уяснил для себя одни важный момент: куда надо перевести родные осмовские координаты, чтобы можно было считать реальную площадь в в метрах квадратных?
Смотря какие у вас площади и на каких широтах.
Ну вот я смотрел чужие примеры, там использовали EPSG 900913 для отображения (если я не ошибкся).
От меня ссылки, формулы, примеры только завтра, когда протон будет род рукой. Куда переводить, уже писал, в метрическую систему координат. А уж как она в коде задаётся с мрбилы с ходу не скажу…
В равноплощадную проекцию, ака плоскую систему координат, очевидно, или мучиться сферической тригонометрией.
Если площади не очень большие, а широты не слишком высокие, можно просто посчитать площадь в квадратных градусах, а потом перевести в метры квадратные, принимая градус за 111319 м и умножив на косинус широты.
Самое простое (но не самое точное) - конвертировать в меркатор, и считать там (не забывая умножать/делить на косинус широты в нужных местах).
Чем больше размеры - тем больше будет погрешность, на размерах в сотни км этот способ лучше не применять.
Предлагаю перевести координаты в трансверсивный меркатор с динамически задаваемым центральным меридианом. Около центрального меридиана искажений нет!
В качестве centerLat и centerLon можно взять центр тяжести многоугольника, т.е. среднее каждой координаты.
Код на питоне здесь
Пример использования
from transverse_mercator import TransverseMercator
projection = TransverseMercator(lat=centerLat, lon=centerLon)
# конвертация из прямоугольной системы координат в географическую систему координат:
(lat, lon) = projection.toGeographic(x, y)
# конвертация из географической системы координат в прямоугольную систему координат:
(x, y) = projection.fromGeographic(lat, lon)
ооо… во я вляпался
значит, чтобы получить более точный результат на больших площадях (тыс км2), мне нужно погрузится в трехмерный тригонометрический ад с эллипсоидами и косинусами
…либо выбирать проекции, которые дадут приемлимую погрешность
vvoovv: я правильно понимаю, что на входе используются осмовские координаты, а на выходе будут метры, т.е. никаких дополнительных преобразований не понадобится?
трехмерный ад оказался проще, чем я ожидал
спасибо всем за помощь
вот так выглядят два варианта варианта расчета: напрямую и через проекцию (если, конечно же, я нигде не напортачил)
#! /usr/bin/python2
from geographiclib.geodesic import Geodesic
from transverse_mercator import TransverseMercator
from shapely.geometry.polygon import Polygon
def p(lat,lon):
return {'lat': lat, 'lon': lon}
print("surface area calculation")
poly = [
[56.7169634463617, 27.552240923220605],
[56.82513609211995, 31.61305255781721],
[58.938431723756175, 31.152451238152313],
[58.821834626837564, 27.420640546173495]
]
center = [57.884769367436945, 29.417133057759603]
area1 = Geodesic.WGS84.Area([p(lat,lon) for (lat,lon) in poly])
print("geographic library {:.2f} km2".format(area1["area"]/1000000))
projection = TransverseMercator(lat=center[0],lon=center[1])
poly_projected = [projection.fromGeographic(lat, lon) for (lat,lon) in poly]
poly_geom = Polygon(poly_projected)
print("transverse mercator {:.2f} km2".format(poly_geom.area/1000000))
результат
surface area calculation
geographic library 54578.87 km2
transverse mercator 54421.12 km2
JOSM (плагин measurment) выдал цифру 54313.134
Учтите что “прямая” (точнее дуга большого круга) на сфере/эллипсоиде и прямая в прямоугольной проекции - это разные линии, соответственно разные контуры, разные площади.
Во-вторых в проекции TM существует масштаб как линейный, так и для площадных объектов. Прямоугольная проекция годится только для небольших объектов, а тут трапеция шириной 4 градуса…
Впрочем эти замечания исключительно с целью пояснить разницу в полученных результатах.
Для ОСМ-задач, думаю, годятся любые способы т.к. избыточная точность ни к чему, учитывая точность входного материала.
Да, правильно.
Мой способ имеет смысл, конечно, для небольших объектов, например, площадь контура здания, придомового участка, квартала в городе. Если объект простирается на километры, то погрешность будет заметной.