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

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

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

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

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

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

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

$f={\frac {a-b}{a}}\qquad b=a(1-f)$

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

$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 и направлена на восток.

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

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

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

${\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 — так называемый радиус кривизны первого вертикала:

$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)

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

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

$\operatorname {tg\,} L={\frac {y}{x}}$

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

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

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

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

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}}\qquad r={\sqrt {x^{2}+y^{2}}}$

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

$\operatorname {tg\,} B={\frac {z+e'^{2}b\sin ^{3}\theta }{p-e^{2}a\cos ^{3}\theta }}$

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

$\operatorname {tg\,} \theta =(1-f)\operatorname {tg\,} B$

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

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

$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)

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

Обратная задача.

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

Прямая задача.

Программы

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