Местная система координат линейного объекта
по адресу http://gis-lab.info/qa/local-cs-linear-object.html
Конструирование проекции для представления системы координат линейного объекта в ГИС
Введение
Система координат линейного объекта строится для обеспечения строительства или эксплуатации протяжённого инженерного сооружения. Целью при этом является минимизация искажений, присущих проекции, в полосе объекта.
Постановка задачи
Пусть ось сооружения задана положением двух крайних точек в глобальной системе координат (ГСК).
В местной системе (МСК) совместим начало координат с точкой в середине отрезка геодезической линии, соединяющей крайние точки. Ось OX направим вдоль оси сооружения. Потребуем, чтобы расстояние между крайними точками равнялось априори заданной величине L.
Требуется подобрать проекцию, подходящую для представления такой МСК в ГИС и в программах, используемых геодезистами.
О проекции
Выбор проекции однозначен. Это косая проекция Меркатора с такими значениями параметров, чтобы так называемая начальная линия (линия наименьшего масштаба) проходила через конечные точки, а расстояние между этими точками равнялось L.
Для косой проекции Меркатора задаются следующие параметры:
- широта и долгота центра проекции φ₀, λ₀
- азимут начальной линии α
- масштаб на начальной линии k₀
- прямоугольные координаты в центре проекции x₀, y₀
- разворот координатных осей γ
Азимут начальной линии должен находиться в диапазоне −90° < α < +90°. Таким образом, если разворот γ равен нулю, ось OY будет направлена вдоль начальной линии в северную полуплоскость, OX в восточную.
Разворот γ обычно приравнивается значению α, чтобы компенсировать начальный разворот осей и вернуть оси OY направление строго на север. Возможность его явного задания позволяет произвольно управлять ориентацией осей МСК. Если задать нулевой разворот γ, ось OY будет направлена вдоль начальной линии в северную полуплоскость, OX перпендикулярно к начальной линии в восточную.
Определение параметров
Приведём данные тестового примера. Осевая линия задана координатами конечных точек на эллипсоиде WGS 84:
| NN | φ | λ |
|---|---|---|
| 1 | 52°00′03.358″N | 23°07′37.837″E |
| 2 | 52°00′46.722″N | 23°10′15.918″E |
Расстояние вдоль оси задано значением L = 3300.000 м.
Рассмотрим последовательность решения задачи с использованием PROJ. Вид строки параметров таков:
+proj=omerc +lat_0=φ₀ +lonc=λ₀ +alpha=α +k_0=k₀ +x_0=x₀ +y_0=y₀ +gamma=γ
Задачу помещения центра проекции в середину линии, соединяющей конечные точки, решим в два этапа. Сначала решим обратную геодезическую задачу, что даст азимут с первой точки на вторую α₁₂, азимут со второй точки на первую α₂₁ и длину отрезка геодезической линии между ними S. Затем решим прямую геодезическую задачу (ПГЗ), чтобы получить координаты средней точки и азимуты направлений с неё на конечные точки.
Решение обратной геодезической задачи
Решим ОГЗ с помощью утилиты geod из пакета PROJ.4:
$ geod -I -f "%.17g" -F "%.17g" +ellps=WGS84 +units=m <<EOF
> 52d00'03.358"N 023d07'37.837"E 52d00'46.722"N 023d10'15.918"E
> EOF
Программа выдаёт решение на эллипсоиде в виде строки значений α₁₂, α₂₁, S₁₂:
66.017759443956336 -113.94763462689073 3299.7360258541303
Решение прямой геодезической задачи
Цель — получить координаты и азимут середины отрезка геодезической линии. Прежде всего вычислим половину длины отрезка:
3299.7360258541303 / 2 = 1649.86801292706515
Для контроля решим ПГЗ дважды, от обоих концов линии. Используем ту же утилиту geod:
$ geod -f "%.17g" +ellps=WGS84 +units=m <<EOF
> 52d00'03.358"N 023d07'37.837"E 66.017759443956336 1649.86801292706515
> 52d00'46.722"N 023d10'15.918"E -113.94763462689073 1649.86801292706515
> EOF
Результатом будут две строки значений: φ₀, λ₀, α:
52.006957604612808 23.14912969166868 -113.96494062446807 52.006957604612793 23.149129691668666 66.035059375531944
Координаты центра проекции практически совпадают, азимуты обратных направлений отличаются на 180°.
Построение проекции
По результатам решения ПГЗ построим проекцию в первом приближении. Параметры lat_0 и lonc примем равными φ₀ и λ₀. Параметр alpha должен быть в диапазоне ±90°, примем для него значение α₀₂ = α₀₁ ± 180°. Чтобы направить ось 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
Выполним координаты конечных точек в проекции:
$ proj -r -f "%.17g" +proj=omerc +lat_0=52.0069576046128 +lonc=23.1491296916687 +alpha=66.0350593755319 +k_0=1 +x_0=0 +y_0=0 +gamma=90 +ellps=WGS84 +units=m <<EOF
> 52d00'03.358"N 023d07'37.837"E
> 52d00'46.722"N 023d10'15.918"E
> EOF
Программа выдаёт координаты первой и второй точек x₁, y₁ и x₂, y₂:
-1649.8680129311654 -1.0102527905258779e-13 1649.8680129226682 1.0631327490867519e-09
Вычислим масштаб k_0 как отношение заданной длины L к разности координат x₂ − x₁: k_0 = 3300 / (1649.8680129226682 + 1649.8680129311654) = 1.0000799985647634. Подставим это значение вместо единицы:
$ proj -r -f "%.17g" +proj=omerc +lat_0=52.0069576046128 +lonc=23.1491296916687 +alpha=66.0350593755319 +k_0=1.00007999856476 +x_0=0 +y_0=0 +gamma=90 +ellps=WGS84 +units=m <<EOF
> 52d00'03.358"N 023d07'37.837"E
> 52d00'46.722"N 023d10'15.918"E
> EOF
Вывод программы:
-1650.0000000040786 -1.0103336092990638e-13 1649.9999999969975 1.0632177981808281e-09
Точки практически лежат на оси OX, расстояние между ними 3300.000 м. Поставленная задача решена, проекция построена.
Пользовательская проекция в QGIS
Создадим пользовательскую проекцию в формате WKT.
В начале введём название системы координат латиницей "Biala Podlaska airdrome".
Таблица соответствия параметров проекции WKT и PROJ:
| Latitude of projection centre | lat_0 |
| Longitude of projection centre | lon_c |
| Azimuth of initial line | alpha |
| Angle from Rectified to Skew Grid | gamma |
| Scale factor on initial line | k_0 |
| Easting at projection centre | x_0 |
| Northing at projection centre | y_0 |
В конце вставим название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате φmin, λmin, φmax, λmax.
Проекция в формате WKT:
PROJCRS["Biala Podlaska airdrome",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Biala Podlaska airdrome",
METHOD["Hotine Oblique Mercator (variant B)",
ID["EPSG",9815]],
PARAMETER["Latitude of projection centre",52.0069576046128,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8811]],
PARAMETER["Longitude of projection centre",23.1491296916687,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8812]],
PARAMETER["Azimuth of initial line",66.0350593755319,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8813]],
PARAMETER["Angle from Rectified to Skew Grid",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8814]],
PARAMETER["Scale factor on initial line",1.00007999856476,
SCALEUNIT["unity",1],
ID["EPSG",8815]],
PARAMETER["Easting at projection centre",0,
LENGTHUNIT["metre",1],
ID["EPSG",8816]],
PARAMETER["Northing at projection centre",0,
LENGTHUNIT["metre",1],
ID["EPSG",8817]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["unknown"],
AREA["Europe - Poland - Biala Podlaska"],
BBOX[51.9,23.0,52.1,23.3]]]
Заключение
Рассмотренный способ построения проекции прост, поскольку позволяет заменить знание математической картографии обращением к утилите geod из библиотеки PROJ, которая используется как чёрный ящик. Этот подход не совсем корректен, поскольку геодезическая линия, соединяющая две точки на эллипсоиде, в косой проекции Меркатора отображается в кривую на апосфере, близкую к дуге большого круга, но не совпадающую с ней. К счастью, это несущественно даже для объектов длиной в сотни и тысячи километров.