Добавление местной координатной системы в GIS
м (→Быстрая реализация) |
(→Создание МСК в QGIS) |
||
(не показаны 78 промежуточных версий 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'' в МСК. Требуется подобрать проекцию, удовлетворительно представляющую МСК в ГИС. При подборе параметров предполагается использовать один из пунктов в качестве центральной точки преобразования. |
− | + | == Подготовка данных == | |
− | + | === Исходные данные === | |
− | + | ||
− | + | ||
− | + | ||
− | == | + | |
− | + | Имеются два каталога. Текстовый файл '''cat_s42z4.tsv''' содержит координаты пунктов в государственной системе (ГСК), а именно в четвёртой зоне СК-42: | |
<pre> | <pre> | ||
Строка 45: | Строка 43: | ||
Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы. | Каждая строка в обоих файлах соответствует одному и тому же пункту. В первой строке центральный пункт системы. | ||
− | + | === Дополнительные данные === | |
− | + | ||
− | + | ||
− | + | ||
− | + | Очень важно помнить, что с точки зрения математической картографии МСК остаётся проекцией Гаусса-Крюгера и наследует её искажения. Поэтому важно знать, на какой именно ГСК основана МСК. Зачастую это заранее неизвестно, и приходится проводить предварительное исследование для выяснения этого вопроса. | |
− | + | В нашем примере мы предполагаем, что базовая ГСК либо СК-42 зона 4, либо СК-63 зона C0. | |
− | + | Каталог в первой системе имеется, нужно подготовить каталог в альтернативной системе. | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | Параметры СК-63 зона C0 известны, это EPSG:3350 "Pulkovo 1942 / CS63 zone C0". | |
+ | Создадим каталог в СК-63 с помощью утилиты '''proj''': | ||
<syntaxhighlight lang="bash"> | <syntaxhighlight lang="bash"> | ||
− | $ | + | $ 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> | ||
− | + | В файл '''lonlat.tsv''' запишутся географические координаты (долготы и широты) в СК-42, а в '''cat_s63c0.tsv''' координаты в СК-63 зона C0: | |
<pre> | <pre> | ||
− | + | 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> | ||
− | + | == Создание МСК == | |
− | + | Для вычислений используем консольную утилиту '''findkey''', о которой сказано ниже в приложении. | |
− | + | ||
− | + | ||
− | == | + | === Определение базовой ГСК === |
+ | |||
+ | Вычислим параметры конформного преобразования координат из СК-42 в МСК. Для этого в командной строке запустим программу '''findkey''' с аргументами '''cat_s42z4.tsv''' и '''cat_local.tsv''': | ||
− | |||
<syntaxhighlight lang="bash"> | <syntaxhighlight lang="bash"> | ||
− | $ | + | $ findkey cat_s42z4.tsv cat_local.tsv |
</syntaxhighlight> | </syntaxhighlight> | ||
− | + | Программа создаст два файла: '''key.txt''' с параметрами конформного преобразования и '''var.csv''' с координатными невязками. | |
+ | Посмотрим на содержимое '''var.csv''': | ||
− | + | {| class="wikitable" | |
− | + | |- | |
− | + | |- align="right" | |
− | + | | -0.006 || 0.007 | |
− | - | + | |- align="right" |
− | + | | 0.182 || 0.046 | |
+ | |- align="right" | ||
+ | | -0.166 || 0.110 | ||
+ | |- align="right" | ||
+ | | 0.019 || -0.185 | ||
+ | |- align="right" | ||
+ | | 0.148 || 0.100 | ||
+ | |- align="right" | ||
+ | | -0.146 || 0.094 | ||
+ | |- align="right" | ||
+ | | -0.031 || -0.171 | ||
+ | |} | ||
− | + | Вычислим параметры конформного преобразования координат из СК-63 в МСК: | |
<syntaxhighlight lang="bash"> | <syntaxhighlight lang="bash"> | ||
− | $ | + | $ findkey cat_s63c0.tsv cat_local.tsv |
</syntaxhighlight> | </syntaxhighlight> | ||
− | + | Теперь содержимое '''var.csv''' будет таким: | |
− | + | ||
− | + | {| class="wikitable" | |
− | + | |- | |
− | + | |- align="right" | |
+ | | 0.000 || -0.002 | ||
+ | |- align="right" | ||
+ | | -0.001 || 0.002 | ||
+ | |- align="right" | ||
+ | | -0.001 || 0.002 | ||
+ | |- align="right" | ||
+ | | 0.004 || 0.000 | ||
+ | |- align="right" | ||
+ | | -0.002 || 0.001 | ||
+ | |- align="right" | ||
+ | | 0.002 || -0.002 | ||
+ | |- align="right" | ||
+ | | -0.002 || -0.001 | ||
+ | |} | ||
− | + | Сравнение невязок позволяет сделать вывод, что базовая ГСК — СК-63 зона C0. | |
− | + | === Полученные параметры === | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | Изучим содержимое файла '''key.txt''', соответствующего СК-63: | |
− | < | + | <pre> |
− | + | 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 | |
− | + | ||
− | + | ||
− | 1. | + | |
− | + | ||
</pre> | </pre> | ||
− | + | Мы видим три группы чисел: WKT, MapInfo и Alternate set. | |
− | + | Обратим внимание на параметр разворота angle из третьей группы. Если он мал, в пределах первых десятых долей градуса, имеет смысл отказаться от использования конформного преобразования и вместо этого смещать осевой меридиан для устранения угла разворота. | |
− | + | === Создание МСК в MapInfo === | |
− | + | Запись базовой СК-63 зона C0 в файле '''MAPINFOW.PRJ''' должна выглядеть так: | |
− | + | ||
− | + | ||
− | + | <pre> | |
+ | "CS63 zone C0", 8, 1001, 7, 21.95, 0.1, 1, 250000, 0 | ||
+ | </pre> | ||
− | + | Используем вторую группу из файла '''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: | |
− | + | <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> | ||
− | === | + | === Создание МСК в QGIS === |
− | + | Возьмём коэффициенты из первой группы в файле '''key.txt''' и создадим МСК в формате WKT как аффинное преобразование на основе проекции "Pulkovo 1942 / CS63 zone C0": | |
− | + | <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]]], | ||
+ | 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> | ||
− | + | В первой строке записали название системы координат латиницей "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=" | + | <syntaxhighlight lang="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); | |
− | + | } | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
</syntaxhighlight> | </syntaxhighlight> | ||
− | + | Сохраним код в файл '''findkey.c'''. Создадим исполняемый модуль компилятором '''gcc''': | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
<syntaxhighlight lang="bash"> | <syntaxhighlight lang="bash"> | ||
− | + | $ gcc -o findkey findkey.c -lm | |
</syntaxhighlight> | </syntaxhighlight> | ||
− | + | Пользователи MS Windows могут загрузить уже скомпилированную [https://wiki.gis-lab.info/images/4/43/Findkey.zip программу]. | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
== Ссылки == | == Ссылки == | ||
* [http://pubs.usgs.gov/pp/1395/report.pdf Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987] | * [http://pubs.usgs.gov/pp/1395/report.pdf Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987] | ||
− | * [http:// | + | * [http://www.epsg.org/Portals/0/373-07-2.pdf Coordinate Conversions and Transformations including Formulas, EPSG Guidance Note 7-2] |
− | * [ | + | * [https://desktop.arcgis.com/ru/arcmap/latest/map/projections/local-cartesian-projection.htm Справка ArcGIS — Локальная проекция Декартовой системы координат] |
− | + | * [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 могут загрузить уже скомпилированную программу.
[править] Ссылки
- Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987
- Coordinate Conversions and Transformations including Formulas, EPSG Guidance Note 7-2
- Справка ArcGIS — Локальная проекция Декартовой системы координат
- Using PROJ
- Конформное преобразование
- Defining a coordinate system in WKT or PROJ format that has an Affine transformaiton and bounds