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

Материал из GIS-Lab
(Различия между версиями)
Перейти к: навигация, поиск
(Быстрая реализация)
(Ссылки)
 
(не показана 21 промежуточная версия 1 участника)
Строка 50: Строка 50:
 
Каталог в первой системе имеется, нужно подготовить каталог в альтернативной системе.
 
Каталог в первой системе имеется, нужно подготовить каталог в альтернативной системе.
  
Параметры СК-63 зона C0 нам известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0".  
+
Параметры СК-63 зона C0 известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0".  
 
Создадим каталог в СК-63 с помощью утилиты '''proj''':
 
Создадим каталог в СК-63 с помощью утилиты '''proj''':
  
Строка 70: Строка 70:
 
</pre>
 
</pre>
  
== Много букв ==
+
== Создание МСК ==
  
=== Выбор исходной проекции ГСК ===
+
Для вычислений используем консольную утилиту '''findkey''', о которой сказано ниже в приложении.
  
Важной особенностью МСК является то, что с точки зрения математической картографии она остаётся проекцией Гаусса-Крюгера и наследует её искажения.
+
=== Определение базовой ГСК ===
  
Удобно, если осевой меридиан ГСК проходит через город при малом угле разворота, как в Пензе. В этом случае МСК строится как проекция Гаусса-Крюгера, отличающаяся от исходной ГСК малым сдвигом осевого меридиана для компенсации угла разворота. Жаль, что это случается редко.
+
Вычислим параметры конформного преобразования координат из СК-42 в МСК. Для этого в командной строке запустим программу '''findkey''' с аргументами '''cat_s42z4.tsv''' и '''cat_local.tsv''':
  
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
+
<syntaxhighlight lang="bash">
 +
$ findkey cat_s42z4.tsv cat_local.tsv
 +
</syntaxhighlight>
  
Хорошо, если исходная ГСК точно известна. Это, в частности, случай, когда вместо каталогов имеется надёжный ключ.
+
Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками.
 
+
Посмотрим на содержимое '''var.csv''':
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:
+
  
 
{| class="wikitable"
 
{| class="wikitable"
Строка 102: Строка 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"
Строка 135: Строка 129:
 
|}
 
|}
  
Сравнение невязок не в пользу СК-42. Исходная ГСК — СК-63 зона C0, и правильная проекция для QGIS такая:
+
Сравнение невязок позволяет сделать вывод, что базовая ГСК — СК-63 зона C0.
 +
 
 +
=== Полученные параметры ===
 +
 
 +
Изучим содержимое файла '''key.txt''', соответствующего СК-63:
  
 
<pre>
 
<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
+
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
 
</pre>
 
</pre>
  
=== Подбор априорных параметров МСК ===
+
Мы видим три группы чисел: WKT, MapInfo и Alternate set.
  
В плане распределения искажений проекция МСК должна имитировать ГСК. В идеале выбор параметров таков, что центр проекции лежит на точке осевого меридана ГСК, ближайшей к середине города, а направление начальной линии совпадает с направлением меридиана. Итак,
+
Обратим внимание на параметр разворота angle из третьей группы. Если он мал, в пределах первых десятых долей градуса, имеет смысл отказаться от использования конформного преобразования и вместо этого смещать осевой меридиан для устранения угла разворота.
* долгота центра проекции равна долготе осевого меридиана ГСК;
+
* широта центра проекции требует вычисления;
+
* азимут начальной линии равен нулю.
+
  
Довольно хорошее приближение для широты центра проекции даёт формула
+
=== Создание МСК в MapInfo ===
  
: tg ''φ''₂ = tg ''φ''₁ / cos (''λ''₂ − ''λ''₁) ,
+
Запись базовой СК-63 зона C0 в файле '''MAPINFOW.PRJ''' должна выглядеть так:
  
где ''φ''₁, ''λ''₁ — координаты середины города; , ''φ''₂, ''λ''₂ — координаты центра проекции. Точное вычисление выполняется на апосфере.
+
<pre>
 +
"CS63 zone C0", 8, 1001, 7, 21.95, 0.1, 1, 250000, 0
 +
</pre>
  
На практике идеал разбивается о нежелание авторов PROJ разрешить нулевой азимут начальной линии. Приходится использовать какое-нибудь малое его значение, но это приводит к смещению точки центра проекции, больше по широте, меньше по долготе.
+
Используем вторую группу из файла '''key.txt'''. Впишем МСК в '''MAPINFOW.PRJ''' как базовую, дополненную аффинным преобразованием:
  
Впрочем, тема вычисления оптимального положения центра проекции выходит за рамки данной статьи. Между тем предложенный выше подход, при котором долгота равна долготе осевого меридиана, а широта равна широте середины территории, даёт вполне удовлетворительные результаты.
+
<pre>
 +
"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
 +
</pre>
  
=== Проверка параметров ===
+
Полноценное определение нуждается в параметрах Bounds:
  
Прежде чем копировать полученную строку параметров в пользовательскую проекцию QGIS, полезно её проверить:
+
<pre>
 +
"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
 +
</pre>
  
<syntaxhighlight lang="bash">
+
=== Создание МСК в QGIS ===
$ proj -f "%.16g" +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 cat_longlat.tsv > cat.tsv
+
$ helmkey cat.tsv cat_local.tsv
+
</syntaxhighlight>
+
  
В выводе программы должна быть практическая единица для масштаба и нули для остальных параметров:
+
Возьмём коэффициенты из первой группы в файле '''key.txt''' и создадим МСК в формате WKT как аффинное преобразование на основе проекции "Pulkovo 1942 / CS63 zone C0":
  
 
<pre>
 
<pre>
+k=0.99999999999999 +x_0=5.956435757244818e-10 +y_0=2.582964953035116e-10 +gamma=1.139084427878181e-13
+
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>
 
</pre>
  
Координаты в файле '''cat.tsv''' будут похожи на исходные координаты МСК из файла '''cat_local.tsv'''.
+
В первой строке записали название системы координат латиницей "Biala Podlaska".
 +
В конце вставили название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате ''φ''<sub>min</sub>, ''λ''<sub>min</sub>, ''φ''<sub>max</sub>, ''λ''<sub>max</sub>.
  
== Заключение ==
+
=== Создание МСК в Global Mapper ===
  
На момент написания статьи для практического применения в QGIS рекомендуется косая проекция Меркатора. Если осевой меридиан близок к центру территории и угол разворота мал, следует использовать проекцию Гаусса-Крюгера.
+
Третья группа параметров в файле '''key.txt''' содержит масштабный коэффициент scale и угол поворота angle. Угол вставляем как есть, а масштабный коэффициент умножим на соответствующий параметр базовой СК.
 +
 
 +
== Заключение ==
  
Возможно, в ближайшем будущем среди атомарных операций PROJ появится аффинное преобразование. Это позволит непосредственно дополнять классическую проекцию Гаусса-Крюгера дополнительным преобразованием прямоугольной системы координат, как это делается в MapInfo Pro.
+
Задачи построения МСК для QGIS и MapInfo выполнены, цель достигнута.
  
 
== Приложение. Утилита findkey ==
 
== Приложение. Утилита findkey ==
Строка 290: Строка 363:
 
   mu = hypot(h[0], h[1]);
 
   mu = hypot(h[0], h[1]);
 
   theta = atan2(h[1], h[0]) / M_PI * 180.;
 
   theta = atan2(h[1], h[0]) / M_PI * 180.;
 
  /* output parameters */
 
  
 
   /* output parameters */
 
   /* output parameters */
Строка 332: Строка 403:
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Пользователи MS Windows могут загрузить уже скомпилированную [http://wiki.gis-lab.info/images/4/49/Helmkey.zip программу].
+
Пользователи MS Windows могут загрузить уже скомпилированную [https://wiki.gis-lab.info/images/4/43/Findkey.zip программу].
 
+
== Примечания ==
+
 
+
<references />
+
  
 
== Ссылки ==
 
== Ссылки ==
Строка 345: Строка 412:
 
* [https://proj4.org/usage/index.html Using PROJ]
 
* [https://proj4.org/usage/index.html Using PROJ]
 
* [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 могут загрузить уже скомпилированную программу.

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

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

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