Местная система координат линейного объекта: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
 
(не показано 69 промежуточных версий этого же участника)
Строка 5: Строка 5:
== Введение ==
== Введение ==


Система координат линейного объекта строится для эксплуатации протяжённого инженерного сооружения.
Система координат линейного объекта строится для обеспечения строительства или эксплуатации протяжённого инженерного сооружения.
Принципы построения проекции сходны с классическим подходом, изложенным в статье [http://gis-lab.info/qa/local-cs.html «Добавление местной координатной системы в GIS»].
Целью при этом является минимизация искажений, присущих проекции, в полосе объекта.
Однако постановка задачи отличается.


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


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


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


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


== О проекции ==
== О проекции ==
Строка 28: Строка 29:
* разворот координатных осей ''γ''
* разворот координатных осей ''γ''


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


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


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


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


Рассмотрим последовательность решения задачи с использованием '''PROJ.4'''. Вид строки параметров таков:
{| class="wikitable"
|-
! NN
! ''φ''
! ''λ''
|-
| 1
| 52°00′03.358″N
| 23°07′37.837″E
|-
| 2
| 52°00′46.722″N
| 23°10′15.918″E
|}


<pre>
Расстояние вдоль оси задано значением ''L'' = 3300.000 м.
+proj=omerc +lat_0=φ₀ +lonc=λ₀ +alpha=α +k_0=k₀ +x_0=x₀ +y_0=y₀ +gamma=γ
</pre>


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


<pre>
<pre>
+lat_0=51 +lonc=22
+proj=omerc +lat_0=φ₀ +lonc=λ₀ +alpha=α +k_0=k₀ +x_0=x₀ +y_0=y₀ +gamma=γ
</pre>
</pre>


Для нахождения остальных параметров нужно решить обратную геодезическую задачу (ОГЗ).
Задачу помещения центра проекции в середину линии, соединяющей конечные точки, решим в два этапа. Сначала решим обратную геодезическую задачу, что даст азимут с первой точки на вторую ''α''₁₂, азимут со второй точки на первую ''α''₂₁ и длину отрезка геодезической линии между ними ''S''. Затем решим прямую геодезическую задачу (ПГЗ), чтобы получить координаты средней точки и азимуты направлений с неё на конечные точки.
 
При значительной длине трассы (точнее, при заметном изменении широты вдоль линии трассы) корректнее поместить центр проекции в середину геодезической линии, соединяющей конечные точки. В это случае придётся сначала решить обратную геодезическую задачу, что даст азимут с первой точки на вторую ''α''₁₂, азимут со второй точки на первую ''α''₂₁ и длину отрезка геодезической линии между ними ''S'', а затем решить прямую геодезическую задачу (ПГЗ), чтобы получить координаты средней точки и азимуты направлений с неё на конечные точки.
<!--
Для определения параметра ''alpha'' нужно решить обратную геодезическую задачу и найти азимут в первой точке на вторую ''α''₁₂. Здесь возможны два случая:
* вторая точка севернее первой, −90° < ''α''₁₂ < +90°; ''alpha'' = ''α''₁₂
* вторая точка южнее первой, 90° < ''α''₁₂ < 270°; ''alpha'' = ''α''₁₂ − 180°


Поскольку поставлена задача направить ось ''OX'' вдоль начальной линии в противоположную сторону от второй точки, параметру ''gamma'' присвоим значение −90° в первом случае и +90° во втором.
Значение параметра ''k_0'' также можно оценить по результатам решения обратной геодезической задачи: ''k_0'' = ''L'' / ''S'', где ''L'' — заданная длина линии, ''S'' — длина геодезической из решения ОГЗ.
-->
=== Решение обратной геодезической задачи ===
=== Решение обратной геодезической задачи ===
<!--
Итак, из решения ОГЗ мы хотим получить азимут ''α''₁₂ и расстояние ''S'', нужные для определения параметров ''alpha'' и ''k_0'' соответственно.


Пикантность ситуации придаёт тот факт, что на эллипсоиде через две точки проходит геодезическая линия, которая в блестящей математике косой проекции сэра Мартина Хотайна отображается в кривую на апосфере, не совпадающую с дугой большого круга. Расчёты показывают, что вследствие этого точность соответствия ''alpha'' = ''α''₁₂ и ''k_0'' = ''L'' / ''S'' при расстояниях в несколько десятков километров перестаёт удовлетворять требованиям геодезии.
Цель — вычисление азимутов прямого и обратного направлений, а также длины отрезка геодезической линии на поверхности эллипсоида.


Впрочем, многие объекты имеют скромные размеры «всего» в несколько километров, и для них этот неприятный момент имеет исключительно теоретическое значение.
Решим ОГЗ с помощью утилиты '''geod''' из пакета '''PROJ.4''':
 
Однако есть ещё одно обстоятельство. Утилита '''geod''' из пакета '''PROJ.4''' версии 4.9.0 будет решать геодезические задачи на эллипсоиде, используя библиотеку '''GeographicLib'''. Но сегодня в ходу '''PROJ.4''' версии 4.8.0, и '''geod''' умеет считать только для сферы.
 
Так или иначе, надо научиться обходиться подручными средствами.
-->
Подготовим файл данных с координатами конечных пунктов '''inv.dat''':
<pre>
51N 22E 50N 20E
</pre>
и решим ОГЗ:
<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
$ geod -I -f "%.10f" -F "%f" +ellps=WGS84 +units=m inv.dat
$ 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
</syntaxhighlight>
</syntaxhighlight>
<!--
Программа выдаёт решение на эллипсоиде в виде строки значений ''α''₁₂, ''α''₂₁, ''S''₁₂:
Возможно, в следующей версии '''PROJ.4''' результатом этой команды будет решение на эллипсоиде. Сегодня же из решения на сфере радиусом, равным экваториальному радиусу эллипсоида WGS 84, получилась такая строка значений ''α''₁₂, ''α''₂₁, ''S'':
<pre>
<pre>
-127.3948484062 51.0617802663  180119.673397
66.017759443956336 -113.94763462689073 3299.7360258541303
</pre>
-->
Программа выдаёт решение на эллипсоиде в виде строки значений ''α''₁₂, ''α''₂₁, ''S'':
<pre>
-127.3190808614 51.1375475317 180292.395229
</pre>
</pre>


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


Цель — получить значения для середины отрезка геодезической линии. Прежде всего вычислим половину длины отрезка:
Цель — получить координаты и азимут середины отрезка геодезической линии. Прежде всего вычислим половину длины отрезка:
<pre>
<pre>
180292.395229 / 2 = 90146.1976145
3299.7360258541303 / 2 = 1649.86801292706515
</pre>
</pre>
Для контроля решим ПГЗ дважды, от обеих концевых точек линии. На основе данных ОГЗ создадим файл '''dir.dat''':
Для контроля решим ПГЗ дважды, от обоих концов линии. Используем ту же утилиту '''geod''':
<syntaxhighlight lang="bash">
$ 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
</syntaxhighlight>
Результатом будут две строки значений: ''φ''₀, ''λ''₀, ''α'':
<pre>
<pre>
51N 22E -127.3190808614 90146.1976145
52.006957604612808 23.14912969166868 -113.96494062446807
50N 20E  51.1375475317 90146.1976145
52.006957604612793 23.149129691668666 66.035059375531944
</pre>
</pre>
Используем ту же утилиту '''geod''' для решения прямой задачи:
Координаты центра проекции практически совпадают, азимуты обратных направлений отличаются на 180°.
<syntaxhighlight lang="bash">
$ geod -f "%.11f" +ellps=WGS84 +units=m dir.dat
</syntaxhighlight>
Результатом будут две строки значений


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


По результатам решения ОГЗ построим черновую проекцию. Поскольку ''φ''₂ < ''φ''₁, имеет место второй вариант; примем ''alpha'' = ''α''₁₂ + 180° = 52.6051515938°, ''gamma'' = +90°. Масштабный коэффициент пока приравняем единице: ''k_0'' = 1. Получен предварительный набор параметров:
По результатам решения ПГЗ построим проекцию в первом приближении. Параметры ''lat_0'' и ''lonc'' примем равными ''φ''₀ и ''λ''₀. Параметр ''alpha'' должен быть в диапазоне ±90°, примем для него значение ''α''₀₂ = ''α''₀₁ ± 180°. Чтобы направить ось ''OX'' вдоль направления 0–1, параметр разворота ''gamma'' примем равным 90°. Вот предварительный набор:
<pre>
<pre>
+proj=omerc +lat_0=51 +lonc=22 +alpha=52.6051515938 +gamma=90 +k_0=1 +x_0=0 +y_0=0
+lat_0=50.504316101 +lonc=20.989441172 +alpha=51.89831026 +k_0=1 +x_0=0 +y_0=0 +gamma=90
</pre>
</pre>
Подготовим файл с координатами конечных точек '''p12.dat''':
Вычислим координаты конечных точек в проекции:
<syntaxhighlight lang="bash">
$ 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
</syntaxhighlight>
Программа выдаёт координаты первой и второй точек ''x''₁, ''y''₁ и ''x''₂, ''y''₂:
<pre>
<pre>
22 51
-1649.8680129311654 -1.0102527905258779e-13
20 50
1649.8680129226682 1.0631327490867519e-09
</pre>
</pre>
Выполним команду:
Вычислим масштаб ''k_0'' как отношение заданной длины ''L'' к разности координат ''x''₂ − ''x''₁: ''k_0'' = 3300 / (1649.8680129226682 + 1649.8680129311654) = 1.0000799985647634. Подставим это значение вместо единицы:
<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
$ proj -f "%f" +proj=omerc +lat_0=51 +lonc=22 +alpha=52.6051515938 +gamma=90 +k_0=1 +x_0=0 +y_0=0 +ellps=WGS84 p12.dat
$ 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
</syntaxhighlight>
</syntaxhighlight>
Программа выдаёт координаты первой и второй точек ''x''₁, ''y''₁ и ''x''₂, ''y''₂ в МСК:
Вывод программы:
<pre>
<pre>
-0.000000      0.000000
-1650.0000000040786 -1.0103336092990638e-13
-180292.238188  238.386305
1649.9999999969975 1.0632177981808281e-09
</pre>
</pre>
Координаты первой точки равны нулю, как и должно быть. Координаты второй точки отличаются от ожидаемых значений −''L'' и 0. Ненулевая величина ''y''₂ говорит о том, что начальная линия проходит мимо второй точки. Значит, нужно подобрать параметр ''alpha'', чтобы исправить промах.
Точки практически лежат на оси ''OX'', расстояние между ними 3300.000 м. Поставленная задача решена, проекция построена.


Численные методы отлично справляются с поиском корня уравнения ''x''₂(''alpha'') = 0. Добавим немного геометрии в задачу и будем улучшать значение параметра по формуле ''alpha''′ = ''alpha'' − arctg [(''y''₂ − ''y''₁) / (''x''₂ − ''x''₁)]. Может, это не очень эффективно, но через пять итераций при ''alpha'' = 52.6809193468 получаем такие координаты:
== Пользовательская проекция в QGIS ==
<pre>
 
-0.000000      0.000000
Создадим пользовательскую систему координат в формате WKT.
-180292.395746  0.000000
 
</pre>
=== Вариант B ===
Вычисляем параметр ''k_0'' по формуле ''k_0'' = −''L'' / (''x''₂ − ''x''₁) = 180300 / 180292.238188 = 1.0000421773.
 
Косая проекция Меркатора может быть задана в вариантах A и B. Начнём со второго, поскольку в нём плоские координаты отсчитываются от центра проекции, и для него всё готово.
 
В начале введём название системы координат латиницей "Biala Podlaska airdrome".


Запуск '''proj''' с окончательным набором параметров:
При вводе параметров проекции поможет следующая таблица соответствия:
<syntaxhighlight lang="bash">
proj -f "%f" +proj=omerc +lat_0=51 +lonc=22 +alpha=52.6809193468 +gamma=90 +k_0=1.0000421773 +x_0=0 +y_0=0 +ellps=WGS84 p12.dat
</syntaxhighlight>
Результаты:
<pre>
-0.000  0.000
-180300.000    0.000
</pre>
Проекция построена.


== Вторая проекция ==
{| class="wikitable"
|-
! WKT !! PROJ
|-
| Latitude of projection centre || lat_0
|-
| Longitude of projection centre || lonc
|-
| Azimuth of initial line || alpha
|-
| Angle from Rectified to Skew Grid || gamma
|-
| Scale factor on initial line || k, k_0
|-
| Easting at projection centre || x_0
|-
| Northing at projection centre || y_0
|}


Нередко требуется вторая проекция, являющаяся зеркальным отражением первой: начало координат МСК-2 во второй точке, ось ''OX'' направлена вдоль оси в сторону, противоположную направлению на первую точку. Таким образом МСК-2 развёрнута по отношению к МСК-1 на 180° и смещена вдоль оси ''OX'' на длину ''L''.
В конце вставим название покрываемой территории "Europe - Poland - Biala Podlaska" и охват в формате ''φ''<sub>min</sub>, ''λ''<sub>min</sub>, ''φ''<sub>max</sub>, ''λ''<sub>max</sub>.


Последнее предложение имеет особый смысл для выбора способа построения МСК-2. Если выбрать способ, изложенный для МСК-1, только с центром проекции во второй точке и опорой на азимут ''α''₂₁, выяснится, что апосфера во втором случае будет не та, что в первом, и большие круги, проходящие через две точки, не совпадают. Правда, разница незаметна, пока расстояние не достигает величин в несколько десятков километров.
Готовая система координат в формате WKT:


Таким образом, если нужна пара взаимоувязанных МСК, вторая система строится на параметрах первой. Параметр ''gamma'' изменяем на 180°, параметру ''x_0'' присваиваем значение ''−L'', всего-то и делов.
<syntaxhighlight lang="xml">
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]]]
</syntaxhighlight>


== Тестирование ==
=== Вариант A ===


Создадим файл с координатами двух точек '''pt34.dat''' на эллипсоиде:
В варианте A прямоугольные координаты отсчитываются от точки пересечения начальной линии с экватором апосферы.
<pre>
Практически он отличается от B числовыми значениями параметров ''x_0'' и ''y_0''.
21 51
 
21 50
Вычислим координаты центра проекции с параметром ''no_off'':
</pre>
Вычислим координаты в МСК-1:
<syntaxhighlight lang="bash">
<syntaxhighlight lang="bash">
proj -f "%.3f" +proj=omerc +lat_0=51 +lonc=22 +alpha=52.6809193468 +gamma=90 +k_0=1.0000421773 +x_0=0 +y_0=0 +ellps=WGS84 p34.dat
$ 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 +no_off +units=m <<EOF
> 52.0069576046128 23.1491296916687
> EOF
</syntaxhighlight>
</syntaxhighlight>
Вывод программы:
<pre>
<pre>
-55539.071      42936.465
8064096.0024511777 -2.1496104200899204e-10
-124171.432    -44612.843
</pre>
</pre>
Вычислим координаты в МСК-2:
Поменяв знаки на противоположные, получим параметры ''x_0'' = -8064096.00245118, ''y_0'' = 0.
<syntaxhighlight lang="bash">
 
proj -f "%.3f" +proj=omerc +lat_0=51 +lonc=22 +alpha=52.6809193468 +gamma=-90 +k_0=1.0000421773 +x_0=-180300 +y_0=0 +ellps=WGS84 p34.dat
Описание проекции в формате WKT по варианту A:
</syntaxhighlight>
 
<pre>
<pre>
-124760.929     -42936.465
PROJCRS["Biala Podlaska airdrome",
-56128.568      44612.843
    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 A)",
            ID["EPSG",9812]],
        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",-8064096.00245118,
            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]]]
</pre>
</pre>
Калькулятор подтверждает, что:
* суммы координат ''x'' соответствующих точек равны ''−L''
* суммы координат ''y'' соответствующих точек равны нулю


== Заключение ==
== Заключение ==


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


== Ссылки ==
== Ссылки ==

Текущая версия от 05:44, 12 августа 2022

Эта страница опубликована в основном списке статей сайта
по адресу 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.

Вариант B

Косая проекция Меркатора может быть задана в вариантах A и B. Начнём со второго, поскольку в нём плоские координаты отсчитываются от центра проекции, и для него всё готово.

В начале введём название системы координат латиницей "Biala Podlaska airdrome".

При вводе параметров проекции поможет следующая таблица соответствия:

WKT PROJ
Latitude of projection centre lat_0
Longitude of projection centre lonc
Azimuth of initial line alpha
Angle from Rectified to Skew Grid gamma
Scale factor on initial line k, 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]]]

Вариант A

В варианте A прямоугольные координаты отсчитываются от точки пересечения начальной линии с экватором апосферы. Практически он отличается от B числовыми значениями параметров x_0 и y_0.

Вычислим координаты центра проекции с параметром no_off: 

$ 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 +no_off +units=m <<EOF
> 52.0069576046128 23.1491296916687
> EOF

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

8064096.0024511777 -2.1496104200899204e-10

Поменяв знаки на противоположные, получим параметры x_0 = -8064096.00245118, y_0 = 0.

Описание проекции в формате WKT по варианту A: 

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 A)",
            ID["EPSG",9812]],
        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",-8064096.00245118,
            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, которая используется как чёрный ящик. Этот подход не совсем корректен, поскольку геодезическая линия, соединяющая две точки на эллипсоиде, в косой проекции Меркатора отображается в кривую на апосфере, близкую к дуге большого круга, но не совпадающую с ней. К счастью, это несущественно даже для объектов длиной в сотни и тысячи километров.

Ссылки

Источник — https://wiki.gis-lab.info/index.php?title=Местная_система_координат_линейного_объекта&oldid=26804