Местная система координат линейного объекта

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

Конструирование проекции для представления системы координат линейного объекта в ГИС

Введение

Система координат линейного объекта строится для эксплуатации протяжённого инженерного сооружения. Принципы построения проекции сходны с классическим подходом, изложенным в статье «Добавление местной координатной системы в GIS». Однако постановка задачи отличается.

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

На оси сооружения задана линия положением двух его конечных точек в глобальной системе координат (ГСК).

Пусть в местной системе (МСК) начало координат совмещено с первой точкой, расстояние между точками задано величиной L, а ось OX направлена вдоль оси сооружения наружу. В такой системе координаты второй точки будут равны X = −L, Y = 0.

Требуется подобрать проекцию, подходящую для представления такой МСК в ГИС.

О проекции

Выбор проекции однозначен. Это косая проекция Меркатора с такими значениями параметров, чтобы так называемая начальная линия (линия наименьшего масштаба) проходила через конечные точки, а расстояние между этими точками равнялось L.

Для косой проекции Меркатора задаются следующие параметры:

• широта и долгота центра проекции φ₀, λ₀
• азимут начальной линии α
• масштаб на начальной линии k₀
• прямоугольные координаты в центре проекции x₀, y₀
• разворот координатных осей γ

Азимут начальной линии должен находиться в диапазоне −90° < α < +90°. Таким образом, если разворот γ равен нулю, ось OY будет направлена вдоль начальной линии в северную полуплоскость, OX в восточную.

Из-за ограничений реализации в библиотеке PROJ.4 азимут α не может равняться 0°. Это не проблема, — если ось направлена вдоль меридиана, можно выбрать проекцию Гаусса-Крюгера. Также α не может принимать значения ±90°. Это тоже не проблема, поскольку в подобных случаях азимут вдоль геодезической линии меняется довольно быстро, и можно выбрать центр проекции на некотором удалении от первоначально выбранной точки. Если же трасса изгибается вдоль параллели, лучше остановить выбор на равноугольной конической проекции Ламберта.

Разворот γ был введён для компенсации начального разворота осей, чтобы вернуть оси OY направление строго на север. Поэтому, если он не задан, подразумевается его численное равенство α. Возможность его явного задания позволяет произвольно управлять ориентацией осей МСК.

Определение параметров

Приведём данные тестового примера. Осевая линия задана положением конечных точек на эллипсоиде WGS 84: φ₁ = 51° с.ш., λ₁ = 22° в.д., φ₂ = 50° с.ш., λ₂ = 20° в.д. Расстояние вдоль оси задано длиной L = 180300 м.

Рассмотрим последовательность решения задачи с использованием PROJ.4. Вид строки параметров таков:

+proj=omerc +lat_0=φ₀ +lonc=λ₀ +alpha=α +k_0=k₀ +x_0=x₀ +y_0=y₀ +gamma=γ

Простой подход состоит в том, чтобы поместить центр проекции в одну из конечных точек. Тогда два параметра можно определить сразу:

+lat_0=51 +lonc=22

Для нахождения остальных параметров нужно решить обратную геодезическую задачу (ОГЗ).

При значительной длине трассы (точнее, при заметном изменении широты вдоль линии трассы) корректнее поместить центр проекции в середину линии, соединяющей конечные точки. В это случае придётся сначала решить обратную геодезическую задачу, что даст азимут с первой точки на вторую α₁₂, азимут со второй точки на первую α₂₁ и длину отрезка геодезической линии между ними S, а затем решить прямую геодезическую задачу (ПГЗ), чтобы получить координаты средней точки и азимуты направлений с неё на конечные точки.

Решение обратной геодезической задачи

Подготовим файл данных с координатами конечных пунктов inv.dat:

51N 22E 50N 20E

и решим ОГЗ с помощью утилиты geod из пакета PROJ.4:

$ geod -I -f "%.8f" -F "%.4f" +ellps=WGS84 +units=m inv.dat

Программа выдаёт решение на эллипсоиде в виде строки значений α₁₂, α₂₁, S:

-127.31908086 51.13754753 180292.3952

Решение прямой геодезической задачи

Цель — получить значения для середины отрезка геодезической линии. Прежде всего вычислим половину длины отрезка:

180292.3952 / 2 = 90146.1976

Для контроля решим ПГЗ дважды, от обоих концов линии. На основе данных ОГЗ создадим файл dir.dat:

51N 22E -127.31908086 90146.1976
50N 20E   51.13754753 90146.1976

Используем ту же утилиту geod для решения прямой задачи:

$ geod -f "%.9f" +ellps=WGS84 +units=m dir.dat

Результатом будут две строки значений: φ₀, λ₀, α:

50.504316101 20.989441172 51.898310258
50.504316101 20.989441172 -128.101689745

Координаты центра проекции совпадают, азимуты обратных направлений отличаются на 180°.

Построение проекции

По результатам решения ПГЗ построим проекцию в первом приближении. Параметры lat_0 и lonc примем равными φ₀ и λ₀. Параметр alpha должен быть в диапазоне ±90°, примем для него значение α₀₁ из первого решения. Чтобы направить ось OX вдоль направления 0–1, параметру разворота gamma присвоим значение +90°. Вот предварительный набор:

+lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +k_0=1 +x_0=0 +y_0=0 +gamma=90

Подготовим файл с координатами конечных точек p12.dat:

22 51
20 50

Выполним команду:

$ proj -f "%.4f" +proj=omerc +lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +gamma=90 +k_0=1 +x_0=0 +y_0=0 +datum=WGS84 p12.dat

Программа выдаёт координаты первой и второй точек x₁, y₁ и x₂, y₂:

90146.1976 0.0000
-90146.1976 0.0000

Вычислим масштаб k_0 как отношение заданной длины L к разности координат x₁ − x₂: k_0 = 180300 / (90146.1976 + 90146.1976) = 1.0000421804. Подставим это значение вместо единицы:

$ proj -f "%.4f" +proj=omerc +lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +gamma=90 +k_0=1.0000421804 +x_0=0 +y_0=0 +datum=WGS84 p12.dat

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

90150.0000 0.0000
-90150.0001 0.0000

Перенесём начало координат в первую точку: x₁ = −90150, y₁ = 0. Выполним команду

$ proj -f "%.4f" +proj=omerc +lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +gamma=90 +k_0=1.0000421804 +x_0=-90150 +y_0=0 +datum=WGS84 p12.dat

и получим координаты в МСК:

-0.0000 0.0000
-180300.0001 0.0000

Координаты первой точки x₁ = 0, y₁ = 0. Координаты второй точки x₂ = −L, y₂ = 0. Поставленная задача решена, проекция построена.

Вторая проекция

Нередко требуется вторая проекция, являющаяся зеркальным отражением первой: начало координат МСК-2 во второй точке, ось OX направлена вдоль оси в сторону, противоположную направлению на первую точку. Таким образом МСК-2 развёрнута по отношению к МСК-1 на 180° и смещена вдоль оси OX на длину L.

Вторая система строится на параметрах первой, только параметр gamma изменяем на 180°.

$ proj -f "%.4f" +proj=omerc +lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +gamma=-90 +k_0=1.0000421804 +x_0=-90150 +y_0=0 +datum=WGS84 p12.dat

-180300.0000 -0.0000
0.0001 -0.0000

Тестирование

Создадим файл с координатами двух точек pt34.dat на эллипсоиде:

21 51
21 50

Вычислим координаты в МСК-1:

$ proj -f "%.4f" +proj=omerc +lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +gamma=90 +k_0=1.0000421804 +x_0=-90150 +y_0=0 +datum=WGS84 p34.dat

-55539.0711 42936.4649
-124171.4321 -44612.8429

Вычислим координаты в МСК-2:

$ proj -f "%.4f" +proj=omerc +lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +gamma=-90 +k_0=1.0000421804 +x_0=-90150 +y_0=0 +datum=WGS84 p34.dat

-124760.9289 -42936.4649
-56128.5679 44612.8429

Калькулятор подтверждает, что:

• суммы координат x соответствующих точек равны −L;
• суммы координат y соответствующих точек равны нулю.

Заключение

Рассмотренный способ построения проекции прост, поскольку позволяет заменить знание математической картографии обращением к утилите geod из библиотеки PROJ.4, которая используется как чёрный ящик. Этот подход не совсем корректен, поскольку геодезическая линия, соединяющая две точки на эллипсоиде, в косой проекции Меркатора отображается в кривую на апосфере, близкую к дуге большого круга, но не совпадающую с ней. К счастью, это несущественно даже для объектов длиной в сотни и тысячи километров.

Ссылки

• Map Projections — A Working Manual, Snyder J. P., USGS Professional Paper 1395, 1987
• Coordinate Conversions and Transformations including Formulas, EPSG Guidance Note 7, 2002
• Hotine Oblique Mercator
• man_proj – PROJ.4
• man_geod – PROJ.4
• GeographicLib
• Добавление местной координатной системы в GIS