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

Эта страница является черновиком статьи.

Введение

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

Характерная особенность МСК — за редкими исключениями малые углы разворота. Для малых населённых пунктов это достигалось ориентированием по магнитной буссоли и введением поправки за магнитное склонение. Погрешность при этом обычно находится в пределах одного–двух градусов. Для крупных городов ориентация могла выполняться на основе астрономических наблюдений. Если в таком случае намеренно или случайно расположить центр преобразования вблизи меридиана астропункта, то угол разворота составит несколько угловых секунд.

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

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

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

В качестве рабочей среды будем использовать командную строку UNIX. Это идеальный инструмент для экспериментирования, позволяющий непринуждённо сочетать PROJ.4, findkey и утилиты обработки текстовых потоков awk, pr.

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

Тестовый пример основан на сгенерированных данных. В таблице X, Y — координаты пункта в государственной системе (ГСК), а именно в седьмой зоне СК-42, x, y — координаты в местной системе (МСК), p — вес пункта.

ID	X	Y	x	y	p
1	7383477.64	6087377.60	1334.71	285.94	1.0
2	7382557.14	6081916.51	563.67	−5197.34	1.0
3	7386610.19	6088160.39	4444.27	1153.79	1.0
4	7381962.05	6090016.34	−252.07	2881.90	1.0

Подготовим файл исходных данных data0.dat:

1 7383477.64 6087377.60 1334.71   285.94 1.0
2 7382557.14 6081916.51  563.67 -5197.34 1.0
3 7386610.19 6088160.39 4444.27  1153.79 1.0
4 7381962.05 6090016.34 -252.07  2881.90 1.0

Убедимся в однородности данных обоих каталогов. Для этого вычислим параметры конформного преобразования:

$ findkey data0.dat key0.dat var0.dat

Сами по себе параметры нас не интересуют. Важны величины невязок, выведенных в файл var0.dat:

1 … -0.002 -0.001
2 … -0.017 0.013
3 … 0.031 0.017
4 … -0.012 -0.029

Малые значения невязок говорят о том, что данные достаточно хороши, чтобы с ними работать.

Пересчитаем координаты из ГСК в долготу/широту:

$ awk '{print $2, " ", $3}' data0.dat | proj -I -f "%.9f" +proj=tmerc +lat_0=0 +lon_0=39 +k=1 +x_0=7500000 +y_0=0 +ellps=krass > pt_longlat.dat

Команда awk извлекает из файла data0.dat вторую и третью колонки и передаёт их утилите proj, которая пересчитывает координаты из седьмой зоны СК-42 в долготу/широту и записывает результаты в файл pt_longlat.dat:

37.183749701 54.896961118
37.171630976 54.847714659
37.232243035 54.904709276
37.159058387 54.920296878

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

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

$ pr -m -t -s\  data0.dat pt_tmerc.dat | awk '{print $1, $7, $8, $4, $5, $6}' > data1.dat

Команда pr построчно объединяет данные из шести колонок 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 | tr -d '\r' > pt_omerc.dat

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

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

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

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

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

В нашем примере разворот оказался равным θ = 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_1:

+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

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