Обсудить в форуме Комментариев 19Редактировать в вики
Для тех, кто хочет разобраться в математике пересчета и\или встроить конвертацию в свою программу.
В этой статье приводятся функции для преобразования геодезических координат из координатной системы Пулково 1942 в координатную систему WGS 1984 и обратно.
В данном примере используются три параметра трансформирования (dx, dy, dz), остальные параметры равны 0. Если вы хотите использовать 7 параметров (дополнительно wx, wy, wz, ms), измените значения угловых элементов трансформирования. Линейные элементы трансформирования также могут быть другими, в данной примере используются линейные элементы трансформирования из ГОСТ 51794-2001 (другие возможные наборы элементов) . Для рассчета используется формулы Бурса-Вольфа (подробнее).
Все угловые значения передаются и возвращаются в десятичных градусах (dd.ddddd), высоты передаются и возвращаются в метрах.
Результаты пересчета проверены с помощью Proj и совпадают с его результатами до 7-го знака после запятой.
Const Pi As Double = 3.14159265358979 ' Число Пи Const ro As Double = 206264.8062 ' Число угловых секунд в радиане ' Эллипсоид Красовского Const aP As Double = 6378245 ' Большая полуось Const alP As Double = 1 / 298.3 ' Сжатие Const e2P As Double = 2 * alP - alP ^ 2 ' Квадрат эксцентриситета ' Эллипсоид WGS84 (GRS80, эти два эллипсоида сходны по большинству параметров) Const aW As Double = 6378137 ' Большая полуось Const alW As Double = 1 / 298.257223563 ' Сжатие Const e2W As Double = 2 * alW - alW ^ 2 ' Квадрат эксцентриситета ' Вспомогательные значения для преобразования эллипсоидов 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 = 23.92 Const dy As Double = -141.27 Const dz As Double = -80.9 ' Угловые элементы трансформирования, в секундах Const wx As Double = 0 Const wy As Double = 0 Const wz As Double = 0 ' Дифференциальное различие масштабов Const ms As Double = 0 Function WGS84_SK42_Lat(Bd, Ld, H) As Double WGS84_SK42_Lat = Bd - dB(Bd, Ld, H) / 3600 End Function Function SK42_WGS84_Lat(Bd, Ld, H) As Double SK42_WGS84_Lat = Bd + dB(Bd, Ld, H) / 3600 End Function Function WGS84_SK42_Long(Bd, Ld, H) As Double WGS84_SK42_Long = Ld - dL(Bd, Ld, H) / 3600 End Function Function SK42_WGS84_Long(Bd, Ld, H) As Double SK42_WGS84_Long = Ld + dL(Bd, Ld, H) / 3600 End Function Function dB(Bd, Ld, H) As Double Dim B, L, M, N As Double B = Bd * Pi / 180 L = Ld * Pi / 180 M = a * (1 - e2) / (1 - e2 * Sin(B) ^ 2) ^ 1.5 N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5 dB = ro / (M + H) * (N / a * e2 * Sin(B) * Cos(B) * da _ + (N ^ 2 / a ^ 2 + 1) * N * Sin(B) * Cos(B) * de2 / 2 _ - (dx * Cos(L) + dy * Sin(L)) * Sin(B) + dz * Cos(B)) _ - wx * Sin(L) * (1 + e2 * Cos(2 * B)) _ + wy * Cos(L) * (1 + e2 * Cos(2 * B)) _ - ro * ms * e2 * Sin(B) * Cos(B) End Function Function dL(Bd, Ld, H) As Double Dim B, L, N As Double B = Bd * Pi / 180 L = Ld * Pi / 180 N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5 dL = ro / ((N + H) * Cos(B)) * (-dx * Sin(L) + dy * Cos(L)) _ + Tan(B) * (1 - e2) * (wx * Cos(L) + wy * Sin(L)) - wz End Function Function WGS84Alt(Bd, Ld, H) As Double Dim B, L, N, dH As Double B = Bd * Pi / 180 L = Ld * Pi / 180 N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5 dH = -a / N * da + N * Sin(B) ^ 2 * de2 / 2 _ + (dx * Cos(L) + dy * Sin(L)) * Cos(B) + dz * Sin(B) _ - N * e2 * Sin(B) * Cos(B) * (wx / ro * Sin(L) - wy / ro * Cos(L)) _ + (a ^ 2 / N + H) * ms WGS84Alt = H + dH End Function
Код пересчетов базируется на примере Olexa Riznyk, www.olexa.com.ua
Приведенный ниже код нужно скопировать и открыть окно редактора Visual Basic программы Microsoft Excel (Сервис\Макрос\Редактор Visual Basic). Далее, необходимо создать новый модуль, для этого нужно выбрать в меню Insert\Module и вставить в появившийся модуль скопированный код.
Так же можно загрузить этот код в виде готового модуля и импортировать его (File\Import File...) в окне редактора Visual Basic.
После установки, в рабочей области Excel станут доступны следующие функции:
Пересчет широты из WGS-84 в СК-42: WGS84_SK42_Lat(Lat,Long,Height)
Пересчет долготы из WGS-84 в СК-42: WGS84_SK42_Long(Lat,Long,Height)
Пересчет широты из СК-42 в WGS-84: SK42_WGS84_Lat(Lat,Long,Height)
Пересчет долготы из СК-42 в WGS-84: SK42_WGS84_Long(Lat,Long,Height)
Для выполнения пересчета, нужно использовать вышепреведенные функции подставляя им в качестве аргументов значения широты и долготы в десятичных градусах и высоты в метрах, например, в ячейку Excel можно ввести:
=WGS84_SK42_Lat(50;50;0)
и получить результат 49.99980414
Если ввод формулы приводит к ошибке, убедитесь, что в документе разрешены Макросы (Сервис\Макрос\Безопасность - Средняя)
Для вычислений также можно использовать готовую таблицу в формате MS Excel (скачать).
Обсудить в форуме Комментариев 19Редактировать в вики
Последнее обновление: 2014-05-14 23:39
Дата создания: 07.12.2004
Автор(ы): Максим Дубинин
© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).