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

Материал из GIS-Lab
(Различия между версиями)
Перейти к: навигация, поиск
м (Подробнее о выборе априорных параметров)
(Ссылки)
 
(не показаны 66 промежуточных версий 1 участника)
Строка 6: Строка 6:
  
 
Под местной системой координат (МСК) будет подразумеваться так называемая «городская» система, построенная независимо от государственной системы (ГСК) и включенная в неё [http://gis-lab.info/qa/helmert2d.html#.D0.92.D0.B2.D0.B5.D0.B4.D0.B5.D0.BD.D0.B8.D0.B5 заданием ключей перехода к СК-42 или СК-63]. МСК крупных территорий, сравнимых с размерами субъектов Федерации, не являются предметом данной статьи, поскольку относятся к классическим картографическим проекциям.
 
Под местной системой координат (МСК) будет подразумеваться так называемая «городская» система, построенная независимо от государственной системы (ГСК) и включенная в неё [http://gis-lab.info/qa/helmert2d.html#.D0.92.D0.B2.D0.B5.D0.B4.D0.B5.D0.BD.D0.B8.D0.B5 заданием ключей перехода к СК-42 или СК-63]. МСК крупных территорий, сравнимых с размерами субъектов Федерации, не являются предметом данной статьи, поскольку относятся к классическим картографическим проекциям.
 +
 +
Многие программы ГИС по примеру геодезических программ позволяют реализовать работу в МСК непосредственно. Так, в QGIS и в MapInfo Pro любая проекция может быть дополнена аффинным преобразованием, а в Global Mapper конформные проекции дополняются разворотом. В данной статье рассматривается создание МСК в программах QGIS и MapInfo.
  
 
== Постановка задачи ==
 
== Постановка задачи ==
  
Имеется множество пунктов, для которых известны координаты ''X'', ''Y'' в ГСК и ''x'', ''y'' в МСК. Требуется подобрать проекцию и вычислить её параметры, удовлетворительно представляющую МСК в ГИС. При подборе параметров предполагается использовать один из пунктов в качестве центральной точки преобразования.
+
Имеется множество пунктов, для которых известны координаты ''X'', ''Y'' в ГСК и ''x'', ''y'' в МСК. Требуется подобрать проекцию, удовлетворительно представляющую МСК в ГИС. При подборе параметров предполагается использовать один из пунктов в качестве центральной точки преобразования.
  
Некоторые программы позволяют реализовать работу в МСК непосредственно. Так, в MapInfo любая проекция может быть дополнена аффинным преобразованием, а в Global Mapper конформным. ArcGIS в качестве МСК предлагает локальную проекцию: аналог ортографической проекции на эллипсоиде, дополненный разворотом и масштабированием.
+
== Подготовка данных ==
  
Другие программы, включая QGIS, работают только с «обыкновенными» проекциями, не допуская дополнительных геометрических преобразований. Задача статьи — показать, как можно конструировать проекции, позволяющие имитировать МСК в QGIS.
+
=== Исходные данные ===
 
+
В качестве рабочей среды будем использовать командную строку UNIX или OSGeo4W.
+
 
+
== Подготовка данных ==
+
  
Нужны два файла исходных данных. Пусть текстовый файл '''cat_s42z4.tsv''' содержит координаты пунктов в государственной системе (ГСК), а именно в четвёртой зоне СК-42:
+
Имеются два каталога. Текстовый файл '''cat_s42z4.tsv''' содержит координаты пунктов в государственной системе (ГСК), а именно в четвёртой зоне СК-42:
  
 
<pre>
 
<pre>
Строка 45: Строка 43:
 
Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы.
 
Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы.
  
Вычислим географические координаты пунктов:
+
=== Дополнительные данные ===
<syntaxhighlight lang="bash">
+
$ proj -I -f "%.16g" +init=epsg:28404 cat_s42z4.tsv > cat_longlat.tsv
+
</syntaxhighlight>
+
 
+
В файл '''cat_longlat.tsv''' запишутся географические долготы и широты в СК-42:
+
 
+
<pre>
+
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
+
</pre>
+
  
== Быстрая реализация ==
+
Очень важно помнить, что с точки зрения математической картографии МСК остаётся проекцией Гаусса-Крюгера и наследует её искажения. Поэтому важно знать, на какой именно ГСК основана МСК. Зачастую это заранее неизвестно, и приходится проводить предварительное исследование для выяснения этого вопроса.
  
Работа заключается в задании трёх параметров косой проекции Меркатора и вычислении оставшихся четырёх. Априорные параметры определяются так: центр проекции должен лежать на осевом меридиане исходной проекции СК-42/СК-63 на широте середины города, а азимут начальной линии должен быть малым. Пусть широта центра проекции равна широте первого пункта из файла '''cat_longlat.tsv''' 52.02642240080064°. Долготу центра проекции приравняем долготе осевого меридиана четвёртой зоны 21°. Азимут начальной линии приравняем произовольно малой величине −0,0001°. Масштабный множитель ''k'' приравняем единице, оставшиеся параметры обнулим, вычислим координаты в этой переходной проекции и запишем их в файл '''cat.tsv''':<ref>
+
В нашем примере мы предполагаем, что базовая ГСК либо СК-42 зона 4, либо СК-63 зона C0.
 +
Каталог в первой системе имеется, нужно подготовить каталог в альтернативной системе.
  
При работе в MSYS под Windows придётся позаботиться об удалении мусорных символов:
+
Параметры СК-63 зона C0 известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0".
 +
Создадим каталог в СК-63 с помощью утилиты '''proj''':
  
 
<syntaxhighlight lang="bash">
 
<syntaxhighlight lang="bash">
$ 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 | tr -d '\r' > cat.tsv
+
$ proj -I -f "%.17g" +init=epsg:28404 cat_s42z4.tsv > lonlat.tsv
 +
$ proj -f "%.17g" +proj=tmerc +lat_0=0.1 +lon_0=21.95 +k=1 +x_0=250000 +y_0=0 +ellps=krass lonlat.tsv > cat_s63c0.tsv
 
</syntaxhighlight>
 
</syntaxhighlight>
  
</ref>
+
В файл '''lonlat.tsv''' запишутся географические координаты (долготы и широты) в СК-42, а в '''cat_s63c0.tsv''' координаты в СК-63 зона C0:
 
+
<syntaxhighlight lang="bash">
+
$ 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
+
</syntaxhighlight>
+
 
+
Найдём оставшиеся параметры с помощью программы '''helmkey''':
+
 
+
<syntaxhighlight lang="bash">
+
$ helmkey cat.tsv cat_local.tsv
+
</syntaxhighlight>
+
 
+
Программа выведет на экран искомую строку параметров:
+
  
 
<pre>
 
<pre>
+k=0.9998372145697554 +x_0=-78707.08698765696 +y_0=32225.08703617829 +gamma=1.677027577390829
+
330797.45370592922 5755981.4751984337
 +
346208.04426388565 5757327.0953416284
 +
335051.73824425979 5770735.1441946141
 +
319180.795365442  5766563.182877643
 +
316233.72446517192 5752221.1970210578
 +
326744.90098082455 5742157.2318198672
 +
340603.31654394628 5744670.1910534762
 
</pre>
 
</pre>
  
Это всё! Осталось добавить в QGIS пользовательскую проекцию, не забыв необходимые стандартные параметры:
+
== Создание МСК ==
  
<pre>
+
Для вычислений используем консольную утилиту '''findkey''', о которой сказано ниже в приложении.
+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
+
</pre>
+
  
== Много букв ==
+
=== Определение базовой ГСК ===
  
=== Выбор исходной проекции ===
+
Вычислим параметры конформного преобразования координат из СК-42 в МСК. Для этого в командной строке запустим программу '''findkey''' с аргументами '''cat_s42z4.tsv''' и '''cat_local.tsv''':
  
Важной особенностью МСК является то, что с точки зрения математической картографии она остаётся проекцией Гаусса-Крюгера и наследует её искажения.
+
<syntaxhighlight lang="bash">
 
+
$ findkey cat_s42z4.tsv cat_local.tsv
Удобно, если осевой меридиан ГСК проходит через город при малом угле разворота, как в Пензе, например. В этом случае МСК строится как проекция Гаусса-Крюгера, отличающаяся от исходной ГСК малым сдвигом осевого меридиана для компенсации угла разворота. Жаль, что это случается редко.
+
</syntaxhighlight>
 
+
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
+
 
+
Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов координат имеется надёжный ключ в виде строки MapInfo.
+
  
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:
+
Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками.
 +
Посмотрим на содержимое '''var.csv''':
  
 
{| class="wikitable"
 
{| class="wikitable"
Строка 128: Строка 103:
 
|}
 
|}
  
В нашем случае второй кандидат на исходную ГСК — СК-63 зона C0 с осевым меридианом 21°57′. Создадим промежуточную проекцию на её основе, вычислим координаты для пунктов в ней и получим остающиеся параметры утилитой '''helmkey:
+
Вычислим параметры конформного преобразования координат из СК-63 в МСК:
  
 
<syntaxhighlight lang="bash">
 
<syntaxhighlight lang="bash">
$ 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
+
$ findkey cat_s63c0.tsv cat_local.tsv
$ helmkey cat.tsv cat_local.tsv
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Вывод программы:
+
Теперь содержимое '''var.csv''' будет таким:
 
+
<pre>
+
+k=1.00001854602098 +x_0=-13532.31228709544 +y_0=30742.75844830525 +gamma=0.9279005081790233
+
</pre>
+
 
+
Содержимое '''var.csv''':
+
  
 
{| class="wikitable"
 
{| class="wikitable"
 
|-
 
|-
 
|- align="right"
 
|- align="right"
| -0.000 || -0.002
+
| 0.000 || -0.002
 
|- align="right"
 
|- align="right"
 
| -0.001 || 0.002
 
| -0.001 || 0.002
Строка 161: Строка 129:
 
|}
 
|}
  
Сравнение невязок не в пользу СК-42. Исходная ГСК — СК-63 зона C0, и правильная проекция для QGIS такая:
+
Сравнение невязок позволяет сделать вывод, что базовая ГСК — СК-63 зона C0.
  
<pre>
+
=== Полученные параметры ===
+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
+
</pre>
+
  
=== Подробнее о выборе априорных параметров ===
+
Изучим содержимое файла '''key.txt''', соответствующего СК-63:
  
В плане распределения искажений проекция МСК должна имитировать ГСК. В идеале выбор параметров таков, что центр проекции лежит на точке осевого меридана ГСК, ближайшей к середине города, а направление начальной линии совпадает с направлением меридиана. Итак,
+
<pre>
* долгота центра проекции равна долготе осевого меридиана ГСК;
+
WKT:
* широта центра проекции требует вычисления;
+
A0 = -356718.938772419
* азимут начальной линии равен нулю.
+
A1 = 0.999887380509183
 +
A2 = 0.0161962611321084
 +
B0 = -5719887.1597502
 +
B1 = -0.0161962611321084
 +
B2 = 0.999887380509183
  
Довольно хорошее приближение для широты центра проекции даёт формула
+
MapInfo:
 +
0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198
  
tg ''φ''₂ = tg ''φ''₁ / cos (''λ''₂ − ''λ''₁) ,
+
Alternate set:
 +
scale = 1.000018546116108
 +
angle = 0.92800077023443683
 +
</pre>
  
где ''φ''₁, ''λ''₁ — координаты середины города; , ''φ''₂, ''λ''₂ — координаты центра проекции. Точное вычисление выполняется на апосфере.
+
Мы видим три группы чисел: WKT, MapInfo и Alternate set.
  
На практике идеал разбивается о нежелание авторов PROJ разрешить нулевой азимут начальной линии. Приходится использовать какое-нибудь малое его значение, но это приводит к смещению точки центра проекции, больше по широте, меньше по долготе.
+
Обратим внимание на параметр разворота angle из третьей группы. Если он мал, в пределах первых десятых долей градуса, имеет смысл отказаться от использования конформного преобразования и вместо этого смещать осевой меридиан для устранения угла разворота.
  
Впрочем, тема вычисления оптимального положения центра проекции выходит за рамки данной статьи. Между тем предложенный выше подход, при котором долгота равна долготе осевого меридиана, а широта равна широте середины территории, даёт вполне удовлетворительные результаты.
+
=== Создание МСК в MapInfo ===
  
== Переходная проекция ==
+
Запись базовой СК-63 зона C0 в файле '''MAPINFOW.PRJ''' должна выглядеть так:
  
Создадим в качестве переходной проекцию Гаусса-Крюгера, осевой меридиан и начальная параллель которой пересекаются в первом пункте, а координаты в этой точке совпадают с координатами МСК этого пункта. Вычислим координаты в переходной проекции и запишем их в файл '''pt_tmerc''':
+
<pre>
<syntaxhighlight lang="bash">
+
"CS63 zone C0", 8, 1001, 7, 21.95, 0.1, 1, 250000, 0
$ 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
+
</pre>
</syntaxhighlight>
+
  
Содержимое '''pt_tmerc.dat''':
+
Используем вторую группу из файла '''key.txt'''. Впишем МСК в '''MAPINFOW.PRJ''' как базовую, дополненную аффинным преобразованием:
  
 
<pre>
 
<pre>
1334.7100 285.9400
+
"Biala Podlaska", 1008, 1001, 7, 21.95, 0.1, 1, 250000, 0, 7, 0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198
556.2353 -5196.2593
+
4445.4005 1149.5698
+
-248.5461 2884.0431
+
 
</pre>
 
</pre>
  
Подготовим входной файл '''data1.dat''' для '''findkey'''. Нужно на место второй и третьей колонок в '''data0.dat''' записать данные '''pt_tmerc.dat'''.
+
Полноценное определение нуждается в параметрах Bounds:
 
+
<syntaxhighlight lang="bash">
+
$ paste data0.dat pt_tmerc.dat | awk '{print $1, $7, $8, $4, $5, $6}' > data1.dat
+
</syntaxhighlight>
+
 
+
Команда '''paste''' построчно объединяет данные из шести колонок '''data0.dat''' и двух колонок '''pt_tmerc.dat''' в общую строку, команда '''awk''' выводит в '''data1.dat''' нужные колонки в нужном же порядке:<ref>
+
Также можно использовать утилиту '''pr''' (после '''-s\''' два пробела!):
+
 
+
<syntaxhighlight lang="bash">
+
$ pr -m -t -s\  data0.dat pt_tmerc.dat | awk '{print $1, $7, $8, $4, $5, $6}' > data1.dat
+
</syntaxhighlight>
+
 
+
</ref>
+
  
 
<pre>
 
<pre>
1 1334.7100 285.9400 1334.71 285.94 1.0
+
"Biala Podlaska", 3008, 1001, 7, 21.95, 0.1, 1, 250000, 0, 7, 0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198, 52000, 15000, 82000, 45000
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
+
 
</pre>
 
</pre>
  
Нетрудно заметить, что координаты первого пункта в переходной проекции совпадают с таковыми в МСК. Для остальных пунктов это не так, что говорит о наличии масштабирования и разворота. Убедимся в этом:
+
=== Создание МСК в QGIS ===
  
<syntaxhighlight lang="bash">
+
Возьмём коэффициенты из первой группы в файле '''key.txt''' и создадим МСК в формате WKT как аффинное преобразование на основе проекции "Pulkovo 1942 / CS63 zone C0":
$ findkey data1.dat key1.dat var1.dat
+
</syntaxhighlight>
+
 
+
Содержимое '''key1.dat''':
+
  
 
<pre>
 
<pre>
0.390
+
DERIVEDPROJCRS["Biala Podlaska",
-1.813
+
    BASEPROJCRS["Pulkovo 1942 / CS63 zone C0",
1.0000052644
+
        BASEGEOGCRS["Pulkovo 1942",
0.0013554914
+
            DATUM["Pulkovo 1942",
1.0000061831
+
                ELLIPSOID["Krassowsky 1940",6378245,298.3,
+0.07766348
+
                    LENGTHUNIT["metre",1]]],
 +
            PRIMEM["Greenwich",0,
 +
                ANGLEUNIT["Degree",0.0174532925199433]],
 +
            ID["EPSG",4284]],
 +
        CONVERSION["CS63 zone C0",
 +
            METHOD["Transverse Mercator",
 +
                ID["EPSG",9807]],
 +
            PARAMETER["Latitude of natural origin",0.1,
 +
                ANGLEUNIT["degree",0.0174532925199433],
 +
                ID["EPSG",8801]],
 +
            PARAMETER["Longitude of natural origin",21.95,
 +
                ANGLEUNIT["degree",0.0174532925199433],
 +
                ID["EPSG",8802]],
 +
            PARAMETER["Scale factor at natural origin",1,
 +
                SCALEUNIT["unity",1],
 +
                ID["EPSG",8805]],
 +
            PARAMETER["False easting",250000,
 +
                LENGTHUNIT["metre",1],
 +
                ID["EPSG",8806]],
 +
            PARAMETER["False northing",0,
 +
                LENGTHUNIT["metre",1],
 +
                ID["EPSG",8807]]]
 +
    ],
 +
    DERIVINGCONVERSION["Affine",
 +
        METHOD["Affine parametric transformation",
 +
            ID["EPSG",9624]],
 +
        PARAMETER["A0",-356718.938772649,
 +
            LENGTHUNIT["metre",1],
 +
            ID["EPSG",8623]],
 +
        PARAMETER["A1",0.999887380509304,
 +
            SCALEUNIT["coefficient",1],
 +
            ID["EPSG",8624]],
 +
        PARAMETER["A2",0.0161962611321413,
 +
            SCALEUNIT["coefficient",1],
 +
            ID["EPSG",8625]],
 +
        PARAMETER["B0",-5719887.15975089,
 +
            LENGTHUNIT["metre",1],
 +
            ID["EPSG",8639]],
 +
        PARAMETER["B1",-0.0161962611321413,
 +
            SCALEUNIT["coefficient",1],
 +
            ID["EPSG",8640]],
 +
        PARAMETER["B2",0.999887380509304,
 +
            SCALEUNIT["coefficient",1],
 +
            ID["EPSG",8641]]],
 +
    CS[Cartesian,2],
 +
        AXIS["northing (X)",north,
 +
            ORDER[1],
 +
            LENGTHUNIT["metre",1]],
 +
        AXIS["easting (Y)",east,
 +
            ORDER[2],
 +
            LENGTHUNIT["metre",1]],
 +
    USAGE[
 +
        SCOPE["unknown"],
 +
        AREA["Europe - Poland - Biala Podlaska"],
 +
        BBOX[51.90,22.90,52.1,23.3]]]
 
</pre>
 
</pre>
  
Первые четыре параметра — ''a''<sub>0</sub>, ''b''<sub>0</sub>, ''a''<sub>1</sub>, ''b''<sub>1</sub>, оставшиеся два — масштабный множитель ''m'' и разворот ''θ'' в градусах.
+
В первой строке записали название системы координат латиницей "Biala Podlaska".
 +
В конце вставили название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате ''φ''<sub>min</sub>, ''λ''<sub>min</sub>, ''φ''<sub>max</sub>, ''λ''<sub>max</sub>.
  
== Косая проекция Меркатора ==
+
=== Создание МСК в Global Mapper ===
  
Всё готово для конструирования косой проекции Меркатора. Зададим, как для переходной проекции, долготу и широту первого пункта в качестве начальных, его координаты в МСК в качестве «фальшивых». Кроме того, перенесём масштаб конформного преобразования в масштаб проекции, а разворот — в азимут линии минимального масштаба. Вычислим координаты в этой проекции и запишем в файл '''pt_omerc.dat''':
+
Третья группа параметров в файле '''key.txt''' содержит масштабный коэффициент scale и угол поворота angle. Угол вставляем как есть, а масштабный коэффициент умножим на соответствующий параметр базовой СК.
  
<syntaxhighlight lang="bash">
+
== Заключение ==
$ 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
+
</syntaxhighlight>
+
  
Вычислим параметры конформного преобразования для '''pt_omerc.dat''': ''m'' = 1.0000000000, ''θ'' = −0.00000022° = −0.0008″. Практически масштаб сведён к единице, разворот к нулю, а вычисленные из значений долготы/широты координаты совпадают с координатами, заданными в МСК. Таким образом, проекция, эквивалентная МСК, построена.
+
Задачи построения МСК для QGIS и MapInfo выполнены, цель достигнута.
  
Следует заметить, что многие программы (например, MapInfo) используют косую проекцию Меркатора в версии Хотина, когда началом координат служит точка пересечения линии минимального масштаба с экватором. Азимут линии и плоские координаты в этой точке, естественно, отличаются от заданных нами, и понадобятся дополнительные вычисления для их определения.
+
== Приложение. Утилита findkey ==
  
Реализация косой проекции Меркатора в PROJ.4 имеет одну особенность: программы просто не работают со значениями параметра '''alpha''' 0 или 90°. Если угол разворота очень мал (скажем, несколько угловых секунд), следует использовать проекцию Гаусса-Крюгера.
+
Программа '''findkey''' вычисляет параметры конформного преобразования. Написана на языке C. Вот листинг:
  
== Проекция Гаусса-Крюгера ==
+
<syntaxhighlight lang="c">
 +
#include <stdio.h>
 +
#include <stdlib.h>
 +
#include <math.h>
  
Встречаются программы, набор проекций в которых ограничен. Проекция Гаусса-Крюгера — она же Transverse Mercator — найдётся в самом бедном наборе. В приёмниках GARMIN пользовательская система координат может быть задана только в этой проекции.
+
#define SEP ';' /* var-file column separator */
  
=== Долгота осевого меридиана ===
+
/* --------------------------------------------------------------------------
 +
* findkey
 +
*
 +
* Program to compute Helmert 2D transformation parameters
 +
*
 +
* Usage: findkey <coord1> <coord2>
 +
*
 +
* Input files: coord1 coord2
 +
*    coord1 - source coordinate 'x1 y1'
 +
*    coord2 - destination coordinate 'x2 y2'
 +
*              a row per a point; 3+ points
 +
*
 +
* Output files:
 +
*    key.txt - transformation parameters
 +
*    var.csv - SEP separated residuals 'dx dy'
 +
* -------------------------------------------------------------------------- */
 +
int main(int argc, char *argv[])
 +
{
 +
  char buf0[1024], buf1[1024];
 +
  double x[2], y[2];
 +
  double xc[2], yc[2];
 +
  double dx[2], dy[2];
 +
  double s[7] = {0., 0., 0., 0., 0.};
 +
  double det, h[6];
 +
  double mu, theta;
 +
  double yh[2];
 +
  int i;
 +
  FILE *fp0, *fp1, *fp2;
  
Если угол разворота не превышает первых десятков угловых минут, можно компенсировать его переносом осевого меридиана в сторону и не бояться потери геодезической точности. Даже первые градусы не должны существенно повлиять на навигационные приложения. Дальнейшее увеличение угла приводит к стремительному росту градиента масштаба в направлении запад-восток и, следовательно, искажению подобия фигур.
+
  if (argc < 3) {
 +
    printf("Usage: findkey <coord1> <coord2>\n");
 +
    exit(EXIT_FAILURE);
 +
  }
  
В нашем примере разворот оказался небольшим: ''θ'' = 0.07766348° = 0°04′39.59″. Значит, в данном случае на основе проекции Гаусса-Крюгера можно создать довольно точный аналог МСК.
+
  if ((fp0 = fopen(argv[1], "r")) == NULL) {
 +
    printf("can't open %s\n", argv[1]);
 +
    exit(EXIT_FAILURE);
 +
  }
  
Предварительное значение перемещения осевого меридиана можно оценить по формуле:
+
  if ((fp1 = fopen(argv[2], "r")) == NULL) {
 +
    printf("can't open %s\n", argv[2]);
 +
    exit(EXIT_FAILURE);
 +
  }
  
<math>\Delta L = -\operatorname{arc\,tg}\frac{\operatorname{tg} \theta}{\sin B}</math>
+
  /* coordinate sums */
 +
  while (fgets(buf0, 1024, fp0) != NULL && fgets(buf1, 1024, fp1) != NULL) {
 +
    sscanf(buf0, "%lf %lf", &x[0], &x[1]);
 +
    sscanf(buf1, "%lf %lf", &y[0], &y[1]);
 +
    s[0] += x[0];
 +
    s[1] += x[1];
 +
    s[2] += y[0];
 +
    s[3] += y[1];
 +
    s[4] += 1.;
 +
  }
 +
  rewind(fp0);
 +
  rewind(fp1);
  
где ∆''L'' — перенос осевого меридиана, ''θ'' — угол разворота, ''B'' — широта первого пункта. Получаем ∆''L'' = −0.0949292655, ''L'' = 37.0888204345. Построим проекцию Гаусса-Крюгера на основе переходной, подставив ''L'' в долготу осевого меридиана и нули в «фальшивые» координаты:
+
  /* centrum gravitatis */
 +
  for (i = 0; i < 2; i++) {
 +
    xc[i] = s[i] / s[4];
 +
    yc[i] = s[2 + i] / s[4];
 +
  }
  
<syntaxhighlight lang="bash">
+
  /* sums of products */
$ 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
+
  for (i = 0; i < 7; i++)
</syntaxhighlight>
+
    s[i] = 0.;
 +
  while (fgets(buf0, 1024, fp0) != NULL && fgets(buf1, 1024, fp1) != NULL) {
 +
    sscanf(buf0, "%lf %lf", &x[0], &x[1]);
 +
    sscanf(buf1, "%lf %lf", &y[0], &y[1]);
 +
    /* coordinate differences */
 +
    dx[0] = x[0] - xc[0];
 +
    dx[1] = x[1] - xc[1];
 +
    dy[0] = y[0] - yc[0];
 +
    dy[1] = y[1] - yc[1];
 +
    /* summation */
 +
    s[0] += dx[0] * dx[0];
 +
    /*s[1] += dx[0] * dx[1];*/
 +
    s[2] += dx[1] * dx[1];
 +
    s[3] += dx[0] * dy[0];
 +
    s[4] += dx[1] * dy[0];
 +
    s[5] += dx[0] * dy[1];
 +
    s[6] += dx[1] * dy[1];
 +
  }
 +
  rewind(fp0);
 +
  rewind(fp1);
  
Последующее вычисление параметров конформного преобразования выдаёт сравнительно большой угол разворота. Однако уже после второй итерации с ''L'' = 37.0888079 получаем ''m'' = 1.0000057306, ''θ'' = +0.00000044.
+
  /* Helmert parameters */
 +
  det = s[0] + s[2];
 +
  h[0] = (s[3] + s[6]) / det;
 +
  h[1] = (s[4] - s[5]) / det;
 +
  h[2] = yc[0] - h[0] * xc[0] - h[1] * xc[1];
 +
  h[3] = -h[1];
 +
  h[4] = h[0];
 +
  h[5] = yc[1] - h[3] * xc[0] - h[4] * xc[1];
  
=== Масштаб проекции ===
+
  /* alternative Helmert parameter set */
 +
  mu = hypot(h[0], h[1]);
 +
  theta = atan2(h[1], h[0]) / M_PI * 180.;
  
Перенесём ''m'' в масштаб проекции:
+
  /* output parameters */
 +
  if ((fp2 = fopen("key.txt", "w")) == NULL) {
 +
    printf("can't create %s\n", "key.txt");
 +
    exit(EXIT_FAILURE);
 +
  }
 +
  fprintf(fp2, "WKT:\nA0 = %.15g\nA1 = %.15g\nA2 = %.15g\nB0 = %.15g\nB1 = %.15g\nB2 = %.15g\n",
 +
  h[2], h[0], h[1], h[5], h[3], h[4]);
 +
  fprintf(fp2, "\nMapInfo:\n%.12f, %.12f, %.17g, %.12f, %.12f, %.17g\n",
 +
  h[0], h[1], h[2], h[3], h[4], h[5]);
 +
  fprintf(fp2, "\nAlternate set:\nscale = %.17g\nangle = %.17g\n", mu, theta);
 +
  fclose(fp2);
  
<syntaxhighlight lang="bash">
+
  /* output residuals */
+proj=tmerc +lat_0=54.8969611 +lon_0=37.0888079 +k=1.00000573 +x_0=0 +y_0=0 +ellps=krass
+
  if ((fp2 = fopen("var.csv", "w")) == NULL) {
</syntaxhighlight>
+
    printf("can't create %s\n", "var.csv");
 +
    exit(EXIT_FAILURE);
 +
  }
 +
  while (fgets(buf0, 1024, fp0) != NULL && fgets(buf1, 1024, fp1) != NULL) {
 +
    sscanf(buf0, "%lf %lf", &x[0], &x[1]);
 +
    sscanf(buf1, "%lf %lf", &y[0], &y[1]);
 +
    /* model y */
 +
    yh[0] = h[0] * x[0] + h[1] * x[1] + h[2];
 +
    yh[1] = h[3] * x[0] + h[4] * x[1] + h[5];
 +
    fprintf(fp2, "%.15g%c%.15g\n", yh[0] - y[0], SEP, yh[1] - y[1]);
 +
  }
 +
  fclose(fp2);
 +
  fclose(fp1);
 +
  fclose(fp0);
  
Вновь вычислим параметры конформного преобразования. Теперь ''m'' = 1.0000000039, т.е. практически единица.
+
  exit(EXIT_SUCCESS);
 
+
}
=== Координаты в начальной точке проекции ===
+
 
+
Наконец обратим внимание на первичные параметры ''a''<sub>0</sub> = −4756.693 и ''a''<sub>1</sub> = 281.807. Перенесём их значения в параметры проекции '''x_0''' и '''y_0''':
+
 
+
<syntaxhighlight lang="bash">
+
+proj=tmerc +lat_0=54.8969611 +lon_0=37.0888079 +k=1.00000573 +x_0=-4756.69 +y_0=281.81 +ellps=krass
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Проекция создана.
+
Сохраним код в файл '''findkey.c'''. Создадим исполняемый модуль компилятором '''gcc''':
 
+
=== Альтернативный набор параметров ===
+
 
+
Пользовательская проекция старых моделей навигаторов GARMIN допускает ограниченный набор параметров. В ней не используется '''lat_0'''. Компенсируем его отсутствие коррекцией '''y_0'''. Напишем и запустим скрипт:
+
  
 
<syntaxhighlight lang="bash">
 
<syntaxhighlight lang="bash">
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
+
$ gcc -o findkey findkey.c -lm
37.0888079 0
+
EOF
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Вывод скрипта даёт координаты пересечения осевого меридиана с экватором:
+
Пользователи MS Windows могут загрузить уже скомпилированную [https://wiki.gis-lab.info/images/4/43/Findkey.zip программу].
 
+
<pre>-4756.6900 -6085619.5031</pre>
+
 
+
Вот эквивалентный набор параметров:
+
 
+
<pre>+proj=tmerc +lat_0=0 +lon_0=37.0888079 +k=1.00000573 +x_0=-4756.69 +y_0=-6085619.50 +ellps=krass</pre>
+
 
+
Может возникнуть вопрос, почему нельзя было сразу начать со значения '''lat_0''' = 0, а не с параллели первого пункта? Ответ прост: устойчивые решения получаются только в окрестности заданных пунктов, и перенос начала координат на экватор должен быть выполнен в самом конце.
+
 
+
== Коническая равноугольная проекция Ламберта ==
+
 
+
Эта проекция по использованию в геодезии занимает второе место после проекции Гаусса-Крюгера. Она, в отличие от последней, предоставляет полную свободу в выборе центрального меридиана, поскольку градиент масштаба вдоль параллели равен нулю.
+
 
+
Редко встречается бочка мёда без ложки дёгтя. Некоторые программы (MapInfo, например) предоставляют возможность задания этой проекции только в варианте с двумя стандартными параллелями, в то время как нам нужна единственная стандартная параллель, проходящая через выбранный пункт. По единственной стандартной параллели и масштабу на ней можно вычислить две стандартные широты тогда и только тогда, когда масштаб на единственной параллели меньше единицы. В общем, возможность имитации МСК таким вариантом с двумя стандартами 50:50.
+
 
+
К счастью, PROJ.4 не имеет этого досадного ограничения, и коническая равноугольная проекция Ламберта является неплохим вариантом создания аналога МСК для QGIS.
+
 
+
=== Долгота центрального меридиана ===
+
 
+
По развороту ''θ'' = 0.07766348° определим величину переноса центрального меридиана:
+
 
+
<math>\Delta L = -\frac{\theta}{\sin B}</math>
+
 
+
Получим ∆''L'' = −0.0949292942, ''L'' = 37.0888204068. Запишем ''L'' в долготу центрального меридиана, нули в '''x_0''' и '''y_0''':
+
 
+
<syntaxhighlight lang="bash">
+
$ 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
+
</syntaxhighlight>
+
 
+
Через пять итераций приходим к значению '''lon_0''' = 37.0824027, при котором ''θ'' практически сводится к нулю.
+
 
+
=== Масштаб проекции ===
+
 
+
Полученное последним значение ''m'' = 1.00002378 запишем в параметр '''k_0''':
+
 
+
<syntaxhighlight lang="bash">
+
+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
+
</syntaxhighlight>
+
 
+
=== Координаты в начале проекции ===
+
 
+
Вновь вычислим параметры конформного преобразования и перенесём ''a''<sub>0</sub> = −5167.777 и ''a''<sub>1</sub> = 281.494 в параметры проекции '''x_0''', '''y_0''':
+
 
+
<syntaxhighlight lang="bash">
+
+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
+
</syntaxhighlight>
+
 
+
Проекция готова. Невязки по результатам конформного преобразования достигают десяти сантиметров, что почти на порядок больше невязок для проекций Гаусса-Крюгера и косой Меркатора. Возможно, этого следовало ожидать, поскольку местная система координат построена на основе проекции Гаусса-Крюгера.
+
 
+
== Косая стереографическая проекция ==
+
 
+
Эта проекция занимает незаслуженное третье место по широте использования геодезистами. Причина этого, возможно, кроется в вычислительных трудностях докомпьютерной эпохи: коническая равноугольная Ламберта выгодно отличается аскетической простотой формул, а проекция Гаусса-Крюгера — возможностью построения координатных таблиц, одинаковых для любой зоны. Между тем, если бы речь шла о создании МСК на пустом месте, а не об имитации уже созданной на основе другой проекции, то именно стереографическая является идеальным воплощением мечты о конформной проекции, искажения в которой минимальны на краях территории произвольной формы.
+
 
+
Можно реализовать имитацию МСК нашего примера на основе стереографической проекции. Если бы PROJ.4 включало параметр разворота, это было бы так же просто, как и реализация на основе косой проекции Меркатора. Следует ожидать, что при этом невязки были бы вдвое меньше, чем в реализации на основе конической равноугольной проекции Ламберта.
+
 
+
Отсутствие параметра разворота заставляет пойти по пути, дважды пройденному выше: перенос осевого меридиана для компенсации разворота, задание масштаба проекции, коррекция плоских координат. Оставим реализацию в качестве упражнения для любителей вычислений.
+
 
+
== Заключение ==
+
 
+
Рассмотрена возоможность имитации МСК с использованием практически всех конформных проекций, широко используемых в геодезии: косой Меркатора, Гаусса-Крюгера, конической равноугольной Ламберта, косой стереографической. Для практического применения в QGIS рекомендуется косая проекция Меркатора. Если программное обеспечение не позволяет реализовать эту проекцию, либо угол разворота мал, предпочтительной является проекция Гаусса-Крюгера.
+
 
+
== Примечания ==
+
 
+
<references />
+
  
 
== Ссылки ==
 
== Ссылки ==
  
 
* [http://pubs.usgs.gov/pp/1395/report.pdf Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987]
 
* [http://pubs.usgs.gov/pp/1395/report.pdf Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987]
* [http://remotesensing.org/geotiff/proj_list/guid7.html Coordinate Conversions and Transformations including Formulas, EPSG Guidance Note 7, 2002]
+
* [http://www.epsg.org/Portals/0/373-07-2.pdf Coordinate Conversions and Transformations including Formulas, EPSG Guidance Note 7-2]
* [http://remotesensing.org/geotiff/proj_list/index.html Projections Transform List]
+
* [https://desktop.arcgis.com/ru/arcmap/latest/map/projections/local-cartesian-projection.htm Справка ArcGIS — Локальная проекция Декартовой системы координат]
* [http://resources.arcgis.com/ru/help/main/10.1/index.html#//003r00000035000000 Справка ArcGIS 10.1 — Локальная проекция Декартовой системы координат]
+
* [https://proj4.org/usage/index.html Using PROJ]
* [http://trac.osgeo.org/proj/wiki/man_proj man_proj – PROJ.4]
+
 
* [http://gis-lab.info/qa/helmert2d.html Конформное преобразование]
 
* [http://gis-lab.info/qa/helmert2d.html Конформное преобразование]
 +
* [https://gis.stackexchange.com/questions/353022/defining-a-coordinate-system-in-wkt-or-proj-format-that-has-an-affine-transforma Defining a coordinate system in WKT or PROJ format that has an Affine transformaiton and bounds]

Текущая версия на 05:44, 13 мая 2020

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


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

Содержание

[править] Введение

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

Многие программы ГИС по примеру геодезических программ позволяют реализовать работу в МСК непосредственно. Так, в QGIS и в MapInfo Pro любая проекция может быть дополнена аффинным преобразованием, а в Global Mapper конформные проекции дополняются разворотом. В данной статье рассматривается создание МСК в программах QGIS и MapInfo.

[править] Постановка задачи

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

[править] Подготовка данных

[править] Исходные данные

Имеются два каталога. Текстовый файл 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

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

[править] Дополнительные данные

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

В нашем примере мы предполагаем, что базовая ГСК либо СК-42 зона 4, либо СК-63 зона C0. Каталог в первой системе имеется, нужно подготовить каталог в альтернативной системе.

Параметры СК-63 зона C0 известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0". Создадим каталог в СК-63 с помощью утилиты proj:

$ proj -I -f "%.17g" +init=epsg:28404 cat_s42z4.tsv > lonlat.tsv
$ proj -f "%.17g" +proj=tmerc +lat_0=0.1 +lon_0=21.95 +k=1 +x_0=250000 +y_0=0 +ellps=krass lonlat.tsv > cat_s63c0.tsv

В файл lonlat.tsv запишутся географические координаты (долготы и широты) в СК-42, а в cat_s63c0.tsv координаты в СК-63 зона C0:

330797.45370592922 5755981.4751984337
346208.04426388565 5757327.0953416284
335051.73824425979 5770735.1441946141
319180.795365442   5766563.182877643
316233.72446517192 5752221.1970210578
326744.90098082455 5742157.2318198672
340603.31654394628 5744670.1910534762

[править] Создание МСК

Для вычислений используем консольную утилиту findkey, о которой сказано ниже в приложении.

[править] Определение базовой ГСК

Вычислим параметры конформного преобразования координат из СК-42 в МСК. Для этого в командной строке запустим программу findkey с аргументами cat_s42z4.tsv и cat_local.tsv:

$ findkey cat_s42z4.tsv cat_local.tsv

Программа создаст два файла: key.txt с параметрами конформного преобразования и var.csv с координатными невязками. Посмотрим на содержимое 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 в МСК:

$ findkey cat_s63c0.tsv cat_local.tsv

Теперь содержимое 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

Сравнение невязок позволяет сделать вывод, что базовая ГСК — СК-63 зона C0.

[править] Полученные параметры

Изучим содержимое файла key.txt, соответствующего СК-63:

WKT:
A0 = -356718.938772419
A1 = 0.999887380509183
A2 = 0.0161962611321084
B0 = -5719887.1597502
B1 = -0.0161962611321084
B2 = 0.999887380509183

MapInfo:
0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198

Alternate set:
scale = 1.000018546116108
angle = 0.92800077023443683

Мы видим три группы чисел: WKT, MapInfo и Alternate set.

Обратим внимание на параметр разворота angle из третьей группы. Если он мал, в пределах первых десятых долей градуса, имеет смысл отказаться от использования конформного преобразования и вместо этого смещать осевой меридиан для устранения угла разворота.

[править] Создание МСК в MapInfo

Запись базовой СК-63 зона C0 в файле MAPINFOW.PRJ должна выглядеть так:

"CS63 zone C0", 8, 1001, 7, 21.95, 0.1, 1, 250000, 0

Используем вторую группу из файла key.txt. Впишем МСК в MAPINFOW.PRJ как базовую, дополненную аффинным преобразованием:

"Biala Podlaska", 1008, 1001, 7, 21.95, 0.1, 1, 250000, 0, 7, 0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198

Полноценное определение нуждается в параметрах Bounds:

"Biala Podlaska", 3008, 1001, 7, 21.95, 0.1, 1, 250000, 0, 7, 0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198, 52000, 15000, 82000, 45000

[править] Создание МСК в QGIS

Возьмём коэффициенты из первой группы в файле key.txt и создадим МСК в формате WKT как аффинное преобразование на основе проекции "Pulkovo 1942 / CS63 zone C0":

DERIVEDPROJCRS["Biala Podlaska",
    BASEPROJCRS["Pulkovo 1942 / CS63 zone C0",
        BASEGEOGCRS["Pulkovo 1942",
            DATUM["Pulkovo 1942",
                ELLIPSOID["Krassowsky 1940",6378245,298.3,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["Degree",0.0174532925199433]],
            ID["EPSG",4284]],
        CONVERSION["CS63 zone C0",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0.1,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",21.95,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",1,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",250000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]]
    ],
    DERIVINGCONVERSION["Affine",
        METHOD["Affine parametric transformation",
            ID["EPSG",9624]],
        PARAMETER["A0",-356718.938772649,
            LENGTHUNIT["metre",1],
            ID["EPSG",8623]],
        PARAMETER["A1",0.999887380509304,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8624]],
        PARAMETER["A2",0.0161962611321413,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8625]],
        PARAMETER["B0",-5719887.15975089,
            LENGTHUNIT["metre",1],
            ID["EPSG",8639]],
        PARAMETER["B1",-0.0161962611321413,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8640]],
        PARAMETER["B2",0.999887380509304,
            SCALEUNIT["coefficient",1],
            ID["EPSG",8641]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Europe - Poland - Biala Podlaska"],
        BBOX[51.90,22.90,52.1,23.3]]]

В первой строке записали название системы координат латиницей "Biala Podlaska". В конце вставили название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате φmin, λmin, φmax, λmax.

[править] Создание МСК в Global Mapper

Третья группа параметров в файле key.txt содержит масштабный коэффициент scale и угол поворота angle. Угол вставляем как есть, а масштабный коэффициент умножим на соответствующий параметр базовой СК.

[править] Заключение

Задачи построения МСК для QGIS и MapInfo выполнены, цель достигнута.

[править] Приложение. Утилита findkey

Программа findkey вычисляет параметры конформного преобразования. Написана на языке C. Вот листинг:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define SEP ';' /* var-file column separator */
 
/* --------------------------------------------------------------------------
 * findkey
 *
 * Program to compute Helmert 2D transformation parameters
 *
 * Usage: findkey <coord1> <coord2>
 *
 * Input files: coord1 coord2
 *     coord1 - source coordinate 'x1 y1'
 *     coord2 - destination coordinate 'x2 y2'
 *              a row per a point; 3+ points
 *
 * Output files:
 *    key.txt - transformation parameters
 *    var.csv - SEP separated residuals 'dx dy'
 * -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
  char buf0[1024], buf1[1024];
  double x[2], y[2];
  double xc[2], yc[2];
  double dx[2], dy[2];
  double s[7] = {0., 0., 0., 0., 0.};
  double det, h[6];
  double mu, theta;
  double yh[2];
  int i;
  FILE *fp0, *fp1, *fp2;
 
  if (argc < 3) {
    printf("Usage: findkey <coord1> <coord2>\n");
    exit(EXIT_FAILURE);
  }
 
  if ((fp0 = fopen(argv[1], "r")) == NULL) {
    printf("can't open %s\n", argv[1]);
    exit(EXIT_FAILURE);
  }
 
  if ((fp1 = fopen(argv[2], "r")) == NULL) {
    printf("can't open %s\n", argv[2]);
    exit(EXIT_FAILURE);
  }
 
  /* coordinate sums */
  while (fgets(buf0, 1024, fp0) != NULL && fgets(buf1, 1024, fp1) != NULL) {
    sscanf(buf0, "%lf %lf", &x[0], &x[1]);
    sscanf(buf1, "%lf %lf", &y[0], &y[1]);
    s[0] += x[0];
    s[1] += x[1];
    s[2] += y[0];
    s[3] += y[1];
    s[4] += 1.;
  }
  rewind(fp0);
  rewind(fp1);
 
  /* centrum gravitatis */
  for (i = 0; i < 2; i++) {
    xc[i] = s[i] / s[4];
    yc[i] = s[2 + i] / s[4];
  }
 
  /* sums of products */
  for (i = 0; i < 7; i++)
    s[i] = 0.;
  while (fgets(buf0, 1024, fp0) != NULL && fgets(buf1, 1024, fp1) != NULL) {
    sscanf(buf0, "%lf %lf", &x[0], &x[1]);
    sscanf(buf1, "%lf %lf", &y[0], &y[1]);
    /* coordinate differences */
    dx[0] = x[0] - xc[0];
    dx[1] = x[1] - xc[1];
    dy[0] = y[0] - yc[0];
    dy[1] = y[1] - yc[1];
    /* summation */
    s[0] += dx[0] * dx[0];
    /*s[1] += dx[0] * dx[1];*/
    s[2] += dx[1] * dx[1];
    s[3] += dx[0] * dy[0];
    s[4] += dx[1] * dy[0];
    s[5] += dx[0] * dy[1];
    s[6] += dx[1] * dy[1];
  }
  rewind(fp0);
  rewind(fp1);
 
  /* Helmert parameters */
  det = s[0] + s[2];
  h[0] = (s[3] + s[6]) / det;
  h[1] = (s[4] - s[5]) / det;
  h[2] = yc[0] - h[0] * xc[0] - h[1] * xc[1];
  h[3] = -h[1];
  h[4] = h[0];
  h[5] = yc[1] - h[3] * xc[0] - h[4] * xc[1];
 
  /* alternative Helmert parameter set */
  mu = hypot(h[0], h[1]);
  theta = atan2(h[1], h[0]) / M_PI * 180.;
 
  /* output parameters */
  if ((fp2 = fopen("key.txt", "w")) == NULL) {
    printf("can't create %s\n", "key.txt");
    exit(EXIT_FAILURE);
  }
  fprintf(fp2, "WKT:\nA0 = %.15g\nA1 = %.15g\nA2 = %.15g\nB0 = %.15g\nB1 = %.15g\nB2 = %.15g\n",
	  h[2], h[0], h[1], h[5], h[3], h[4]);
  fprintf(fp2, "\nMapInfo:\n%.12f, %.12f, %.17g, %.12f, %.12f, %.17g\n",
	  h[0], h[1], h[2], h[3], h[4], h[5]);
  fprintf(fp2, "\nAlternate set:\nscale = %.17g\nangle = %.17g\n", mu, theta);
  fclose(fp2);
 
  /* output residuals */
  if ((fp2 = fopen("var.csv", "w")) == NULL) {
    printf("can't create %s\n", "var.csv");
    exit(EXIT_FAILURE);
  }
  while (fgets(buf0, 1024, fp0) != NULL && fgets(buf1, 1024, fp1) != NULL) {
    sscanf(buf0, "%lf %lf", &x[0], &x[1]);
    sscanf(buf1, "%lf %lf", &y[0], &y[1]);
    /* model y */
    yh[0] = h[0] * x[0] + h[1] * x[1] + h[2];
    yh[1] = h[3] * x[0] + h[4] * x[1] + h[5];
    fprintf(fp2, "%.15g%c%.15g\n", yh[0] - y[0], SEP, yh[1] - y[1]);
  }
  fclose(fp2);
  fclose(fp1);
  fclose(fp0);
 
  exit(EXIT_SUCCESS);
}

Сохраним код в файл findkey.c. Создадим исполняемый модуль компилятором gcc:

$ gcc -o findkey findkey.c -lm

Пользователи MS Windows могут загрузить уже скомпилированную программу.

[править] Ссылки

Персональные инструменты
Пространства имён

Варианты
Действия
Статьи
Спецпроекты
Инструменты