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

Материал из GIS-Lab
(Различия между версиями)
Перейти к: навигация, поиск
м (Выбор исходной проекции)
(Ссылки)
 
(не показаны 72 промежуточные версии 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 зона 4, либо СК-63 зона C0.
 +
Каталог в первой системе имеется, нужно подготовить каталог в альтернативной системе.
  
Работа заключается в задании трёх параметров косой проекции Меркатора и вычислении оставшихся четырёх. Априорные параметры определяются так: центр проекции должен лежать на осевом меридиане исходной проекции СК-42/СК-63 на широте середины города, а азимут начальной линии должен быть малым. Пусть широта центра проекции равна широте первого пункта из файла '''cat_longlat.tsv''' 52.02642240080064°. Долготу центра проекции приравняем долготе осевого меридиана четвёртой зоны 21°. Азимут начальной линии приравняем произовольно малой величине −0,0001°. Масштабный множитель ''k'' приравняем единице, оставшиеся параметры обнулим, вычислим координаты в этой переходной проекции и запишем их в файл '''cat.tsv''':
+
Параметры СК-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 > 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>
  
Найдём оставшиеся параметры с помощью программы '''helmkey''':
+
В файл '''lonlat.tsv''' запишутся географические координаты (долготы и широты) в СК-42, а в '''cat_s63c0.tsv''' координаты в СК-63 зона C0:
 
+
<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>
 
+
Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов координат имеется надёжный ключ в виде строки MapInfo.
+
  
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:
+
Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками.
 +
Посмотрим на содержимое '''var.csv''':
  
 
{| class="wikitable"
 
{| class="wikitable"
Строка 118: Строка 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
+
+k=1.00001854602098 +x_0=-13532.31228709544 +y_0=30742.75844830525 +gamma=0.9279005081790233
+
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Посмотрим на '''var.csv''':
+
Теперь содержимое '''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
Строка 146: Строка 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:
  
Создадим в качестве переходной проекцию Гаусса-Крюгера, осевой меридиан и начальная параллель которой пересекаются в первом пункте, а координаты в этой точке совпадают с координатами МСК этого пункта. Вычислим координаты в переходной проекции и запишем их в файл '''pt_tmerc''':
+
<pre>
<syntaxhighlight lang="bash">
+
WKT:
$ 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
+
A0 = -356718.938772419
</syntaxhighlight>
+
A1 = 0.999887380509183
 +
A2 = 0.0161962611321084
 +
B0 = -5719887.1597502
 +
B1 = -0.0161962611321084
 +
B2 = 0.999887380509183
  
Содержимое '''pt_tmerc.dat''':
+
MapInfo:
 +
0.999887380509, 0.016196261132, -356718.93877241889, -0.016196261132, 0.999887380509, -5719887.159750198
  
<pre>
+
Alternate set:
1334.7100 285.9400
+
scale = 1.000018546116108
556.2353 -5196.2593
+
angle = 0.92800077023443683
4445.4005 1149.5698
+
-248.5461 2884.0431
+
 
</pre>
 
</pre>
  
Подготовим входной файл '''data1.dat''' для '''findkey'''. Нужно на место второй и третьей колонок в '''data0.dat''' записать данные '''pt_tmerc.dat'''.
+
Мы видим три группы чисел: WKT, MapInfo и Alternate set.
  
<syntaxhighlight lang="bash">
+
Обратим внимание на параметр разворота angle из третьей группы. Если он мал, в пределах первых десятых долей градуса, имеет смысл отказаться от использования конформного преобразования и вместо этого смещать осевой меридиан для устранения угла разворота.
$ 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>
+
=== Создание МСК в MapInfo ===
Также можно использовать утилиту '''pr''' (после '''-s\''' два пробела!):
+
  
<syntaxhighlight lang="bash">
+
Запись базовой СК-63 зона C0 в файле '''MAPINFOW.PRJ''' должна выглядеть так:
$ 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
+
"CS63 zone C0", 8, 1001, 7, 21.95, 0.1, 1, 250000, 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
+
 
</pre>
 
</pre>
  
Нетрудно заметить, что координаты первого пункта в переходной проекции совпадают с таковыми в МСК. Для остальных пунктов это не так, что говорит о наличии масштабирования и разворота. Убедимся в этом:
+
Используем вторую группу из файла '''key.txt'''. Впишем МСК в '''MAPINFOW.PRJ''' как базовую, дополненную аффинным преобразованием:
  
<syntaxhighlight lang="bash">
+
<pre>
$ findkey data1.dat key1.dat var1.dat
+
"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
</syntaxhighlight>
+
</pre>
  
Содержимое '''key1.dat''':
+
Полноценное определение нуждается в параметрах Bounds:
  
 
<pre>
 
<pre>
0.390
+
"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
-1.813
+
1.0000052644
+
0.0013554914
+
1.0000061831
+
+0.07766348
+
 
</pre>
 
</pre>
  
Первые четыре параметра — ''a''<sub>0</sub>, ''b''<sub>0</sub>, ''a''<sub>1</sub>, ''b''<sub>1</sub>, оставшиеся два — масштабный множитель ''m'' и разворот ''θ'' в градусах.
+
=== Создание МСК в QGIS ===
  
== Косая проекция Меркатора ==
+
Возьмём коэффициенты из первой группы в файле '''key.txt''' и создадим МСК в формате WKT как аффинное преобразование на основе проекции "Pulkovo 1942 / CS63 zone C0":
  
Всё готово для конструирования косой проекции Меркатора. Зададим, как для переходной проекции, долготу и широту первого пункта в качестве начальных, его координаты в МСК в качестве «фальшивых». Кроме того, перенесём масштаб конформного преобразования в масштаб проекции, а разворот — в азимут линии минимального масштаба. Вычислим координаты в этой проекции и запишем в файл '''pt_omerc.dat''':
+
<pre>
 +
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]]]
 +
</pre>
  
<syntaxhighlight lang="bash">
+
В первой строке записали название системы координат латиницей "Biala Podlaska".
$ 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
+
В конце вставили название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате ''φ''<sub>min</sub>, ''λ''<sub>min</sub>, ''φ''<sub>max</sub>, ''λ''<sub>max</sub>.
</syntaxhighlight>
+
  
Вычислим параметры конформного преобразования для '''pt_omerc.dat''': ''m'' = 1.0000000000, ''θ'' = −0.00000022° = −0.0008″. Практически масштаб сведён к единице, разворот к нулю, а вычисленные из значений долготы/широты координаты совпадают с координатами, заданными в МСК. Таким образом, проекция, эквивалентная МСК, построена.
+
=== Создание МСК в Global Mapper ===
  
Следует заметить, что многие программы (например, MapInfo) используют косую проекцию Меркатора в версии Хотина, когда началом координат служит точка пересечения линии минимального масштаба с экватором. Азимут линии и плоские координаты в этой точке, естественно, отличаются от заданных нами, и понадобятся дополнительные вычисления для их определения.
+
Третья группа параметров в файле '''key.txt''' содержит масштабный коэффициент scale и угол поворота angle. Угол вставляем как есть, а масштабный коэффициент умножим на соответствующий параметр базовой СК.
  
Реализация косой проекции Меркатора в PROJ.4 имеет одну особенность: программы просто не работают со значениями параметра '''alpha''' 0 или 90°. Если угол разворота очень мал (скажем, несколько угловых секунд), следует использовать проекцию Гаусса-Крюгера.
+
== Заключение ==
  
== Проекция Гаусса-Крюгера ==
+
Задачи построения МСК для QGIS и MapInfo выполнены, цель достигнута.
  
Встречаются программы, набор проекций в которых ограничен. Проекция Гаусса-Крюгера — она же Transverse Mercator — найдётся в самом бедном наборе. В приёмниках GARMIN пользовательская система координат может быть задана только в этой проекции.
+
== Приложение. Утилита findkey ==
  
=== Долгота осевого меридиана ===
+
Программа '''findkey''' вычисляет параметры конформного преобразования. Написана на языке C. Вот листинг:
  
Если угол разворота не превышает первых десятков угловых минут, можно компенсировать его переносом осевого меридиана в сторону и не бояться потери геодезической точности. Даже первые градусы не должны существенно повлиять на навигационные приложения. Дальнейшее увеличение угла приводит к стремительному росту градиента масштаба в направлении запад-восток и, следовательно, искажению подобия фигур.
+
<syntaxhighlight lang="c">
 +
#include <stdio.h>
 +
#include <stdlib.h>
 +
#include <math.h>
  
В нашем примере разворот оказался небольшим: ''θ'' = 0.07766348° = 0°04′39.59″. Значит, в данном случае на основе проекции Гаусса-Крюгера можно создать довольно точный аналог МСК.
+
#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;
  
<math>\Delta L = -\operatorname{arc\,tg}\frac{\operatorname{tg} \theta}{\sin B}</math>
+
  if (argc < 3) {
 +
    printf("Usage: findkey <coord1> <coord2>\n");
 +
    exit(EXIT_FAILURE);
 +
  }
  
где ∆''L'' — перенос осевого меридиана, ''θ'' — угол разворота, ''B'' — широта первого пункта. Получаем ∆''L'' = −0.0949292655, ''L'' = 37.0888204345. Построим проекцию Гаусса-Крюгера на основе переходной, подставив ''L'' в долготу осевого меридиана и нули в «фальшивые» координаты:
+
  if ((fp0 = fopen(argv[1], "r")) == NULL) {
 +
    printf("can't open %s\n", argv[1]);
 +
    exit(EXIT_FAILURE);
 +
  }
  
<syntaxhighlight lang="bash">
+
  if ((fp1 = fopen(argv[2], "r")) == NULL) {
$ 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
+
    printf("can't open %s\n", argv[2]);
</syntaxhighlight>
+
    exit(EXIT_FAILURE);
 +
  }
  
Последующее вычисление параметров конформного преобразования выдаёт сравнительно большой угол разворота. Однако уже после второй итерации с ''L'' = 37.0888079 получаем ''m'' = 1.0000057306, ''θ'' = +0.00000044.
+
  /* 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];
 +
  }
  
Перенесём ''m'' в масштаб проекции:
+
  /* 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);
  
<syntaxhighlight lang="bash">
+
  /* Helmert parameters */
+proj=tmerc +lat_0=54.8969611 +lon_0=37.0888079 +k=1.00000573 +x_0=0 +y_0=0 +ellps=krass
+
  det = s[0] + s[2];
</syntaxhighlight>
+
  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];
  
Вновь вычислим параметры конформного преобразования. Теперь ''m'' = 1.0000000039, т.е. практически единица.
+
  /* 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);
  
Наконец обратим внимание на первичные параметры ''a''<sub>0</sub> = −4756.693 и ''a''<sub>1</sub> = 281.807. Перенесём их значения в параметры проекции '''x_0''' и '''y_0''':
+
  /* 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);
  
<syntaxhighlight lang="bash">
+
  exit(EXIT_SUCCESS);
+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 могут загрузить уже скомпилированную программу.

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

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

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