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

Материал из GIS-Lab
(Различия между версиями)
Перейти к: навигация, поиск
м (Быстрая реализация)
(Ссылки)
 
(не показаны 60 промежуточных версий 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>
+
В нашем примере мы предполагаем, что базовая ГСК либо СК-42 зона 4, либо СК-63 зона C0.
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>
+
  
== Быстрая реализация ==
+
Параметры СК-63 зона C0 известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0".
 
+
Создадим каталог в СК-63 с помощью утилиты '''proj''':
Работа заключается в задании трёх параметров косой проекции Меркатора и вычислении оставшихся четырёх. Априорные параметры определяются так: центр проекции должен лежать на осевом меридиане исходной проекции СК-42/СК-63 на широте середины города, а азимут начальной линии должен быть малым. Пусть широта центра проекции равна широте первого пункта из файла '''cat_longlat.tsv''' 52.02642240080064°. Долготу центра проекции приравняем долготе осевого меридиана четвёртой зоны 21°. Азимут начальной линии приравняем произовольно малой величине −0,0001°. Масштабный множитель ''k'' приравняем единице, оставшиеся параметры обнулим, вычислим координаты в этой переходной проекции и запишем их в файл '''cat.tsv''':
+
  
 
<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>
 
+
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
+
  
Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов имеется надёжный ключ.
+
Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками.
 
+
Посмотрим на содержимое '''var.csv''':
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:
+
  
 
{| class="wikitable"
 
{| class="wikitable"
Строка 120: Строка 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
Строка 153: Строка 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>
  
Возможно, в ближайшем будущем среди атомарных операций PROJ появится аффинное преобразование. Это позволит непосредственно дополнять классическую проекцию Гаусса-Крюгера конформным преобразованием, как это делается в MapInfo Pro.
+
=== Создание МСК в QGIS ===
  
== Утилита конформного преобразования helmkey ==
+
Возьмём коэффициенты из первой группы в файле '''key.txt''' и создадим МСК в формате WKT как аффинное преобразование на основе проекции "Pulkovo 1942 / CS63 zone C0":
  
Вот листинг программы на языке C:
+
<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>
 +
 
 +
В первой строке записали название системы координат латиницей "Biala Podlaska".
 +
В конце вставили название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате ''φ''<sub>min</sub>, ''λ''<sub>min</sub>, ''φ''<sub>max</sub>, ''λ''<sub>max</sub>.
 +
 
 +
=== Создание МСК в Global Mapper ===
 +
 
 +
Третья группа параметров в файле '''key.txt''' содержит масштабный коэффициент scale и угол поворота angle. Угол вставляем как есть, а масштабный коэффициент умножим на соответствующий параметр базовой СК.
 +
 
 +
== Заключение ==
 +
 
 +
Задачи построения МСК для QGIS и MapInfo выполнены, цель достигнута.
 +
 
 +
== Приложение. Утилита findkey ==
 +
 
 +
Программа '''findkey''' вычисляет параметры конформного преобразования. Написана на языке C. Вот листинг:
  
 
<syntaxhighlight lang="c">
 
<syntaxhighlight lang="c">
Строка 194: Строка 266:
  
 
/* --------------------------------------------------------------------------
 
/* --------------------------------------------------------------------------
  * helmkey
+
  * findkey
 
  *
 
  *
  * Program to compute for Helmert & affine planar transformation parameters
+
  * Program to compute Helmert 2D transformation parameters
 
  *
 
  *
  * Usage: helmkey <coord1> <coord2>
+
  * Usage: findkey <coord1> <coord2>
 
  *
 
  *
 
  * Input files: coord1 coord2
 
  * Input files: coord1 coord2
Строка 205: Строка 277:
 
  *              a row per a point; 3+ points
 
  *              a row per a point; 3+ points
 
  *
 
  *
  * Output: omerc parameter string
+
  * Output files:
  *    k    - Scale factor on initial line
+
  *    key.txt - transformation parameters
*    gamma - Angle from Rectified to Skew Grid
+
  *    var.csv - SEP separated residuals 'dx dy'
*    x_0  - Easting at projection center
+
*    y_0  - Northing at projection center
+
*
+
* Output file:
+
  *    var.csv - output SEP separated variances 'dx dy'
+
 
  * -------------------------------------------------------------------------- */
 
  * -------------------------------------------------------------------------- */
 
int main(int argc, char *argv[])
 
int main(int argc, char *argv[])
Строка 228: Строка 295:
  
 
   if (argc < 3) {
 
   if (argc < 3) {
     printf("Usage: helmkey <coord1> <coord2>\n");
+
     printf("Usage: findkey <coord1> <coord2>\n");
 
     exit(EXIT_FAILURE);
 
     exit(EXIT_FAILURE);
 
   }
 
   }
Строка 295: Строка 362:
 
   /* alternative Helmert parameter set */
 
   /* alternative Helmert parameter set */
 
   mu = hypot(h[0], h[1]);
 
   mu = hypot(h[0], h[1]);
   theta = atan2(h[1], h[0]);
+
   theta = atan2(h[1], h[0]) / M_PI * 180.;
  
 
   /* output parameters */
 
   /* output parameters */
   printf("+k=%.16g +x_0=%.16g +y_0=%.16g +gamma=%.16g\n",
+
   if ((fp2 = fopen("key.txt", "w")) == NULL) {
mu, h[2], h[5], theta / M_PI * 180.);
+
    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 */
 
   /* output residuals */
Строка 312: Строка 387:
 
     yh[0] = h[0] * x[0] + h[1] * x[1] + h[2];
 
     yh[0] = h[0] * x[0] + h[1] * x[1] + h[2];
 
     yh[1] = h[3] * x[0] + h[4] * x[1] + h[5];
 
     yh[1] = h[3] * x[0] + h[4] * x[1] + h[5];
     fprintf(fp2, "%.3f%c%.3f\n", yh[0] - y[0], SEP, yh[1] - y[1]);
+
     fprintf(fp2, "%.15g%c%.15g\n", yh[0] - y[0], SEP, yh[1] - y[1]);
 
   }
 
   }
 
   fclose(fp2);
 
   fclose(fp2);
Строка 318: Строка 393:
 
   fclose(fp0);
 
   fclose(fp0);
  
   return 0;
+
   exit(EXIT_SUCCESS);
 
}
 
}
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Сохраним код в файл '''helmkey.c'''. Исполняемый модуль создадим компилятором '''gcc''':
+
Сохраним код в файл '''findkey.c'''. Создадим исполняемый модуль компилятором '''gcc''':
  
 
<syntaxhighlight lang="bash">
 
<syntaxhighlight lang="bash">
$ gcc -o helmkey helmkey.c -lm
+
$ gcc -o findkey findkey.c -lm
 
</syntaxhighlight>
 
</syntaxhighlight>
  
Пользователи MS Windows могут загрузить уже скомпилированную [http://gis-lab.info/other/helmkey.zip программу].
+
Пользователи MS Windows могут загрузить уже скомпилированную [https://wiki.gis-lab.info/images/4/43/Findkey.zip программу].
 
+
== Примечания ==
+
 
+
<references />
+
  
 
== Ссылки ==
 
== Ссылки ==
Строка 341: Строка 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 могут загрузить уже скомпилированную программу.

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

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

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