June 15, 2018

ppGeoCoord

Пишу больше для себя и непонятно, лучше не читайте. Есть слабая вероятность, что кому-то еще пригодится. Или самому понадобится вспомнить.

Итак, задача написать модуль топопривязки изображения к системам координат. Что это? Есть карта, скачаная, например, из интернета, картинка примерно 3000 на 3000 пикселей. Как узнать координаты точки на карте? Для этого и нужна топопривязка. Пользователь щелкает мышкой по углам карты, в местах, где координаты известны (в узлах сетки), и вводит в модуль координаты. Из этих точек (достаточно три, но лучше больше - так мы узнаем RMS - среднеквадратичное отклонение, т.е. насколько точной получилась привязка) мы вычислим коэффициенты аффинного преобразования, решая СЛАУ методом Гаусса или Крамера (привет Цешковскому!). Крамер сложнее, но лучше. Получив нулевой главный детерминант мы узнаем, если матрица вырождена. Что делать, если сетка в одной системе координат, а на выходе нужна другая? Пишем функцию преобразования координат, в гугле есть все формулы, вроде всё просто.

Но приключения только начинаются.

На сайте Науково-дослідного інституту геодезії і картографії есть онлайн-конвертер, его используем для проверки правильности вычислений:
http://dgm.gki.com.ua/ua/transkord
Координаты он принимает только по территории Украины, но для целей проверки этого достаточно.

Мне пришлось работать с картами генштаба, там советская система координат СК-42, и сетка в прямоугольной системе координат Гаусса-Крюгера. Весь мир работает в WGS-84, значит, наша задача привязать карту в СК-42, а координаты получить в WGS-84. Работаем.

Гуглом находим формулы:
http://gis-lab.info/qa/wgs84-sk42-wgs84-formula.html
На форуме есть даже даже перевод из VBA на C:
https://gis-lab.info/forum/viewtopic.php?t=7281
Работает только прямое преобразование, обратное глючит, правим формулы, сверяясь с российским ГОСТом 51794-2001:
http://gis-lab.info/docs/legislation.html

Гаусс-Крюгер в СК-42 и обратно работает, что подтверждается конвертером института геодезии. Но простая часть (СК-42 в WGS-84 и обратно) даёт ошибку на 4-м знаке после запятой, что может давать сотни метров погрешности.

После муторных проверок понимаем, что перевод из СК-42 в WGS-84 на сайте института геодезии работает не по ГОСТу.

Намёк на возможную причину даёт следующий документ:
http://gis-lab.info/qa/datum-transform-sets.html

У каждой страны имеется своя местная система координат, Паганель доступно об этом пишет:
http://www.hllab.dp.ua/Store/texts/gps/datums.htm

Набор параметров смещения (датумы) для разных стран можно посмотреть здесь:
https://www.colorado.edu/geography/gcraft/notes/datum/edlist.html

Украины здесь нет, на сайте института написано, почему:
"...інформацією з обмеженим доступом. Зокрема це стосується, наприклад, параметрів переходу від просторових прямокутних координат X, Y, Z у системах координат ITRS, WGS-84, ETRS89 до просторових прямокутних координат X, Y, Z..."
http://dgm.gki.com.ua/ua/faq

Стало понятно, почему в подробной документации по онлайн-конвертеру отсутствует датум для WGS-84:
http://dgm.gki.com.ua/files/calc/Transkor_help.pdf

Из списка датумов разных стран (см ссылку выше) ближе всего подошла Польша, она дала минимальную ошибку:
dx=23. dy=-124, dz=-82.

Путём подбора смещений удалось загнать отличия от официального конвертера за 5-тый знак после запятой, что уже в пределах угловой секунды и практическая ошибка не должна превышать 20 метров:
dx = 22.8788; dy = -122.977; dz = -80.113;

Паганель пишет про голландскую утилиту invmol.exe - обратное преобразование Молодецкого, с помощью него можно найти неизвестный датум по известным точкам, но даже с ней точность мне улучшить не удалось. Если кто захочет попробовать, утилита переехала на новый адрес:
https://www.itc.nl/ilwis/download/geodetic-tools/

Итого, написал класс ppGeoCoord под c++ для преобразования Гаусс-Крюгер - СК42 - WGS84 и обратно, для скачки доступен здесь:
https://github.com/podoroges/pplib

Всем спасибо за внимание.

No comments: