Добавление местной координатной системы в GIS: различия между версиями
ErnieBoyd (обсуждение | вклад) Нет описания правки |
ErnieBoyd (обсуждение | вклад) Нет описания правки |
||
| Строка 19: | Строка 19: | ||
Тестовый пример основан на сгенерированных данных. В таблице ''X'', ''Y'' — координаты пункта в государственной системе (ГСК), а именно в седьмой зоне СК-42, ''x'', ''y'' — координаты в местной системе (МСК), ''p'' — вес пункта. | Тестовый пример основан на сгенерированных данных. В таблице ''X'', ''Y'' — координаты пункта в государственной системе (ГСК), а именно в седьмой зоне СК-42, ''x'', ''y'' — координаты в местной системе (МСК), ''p'' — вес пункта. | ||
{| class="wikitable" | {| class="wikitable" | ||
|- | |- | ||
| Строка 33: | Строка 34: | ||
Подготовим файл исходных данных ''data0.dat'': | Подготовим файл исходных данных ''data0.dat'': | ||
<pre> | <pre> | ||
1 7383477.64 6087377.60 1334.71 285.94 1.0 | 1 7383477.64 6087377.60 1334.71 285.94 1.0 | ||
| Строка 41: | Строка 43: | ||
Убедимся в однородности данных обоих каталогов. Для этого вычислим параметры конформного преобразования: | Убедимся в однородности данных обоих каталогов. Для этого вычислим параметры конформного преобразования: | ||
<syntaxhighlight lang="bash"> | <syntaxhighlight lang="bash"> | ||
$ findkey data0.dat key0.dat var0.dat | $ findkey data0.dat key0.dat var0.dat | ||
| Строка 46: | Строка 49: | ||
Сами по себе параметры нас не интересуют. Гораздо интереснее невязки, выведенные в файл '''var0.dat''': | Сами по себе параметры нас не интересуют. Гораздо интереснее невязки, выведенные в файл '''var0.dat''': | ||
<pre> | <pre> | ||
1 … -0.002 -0.001 | 1 … -0.002 -0.001 | ||
| Строка 56: | Строка 60: | ||
Пересчитаем координаты из ГСК в долготу/широту: | Пересчитаем координаты из ГСК в долготу/широту: | ||
<syntaxhighlight lang="bash"> | <syntaxhighlight lang="bash"> | ||
$ awk '{print $2, " ", $3}' data0.dat | proj -I -f "%.9f" +proj=tmerc +lat_0=0 +lon_0=39 +k=1 +x_0=7500000 +y_0=0 +ellps=krass > pt_longlat.dat | $ awk '{print $2, " ", $3}' data0.dat | proj -I -f "%.9f" +proj=tmerc +lat_0=0 +lon_0=39 +k=1 +x_0=7500000 +y_0=0 +ellps=krass > pt_longlat.dat | ||
| Строка 61: | Строка 66: | ||
Команда '''awk''' извлекает из файла '''data0.dat''' вторую и третью колонки и передаёт их утилите '''proj''', которая пересчитывает координаты из седьмой зоны СК-42 в долготу/широту и записывает результаты в файл '''pt_longlat.dat''': | Команда '''awk''' извлекает из файла '''data0.dat''' вторую и третью колонки и передаёт их утилите '''proj''', которая пересчитывает координаты из седьмой зоны СК-42 в долготу/широту и записывает результаты в файл '''pt_longlat.dat''': | ||
<pre> | <pre> | ||
37.183749701 54.896961118 | 37.183749701 54.896961118 | ||
| Строка 74: | Строка 80: | ||
$ proj -f "%.4f" +proj=tmerc +lat_0=54.896961118 +lon_0=37.183749701 +k=1 +x_0=1334.71 +y_0=285.94 +ellps=krass pt_longlat.dat > pt_tmerc.dat | $ proj -f "%.4f" +proj=tmerc +lat_0=54.896961118 +lon_0=37.183749701 +k=1 +x_0=1334.71 +y_0=285.94 +ellps=krass pt_longlat.dat > pt_tmerc.dat | ||
</syntaxhighlight> | </syntaxhighlight> | ||
Содержимое '''pt_tmerc.dat''': | Содержимое '''pt_tmerc.dat''': | ||
<pre> | <pre> | ||
1334.7100 285.9400 | 1334.7100 285.9400 | ||
| Строка 81: | Строка 89: | ||
-248.5461 2884.0431 | -248.5461 2884.0431 | ||
</pre> | </pre> | ||
Подготовим входной файл '''data1.dat''' для '''findkey'''. Нужно на место второй и третьей колонок в '''data0.dat''' записать данные '''pt_tmerc.dat'''. | |||
<syntaxhighlight lang="bash"> | |||
$ pr -m -t -s\ data0.dat pt_tmerc.dat | awk '{print $1, $7, $8, $4, $5, $6}' > data1.dat | |||
</syntaxhighlight> | |||
Команда '''pr''' построчно объединяет данные из шести колонок '''data0.dat''' и двух колонок '''pt_tmerc.dat''' в общую строку, команда '''awk''' выводит в '''data1.dat''' нужные колонки в нужном же порядке: | |||
<pre> | |||
1 1334.7100 285.9400 1334.71 285.94 1.0 | |||
2 556.2353 -5196.2593 563.67 -5197.34 1.0 | |||
3 4445.4005 1149.5698 4444.27 1153.79 1.0 | |||
4 -248.5461 2884.0431 -252.07 2881.90 1.0 | |||
</pre> | |||
Нетрудно заметить, что координаты первого пункта в переходной проекции совпадают с таковыми в МСК. Для остальных пунктов это не так, что говорит о наличии масштабирования и разворота. Убедимся в этом: | |||
<syntaxhighlight lang="bash"> | |||
$ findkey data1.dat key1.dat var1.dat | |||
</syntaxhighlight> | |||
Содержимое '''key1.dat''': | |||
<pre> | |||
0.390 | |||
-1.813 | |||
1.0000052644 | |||
0.0013554914 | |||
1.0000061831 | |||
+0.07766348 | |||
</pre> | |||
Первые четыре параметра — ''a''<sub>0</sub>, ''b''<sub>0</sub>, ''a''<sub>1</sub>, ''b''<sub>1</sub>, оставшиеся два — масштабный множитель ''m'' и разворот ''θ'' в градусах. | |||
<references /> | <references /> | ||
Версия от 15:36, 10 марта 2013
Введение
Под местной системой координат (МСК) будет подразумеваться так называемая «городская» система, построенная независимо от государственной системы (ГСК) и включенная в неё заданием ключей перехода к СК-42 или СК-63. МСК крупных территорий, сравнимых с размерами субъектов Федерации, не являются предметом данной статьи, поскольку относятся к классическим картографическим проекциям.
Постановка задачи
Некоторые программы позволяют реализовать работу в МСК непосредственно. Так, в MapInfo любая проекция может быть дополнена аффинным преобразованием. ArcGIS в качестве МСК предлагает локальную проекцию: аналог ортометрической проекции на эллипсоиде, дополненный разворотом и масштабированием.
Другие программы, включая QGIS, работают только с «обыкновенными» проекциями, не допуская дополнительных геометрических преобразований. Задача статьи — показать, как сконструировать проекцию, позволяющую работать в МСК в таких средах, как QGIS или, скажем, бортовой софт приёмников GARMIN.
В качестве рабочей среды будем использовать командную строку UNIX. Это идеальный инструмент для экспериментирования, позволяющий непринуждённо сочетать PROJ.4, findkey и утилиты обработки текстовых потоков awk, pr.[1]
Подготовка данных
Тестовый пример основан на сгенерированных данных. В таблице X, Y — координаты пункта в государственной системе (ГСК), а именно в седьмой зоне СК-42, x, y — координаты в местной системе (МСК), p — вес пункта.
| ID | X | Y | x | y | p |
|---|---|---|---|---|---|
| 1 | 7383477.64 | 6087377.60 | 1334.71 | 285.94 | 1.0 |
| 2 | 7382557.14 | 6081916.51 | 563.67 | −5197.34 | 1.0 |
| 3 | 7386610.19 | 6088160.39 | 4444.27 | 1153.79 | 1.0 |
| 4 | 7381962.05 | 6090016.34 | −252.07 | 2881.90 | 1.0 |
Подготовим файл исходных данных data0.dat:
1 7383477.64 6087377.60 1334.71 285.94 1.0 2 7382557.14 6081916.51 563.67 -5197.34 1.0 3 7386610.19 6088160.39 4444.27 1153.79 1.0 4 7381962.05 6090016.34 -252.07 2881.90 1.0
Убедимся в однородности данных обоих каталогов. Для этого вычислим параметры конформного преобразования:
$ findkey data0.dat key0.dat var0.dat
Сами по себе параметры нас не интересуют. Гораздо интереснее невязки, выведенные в файл var0.dat:
1 … -0.002 -0.001 2 … -0.017 0.013 3 … 0.031 0.017 4 … -0.012 -0.029
Малые значения невязок говорят о том, что данные достаточно хороши, чтобы с ними работать.
Пересчитаем координаты из ГСК в долготу/широту:
$ awk '{print $2, " ", $3}' data0.dat | proj -I -f "%.9f" +proj=tmerc +lat_0=0 +lon_0=39 +k=1 +x_0=7500000 +y_0=0 +ellps=krass > pt_longlat.dat
Команда awk извлекает из файла data0.dat вторую и третью колонки и передаёт их утилите proj, которая пересчитывает координаты из седьмой зоны СК-42 в долготу/широту и записывает результаты в файл pt_longlat.dat:
37.183749701 54.896961118 37.171630976 54.847714659 37.232243035 54.904709276 37.159058387 54.920296878
Переходная проекция
Создадим в качестве переходной проекцию Гаусса-Крюгера, осевой меридиан и начальная параллель которой пересекаются в первом пункте, а координаты в этой точке совпадают с координатами МСК этого пункта. Вычислим координаты в переходной проекции и запишем их в файл pt_tmerc:
$ proj -f "%.4f" +proj=tmerc +lat_0=54.896961118 +lon_0=37.183749701 +k=1 +x_0=1334.71 +y_0=285.94 +ellps=krass pt_longlat.dat > pt_tmerc.dat
Содержимое pt_tmerc.dat:
1334.7100 285.9400 556.2353 -5196.2593 4445.4005 1149.5698 -248.5461 2884.0431
Подготовим входной файл data1.dat для findkey. Нужно на место второй и третьей колонок в data0.dat записать данные pt_tmerc.dat.
$ pr -m -t -s\ data0.dat pt_tmerc.dat | awk '{print $1, $7, $8, $4, $5, $6}' > data1.dat
Команда pr построчно объединяет данные из шести колонок data0.dat и двух колонок pt_tmerc.dat в общую строку, команда awk выводит в data1.dat нужные колонки в нужном же порядке:
1 1334.7100 285.9400 1334.71 285.94 1.0 2 556.2353 -5196.2593 563.67 -5197.34 1.0 3 4445.4005 1149.5698 4444.27 1153.79 1.0 4 -248.5461 2884.0431 -252.07 2881.90 1.0
Нетрудно заметить, что координаты первого пункта в переходной проекции совпадают с таковыми в МСК. Для остальных пунктов это не так, что говорит о наличии масштабирования и разворота. Убедимся в этом:
$ findkey data1.dat key1.dat var1.dat
Содержимое key1.dat:
0.390 -1.813 1.0000052644 0.0013554914 1.0000061831 +0.07766348
Первые четыре параметра — a0, b0, a1, b1, оставшиеся два — масштабный множитель m и разворот θ в градусах.
- ↑ Команды примера воспроизводятся в среде MinGW, что входит в состав QGIS под MS Windows. Две особенности:
- нужно заменить команду awk на gawk, или лучше в системе определить awk как синоним gawk;
- после некоторых команд (proj, pr) придётся добавить в пайп команду удаления лишних символов CR:
| tr -d '\r'