Добавление местной координатной системы в GIS

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/local-cs.html

Конструирование проекций, имитирующих местные координатные системы, в QGIS

Введение

Под местной системой координат (МСК) будет подразумеваться так называемая «городская» система, построенная независимо от государственной системы (ГСК) и включенная в неё заданием ключей перехода к СК-42 или СК-63. МСК крупных территорий, сравнимых с размерами субъектов Федерации, не являются предметом данной статьи, поскольку относятся к классическим картографическим проекциям.

Постановка задачи

Имеется множество пунктов, для которых известны координаты X, Y в ГСК и x, y в МСК. Требуется подобрать проекцию и вычислить её параметры, удовлетворительно представляющую МСК в ГИС. При подборе параметров предполагается использовать один из пунктов в качестве центральной точки преобразования.

Некоторые программы позволяют реализовать работу в МСК непосредственно. Так, в MapInfo любая проекция может быть дополнена аффинным преобразованием, а в Global Mapper конформным. ArcGIS в качестве МСК предлагает локальную проекцию: аналог ортографической проекции на эллипсоиде, дополненный разворотом и масштабированием.

Другие программы, включая QGIS, работают только с «обыкновенными» проекциями, не допуская дополнительных геометрических преобразований. Задача статьи — показать, как можно конструировать проекции, позволяющие имитировать МСК в QGIS.

В качестве рабочей среды будем использовать командную строку UNIX или OSGeo4W.

Подготовка данных

Нужны два файла исходных данных. Пусть текстовый файл cat_s42z4.tsv содержит координаты пунктов в государственной системе (ГСК), а именно в четвёртой зоне СК-42:

4645997.49 5768521.60
4661392.15 5770068.91
4650059.09 5783332.41
4634241.37 5778952.22
4631481.69 5764570.61
4642125.18 5754643.12
4655952.19 5757337.28

В файле cat_local.tsv — координаты в местной системе (МСК):

67266.64 30088.40
82697.29 31184.27
71759.40 44771.50
55822.67 40857.06
52643.65 26564.42
62990.64 16331.35
76888.20 18619.57

Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы.

Вычислим географические координаты пунктов:

$ proj -I -f "%.16g" +init=epsg:28404 cat_s42z4.tsv > cat_longlat.tsv

В файл cat_longlat.tsv запишутся географические долготы и широты в СК-42:

23.12716113273887 52.02642240080064
23.35199369039325 52.03605540518488
23.19280603770445 52.15834802315441
22.96008928353681 52.12307568371703
22.91428747596251 51.99456156589881
23.06504376727423 51.90277894131356
23.26700321114454 51.92327827668979

Быстрая реализация

Работа заключается в задании трёх параметров косой проекции Меркатора и вычислении оставшихся четырёх. Априорные параметры определяются так: центр проекции должен лежать на осевом меридиане исходной проекции СК-42/СК-63 на широте середины города, а азимут начальной линии должен быть малым. Пусть широта центра проекции равна широте первого пункта из файла cat_longlat.tsv 52.02642240080064°. Долготу центра проекции приравняем долготе осевого меридиана четвёртой зоны 21°. Азимут начальной линии приравняем произовольно малой величине −0,0001°. Масштабный множитель k приравняем единице, оставшиеся параметры обнулим, вычислим координаты в этой переходной проекции и запишем их в файл cat.tsv:

$ proj -f "%.16g" +proj=omerc +lat_0=52.02642240080064 +lonc=21 +alpha=-0.0001 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=krass cat_longlat.tsv > cat.tsv

Найдём оставшиеся параметры с помощью программы helmkey:

$ helmkey cat.tsv cat_local.tsv

Программа выведет на экран искомую строку параметров:

+k=0.9998372145697554 +x_0=-78707.08698765696 +y_0=32225.08703617829 +gamma=1.677027577390829

Это всё! Осталось добавить в QGIS пользовательскую проекцию, не забыв необходимые стандартные параметры:

+proj=omerc +lat_0=52.02642240080064 +lonc=21 +alpha=-0.0001 +k=0.9998372145697554 +x_0=-78707.08698765696 +y_0=32225.08703617829 +gamma=1.677027577390829 +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs

Много букв

Выбор исходной проекции

Важной особенностью МСК является то, что с точки зрения математической картографии она остаётся проекцией Гаусса-Крюгера и наследует её искажения.

Удобно, если осевой меридиан ГСК проходит через город при малом угле разворота, как в Пензе, например. В этом случае МСК строится как проекция Гаусса-Крюгера, отличающаяся от исходной ГСК малым сдвигом осевого меридиана для компенсации угла разворота. Жаль, что это случается редко.

В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.

Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов координат имеется надёжный ключ в виде строки MapInfo.

Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита helmkey. Заглянем в файл var.csv, созданный этой утилитой. Он содержит невязки конформного преобразования:

-0.006	0.007
0.182	0.046
-0.166	0.110
0.019	-0.185
0.148	0.100
-0.146	0.094
-0.031	-0.171

В нашем случае второй кандидат на исходную ГСК — СК-63 зона C0 с осевым меридианом 21°57′. Создадим промежуточную проекцию на её основе, вычислим координаты для пунктов в ней и получим остающиеся параметры утилитой helmkey:

$ proj -f "%.16g" +proj=omerc +lat_0=52.02642240080064 +lonc=21.95 +alpha=-0.0001 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=krass cat_longlat.tsv > cat.tsv
$ helmkey cat.tsv cat_local.tsv

Вывод программы:

+k=1.00001854602098 +x_0=-13532.31228709544 +y_0=30742.75844830525 +gamma=0.9279005081790233

Содержимое var.csv:

-0.000	-0.002
-0.001	0.002
-0.001	0.002
0.004	0.000
-0.002	0.001
0.002	-0.002
-0.002	-0.001

Сравнение невязок не в пользу СК-42. Исходная ГСК — СК-63 зона C0, и правильная проекция для QGIS такая:

+proj=omerc +lat_0=52.02642240080064 +lonc=21.95 +alpha=-0.0001 +k=1.00001854602098 +x_0=-13532.31228709544 +y_0=30742.75844830525 +gamma=0.9279005081790233 +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs

Подробнее о выборе априорных параметров

В плане распределения искажений проекция МСК должна имитировать ГСК. В идеале выбор параметров таков, что центр проекции лежит на точке осевого меридана ГСК, ближайшей к середине города, а направление начальной линии совпадает с направлением меридиана. Итак,

• долгота центра проекции равна долготе осевого меридиана ГСК;
• широта центра проекции требует вычисления;
• азимут начальной линии равен нулю.

Достатоно хорошее приближение для широты центра проекции даёт формула

$\operatorname {tg} \phi _{2}={\frac {\operatorname {tg} \phi _{1}}{\cos(\lambda _{2}-\lambda _{1})}}$

Переходная проекция

Создадим в качестве переходной проекцию Гаусса-Крюгера, осевой меридиан и начальная параллель которой пересекаются в первом пункте, а координаты в этой точке совпадают с координатами МСК этого пункта. Вычислим координаты в переходной проекции и запишем их в файл pt_tmerc:

$ proj -f "%.4f" +proj=tmerc +lat_0=54.896961118 +lon_0=37.183749701 +k=1 +x_0=1334.71 +y_0=285.94 +ellps=krass pt_longlat.dat > pt_tmerc.dat

Содержимое pt_tmerc.dat:

1334.7100 285.9400
556.2353 -5196.2593
4445.4005 1149.5698
-248.5461 2884.0431

Подготовим входной файл data1.dat для findkey. Нужно на место второй и третьей колонок в data0.dat записать данные pt_tmerc.dat.

$ paste data0.dat pt_tmerc.dat | awk '{print $1, $7, $8, $4, $5, $6}' > data1.dat

Команда paste построчно объединяет данные из шести колонок data0.dat и двух колонок pt_tmerc.dat в общую строку, команда awk выводит в data1.dat нужные колонки в нужном же порядке:

1 1334.7100 285.9400 1334.71 285.94 1.0
2 556.2353 -5196.2593 563.67 -5197.34 1.0
3 4445.4005 1149.5698 4444.27 1153.79 1.0
4 -248.5461 2884.0431 -252.07 2881.90 1.0

Нетрудно заметить, что координаты первого пункта в переходной проекции совпадают с таковыми в МСК. Для остальных пунктов это не так, что говорит о наличии масштабирования и разворота. Убедимся в этом:

$ findkey data1.dat key1.dat var1.dat

Содержимое key1.dat:

0.390
-1.813
1.0000052644
0.0013554914
1.0000061831
+0.07766348

Первые четыре параметра — a₀, b₀, a₁, b₁, оставшиеся два — масштабный множитель m и разворот θ в градусах.

Косая проекция Меркатора

Всё готово для конструирования косой проекции Меркатора. Зададим, как для переходной проекции, долготу и широту первого пункта в качестве начальных, его координаты в МСК в качестве «фальшивых». Кроме того, перенесём масштаб конформного преобразования в масштаб проекции, а разворот — в азимут линии минимального масштаба. Вычислим координаты в этой проекции и запишем в файл pt_omerc.dat:

$ proj -f "%.4f" +proj=omerc +lat_0=54.896961118 +lonc=37.183749701 +alpha=0.07766348 +gamma=0 +k=1.0000061831 +x_0=1334.71 +y_0=285.94 +ellps=krass pt_longlat.dat > pt_omerc.dat

Вычислим параметры конформного преобразования для pt_omerc.dat: m = 1.0000000000, θ = −0.00000022° = −0.0008″. Практически масштаб сведён к единице, разворот к нулю, а вычисленные из значений долготы/широты координаты совпадают с координатами, заданными в МСК. Таким образом, проекция, эквивалентная МСК, построена.

Следует заметить, что многие программы (например, MapInfo) используют косую проекцию Меркатора в версии Хотина, когда началом координат служит точка пересечения линии минимального масштаба с экватором. Азимут линии и плоские координаты в этой точке, естественно, отличаются от заданных нами, и понадобятся дополнительные вычисления для их определения.

Реализация косой проекции Меркатора в PROJ.4 имеет одну особенность: программы просто не работают со значениями параметра alpha 0 или 90°. Если угол разворота очень мал (скажем, несколько угловых секунд), следует использовать проекцию Гаусса-Крюгера.

Проекция Гаусса-Крюгера

Встречаются программы, набор проекций в которых ограничен. Проекция Гаусса-Крюгера — она же Transverse Mercator — найдётся в самом бедном наборе. В приёмниках GARMIN пользовательская система координат может быть задана только в этой проекции.

Долгота осевого меридиана

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

В нашем примере разворот оказался небольшим: θ = 0.07766348° = 0°04′39.59″. Значит, в данном случае на основе проекции Гаусса-Крюгера можно создать довольно точный аналог МСК.

Предварительное значение перемещения осевого меридиана можно оценить по формуле:

$\Delta L=-\operatorname {arc\,tg} {\frac {\operatorname {tg} \theta }{\sin B}}$

где ∆L — перенос осевого меридиана, θ — угол разворота, B — широта первого пункта. Получаем ∆L = −0.0949292655, L = 37.0888204345. Построим проекцию Гаусса-Крюгера на основе переходной, подставив L в долготу осевого меридиана и нули в «фальшивые» координаты:

$ proj -f "%.4f" +proj=tmerc +lat_0=54.8969611 +lon_0=37.0888204 +k=1 +x_0=0 +y_0=0 +ellps=krass pt_longlat.dat > pt_tmerc.dat

Последующее вычисление параметров конформного преобразования выдаёт сравнительно большой угол разворота. Однако уже после второй итерации с L = 37.0888079 получаем m = 1.0000057306, θ = +0.00000044.

Масштаб проекции

Перенесём m в масштаб проекции:

+proj=tmerc +lat_0=54.8969611 +lon_0=37.0888079 +k=1.00000573 +x_0=0 +y_0=0 +ellps=krass

Вновь вычислим параметры конформного преобразования. Теперь m = 1.0000000039, т.е. практически единица.

Координаты в начальной точке проекции

Наконец обратим внимание на первичные параметры a₀ = −4756.693 и a₁ = 281.807. Перенесём их значения в параметры проекции x_0 и y_0:

+proj=tmerc +lat_0=54.8969611 +lon_0=37.0888079 +k=1.00000573 +x_0=-4756.69 +y_0=281.81 +ellps=krass

Проекция создана.

Альтернативный набор параметров

Пользовательская проекция старых моделей навигаторов GARMIN допускает ограниченный набор параметров. В ней не используется lat_0. Компенсируем его отсутствие коррекцией y_0. Напишем и запустим скрипт:

proj -f "%.4f" +proj=tmerc +lat_0=54.8969611 +lon_0=37.0888079 +k=1.00000573 +x_0=-4756.69 +y_0=281.81 +ellps=krass <<EOF
37.0888079 0
EOF

Вывод скрипта даёт координаты пересечения осевого меридиана с экватором:

-4756.6900 -6085619.5031

Вот эквивалентный набор параметров:

+proj=tmerc +lat_0=0 +lon_0=37.0888079 +k=1.00000573 +x_0=-4756.69 +y_0=-6085619.50 +ellps=krass

Может возникнуть вопрос, почему нельзя было сразу начать со значения lat_0 = 0, а не с параллели первого пункта? Ответ прост: устойчивые решения получаются только в окрестности заданных пунктов, и перенос начала координат на экватор должен быть выполнен в самом конце.

Коническая равноугольная проекция Ламберта

Эта проекция по использованию в геодезии занимает второе место после проекции Гаусса-Крюгера. Она, в отличие от последней, предоставляет полную свободу в выборе центрального меридиана, поскольку градиент масштаба вдоль параллели равен нулю.

Редко встречается бочка мёда без ложки дёгтя. Некоторые программы (MapInfo, например) предоставляют возможность задания этой проекции только в варианте с двумя стандартными параллелями, в то время как нам нужна единственная стандартная параллель, проходящая через выбранный пункт. По единственной стандартной параллели и масштабу на ней можно вычислить две стандартные широты тогда и только тогда, когда масштаб на единственной параллели меньше единицы. В общем, возможность имитации МСК таким вариантом с двумя стандартами 50:50.

К счастью, PROJ.4 не имеет этого досадного ограничения, и коническая равноугольная проекция Ламберта является неплохим вариантом создания аналога МСК для QGIS.

Долгота центрального меридиана

По развороту θ = 0.07766348° определим величину переноса центрального меридиана:

$\Delta L=-{\frac {\theta }{\sin B}}$

Получим ∆L = −0.0949292942, L = 37.0888204068. Запишем L в долготу центрального меридиана, нули в x_0 и y_0:

$ proj -f "%.4f" +proj=lcc +lat_0=54.8969611 +lon_0=37.0888204 +lat_1=54.8969611 +k_0=1 +x_0=0 +y_0=0 +ellps=krass pt_longlat.dat > pt_lcc.dat

Через пять итераций приходим к значению lon_0 = 37.0824027, при котором θ практически сводится к нулю.

Масштаб проекции

Полученное последним значение m = 1.00002378 запишем в параметр k_0:

+proj=lcc +lat_0=54.8969611 +lon_0=37.0824027 +lat_1=54.8969611 +k_0=1.00002378 +x_0=0 +y_0=0 +ellps=krass

Координаты в начале проекции

Вновь вычислим параметры конформного преобразования и перенесём a₀ = −5167.777 и a₁ = 281.494 в параметры проекции x_0, y_0:

+proj=lcc +lat_0=54.8969611 +lon_0=37.0824027 +lat_1=54.8969611 +k_0=1.00002378 +x_0=-5167.78 +y_0=281.49 +ellps=krass

Проекция готова. Невязки по результатам конформного преобразования достигают десяти сантиметров, что почти на порядок больше невязок для проекций Гаусса-Крюгера и косой Меркатора. Возможно, этого следовало ожидать, поскольку местная система координат построена на основе проекции Гаусса-Крюгера.

Косая стереографическая проекция

Эта проекция занимает незаслуженное третье место по широте использования геодезистами. Причина этого, возможно, кроется в вычислительных трудностях докомпьютерной эпохи: коническая равноугольная Ламберта выгодно отличается аскетической простотой формул, а проекция Гаусса-Крюгера — возможностью построения координатных таблиц, одинаковых для любой зоны. Между тем, если бы речь шла о создании МСК на пустом месте, а не об имитации уже созданной на основе другой проекции, то именно стереографическая является идеальным воплощением мечты о конформной проекции, искажения в которой минимальны на краях территории произвольной формы.

Можно реализовать имитацию МСК нашего примера на основе стереографической проекции. Если бы PROJ.4 включало параметр разворота, это было бы так же просто, как и реализация на основе косой проекции Меркатора. Следует ожидать, что при этом невязки были бы вдвое меньше, чем в реализации на основе конической равноугольной проекции Ламберта.

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

Заключение

Рассмотрена возоможность имитации МСК с использованием практически всех конформных проекций, широко используемых в геодезии: косой Меркатора, Гаусса-Крюгера, конической равноугольной Ламберта, косой стереографической. Для практического применения в QGIS рекомендуется косая проекция Меркатора. Если программное обеспечение не позволяет реализовать эту проекцию, либо угол разворота мал, предпочтительной является проекция Гаусса-Крюгера.

Ссылки

• Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987
• Coordinate Conversions and Transformations including Formulas, EPSG Guidance Note 7, 2002
• Projections Transform List
• Справка ArcGIS 10.1 — Локальная проекция Декартовой системы координат
• man_proj – PROJ.4
• Конформное преобразование