Здравствуйте, столкнулся с довольно простой задачей, но быстрого решения не нашел.
Поэтому возникла идея сделать тему… для решения подобных вопросов.
Вобщем суть в том, что зная все географические координаты нодов(точек) вея (полигона) мне хотелось бы получить его площадь в метрах квадратных.
На данный момент я нашел формулу для расчета кусочка поверхности сферы для “прямоугольного” случая.
Вот код… на питоне. Возможно оно глючное, пока просто набросок. Тут статейка про геометрию сферы… вроде бы все довольно просто.
Используя библиотеку Shapely можно получить прямоугольник ограничиващий любой вей и затем посчитать его площадь… числа вроде бы адекватные)
Но хотелось бы найти нормальную библиотеку для питона, дающие корректные вычисления в метрах квадратных.
Shapely считает площадь полигона, но т.к. координаты заданы в углах, то и площадь получается фиг знает в чем… углы квадратные))
Как их можно быстро и правильно привести к метрам? Я так полагаю нужна другая система координат, но какая?
Попытаюсь задать координаты в метрах от (0,0). Дистанцию умеет считать geopy в разных моделях Земли (нужная WGS-84?).
Именно на питоне? У Лёши в перловском osm2mp площадь волшебным образом вычисляется с использованием стандартных перловских же библиотек. Хотя что там может получиться с проекциями и прочим, фиг знает
Подозреваю, что все мыслимые географические алгоритмы давным-давно реализованы, проверены, отлажены и упакованы в библиотеки.
ой, чегой-то мне говорит что там все волшебство- сперва площадь считается в квадратных градусах (что элементарно), а потом пересчитывается в кв. км. считая землю шаром с длиной большой окружности в 40000 км ровно
Для распределения объектов по масштабным уровням - вполне годный алгоритм)
Да-да! Там (где-то в недрах Перла) даже алгоритм пересечения двух отрезков был кривой. Так что верны только первый и четвертый постулат, но не второй и третий
мда, как оказалось задачка не такая простая…
по большей части нужно как раз для распределения объектов по масштабным уровням.
Но хотелось бы и просто иметь такой алгоритм под рукой с необходимой точностью)
gps-Max, не обязательно на питоне, пойдет и на С++, и на яве… просто может где описание какое есть нормальное. Сам ничего токового нагуглить не смог(
Просто импорт, экспорт я замутил питоновскими скриптами… вдоволь намучавшись с SQL базами я переполз на NoSQL MongoDB. Пытаюсь сделать уровни детализации в векторе.
Zkir, как именно пересчитываются квадратные градусы в метры квадратные? у них же разная длина в метрах… если от экватора подниматься длина градуса одного угла все время уменьшается.
А если еще и несферичность Земли участь…
В osm2mp площадь объекта считается в обычной географической проекции, и умножается на коэффициент - косинус широты. Для первого приближения этого вполне достаточно.
Для несферической земли можно так:
S = lat / 360 * <длина экватора> x lon*cos(lat) / 360 * <длина меридиана>
/* Сразу привести широту и долготу точек полигона в метры, т.е. в очень приблизительную прямоуголку. Можно даже не заморачиваться с центральным меридианом, т.к. нас не волнуют абсолютные значения. / R = 40000000.0 / 2.0 / 3.14159265359 / радиус сферы Земли, получается что-то среднее между большой и малой полуосями WGS-84 / y[n] = radians(lat[n]) * R
x[n] = radians(lon[n]) * cos(radians(lat[n])) * R
/ В цикле, перебирая пары координат, по одному отрезку из полигона, считаем площадь трапеции, ограниченной текущей и последующей точками, их меридианами и экватором (проще нарисовать, чем описать). */ dS = (x[n+1] - x[n]) * (y[n+1] + y[n]) / 2.0
S += dS
n++
ПРИМЕЧАНИЯ:
площадь будет иметь знак;
не забудьте замкнуть полигон, в противном случае значение будет неверное;
для больших площадей надо считать площадь трапеции на сфере.
ИМХО, достоинство этого метода:
считается площадь всего полигона сразу;
получившаяся площадь имеет положительное значение, если полигон по часовой стрелке, и отрицательная, если против часовой;
Недостатки:
площадь самопересекающегося полигона может принять неверное значение, включая 0 (например для фигуры “восьмерки”);
при больших площадях будут сильные искажения (например площадь России так не посчитаешь) (см. новый алгоритм ниже);
надо отслеживать работу на Чукотке при переходе долготы от 180 к -180, в этом случае положительную восточную долготу не трогаем, а к западной долготе прибавляем 360.0 градусов;
при переходе через экватор возможны глюки.
Добавление:
Навскидку получилась следующая формула для больших площадей, в этом случае переводить широту-долготу в метры не надо: R = 40000000.0 / 2.0 / 3.14159265359
/* в цикле, перебираем все точки, не забывая замкнуть полигон */ dS = R * R * radians(lon[n+1] - lon[n]) * sin( radians(lat[n+1] + lat[n]) / 2.0 )
S += dS
n++
Пока потестил вот этот исходничок… почти такой же как и этот
Обе реализации показывают довольно внятные числа… в смысле если и ожидаю результат в сотни квадратных метров, то результат хотя бы имеет тот же порядок)
Выбрал на карте прямоугольник, взял его углы и пропустил через ту функцию… выдает 10270 метров квадратных.
Участок был со сторонами 110м х 207м… что по площади в два раза больше! В принципе уже радует, что хоть такая оценка есть!
Насколько я понимаю с невыпуклыми полигонами алгоритмы могут и не справится… наверно триангуляция нужна)))
сделаю реализацию на питоне - выложу)
Зря значит я тут интегрированием занимался Попробуйте второй метод из моего сообщения, должен работать на любых полигонах без самопересечения, полностью находящихся в одном полушарии (либо северном, либо южном). Точность будет доли процентов, а не в два раза…
chnav, Труд Ваш не пропадет в любом случае!)
просто там чего-то готовенькое… скомпилил по-быстрому да и все…
Вот реализация Вашего алгоритма:
import math
R = 40000000.0 / 2.0 / 3.14159265359
lat = [55.728634, 55.729489, 55.729440, 55.728565, 55.728634];
lon = [37.826026, 37.826097, 37.827862, 37.827780, 37.826026];
S = 0
for n in range(len(lat)-1):
dS = R * R * math.radians(lon[n+1] - lon[n]) * math.sin( math.radians(lat[n+1] + lat[n]) / 2.0 )
S = S+dS
print S
Все правильно?
Результат 10612.008839 я так понимаю в метрах)
немного расходится с моими прошлыми результатами для этого же участка… должно быть около 22880, если я ничего не напутал)
мерил по Яндексовской карте их рулеткой, проверял на забугорных сайтах, которые дистанцию по координатам считают…
надеюсь в цифрах я не опечатался, когда координаты вбивал
Komяpa, в радианы переводил все прально!
на счет косинусов и прочей математики это я уже не соображу так просто… я готовые исходнички проверил особо не разбираясь)
пересчитать теми же треугольниками площадь, взять по модулю
помножить на косинус средней широты. (математики, поправьте: с длинами-ширинами прокатывает, не становится ли косинус косинусом-квадратом при уходе в площадь?)
x10kHz
Загрузил ваш пример в GPSMapEdit, он дает площадь 10605 кв.м, т.е. мой алгоритм кривоват 0.06%. ИМХО это много для такого маленького объекта. Скорее всего сыграло роль что мы считаем на сфере очень примерного радиуса, а не на эллипсоиде как мапэдит.
