Добавление местной координатной системы в 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

Много букв

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

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

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

Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов координат имеется надёжный ключ в виде строки 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

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

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