Площадь по геокоординатам

[Merkator]

Я конечно не студент, но нужна помощь.
Нужно найти площадь многоугольника, заданного вершинами (longitude, latitude)
(Англоязычный гугл исхожен вдоль и поперек, в русскоязычном ничего адекватного не нашел)
Мое решение на питоне pastebin.

Сразу поясню:
В начале мы имеем некий набор координат в формате degrees, minutes
Я перевожу их прям там же в стандартный формат (градусы с десятичной частью).

Дальше по формулам из wiki перевожу их в формат ECEF (earth-centred, earth-fixed).

А вот в дальнейшем коде я не уверен. Я считаю площадь зоны, ограниченной этими точками, методом векторного произведения.
Я сомневаюсь, что это верно.

Нужна Ваша помощь в подсчете площади.

PS для любопытных: координаты в программе - зона в Финском заливе.

[Merkator]

Ладно... видимо задача несколько отличного от форумского профиля.

[l1pton17]

терпеливее будьте или вы думаете, что после размещения вашей проблемы, сразу же весь форум полезет её решать?

[Merkator]

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

Но раз опубликовал, то постараюсь довести решение до логического конца:

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

Ищу теперь иной метод для подсчета.

[Vago]

Сообщение от Merkator

Я понял в чем ошибка в моем коде. Тем методом, которым я считаю площадь, я получаю не площадь нужной мне фигуры, а площадь боковых граней "пирамиды" с центром в (0, 0) и основанием - тем самым полигоном.

Там, у Вас, на самом деле, по-моему, тоже немного не так, но для начала скажите, пожалуйста: а Вам от кординаты Z отказаться - очень тяжело? Ну вот у Вас в примере: 5502.7, 5507.3, 5493.1, 5489.1. Относительная величина "плавания" Z - 3 десятых процента. Может, ну её, эту выпуклость Земли?..

[Merkator]

Vago, мне кажется, что Z - не просто выпуклость Земли. Насколько я понял систему ECEF, то просто из центра Земли определенным образом отложена декартова система координат. И если я откажусь от Z, то скажем не в Финском заливе, а где-нибудь в Атлантике программа перестанет работать. (координата Z будет сильно меняться, а какая-нибудь Y наоборот не будет)

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

Тогда появляется новый вопрос: как правильно нормализовать координаты?

[Vago]

Merkator, я - не специалист по ГИС и, к сожалению, ничего не знаю о картографических проекциях и нормализации. Меня пока смущает другое. У Вас там ответственность за то, чтобы вершины были введены "по порядку", на программисте лежит, да?

[alexBlack]

По-моему нужно просто переместить все точки до совмещения, например, ndata[0] с центром координат. Тогда Ваш цикл будет суммировать площади треугольников:

area = 0;
for i := 1 to len(ndata)-2   // не знаю как это в Питоне
	area += 0.5*sqrt(dot(cross(ndata[i], ndata[i+1]), cross(ndata[i], ndata[i+1])))

[Merkator]

Сообщение от Vago

Merkator, я - не специалист по ГИС и, к сожалению, ничего не знаю о картографических проекциях и нормализации. Меня пока смущает другое. У Вас там ответственность за то, чтобы вершины были введены "по порядку", на программисте лежит, да?

Нет, они уже даны в порядке обхода. А что?

[Merkator]

Сообщение от alexBlack

По-моему нужно просто переместить все точки до совмещения, например, ndata[0] с центром координат. Тогда Ваш цикл будет суммировать площади треугольников:

area = 0;
for i := 1 to len(ndata)-2   // не знаю как это в Питоне
	area += 0.5*sqrt(dot(cross(ndata[i], ndata[i+1]), cross(ndata[i], ndata[i+1])))

Возможно... сейчас попробую.
Спасибо!