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

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

Первые четыре параметра — a₀, b₀, a₁, b₁, оставшиеся два — масштабный множитель m и разворот θ в градусах.

↑
Команды примера воспроизводятся в среде MinGW, что входит в состав QGIS под MS Windows. Две особенности:
- нужно заменить команду awk на gawk, или лучше в системе определить awk как синоним gawk;
- после некоторых команд (proj, pr) придётся добавить в пайп команду удаления лишних символов CR:
```
| tr -d '\r'
```

[1] Команды примера воспроизводятся в среде MinGW, что входит в состав QGIS под MS Windows. Две особенности:
нужно заменить команду awk на gawk, или лучше в системе определить awk как синоним gawk;

после некоторых команд (proj, pr) придётся добавить в пайп команду удаления лишних символов CR:
| tr -d '\r'

[2] нужно заменить команду awk на gawk, или лучше в системе определить awk как синоним gawk;

[3] после некоторых команд (proj, pr) придётся добавить в пайп команду удаления лишних символов CR:
| tr -d '\r'

[1]

@@ Строка 19: / Строка 19: @@
 Тестовый пример основан на сгенерированных данных. В таблице ''X'', ''Y'' — координаты пункта в государственной системе (ГСК), а именно в седьмой зоне СК-42, ''x'', ''y'' — координаты в местной системе (МСК), ''p'' — вес пункта.
 {| class="wikitable"
 |-
@@ Строка 33: / Строка 34: @@
 Подготовим файл исходных данных ''data0.dat'':
 <pre>
 7383477.64 6087377.60 1334.71   285.94 1.0
@@ Строка 41: / Строка 43: @@
 Убедимся в однородности данных обоих каталогов. Для этого вычислим параметры конформного преобразования:
 <syntaxhighlight lang="bash">
 $ findkey data0.dat key0.dat var0.dat
@@ Строка 46: / Строка 49: @@
 Сами по себе параметры нас не интересуют. Гораздо интереснее невязки, выведенные в файл '''var0.dat''':
 <pre>
 … -0.002 -0.001
@@ Строка 56: / Строка 60: @@
 Пересчитаем координаты из ГСК в долготу/широту:
 <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
@@ Строка 61: / Строка 66: @@
 Команда '''awk''' извлекает из файла '''data0.dat''' вторую и третью колонки и передаёт их утилите '''proj''', которая пересчитывает координаты из седьмой зоны СК-42 в долготу/широту и записывает результаты в файл '''pt_longlat.dat''':
 <pre>
 .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
 </syntaxhighlight>
 Содержимое '''pt_tmerc.dat''':
 <pre>
 .7100 285.9400
@@ Строка 81: / Строка 89: @@
 -248.5461 2884.0431
 </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>
+1334.7100 285.9400 1334.71 285.94 1.0
+556.2353 -5196.2593 563.67 -5197.34 1.0
+4445.4005 1149.5698 4444.27 1153.79 1.0
+-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>
+.390
+-1.813
+.0000052644
+.0013554914
+.0000061831
++0.07766348
+</pre>
+Первые четыре параметра — ''a''<sub>0</sub>, ''b''<sub>0</sub>, ''a''<sub>1</sub>, ''b''<sub>1</sub>, оставшиеся два — масштабный множитель ''m'' и разворот ''θ'' в градусах.
 <references />

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

Версия от 14:36, 10 марта 2013

Содержание

Введение

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

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

Переходная проекция

Навигация

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

Версия от 14:36, 10 марта 2013

Введение

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

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

Переходная проекция

Навигация

Поиск