Добавление местной координатной системы в GIS: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
 
(не показано 38 промежуточных версий этого же участника)
Строка 7: Строка 7:
Под местной системой координат (МСК) будет подразумеваться так называемая «городская» система, построенная независимо от государственной системы (ГСК) и включенная в неё [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 и Global Mapper.
Многие программы ГИС по примеру геодезических программ позволяют реализовать работу в МСК непосредственно. Так, в QGIS и в MapInfo Pro любая проекция может быть дополнена аффинным преобразованием, а в Global Mapper конформные проекции дополняются разворотом. В данной статье рассматривается создание МСК в программах QGIS и MapInfo.


== Постановка задачи ==
== Постановка задачи ==
Строка 43: Строка 43:
Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы.
Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы.


=== Выбор базовой ГСК ===
=== Дополнительные данные ===


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


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


Хорошо, если исходная ГСК точно известна. Однако очень часто исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита helmkey.
Параметры СК-63 зона C0 известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0".  
 
Создадим каталог в СК-63 с помощью утилиты '''proj''':
Вычислим географические координаты пунктов:<ref>
 
При запуске '''proj''' в MSYS под Windows придётся позаботиться об удалении мусорных символов:
 
<syntaxhighlight lang="bash">
$ proj -f "%.16g" +proj=omerc +lat_0=52.02642240080064 +lonc=21 +alpha=-0.0001 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=krass cat_longlat.tsv | tr -d '\r' > cat.tsv
</syntaxhighlight>
 
</ref>


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
$ proj -I -f "%.16g" +init=epsg:28404 cat_s42z4.tsv > cat_longlat.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>


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


<pre>
<pre>
23.12716113273887 52.02642240080064
330797.45370592922 5755981.4751984337
23.35199369039325 52.03605540518488
346208.04426388565 5757327.0953416284
23.19280603770445 52.15834802315441
335051.73824425979 5770735.1441946141
22.96008928353681 52.12307568371703
319180.795365442  5766563.182877643
22.91428747596251 51.99456156589881
316233.72446517192 5752221.1970210578
23.06504376727423 51.90277894131356
326744.90098082455 5742157.2318198672
23.26700321114454 51.92327827668979
340603.31654394628 5744670.1910534762
</pre>
</pre>


== Быстрая реализация ==
== Создание МСК ==


Работа заключается в задании трёх параметров косой проекции Меркатора и вычислении оставшихся четырёх. Априорные параметры определяются так: центр проекции должен лежать на осевом меридиане исходной проекции СК-42/СК-63 на широте середины города, а азимут начальной линии должен быть малым. Пусть широта центра проекции равна широте первого пункта из файла '''cat_longlat.tsv''' 52.02642240080064°. Долгота центра проекции будет совпадать с долготой осевого меридиана четвёртой зоны 21°. Азимут начальной линии положим равным произовольно малой величине −0,0001°. Масштабный множитель ''k'' приравняем единице, оставшиеся параметры обнулим, вычислим координаты в этой переходной проекции и запишем их в файл '''cat.tsv''':
Для вычислений используем консольную утилиту '''findkey''', о которой сказано ниже в приложении.


<syntaxhighlight lang="bash">
=== Определение базовой ГСК ===
$ proj -f "%.16g" +proj=omerc +lat_0=52.02642240080064 +lonc=21 +alpha=-0.0001 +k=1 +x_0=0 +y_0=0 +gamma=0 +ellps=krass cat_longlat.tsv > cat.tsv
</syntaxhighlight>


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


<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
$ helmkey cat.tsv cat_local.tsv
$ findkey cat_s42z4.tsv cat_local.tsv
</syntaxhighlight>
</syntaxhighlight>


Программа выведет на экран искомую строку параметров:
Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками.
 
Посмотрим на содержимое '''var.csv''':
<pre>
+k=0.9998372145697554 +x_0=-78707.08698765696 +y_0=32225.08703617829 +gamma=1.677027577390829
</pre>
 
Это всё! Осталось добавить в QGIS пользовательскую проекцию, не забыв необходимые стандартные параметры:
 
<pre>
+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>
 
== Много букв ==
 
=== Выбор исходной проекции ГСК ===
 
Важной особенностью МСК является то, что с точки зрения математической картографии она остаётся проекцией Гаусса-Крюгера и наследует её искажения.
 
Удобно, если осевой меридиан ГСК проходит через город при малом угле разворота, как в Пензе. В этом случае МСК строится как проекция Гаусса-Крюгера, отличающаяся от исходной ГСК малым сдвигом осевого меридиана для компенсации угла разворота. Жаль, что это случается редко.
 
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
 
Хорошо, если исходная ГСК точно известна. Это, в частности, случай, когда вместо каталогов имеется надёжный ключ.
 
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:


{| class="wikitable"
{| class="wikitable"
Строка 135: Строка 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
Строка 168: Строка 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]]],
        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["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Europe - Poland - Biala Podlaska"],
        BBOX[51.9,22.9,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 ===
 
Третья группа параметров в файле '''key.txt''' содержит масштабный коэффициент scale и угол поворота angle. Угол вставляем как есть, а масштабный коэффициент умножим на соответствующий параметр базовой СК.


== Заключение ==
== Заключение ==


На момент написания статьи для практического применения в QGIS рекомендуется косая проекция Меркатора. Если осевой меридиан близок к центру территории и угол разворота мал, следует использовать проекцию Гаусса-Крюгера.
Задачи построения МСК для QGIS и MapInfo выполнены, цель достигнута.


Возможно, в ближайшем будущем среди атомарных операций PROJ появится аффинное преобразование. Это позволит непосредственно дополнять классическую проекцию Гаусса-Крюгера дополнительным преобразованием прямоугольной системы координат, как это делается в MapInfo Pro.
== Приложение. Утилита findkey ==


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


<syntaxhighlight lang="c">
<syntaxhighlight lang="c">
Строка 226: Строка 264:


/* --------------------------------------------------------------------------
/* --------------------------------------------------------------------------
  * 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
Строка 237: Строка 275:
  *              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
*    x_0  - Easting at projection center
  *    var.csv - SEP separated residuals 'dx dy'
*    y_0  - Northing at projection center
*    gamma - Angle from Rectified to Skew Grid
*
* Output file:
  *    var.csv - output SEP separated variances 'dx dy'
  * -------------------------------------------------------------------------- */
  * -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
int main(int argc, char *argv[])
Строка 260: Строка 293:


   if (argc < 3) {
   if (argc < 3) {
     printf("Usage: helmkey <coord1> <coord2>\n");
     printf("Usage: findkey <coord1> <coord2>\n");
     exit(EXIT_FAILURE);
     exit(EXIT_FAILURE);
   }
   }
Строка 327: Строка 360:
   /* 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 */
Строка 344: Строка 385:
     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);
Строка 354: Строка 395:
</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://wiki.gis-lab.info/images/4/49/Helmkey.zip программу].
Пользователи MS Windows могут загрузить уже скомпилированную [https://wiki.gis-lab.info/images/4/43/Findkey.zip программу].
 
== Примечания ==
 
<references />


== Ссылки ==
== Ссылки ==
Строка 373: Строка 410:
* [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]

Текущая версия от 16:56, 3 сентября 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]]],
        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["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Europe - Poland - Biala Podlaska"],
        BBOX[51.9,22.9,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 могут загрузить уже скомпилированную программу.

Ссылки