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

Эта страница является черновиком статьи.

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

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

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

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

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle f = \frac{a - b}{a} \qquad b = a (1 - f)}

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle c = \frac{a}{1 - f} \qquad e^2 = f (2 - f) \qquad e'^2 = \frac{e^2}{1 - e^2}}

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

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)

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

Рассмотрим следующие системы координат.

1. Геоцентрические декартовы прямоугольные координаты:
   • начало координат находится в центре эллипсоида,
   • ось z расположена вдоль оси вращения эллипсоида и направлена в северный полюс,
   • ось x лежит в пересечении экватора и начального меридиана,
   • ось y лежит в пересечении экватора и меридиана с долготой L = 90°.
2. Система геодезических координат:

   геодезическая широта B

   угол между нормалью к поверхности эллипсоида и плоскостью экватора,

   геодезическая долгота L

   двугранный угол между плоскостями данного и начального меридианов,

   геодезическая высота H

   кратчайшее расстояние до поверхности эллипсоида.

3. Топоцентрические декартовы прямоугольные координаты:
   • начало координат находится в некоторой точке Q₀ (B₀, L₀, H₀) над эллипсоидом,
   • ось z расположена вдоль нормали к поверхности эллипсоида и направлена вверх,
   • ось x расположена в плоскости меридиана и направлена на север,
   • ось y перпендикулярна к осям x и z и направлена на восток.

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

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

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \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}}

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle N = \frac{a}{1 - e^2 \sin^2 B} = \frac{c}{1 + e'^2 \cos^2 B}}

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

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)

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

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \operatorname{tg\,} L = \frac{y}{x}}

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

В начале находится предварительная оценка параметра θ:

$\operatorname {tg\,} \theta ={\frac {z}{r}}\left(1+e'^{2}{\frac {b}{\rho }}\right)(1-f)$

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \rho = \sqrt{x^2 + y^2 + z^2} \qquad r = \sqrt{x^2 + y^2}}

Затем вычисляется тангенс широты:

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \operatorname{tg\,} B = \frac{z + e'^2 b \sin^3 \theta}{p - e^2 a \cos^3 \theta}}

и параметр θ уточняется:

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \operatorname{tg\,} \theta = (1 - f) \operatorname{tg\,} B}

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

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

Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle H = \frac{r}{\cos B} - N = \frac{z}{\sin B} - N (1 - e^2)}

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)

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

Конформное преобразование между двумя декартовыми прямоугольными системами координат всегда может быть представлено последовательностью сдвигов и вращений координатной системы. Д

def toTopo(lat, lon, h, x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    z = z + e2 * n * sin_lat
    x, y = rotate(x, y, lon)
    z, x = rotate(z, x, math.pi / 2. - lat)
    z = z - (n + h)
    x = -x
    return (x, y, z)

Для реализации алгоритма понадобится функция вращения:

def rotate(x, y, a):
    c, s = math.cos(a), math.sin(a)
    return (x * c + y * s, -x * s + y * c)

Программы

import math

def rotate(x, y, a):
    c, s = math.cos(a), math.sin(a)
    return (x * c + y * s, -x * s + y * c)

def forward3d(lat1, lon1, h1, x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat1)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    x = -x
    z = z + (n + h1)
    z, x = rotate(z, x, lat1 - math.pi / 2.)
    x, y = rotate(x, y, -lon1)
    z = z - e2 * n * sin_lat
    return toLatLong(x, y, z, a, f)

def inverse3d(lat1, lon1, h1, lat2, lon2, h2, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat1)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    x, y, z = fromLatLong(lat2, lon2, h2, a, f)
    z = z + e2 * n * sin_lat
    x, y = rotate(x, y, lon1)
    z, x = rotate(z, x, math.pi / 2. - lat1)
    z = z - (n + h1)
    x = -x
    return (x, y, z)

# -*- coding: utf-8 -*-
from sys import argv
import math
import latlong

script, fn = argv

a, f = 6378136., 1./298.25784 # ПЗ-90

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 = latlong.forward3d(lat0, lon0, hgt0, x, y, z, a, f)
    print "%.8f %.8f %.3f" % (math.degrees(lat), math.degrees(lon), hgt)
fp.close()

# -*- coding: utf-8 -*-
from sys import argv
import math
import latlong

script, fn = argv

a, f = 6378136., 1./298.25784 # ПЗ-90

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 = latlong.inverse3d(lat0, lon0, hgt0, lat, lon, hgt, a, f)
    print "%.3f %.3f %.3f" % (x, y, z)
fp.close()

Расчёты

-40000 30000 0

64.63992461 45.62743323 695.578

-40000.000 30000.000 0.000

Ссылки

• Robert Burtch, A Comparison of Methods Used in Rectangular to Geodetic Coordinate Transformations, 2006