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

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

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

Программы

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