Просмотр кода страницы Геодезические системы пространственных координат

{{Статья|Черновик}}

{{Аннотация|Рассматриваются преобразования между пространственными координатными системами. Приводится пример программной реализации на языке Питон.}}

== Земной эллипсоид ==

Земным эллипсоидом называется эллипсоид вращения, поверхность которого по форме и размерам довольно близка к поверхности геоида.

Поверхность эллипсоида образуется вращением эллипса вокруг его малой оси.

Эллипс обычно определяется размером его большой полуоси ''a'' и сжатием ''f''. Реже вместо сжатия задаётся размер малой полуоси ''b'':

: <math>f = \frac{a - b}{a} \qquad b = a (1 - f)</math>

В теории и практике вычислений широко используются такие параметры, как полярный радиус кривизны поверхности ''c'', первый эксцентриситет ''e'' и второй эксцентриситет ''e′'':

: <math>c = \frac{a}{1 - f} \qquad e^2 = f (2 - f) \qquad e'^2 = \frac{e^2}{1 - e^2}</math>

Пример функции Питона, вычисляющей по ''a'' и ''f'' параметры ''b'', ''c'', ''e'' и ''e′'':

<syntaxhighlight lang="python">
def initSpher(a, f):
    b = a * (1. - f)
    c = a / (1. - f)
    e2 = f * (2. - f)
    e12 = e2 / (1. - e2)
    return (b, c, e2, e12)
</syntaxhighlight>

== Системы координат ==

Рассмотрим следующие системы координат.
# Геоцентрические декартовы прямоугольные координаты:
#* начало координат находится в центре эллипсоида,
#* ось ''z'' расположена вдоль оси вращения эллипсоида и направлена в северный полюс,
#* ось ''x'' лежит в пересечении экватора и начального меридиана,
#* ось ''y'' лежит в пересечении экватора и меридиана с долготой ''L'' = 90°.
# Система геодезических координат:
#; геодезическая широта ''B'' : угол между нормалью к поверхности эллипсоида и плоскостью экватора,
#; геодезическая долгота ''L'' : двугранный угол между плоскостями данного и начального меридианов,
#; геодезическая высота ''H'' : кратчайшее расстояние до поверхности эллипсоида.
# Топоцентрические декартовы прямоугольные координаты:
#* начало координат находится в некоторой точке ''Q''₀ (''B''₀, ''L''₀, ''H''₀) над эллипсоидом,
#* ось ''z'' расположена вдоль нормали к поверхности эллипсоида и направлена вверх,
#* ось ''x'' расположена в плоскости меридиана и направлена на север,
#* ось ''y'' перпендикулярна к осям ''x'' и ''z'' и направлена на восток.

== Преобразования координат ==

=== Переход от геодезических координат к геоцентрическим ===

Это преобразование выполняется по следующим формулам:

: <math>\begin{array}{rcl}
x & = & (N + H) \cos B \cos L \\
y & = & (N + H) \cos B \sin L \\
z & = & (N + H - e^2 N) \sin B
\end{array}</math>

Здесь ''N'' — так называемый радиус кривизны первого вертикала:

: <math>N = \frac{a}{1 - e^2 \sin^2 B} = \frac{c}{1 + e'^2 \cos^2 B}</math>

Реализация на Питоне:

<syntaxhighlight lang="python">
def fromLatLong(lat, lon, h, a, f):
    b, c, e2, e12 = initSpher(a, f)
    cos_lat = math.cos(lat)
    n = c / math.sqrt(1. + e12 * cos_lat ** 2)
    p = (n + h) * cos_lat
    x = p * math.cos(lon)
    y = p * math.sin(lon)
    z = (n + h - e2 * n) * math.sin(lat)
    return (x, y, z)
</syntaxhighlight>

=== Переход от геоцентрических координат к геодезическим ===

Проще всего вычисляется долгота:

: <math>\operatorname{tg\,} L = \frac{y}{x}</math>

Сложнее с определением широты и высоты. Существует множество способов решения этой задачи. Воспользуемся итеративным методом Боуринга.

В начале находится предварительная оценка широты ''B'':

: <math>\operatorname{tg\,} B = \frac{z}{r} \left( 1 + e'^2 \frac{b}{\rho} \right)</math>

Здесь ''ρ'' — геоцентрический радиус-вектор, ''r'' — расстояние от оси вращения эллипсоида:

: <math>\rho = \sqrt{x^2 + y^2 + z^2} \qquad r = \sqrt{x^2 + y^2}</math>

Затем вычисляется параметр ''θ'' и получается уточнённое значение широты:

: <math>\begin{array}{rcl}
\operatorname{tg\,} \theta & = & (1 - f) \operatorname{tg\,} B \\
\operatorname{tg\,} B & = & \dfrac{z + e'^2 b \sin^3 \theta}{p - e^2 a \cos^3 \theta}
\end{array}</math>

Действия по последним двум формулам предполагается повторять до сходимости к требуемой точности. Как правило, бывает достаточно одной итерации. В примере реализации метода Боуринга, приведённом ниже, запрограммировано две итерации.

В конце определяется высота:

: <math>H = \frac{r}{\cos B} - N = \frac{z}{\sin B} - N (1 - e^2)</math>

<syntaxhighlight lang="python">
def toLatLong(x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    p = math.hypot(x, y)
    if p == 0.:
        lat = math.copysign(math.pi / 2., z)
        lon = 0.
        h = math.fabs(z) - b
    else:
        t = z / p * (1. + e12 * b / math.hypot(p, z))
        for i in range(2):
            t = t * (1. - f)
            lat = math.atan(t)
            cos_lat = math.cos(lat)
            sin_lat = math.sin(lat)
            t = (z + e12 * b * sin_lat ** 3) / (p - e2 * a * cos_lat ** 3)
        lon = math.atan2(y, x)
        lat = math.atan(t)
        cos_lat = math.cos(lat)
        n = c / math.sqrt(1. + e12 * cos_lat ** 2)
        if math.fabs(t) <= 1.:
            h = p / cos_lat - n
        else:
            h = z / math.sin(lat) - n * (1. - e2)
    return (lat, lon, h)
</syntaxhighlight>

=== Переход от геоцентрических координат к топоцентрическим ===

Постановка задачи: начало топоцентрической системы координат задано точкой ''Q''₀ (''B''₀, ''L''₀, ''H''₀); по геоцентрическим координатам точки ''Q'' (''x'', ''y'', ''z'') вычислить её топоцентрические координаты.

Конформное преобразование между двумя декартовыми прямоугольными системами координат всегда может быть представлено последовательностью сдвигов и вращений координатной системы. Данное преобразование можно реализовать по следующему алгоритму:
* сместить начало координат вдоль оси ''z'' на величину ''e''² ''N''₀ sin ''B''₀ до вершины конуса, образованного нормалями, лежащими на параллели с широтой ''B''₀,
* повернуть систему координат вокруг оси ''z'' на угол ''L''₀, чтобы ось ''x'' оказалась в плоскости меридиана точки ''Q''₀,
* повернуть систему координат вокруг оси ''y'' на угол 90° − ''B''₀, чтобы ось ''z'' совпала с нормалью к поверхности эллипсоида в точке ''Q''₀,
* сместить начало координат вдоль оси ''z'' на величину ''N''₀ + ''H''₀ в точку ''Q''₀,
* изменить знак ''x'' на противоположный.

Реализация алгоритма:

<syntaxhighlight lang="python">
def toTopo(lat0, lon0, h0, x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat0)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    z = z + e2 * n * sin_lat
    x, y = rotate(x, y, lon0)
    z, x = rotate(z, x, math.pi / 2. - lat0)
    z = z - (n + h0)
    x = -x
    return (x, y, z)
</syntaxhighlight>

Функция '''toTopo()''' содержит обращения к функции вращения '''rotate()''':

<syntaxhighlight lang="python">
def rotate(x, y, a):
    c, s = math.cos(a), math.sin(a)
    return (x * c + y * s, -x * s + y * c)
</syntaxhighlight>

=== Переход от топоцентрических координат к геоцентрическим ===

Постановка задачи: начало топоцентрической системы координат задано точкой ''Q''₀ (''B''₀, ''L''₀, ''H''₀); по топоцентрическим координатам точки ''Q'' (''x'', ''y'', ''z'') вычислить её геоцентрические координаты.

Алгоритм решения получается обращением алгоритма обратной задачи:
* изменить знак ''x'' на противоположный,
* сместить начало координат вдоль оси ''z'' на величину ''N''₀ + ''H''₀ в точку пересечения с осью вращения эллипсоида,
* повернуть систему координат вокруг оси ''y'' на угол ''B''₀ − 90°, чтобы ось ''z'' совпала с осью вращения эллипсоида,
* повернуть систему координат вокруг оси ''z'' на угол −''L''₀, чтобы ось ''x'' оказалась в плоскости начального меридиана,
* сместить начало координат вдоль оси ''z'' на величину ''e''² ''N''₀ sin ''B''₀ в центр эллипсоида.

Реализация алгоритма:

<syntaxhighlight lang="python">
def fromTopo(lat0, lon0, h0, x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat0)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    x = -x
    z = z + (n + h0)
    z, x = rotate(z, x, lat0 - math.pi / 2.)
    x, y = rotate(x, y, -lon0)
    z = z - e2 * n * sin_lat
    return (x, y, z)
</syntaxhighlight>

=== Переход от геодезических координат к топоцентрическим. Обратная пространственная задача ===

Постановка задачи: начало топоцентрической системы координат задано точкой ''Q''₀ (''B''₀, ''L''₀, ''H''₀); по геодезическим координатам точки ''Q'' (''B'', ''L'', ''H'') вычислить её топоцентрические координаты ''x'', ''y'', ''z''.

Задача решается последовательным применением готовых алгоритмов:
* по геодезическим координатам точки ''B'', ''L'', ''H'' вычислить её геоцентрические координаты ''x'', ''y'', ''z'',
* по геоцентрическим координатам точки вычислить её топоцентрические координаты ''x'', ''y'', ''z''.

Реализация алгоритма:

<syntaxhighlight lang="python">
def inverse3d(lat0, lon0, h0, lat, lon, h, a, f):
    x, y, z = fromLatLong(lat, lon, h, a, f)
    return toTopo(lat0, lon0, h0, x, y, z, a, f)
</syntaxhighlight>

Рассмотренная задача является разновидностью обратной геодезической задачи в пространстве. Вместо декартовых прямоугольных топоцентрических координат может требоваться вычисление каких-то других связанных с ними величин, например, полярных координат «дальность-азимут-зенитное расстояние», варианты могут быть разные. Однако в большинстве случаев сначала находятся топоцентрические ''x'', ''y'', ''z'', по которым и выводятся искомые значения.

=== Переход от топоцентрических координат к геодезическим. Прямая пространственная задача ===

Постановка задачи: начало топоцентрической системы координат задано точкой ''Q''₀ (''B''₀, ''L''₀, ''H''₀); по топоцентрическим координатам точки ''Q'' (''x'', ''y'', ''z'') вычислить её геодезические координаты ''B'', ''L'', ''H''.

Задача решается через вычисление геоцентрических координат:
* по тороцентрическим координатам точки ''x'', ''y'', ''z'' вычислить её геоцентрические координаты,
* по геоцентрическим координатам точки ''x'', ''y'', ''z'' вычислить её геодезические координаты ''B'', ''L'', ''H''.

Реализация алгоритма:

<syntaxhighlight lang="python">
def forward3d(lat0, lon0, h0, x, y, z, a, f):
    x, y, z = fromTopo(lat0, lon0, h0, x, y, z, a, f)
    return toLatLong(x, y, z, a, f)
</syntaxhighlight>

Эта задача является разновидностью прямой геодезической задачи в пространстве. Вместо декартовых прямоугольных топоцентрических координат могут задаваться какие-то другие связанные с ними величины, например, полярные координаты «дальность-азимут-зенитное расстояние», варианты могут быть разные. Однако в большинстве случаев сначала находятся топоцентрические ''x'', ''y'', ''z'', по которым и решается задача.

== Пример программной реализации ==

Коды вышеприведённых функций находятся в архиве [[Медиа:Spheroid.zip|Spheroid.zip]] в файле '''spheroid.py'''. Напишем программы, которые используют их для преобразования координат.

=== Пересчёт топоцентрических координат в геодезические ===

В этом примере программы явно задаются параметры эллипсоида ''a'', ''f'' и геодезические координаты начала топоцентрической системы ''B''₀, ''L''₀, ''H''₀. Координаты точек ''x'', ''y'', ''z'' читаются из файла данных и пересчитанные значения ''B'', ''L'', ''H'' выводятся в консоль.

<syntaxhighlight lang="python">
from sys import argv
import math
import spheroid

script, fn = argv

a, f = 6378137., 1./298.257223563 # WGS 84

lat0, lon0, hgt0 = math.radians(65.), math.radians(45.), 500.

fp = open(fn, 'r')
for line in fp:
    x, y, z = map(float, line.split(" "))
    lat, lon, hgt = spheroid.forward3d(lat0, lon0, hgt0, x, y, z, a, f)
    print "%.8f %.8f %.3f" % (math.degrees(lat), math.degrees(lon), hgt)
fp.close()
</syntaxhighlight>

Этот скрипт находится в архиве [[Медиа:Spheroid.zip|Spheroid.zip]] в файле '''forwrd3d.py'''.

Файл данных должен содержать в каждой строке координаты одной точки ''x'', ''y'', ''z'', разделённые пробелом. Создадим файл данных '''fwd3d.dat''':

<pre>-40000 30000 0</pre>

Выполним скрипт в командной строке:

<pre>$ python forwrd3d.py fwd3d.dat</pre>

Координаты на выходе:

<pre>64.63992461 45.62743323 695.578</pre>

Запишем полученные координаты в файл результатов '''inv3d.dat''':

<pre>$ python forwrd3d.py fwd3d.dat > inv3d.dat</pre>

=== Пересчёт геодезических координат в топоцентрические ===

В этом примере программы явно задаются параметры эллипсоида ''a'', ''f'' и геодезические координаты начала топоцентрической системы ''B''₀, ''L''₀, ''H''₀. Координаты точек ''B'', ''L'', ''H'' читаются из файла данных и пересчитанные значения ''x'', ''y'', ''z'' выводятся в консоль.

<syntaxhighlight lang="python">
from sys import argv
import math
import spheroid

script, fn = argv

a, f = 6378137., 1./298.257223563 # WGS 84

lat0, lon0, hgt0 = math.radians(65.), math.radians(45.), 500.

fp = open(fn, 'r')
for line in fp:
    lat, lon, hgt = map(float, line.split(" "))
    lat = math.radians(lat)
    lon = math.radians(lon)
    x, y, z = spheroid.inverse3d(lat0, lon0, hgt0, lat, lon, hgt, a, f)
    print "%.3f %.3f %.3f" % (x, y, z)
fp.close()
</syntaxhighlight>

Этот скрипт находится в архиве [[Медиа:Spheroid.zip|Spheroid.zip]] в файле '''invers3d.py'''.

Файл данных должен содержать в каждой строке координаты одной точки ''B'', ''L'', ''H'', разделённые пробелом. Используем в качестве файла данных созданный выше '''inv3d.dat''':

<pre>64.63992461 45.62743323 695.578</pre>

Выполним скрипт в командной строке:

<pre>$ python invers3d.py inv3d.dat</pre>

Координаты на выходе:

<pre>-40000.000 30000.000 0.000</pre>

== Ссылки ==
* [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.139.7504&rep=rep1&type=pdf Robert Burtch, A Comparison of Methods Used in Rectangular to Geodetic Coordinate Transformations, 2006]