Географические алгоритмы

x10kHz · 2010-09-06 19:58:57

Здравствуйте, столкнулся с довольно простой задачей, но быстрого решения не нашел.
Поэтому возникла идея сделать тему... для решения подобных вопросов.

Вобщем суть в том, что зная все географические координаты нодов(точек) вея (полигона) мне хотелось бы получить его площадь в метрах квадратных.

На данный момент я нашел формулу для расчета кусочка поверхности сферы для "прямоугольного" случая.
Вот код... на питоне. Возможно оно глючное, пока просто набросок. Тут статейка про геометрию сферы... вроде бы все довольно просто.

import math

def boundSquare(bounds):
    EARTH_RADIUS = 6370000.
    S1 = math.copysign(1-math.cos(math.radians(bounds[0])), bounds[0])
    S2 = math.copysign(1-math.cos(math.radians(bounds[2])), bounds[2])
    return math.fabs(2*math.pi*pow(EARTH_RADIUS,2)*(S1-S2)*(bounds[1]-bounds[3])/360.0)

Используя библиотеку Shapely можно получить прямоугольник ограничиващий любой вей и затем посчитать его площадь... числа вроде бы адекватные)

Но хотелось бы найти нормальную библиотеку для питона, дающие корректные вычисления в метрах квадратных.
Shapely считает площадь полигона, но т.к. координаты заданы в углах, то и площадь получается фиг знает в чем... углы квадратные))
Как их можно быстро и правильно привести к метрам? Я так полагаю нужна другая система координат, но какая?
Попытаюсь задать координаты в метрах от (0,0). Дистанцию умеет считать geopy в разных моделях Земли (нужная WGS-84?).

Last edited by x10kHz (2010-09-06 20:01:08)

gps-Max · 2010-09-06 22:29:17

Именно на питоне? У Лёши в перловском osm2mp площадь волшебным образом вычисляется с использованием стандартных перловских же библиотек. Хотя что там может получиться с проекциями и прочим, фиг знает :-)

Подозреваю, что все мыслимые географические алгоритмы давным-давно реализованы, проверены, отлажены и упакованы в библиотеки.

Zkir · 2010-09-06 22:41:03

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

Для распределения объектов по масштабным уровням - вполне годный алгоритм)

Last edited by Zkir (2010-09-06 22:42:47)

Ivan Komarov · 2010-09-07 01:00:29

gps-Max wrote:

Подозреваю, что все мыслимые географические алгоритмы давным-давно реализованы, проверены, отлажены и упакованы в библиотеки.

Да-да! Там (где-то в недрах Перла) даже алгоритм пересечения двух отрезков был кривой. Так что верны только первый и четвертый постулат, но не второй и третий

x10kHz · 2010-09-07 08:29:16

мда, как оказалось задачка не такая простая....
по большей части нужно как раз для распределения объектов по масштабным уровням.
Но хотелось бы и просто иметь такой алгоритм под рукой с необходимой точностью)

gps-Max, не обязательно на питоне, пойдет и на С++, и на яве... просто может где описание какое есть нормальное. Сам ничего токового нагуглить не смог(
Просто импорт, экспорт я замутил питоновскими скриптами... вдоволь намучавшись с SQL базами я переполз на NoSQL MongoDB. Пытаюсь сделать уровни детализации в векторе.

Zkir, как именно пересчитываются квадратные градусы в метры квадратные? у них же разная длина в метрах.... если от экватора подниматься длина градуса одного угла все время уменьшается.
А если еще и несферичность Земли участь...

Линки что нашел... пишу скорее для себя, чтобы потом не забыть, но возможно кому-нибудь пригодится:
The Geometry of the Sphere - полезные статейки по теме
areaint(lat,lon,ellipsoid,units) Функция Mapping тулбокса матлаба... там вроде в *.m файликах исходнички, если тело функции не в какой-нибудь dll'ке
Measurement of Areas on a Sphere Using Fibonacci and Latitude–Longitude Lattices - какая-то умная пдфка...
A method to compute the area of a spherical polygon - непроверенные кусочки кода)
Дока к либе GeoSphere под GPL есть функция areaPolygon... Compute the area of a polygon on a sphere. Код на основе статейки Computing the area of a spherical polygon of arbitrary shape, которую я достать не смог...
Spherical Polygon что-то полезное от Wolfram)
Equal-area Earth map utility in Java - что-то связанное с проекциями что ли... на которых площадь правильная... по ней бы можно было и померить

liosha · 2010-09-07 08:37:41

Ivan Komarov wrote:

Да-да! Там (где-то в недрах Перла) даже алгоритм пересечения двух отрезков был кривой.

Не в недрах перла, а в одной кривой библиотеке.

В osm2mp площадь объекта считается в обычной географической проекции, и умножается на коэффициент - косинус широты. Для первого приближения этого вполне достаточно.

Для несферической земли можно так:
S = lat / 360 * <длина экватора> x lon*cos(lat) / 360 * <длина меридиана>

Hind · 2010-09-07 08:39:26

Квадратные градусы? Хмм... Может, он стерадианы имел в виду? :3

liosha · 2010-09-07 08:59:12

Из гугла: http://ru-patent.info/21/65-69/2166731.html

chnav · 2010-09-07 09:23:47

Я бы проверил следующий алгоритм

/* Сразу привести широту и долготу точек полигона в метры, т.е. в очень приблизительную прямоуголку. Можно даже не заморачиваться с центральным меридианом, т.к. нас не волнуют абсолютные значения. */
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++

ПРИМЕЧАНИЯ:
1. площадь будет иметь знак;
2. не забудьте замкнуть полигон, в противном случае значение будет неверное;
3. для больших площадей надо считать площадь трапеции на сфере.

ИМХО, достоинство этого метода:
+ считается площадь всего полигона сразу;
+ получившаяся площадь имеет положительное значение, если полигон по часовой стрелке, и отрицательная, если против часовой;
Недостатки:
- площадь самопересекающегося полигона может принять неверное значение, включая 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++

Last edited by chnav (2010-09-07 10:26:31)

x10kHz · 2010-09-07 11:43:18

Пока потестил вот этот исходничок... почти такой же как и этот
Обе реализации показывают довольно внятные числа.... в смысле если и ожидаю результат в сотни квадратных метров, то результат хотя бы имеет тот же порядок)

Выбрал на карте прямоугольник, взял его углы и пропустил через ту функцию... выдает 10270 метров квадратных.
Участок был со сторонами 110м х 207м... что по площади в два раза больше! В принципе уже радует, что хоть такая оценка есть!

Насколько я понимаю с невыпуклыми полигонами алгоритмы могут и не справится... наверно триангуляция нужна)))
сделаю реализацию на питоне - выложу)

chnav · 2010-09-07 11:51:39

Зря значит я тут интегрированием занимался Попробуйте второй метод из моего сообщения, должен работать на любых полигонах без самопересечения, полностью находящихся в одном полушарии (либо северном, либо южном). Точность будет доли процентов, а не в два раза...

Komяpa · 2010-09-07 11:56:58

x10kHz, ошибка в два раза - скорее всего, где-то потеряно домножение на косинус широты. Или забыт перевод из градусов в радианы.

x10kHz · 2010-09-07 12:05:57

chnav wrote:

Зря значит я тут интегрированием занимался Попробуйте второй метод из моего сообщения, должен работать на любых полигонах без самопересечения, полностью находящихся в одном полушарии (либо северном, либо южном). Точность будет доли процентов, а не в два раза...

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, в радианы переводил все прально!
на счет косинусов и прочей математики это я уже не соображу так просто... я готовые исходнички проверил особо не разбираясь)

Last edited by x10kHz (2010-09-07 12:09:00)

Komяpa · 2010-09-07 12:21:29

x10kHz, рекомендую:
1) хранить всё в парах туплов (lon, lat).
2) перед подсчётом площади всё перевести в меркаторовскую проекцию (EPSG:3857).
Функцией, например, from4326() отсюда: http://code.google.com/p/twms/source/br … ections.py
3) пересчитать теми же треугольниками площадь, взять по модулю
4) помножить на косинус средней широты. (математики, поправьте: с длинами-ширинами прокатывает, не становится ли косинус косинусом-квадратом при уходе в площадь?)

liosha · 2010-09-07 12:23:50

Komяpa wrote:

2) перед подсчётом площади всё перевести в меркаторовскую проекцию (EPSG:3857)

8-O
Это ещё зачем??

Komяpa · 2010-09-07 12:29:06

liosha, чтобы вылезли метры. Тогда можно не заморачиваться с радиусами планеты и тому подобным

liosha · 2010-09-07 12:31:00

Komяpa, в метры лучше как угодно переводить, только не через меркатор.

chnav · 2010-09-07 12:31:23

x10kHz
Загрузил ваш пример в GPSMapEdit, он дает площадь 10605 кв.м, т.е. мой алгоритм кривоват 0.06%. ИМХО это много для такого маленького объекта. Скорее всего сыграло роль что мы считаем на сфере очень примерного радиуса, а не на эллипсоиде как мапэдит.

Last edited by chnav (2010-09-07 12:33:02)

liosha · 2010-09-07 12:34:02

Если нужна хорошая точность, надо взять в билиотеке какой-нть учебник по геодезическим вычислениям

liosha · 2010-09-07 12:36:33

Ещё раз, мой способ:

1) Преобразуем координаты
mlat = lat * (длина экватора) /360
mlon = lon * cos(lat) * (длина меридиана) /360

2) Используем любую функцию из любой библиотеки для площади многоугольника

x10kHz · 2010-09-07 21:22:13

chnav wrote:

x10kHzЗагрузил ваш пример в GPSMapEdit, он дает площадь 10605 кв.м, т.е. мой алгоритм кривоват 0.06%.

Линейка Яндекса показывает что стороны именно такие и площадь должна быть большей... но координаты я с гугла брал, проверял ставил точку... все совпадает. Мне кажется GPSMapEdit тоже наврал немного... Поставил его демку, но что-то не подружились мы)

liosha, вот реализация Вашей идеи...

import math

from shapely.geometry import Polygon

Lekv = 40075700
Lmer = 19980000

lat = [55.728634, 55.729489, 55.729440, 55.728565, 55.728634]
lon = [37.826026, 37.826097, 37.827862, 37.827780, 37.826026]

nds = [];
for i in range(len(lat)):
    mlat = lat[i]*Lekv/360.0
    mlon = lon[i]*math.cos(math.radians(lat[i]))*Lmer/360.0
    nds.append((mlat, mlon))

polygon = Polygon(nds)
print polygon.area

Результат 5310.72967046

Я уже сомневаюсь в правильности исходных координат, в реальной площади... есть ли какая-нибудь продвинутая считалка у кого-нибудь, может САПР какой?
Получить бы точные исходные данные и площадь для них заведомо правильно посчитанную...

Zkir · 2010-09-07 22:08:09

У вас длина меридиана вполовину меньше чем надо.

x10kHz · 2010-09-07 23:50:28

Zkir, спасибо!
Длину меридиана брал из статьи "Эллипсоид Красовского"... гугл выдал

Поискал еще и нашел значение 40008.55 км

Lmer = 40008550

Теперь результат 10634.3640418
Это очень похоже на все остальные и особенно на результаты что получил chnav в GPSMapEdit
качаю Civil от автодеска... вроде бы такая продвинутая прога для картографии и всяких там инженерных расчетов... посмотрю в ней можно ли посчитать и что покажет

Announcement

#1 2010-09-06 19:58:57

Географические алгоритмы

#2 2010-09-06 22:29:17

Re: Географические алгоритмы

#3 2010-09-06 22:41:03

Re: Географические алгоритмы

#4 2010-09-07 01:00:29

Re: Географические алгоритмы

#5 2010-09-07 08:29:16

Re: Географические алгоритмы

#6 2010-09-07 08:37:41

Re: Географические алгоритмы

#7 2010-09-07 08:39:26

Re: Географические алгоритмы

#8 2010-09-07 08:59:12

Re: Географические алгоритмы

#9 2010-09-07 09:23:47

Re: Географические алгоритмы

#10 2010-09-07 11:43:18

Re: Географические алгоритмы

#11 2010-09-07 11:51:39

Re: Географические алгоритмы

#12 2010-09-07 11:56:58

Re: Географические алгоритмы

#13 2010-09-07 12:05:57

Re: Географические алгоритмы

#14 2010-09-07 12:21:29

Re: Географические алгоритмы

#15 2010-09-07 12:23:50

Re: Географические алгоритмы

#16 2010-09-07 12:29:06

Re: Географические алгоритмы

#17 2010-09-07 12:31:00

Re: Географические алгоритмы

#18 2010-09-07 12:31:23

Re: Географические алгоритмы

#19 2010-09-07 12:34:02

Re: Географические алгоритмы

#20 2010-09-07 12:36:33

Re: Географические алгоритмы

#21 2010-09-07 21:22:13

Re: Географические алгоритмы

#22 2010-09-07 22:08:09

Re: Географические алгоритмы

#23 2010-09-07 23:50:28

Re: Географические алгоритмы

Board footer