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

Материал из GIS-Lab
(Различия между версиями)
Перейти к: навигация, поиск
м
м (Выбор исходной проекции ГСК)
Строка 106: Строка 106:
 
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
 
В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.
  
Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов координат имеется надёжный ключ в виде строки MapInfo.
+
Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов имеется надёжный ключ в виде строки MapInfo.
  
 
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:
 
Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита '''helmkey'''. Заглянем в файл '''var.csv''', созданный этой утилитой. Он содержит невязки конформного преобразования:

Версия 10:49, 20 мая 2018

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/local-cs.html


Конструирование проекций, имитирующих местные координатные системы, в QGIS

Содержание

Введение

Под местной системой координат (МСК) будет подразумеваться так называемая «городская» система, построенная независимо от государственной системы (ГСК) и включенная в неё заданием ключей перехода к СК-42 или СК-63. МСК крупных территорий, сравнимых с размерами субъектов Федерации, не являются предметом данной статьи, поскольку относятся к классическим картографическим проекциям.

Постановка задачи

Имеется множество пунктов, для которых известны координаты X, Y в ГСК и x, y в МСК. Требуется подобрать проекцию и вычислить её параметры, удовлетворительно представляющую МСК в ГИС. При подборе параметров предполагается использовать один из пунктов в качестве центральной точки преобразования.

Некоторые программы позволяют реализовать работу в МСК непосредственно. Так, в MapInfo любая проекция может быть дополнена аффинным преобразованием, а в Global Mapper конформным. ArcGIS в качестве МСК предлагает локальную проекцию: аналог ортографической проекции на эллипсоиде, дополненный разворотом и масштабированием.

Другие программы, включая QGIS, работают только с «обыкновенными» проекциями, не допуская дополнительных геометрических преобразований. Задача статьи — показать, как можно конструировать проекции, позволяющие имитировать МСК в QGIS.

В качестве рабочей среды будем использовать командную строку UNIX или OSGeo4W.

Подготовка данных

Нужны два файла исходных данных. Пусть текстовый файл 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

Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы.

Вычислим географические координаты пунктов:

$ proj -I -f "%.16g" +init=epsg:28404 cat_s42z4.tsv > cat_longlat.tsv

В файл cat_longlat.tsv запишутся географические долготы и широты в СК-42:

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

Быстрая реализация

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

$ 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

Найдём оставшиеся параметры с помощью программы helmkey:

$ helmkey cat.tsv cat_local.tsv

Программа выведет на экран искомую строку параметров:

+k=0.9998372145697554 +x_0=-78707.08698765696 +y_0=32225.08703617829 +gamma=1.677027577390829

Это всё! Осталось добавить в QGIS пользовательскую проекцию, не забыв необходимые стандартные параметры:

+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

Много букв

Выбор исходной проекции ГСК

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

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

В остальных случаях следует использовать косую проекцию Меркатора, так как она, во-первых, при правильном подборе параметров близка к проекции Гаусса-Крюгера в окрестности заданного центра проекции и, во-вторых, в PROJ.4 позволяет произвольно задавать угол разворота МСК.

Хорошо, если исходная ГСК точно известна. Это, например, случай, когда вместо каталогов имеется надёжный ключ в виде строки MapInfo.

Однако встречаются ситуации, когда исходная ГСК неизвестна, и приходится выбирать между СК-42 и СК-63, или даже между одной зоной одной из них и двумя зонами другой. Ответ подскажут свойства самих проекций и утилита helmkey. Заглянем в файл 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 зона C0 с осевым меридианом 21°57′. Создадим промежуточную проекцию на её основе, вычислим координаты для пунктов в ней и получим остающиеся параметры утилитой helmkey:

$ 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
$ helmkey cat.tsv cat_local.tsv

Вывод программы:

+k=1.00001854602098 +x_0=-13532.31228709544 +y_0=30742.75844830525 +gamma=0.9279005081790233

Содержимое 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

Сравнение невязок не в пользу СК-42. Исходная ГСК — СК-63 зона C0, и правильная проекция для QGIS такая:

+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

О выборе априорных параметров МСК

В плане распределения искажений проекция МСК должна имитировать ГСК. В идеале выбор параметров таков, что центр проекции лежит на точке осевого меридана ГСК, ближайшей к середине города, а направление начальной линии совпадает с направлением меридиана. Итак,

  • долгота центра проекции равна долготе осевого меридиана ГСК;
  • широта центра проекции требует вычисления;
  • азимут начальной линии равен нулю.

Довольно хорошее приближение для широты центра проекции даёт формула

tg φ₂ = tg φ₁ / cos (λ₂ − λ₁) ,

где φ₁, λ₁ — координаты середины города; , φ₂, λ₂ — координаты центра проекции. Точное вычисление выполняется на апосфере.

На практике идеал разбивается о нежелание авторов PROJ разрешить нулевой азимут начальной линии. Приходится использовать какое-нибудь малое его значение, но это приводит к смещению точки центра проекции, больше по широте, меньше по долготе.

Впрочем, тема вычисления оптимального положения центра проекции выходит за рамки данной статьи. Между тем предложенный выше подход, при котором долгота равна долготе осевого меридиана, а широта равна широте середины территории, даёт вполне удовлетворительные результаты.

Заключение

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

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

Утилита конформного преобразования helmkey

Вот листинг программы на языке C:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define SEP ';' /* var-file column separator */
 
/* --------------------------------------------------------------------------
 * helmkey
 *
 * Program to compute for Helmert & affine planar transformation parameters
 *
 * Usage: helmkey <coord1> <coord2>
 *
 * Input files: coord1 coord2
 *     coord1 - source coordinate 'x1 y1'
 *     coord2 - destination coordinate 'x2 y2'
 *              a row per a point; 3+ points
 *
 * Output: omerc parameter string
 *    k     - Scale factor on initial line
 *    gamma - Angle from Rectified to Skew Grid
 *    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[])
{
  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: helmkey <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]);
 
  /* output parameters */
  printf("+k=%.16g +x_0=%.16g +y_0=%.16g +gamma=%.16g\n",
	 mu, h[2], h[5], theta / M_PI * 180.);
 
  /* 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, "%.3f%c%.3f\n", yh[0] - y[0], SEP, yh[1] - y[1]);
  }
  fclose(fp2);
  fclose(fp1);
  fclose(fp0);
 
  return 0;
}

Сохраним код в файл helmkey.c. Исполняемый модуль создадим компилятором gcc:

$ gcc -o helmkey helmkey.c -lm

Пользователи MS Windows могут загрузить уже скомпилированную программу.

Примечания

  1. При работе в MSYS под Windows придётся позаботиться об удалении мусорных символов:
    $ 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

Ссылки

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

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