Тип плагинов: Параметры полигона

[vdemidov]

Раз уж на форуме зашла речь о вычислении площади полигона, а мне ее точное вычисление делать некогда, созрела мысль о таком типе плагинов.

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

Код: Выделить всё

  ILonLatPolygonGetSimpleInfo = interface
    function GetPolygonInfo(APointsCount: Cardinal; APoints: PDoublePoint): WideString;
  end;


APointsCount - количество вершин полигона
APoints - указатель на массив с географическими координатами вершин полигона
Возвращать можно просто текст, или html

[DJ VK]

Пока только в голове крутится решение вида разбить полигон на треугольники и вычислить по формуле Жирара площадь каждого. Насколько на локальном участке возможно считать Землю за сферу?
[url]http://ru.wikipedia.org/wiki/Сферическая_геометрия[/url]

[PavelML]

Пока только в голове крутится решение вида разбить полигон на треугольники и вычислить по формуле Жирара площадь каждого. Насколько на локальном участке возможно считать Землю за сферу?
[url]http://ru.wikipedia.org/wiki/Сферическая_геометрия[/url]

Извините, что на васике, но полагаю формула понятна...

Код: Выделить всё

Public Function AreaPol(ArrXY As Variant) As Double
'функция вычисления площади контура по двумерному массиву координат (координата X в столбце 0, Y - в столбце 1)
Dim S1 As Double
Dim S2 As Double
Dim MM As Long
Dim I As Long
MM = UBound(ArrXY, 2)
If MM > 1 Then
S1 = 0
S2 = 0
S1 = S1 + Int(0.5 + (100 * ArrXY(1, 0))) * (Int(0.5 + (100 * ArrXY(0, MM))) - Int(0.5 + (100 * ArrXY(0, 1))))
S2 = S2 + Int(0.5 + (100 * ArrXY(0, 0))) * (Int(0.5 + (100 * ArrXY(1, 1))) - Int(0.5 + (100 * ArrXY(1, MM))))
For I = 1 To MM - 1
    S1 = S1 + Int(0.5 + (100 * ArrXY(1, I))) * (Int(0.5 + (100 * ArrXY(0, I - 1))) - Int(0.5 + (100 * ArrXY(0, I + 1))))
    S2 = S2 + Int(0.5 + (100 * ArrXY(0, I))) * (Int(0.5 + (100 * ArrXY(1, I + 1))) - Int(0.5 + (100 * ArrXY(1, I - 1))))
Next I
S1 = S1 + Int(0.5 + (100 * ArrXY(1, MM))) * (Int(0.5 + (100 * ArrXY(0, MM - 1))) - Int(0.5 + (100 * ArrXY(0, 0))))
S2 = S2 + Int(0.5 + (100 * ArrXY(0, MM))) * (Int(0.5 + (100 * ArrXY(1, 0))) - Int(0.5 + (100 * ArrXY(1, MM - 1))))
AreaPol = Abs((S1 + S2) / 40000)
Else
AreaPol = 0
End If
End Function

В данном коде много лишнего, он используется для вычисления прощади в кадастровой программе. Используются округления координат до второго знака (до сантиметров). Для простоты вычислений все вычисляется в квадратных сантиметрах а потом делится на 10000 для получения квадратных метров. Точнее - делится на 40000, поскольку S1 и S2 - каждая удвоенные значения площадей. Там вычисляется также модуль площади - она может получиться либо положительной либо отрицательной в зависимости от направления обхода границы контура (по часовой или против часовой стрелки). Кстати - знак иногда полезен именно для определения направления нумерации точек.
Что касается массива координат X и Y - потребуется преобразовать широту и долготу в метрические прямоугольные координаты. Поскольку абсолютная точность позиционирования не важна - можно использовать приблизительные формулы с условной точкой отсчета. По широте каждый градус достаточно четко выражен в метрах, а по долготе он зависит от широты и вычисляется примерно по той же формуле, по которой вычисляется местоположение квадратиков космических снимков. Либо используем сложные формулы из ГОСТ Р 51794-2001.
Само собой будет вычисляться площадь горизонтальной проекции. Участок на склоне горы будет по поверхности намного больше, чем по горизонтальной проекции.

[i] слегка подправил тэги

[vdemidov]

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

[PavelML]

Ну вообще-то я здесь не спрашивал форму вычисления площади полигона на плоскости. Ее я слава богу еще со школьного курса знаю.

Господин Демидов, я совершенно неспособен узнать дистанционно - что Вы изволите знать со "школьного курса" - а что только предполагаете узнать...
Было спрошено "зашла речь о вычислении площади полигона, а мне ее точное вычисление делать некогда"
и было отвечено формулой вычисления площади полигона.
Если бы было спрошено преобразование географических координат в координаты прямоугольной проекции - было бы отвечено формулой вычисления координат прямоугольной проекции. Вообще конечно рулит поиск, для чего я указал - что именно нужно искать: "ГОСТ Р 51794-2001".
Там все эти формулы указаны с условной точностью вычисления в миллиметры в пределах полигонов стокилометровой размерности. Конечно - речь идет о проекции гаусса-крюгера. Потому что площади более-менее крупных полигонов в проекции меркатора (испольтзуемой в навигации) - вычислять бессмысленно, как бессмысленно приравнивать поверхность цилиндра к поверхности шара того же диаметра. А Вы кажется точность хотели?

для гаусса крюгера без привязки к местной системе координат (то есть совпадения с кадастровой площадью НЕ будет):

Const Pi As Double = 3.14159265358979 ' Число Пи
Const ro As Double = 206264.8062 ' Число угловых секунд в радиане
' Эллипсоид Красовского (Пулково 1942)
Const aP As Double = 6378245 ' Большая полуось
Const alP As Double = 1 / 298.3 ' Сжатие
Const e2P As Double = 2 * alP - alP ^ 2 ' Квадрат эксцентриситета
' Эллипсоид GRS80 (WGS84)
Const aW As Double = 6378137 ' Большая полуось
Const alW As Double = 1 / 298.257223563 ' Сжатие
Const e2W As Double = 2 * alW - alW ^ 2 ' Квадрат эксцентриситета
Const e As Double = e2W ^ 0.5 ' Эксцентриситет
Const e4 As Double = e ^ 4 ' Эксцентриситет в 4 степени
Const e6 As Double = e ^ 6 ' Эксцентриситет в 6 степени
' Вспомогательные значения для преобразования эллипсоидов
Const a As Double = (aP + aW) / 2
Const e2 As Double = (e2P + e2W) / 2
Const da As Double = aW - aP
Const de2 As Double = e2W - e2P
' Линейные элементы трансформирования, в метрах
Const dx As Double = 25
Const dy As Double = -141
Const dz As Double = -80
' Угловые элементы трансформирования, в секундах
Const wx As Double = 0
Const wy As Double = 0
Const wz As Double = 0
' Дифференциальное различие масштабов
Const ms As Double = 0

Function G2X(LatP, LonP, CMer) As Double
' функция вычисления координаты X по широте-долготе WGS-84
' входные значения - в десятичных градусах широты LatP, долготы LonP и
' центрального меридиана CMer (принимаем для вычисления площади полигона долготу первой точки полигона)
Dim LatR, LonR, Lat0, C, nu, Aa, T, CMerR, M, M0 As Double
LatR = LatP * (Pi / 180) 'преобразуем градусы в радианы
LonR = LonP * (Pi / 180) '-"
Lat0 = 0 * (Pi / 180) '
CMerR = CMer * (Pi / 180) '
C = e2W * Cos(LatR) ^ 2 / (1 - e2W)
nu = aW / Sqr(1 - e2W * Sin(LatR) ^ 2)
Aa = (LonR - CMerR) * Cos(LatR)
T = Tan(LatR) ^ 2
M = aW * ((1 - e2W / 4 - 3 * e4 / 64 - 5 * e6 / 256) * LatR - (3 * e2W / 8 + 3 * e4 / 32 + 45 * e6 / 1024) * Sin(2 * LatR) + _
(15 * e4 / 256 + 45 * e6 / 1024) * Sin(4 * LatR) - (35 * e6 / 3072) * Sin(6 * LatR))
M0 = aW * ((1 - e2W / 4 - 3 * e4 / 64 - 5 * e6 / 256) * Lat0 - (3 * e2W / 8 + 3 * e4 / 32 + 45 * e6 / 1024) * Sin(2 * Lat0) + _
(15 * e4 / 256 + 45 * e6 / 1024) * Sin(4 * Lat0) - (35 * e6 / 3072) * Sin(6 * Lat0))
G2X = M - M0 + nu * Tan(LatR) * Aa ^ 2 / 2 + (5 - T + 9 * C + 4 * C ^ 2) * Aa ^ 4 / 24 + (61 - 58 * T + T ^ 2 + 600 * C - 330 * e2W) * Aa ^ 6 / 720
End Function

Function G2Y(LatP, LonP, CMer) As Double
' функция вычисления координаты Y по широте-долготе WGS-84
' входные значения - в десятичных градусах широты LatP, долготы LonP и
' центрального меридиана CMer (принимаем для вычисления площади полигона долготу первой точки полигона)
Dim LatR, LonR, Lat0, CMerR, C, nu, Aa, T As Double
LatR = LatP * (Pi / 180)
LonR = LonP * (Pi / 180)
Lat0 = 0 * (Pi / 180)
CMerR = CMer * (Pi / 180)
C = e2W * Cos(LatR) ^ 2 / (1 - e2W)
nu = aW / Sqr(1 - e2W * Sin(LatR) ^ 2)
Aa = (LonR - CMerR) * Cos(LatR)
T = Tan(LatR) ^ 2
G2Y = nu * (Aa + (1 - T + C) * Aa ^ 3 / 6 + (5 - 18 * T + T ^ 2 + 72 * C - 58 * e2W) * Aa ^ 5 / 120)
End Function

[PavelML]

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

[bk99]

Не знаю, будет ли полезно - плагин "Расчет площади произвольных многоугольных фигур":
http://plugins.2gis.ru/download/

[vdemidov]

Не знаю, будет ли полезно - плагин "Расчет площади произвольных многоугольных фигур":
http://plugins.2gis.ru/download/

Мне точно неинтересно, меня устраивает текущая точность. Но если, кто-то возьмется реализовывать, то может и сгодится.

[Guam]

Сейчас программа вычисляет площадь совсем не верно, и это не погрешность из-за сферичности земли. Вот, к примеру, последняя версия программы площадь квадратной площадки ок. 1.07х1.07 км она вычисляет как 0.69 км2, но такая ошибка по-моему была не всегда. Реальная же ее площадь с учетом сферичности земли 1139447,89 м2, а без учета 1139447,88 м2. На малых площадях расчет площади с учетом сферичности или без учета дает почти одинаковый результат. Расхождения заметны на площадях порядка административных областей. Я могу выложить процедуру по вычислению площади с учетом сферичности, но как оформлять плагин я не в курсе.

[vdemidov]

Оформите процедуру как объект с инетрфейсом описанным в первом посте этой темы. Если еще не догадались, то язык Delphi.