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

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

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

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

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

1. Геоцентрические декартовы прямоугольные координаты:

* Начало координат находится в центре эллипсоида,
* ось z расположена вдоль оси вращения эллипсоида и направлена в северный полюс,
* ось x лежит в пересечении экватора и начального меридиана,
* ось y лежит в пересечении экватора и меридиана с долготой L = 0.

• B
• L
• 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 initSpher(a, f):
    b = a * (1. - f)
    c = a / (1. - f)
    e2 = f * (2. - f)
    e12 = e2 / (1. - e2)
    return (b, c, e2, e12)

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)

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 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.63992455 45.62743333 695.578

-40000.000 30000.000 0.000

Ссылки

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