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

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Строка 72: Строка 72:
== Вычисление параметров ==
== Вычисление параметров ==


Для вычислений используем утилиту '''findkey''', о которой сказано ниже в приложении. Она запускается в командной строке вот так:
Для вычислений используем консольную утилиту '''findkey''', о которой сказано ниже в приложении.
 
<pre></pre>


=== Определение базовой ГСК ===
=== Определение базовой ГСК ===


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


Удобно, если осевой меридиан ГСК проходит через город при малом угле разворота, как в Пензе. В этом случае МСК строится как проекция Гаусса-Крюгера, отличающаяся от исходной ГСК малым сдвигом осевого меридиана для компенсации угла разворота. Жаль, что это случается редко.
<syntaxhighlight lang="bash">
 
$ findkey cat_s42z4.tsv cat_local.tsv
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
</syntaxhighlight>


Хорошо, если исходная ГСК точно известна. Это, в частности, случай, когда вместо каталогов имеется надёжный ключ.
Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками.
 
Посмотрим на содержимое '''var.csv''':
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:


{| class="wikitable"
{| class="wikitable"
Строка 106: Строка 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"
Строка 139: Строка 129:
|}
|}


Сравнение невязок не в пользу СК-42. Исходная ГСК — СК-63 зона C0, и правильная проекция для QGIS такая:
Сравнение невязок позволяет прийти к выводу, что базовая ГСК — СК-63 зона C0.
Изучим содержимое соответствующего файла '''key.txt''':


<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:
</pre>
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
* азимут начальной линии равен нулю.
 
Довольно хорошее приближение для широты центра проекции даёт формула
 
: tg ''φ''₂ = tg ''φ''₁ / cos (''λ''₂ − ''λ''₁) ,
 
где ''φ''₁, ''λ''₁ — координаты середины города; , ''φ''₂, ''λ''₂ — координаты центра проекции. Точное вычисление выполняется на апосфере.
 
На практике идеал разбивается о нежелание авторов PROJ разрешить нулевой азимут начальной линии. Приходится использовать какое-нибудь малое его значение, но это приводит к смещению точки центра проекции, больше по широте, меньше по долготе.
 
Впрочем, тема вычисления оптимального положения центра проекции выходит за рамки данной статьи. Между тем предложенный выше подход, при котором долгота равна долготе осевого меридиана, а широта равна широте середины территории, даёт вполне удовлетворительные результаты.
 
=== Проверка параметров ===
 
Прежде чем копировать полученную строку параметров в пользовательскую проекцию QGIS, полезно её проверить:
 
<syntaxhighlight lang="bash">
$ 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>
 
В выводе программы должна быть практическая единица для масштаба и нули для остальных параметров:
 
<pre>
+k=0.99999999999999 +x_0=5.956435757244818e-10 +y_0=2.582964953035116e-10 +gamma=1.139084427878181e-13
</pre>
</pre>


Координаты в файле '''cat.tsv''' будут похожи на исходные координаты МСК из файла '''cat_local.tsv'''.
Мы видим три группы чисел: WKT, MapInfo и Alternate set.


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

Версия от 20:11, 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:

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.

Заключение

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

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

Приложение. Утилита 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 */

  /* 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 могут загрузить уже скомпилированную программу.

Примечания


Ссылки