Добавление местной координатной системы в GIS
по адресу 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 могут загрузить уже скомпилированную программу.