TL; DR: пожалуйста, помогите повысить точность расчета между двумя Geodeti c Datums
Мне нужна помощь в решении проблемы с калькулятором Geodeti c, который я построил. Я немного занимаюсь кодированием (для развлечения), но это одна из моих первых "правильных" программ. Это началось как практика для A) алгебры и алгоритмов кодирования и B) я геодезист и хотел вспомнить некоторые теории, которые я изучил в Университете (долгое время go!)
Алгоритм Я использую трехэтапный процесс:
- Преобразование широты / долготы в декартовы координаты Преобразование декартовых координат
- Координаты в новую точку отсчета
- Преобразование декартовых координат обратно в Широта / Долгота
Проблема в том, что полученный результат находится на расстоянии около 65 метров, причем ошибка почти полностью связана с широтой. Я думаю, что ошибка находится где-то в моей математике (хотя я не сбрасываю со счетов ошибку в моем коде), вероятно, в преобразовании между эллипсоидальными и декартовыми координатами. Однако я долго смотрел на это и не могу получить более точную информацию.
Программа разделена на три файла:
- GeodeticConverter - Калькулятор
- ECinPy - Функции для трех шагов
- GenGeo - Методы, которые использовались несколько раз в этом коде и будут использоваться в будущих проектах
Уравнения для шагов 1 и 3 из:
https://gssc.esa.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion
и шаг 2 из:
https://proj.org/operations/transformations/helmert.html
Текущий пример настроен для использования следующей контрольной точки:
- От: WGS 84 - 58 ° 00 '00.31254 "N 000 ° 43' 24.48956" E
- Кому: ED50 - 58 ° 00 '02.59433 "N 000 ° 43' 30.28666" E
Текущий результат:
- 58 ° 00 '00.497437135732071 "N 000 ° 43' 30.286701255138723" E
- C -O: 65 м @ 181,93 °
Код:
Модуль 1 - GeodeticConverter:
#! Python3
"""
Geodetic Calculator that transforms between coordinate systems
"""
from ECinPy import Geo2Cart, Cart2Geo, Transform
from GenGeo import dms2dd, dd2dms, arcs2rad
from math import radians, degrees
# Input Parameters
# From Ellipsoid Parameters
datum_in = "WGS84"
semi_major_in = 6378137.000
flattening_in = 1 / 298.257223563
# To Ellipsoid Parameters
datum_out = "ED50"
semi_major_out = 6378388.000
semi_minor_out = 6356911.946
flattening_out = 1/297.000
# Translation
Tx = 89.5 # m
Ty = 93.8
Tz = 123.1
# Rotation
Rx_in = 0 # acrsec
Ry_in = 0
Rz_in = 0.156
# Scale
s = -1.2
# Latitude
lat_d_in = 58
lat_m_in = 0
lat_s_in = 00.312554
# Longitude
lon_d_in = 0
lon_m_in = 43
lon_s_in = 24.48956
# Height
h_in = 0
# Initial Conversions: Convert DMS to D.ddd Ellipsoidal Coordinates
# Lat & Long from DMS to D.dd
latDD_in = dms2dd(lat_d_in, lat_m_in, lat_s_in)
lonDD_in = dms2dd(lon_d_in, lon_m_in, lon_s_in)
# Lat & Long from Degrees to Radians
latDD = radians(latDD_in)
lonDD = radians(lonDD_in)
# Convert R values from arcsec to Radians
Rx_in = arcs2rad(Rx_in)
Ry_in = arcs2rad(Ry_in)
Rz_in = arcs2rad(Rz_in)
# Step 1: Convert Ellipsoidal to Cartesian Class
g2c = Geo2Cart(semi_major_in, flattening_in, latDD, lonDD, h_in)
# Generate input X, Y, Z
Vax = g2c.x_in()
Vay = g2c.y_in()
Vaz = g2c.z_in()
# Step 2: Transform Cartesian
trans = Transform(Tx, Ty, Tz, Rx_in, Ry_in, Rz_in, s, Vax, Vay, Vaz)
# Run 7 Parameter Helmert Transformation
transformed_coords = trans.trans()
# Cartesian Coordinates in new Datum
x_out = transformed_coords[0]
y_out = transformed_coords[1]
z_out = transformed_coords[2]
# Step 3: Convert Cartesian to Ellipsoidal
c2g = Cart2Geo(semi_major_out, flattening_out, x_out, y_out, z_out)
# Calculate the Lat & Long from Cartesian Coordinates
lat_out, iter_count = c2g.lat_out()
lon_out = c2g.lon_out()
# Convert from Radians to Degrees and then convert D.dd to DMS
lat_out_D, lat_out_M, lat_out_S = dd2dms(degrees(lat_out))
lon_out_D, lon_out_M, lon_out_S = dd2dms(degrees(lon_out))
# Print Results
print(f"lat - {lat_out_D}° {lat_out_M}' {lat_out_S}\"")
print(f"lon - {lon_out_D}° {lon_out_M}' {lon_out_S}\"")
print(f"Latitude function iterated {iter_count} time[s]")
Модуль 2 - ECinPy
"""
Lib to convert Elliptical Coordinates to Cartesian Coordinates and Back
Using the Formulas as provided by:
https://gssc.esa.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion
and
https://proj.org/operations/transformations/helmert.html
"""
from math import sin, cos, atan
from GenGeo import n_fn, e_2_fn, p_fn, h_fn
import numpy as np
class Geo2Cart:
"""
Equations from Ellopsoidal to are:
x = ( N + h ) cos φ cos λ
y = ( N + h ) cos φ sin λ
z = ( (1 − e^2) N + h ) sin φ
"""
def __init__(self, semi_major, flattening, lat, lon, h):
"""
:param semi_major: Of Ellipsoid being converted from
:param flattening: Of Ellipsoid being converted from
:param lat: Of Ellipsoid being converted from
:param lon: Of Ellipsoid being converted from
:param h: Of Ellipsoid being converted from
"""
self.semi_major = semi_major # 'a' in the equations
self.flattening = flattening #
self.lat = lat
self.lon = lon
self.h = h
self.e_2 = e_2_fn(self.flattening)
def x_in(self):
"""
x = ( N + h ) cos φ cos λ
:return: Cartesian X
"""
x = (n_fn(self.semi_major, self.e_2, self.lat) + self.h) \
* cos(self.lat) * cos(self.lon)
return x
def y_in(self):
"""
y = ( N + h ) cos φ sin λ
:return: Cartesian Y
"""
y = (n_fn(self.semi_major, self.e_2, self.lat) + self.h) \
* cos(self.lat) * sin(self.lon)
return y
def z_in(self):
"""
z = ( (1 − e^2) N + h ) sin φ
:return Cartesian Z
"""
z = ((1 - self.e_2) * n_fn(self.semi_major, self.e_2, self.lat)
+ self.h) * sin(self.lat)
return z
class Transform:
"""
Transformation of Cartesian Coordinates into new Datum
7 Parameter (3D) Helmert Transformation:
V(b) = T + (1 + s × 10^-6) R V(a)
"""
def __init__(self, tx, ty, tz, rx, ry, rz, s, vax, vay, vaz):
"""
:param tx: Transformation - X
:param ty: Transformation - Y
:param tz: Transformation - Z
:param rx: Rotation - x
:param ry: Rotation - Y
:param rz: Rotation - Z
:param s: Scale Factor
:param vax: Initial Vector - X
:param vay: Initial Vector - Y
:param vaz: Initial Vector - Z
"""
self.tx = tx
self.ty = ty
self.tz = tz
self.rx = rx
self.ry = ry
self.rz = rz
self.s = s
self.vax = vax
self.vay = vay
self.vaz = vaz
# Arrays
def t_array(self):
"""Transformation Values Array"""
t = np.array([[self.tx],
[self.ty],
[self.tz]])
return t
def r_array(self):
"""Rotation values Array"""
r = np.array([[1, -self.rz, self.ry],
[self.rz, 1, -self.rx],
[-self.ry, self.rx, 1]])
return r
def va_array(self):
"""Initial Values"""
va = np.array([[self.vax],
[self.vay],
[self.vaz]])
return va
def ra(self):
"""Multiply Rotation with initial Values"""
rva = self.r_array().dot(self.va_array()) # array1.dot(array2) is used to multiply array
return rva
def sf(self):
"""Calculate scale factor"""
s = 1 + self.s * pow(10, -6)
return s
def srva(self):
"""Multiply Scale Factor to Rotation Applied Input Vector"""
srva_fn = self.sf() * self.ra()
return srva_fn
def trans(self):
"""Transform Scaled & Rotation Coordinates"""
vb = self.t_array() + self.srva()
return vb
class Cart2Geo:
"""
The ellipsoidal coordinates of a point (φ,λ,h) can be obtained from the cartesian coordinates (x,y,z) as follows:
The longitude λ is given by:
λ = arctan(y/x)
The latitude is computed by an iterative procedure.
1. The initial value is given by:
φ(0) = arctan[ z / ( (1−e^2) p ) ]
with
p = sqrt(x^2 + y^2)
2. Improved values of φ, as well as the height h, are computed iterating in the equations:
N(i) = a / sqrt( (1 − e^2 sin^2 φ(i−1) )
h(i) = ( p / (cos φ(i−1)) ) − N(i)
φ(i) = arctan[ z / ( 1 − e^2 ( N(i) / N(i) + h(i) )p ) ]
"""
def __init__(self, semi_major, flattening, x, y, z):
"""
:param semi_major: Of Ellipsoid being converted to
:param flattening: Of Ellipsoid being converted to
:param x: Cartesian Calculated in Transformation
:param y: Cartesian Calculated in Transformation
:param z: Cartesian Calculated in Transformation
"""
self.semi_major = semi_major # 'a' in the equations
self.flattening = flattening #
self.x = x
self.y = y
self.z = z
self.e_2 = e_2_fn(self.flattening)
self.p = p_fn(self.x, self.y)
def lon_out(self):
"""
λ = arctan(y/x)
:return: λ
"""
lo = atan(self.y / self.x)
return lo
def lat_initial(self):
"""
φ(0) = arctan[ z / ( (1−e^2) p ) ]
p = sqrt(x^2 + y^2)
:return: φ(0)
"""
phi_init = atan(self.z / ((1 - self.e_2) * self.p))
return phi_init
def lat_out(self):
"""
N(i) = a / sqrt( (1 − e^2 sin^2 φ(i−1) )
h(i) = (p / cos φ(i−1)) − N(i)
φ(i) = arctan[ z / ( 1 − e^2 ( N(i) / N(i) + h(i) )p ) ]
Iterates until 5 decimal places of accuracy
:return: φ
"""
# ToDo: List with previous and current iteration value
iter_values = [0, self.lat_initial()]
count = 0
# Run iteration until [n-1] == [n] or the iteration has run 500 times
while iter_values[0] != iter_values[1] and count < 500:
# for a in range(500):
i_1 = iter_values[0]
# Set up N(i) import
n_i = n_fn(self.semi_major, self.e_2, i_1)
# Set up h(i) import
h_i = h_fn(p_fn(self.x, self.y), i_1, n_i)
# Phi formula
phi_i = atan(self.z / (((1 - self.e_2) * (n_i / (n_i + h_i))) * self.p))
# Update iter_values
iter_values[0] = iter_values[1]
iter_values[1] = phi_i
# Format the list to 5 decimal places
iter_values = [float('%.10f' % elem) for elem in iter_values]
# print(f"count = {count + 1} : n-1 = {iter_values[0]} : n = {iter_values[1]}")
count += 1
return iter_values[1], count
Модуль 3 - GenGeo
"""
General equations that are used in Geodetic Calculations
"""
from math import sqrt, sin, cos, radians, degrees
def dms2dd(d_in, m_in, s_in):
"""Convert DMS Coordinate inputs to D.DD"""
dd = d_in + (m_in / 60) + (s_in / 3600)
return dd
def dd2dms(dd_in):
"""Convert D.DD to DMS"""
d_out = int(dd_in)
m_out = int((dd_in - d_out) * 60)
# s_out = (((dd_in - d_out) * 60) - int((dd_in - d_out) * 60)) * 60
s_out = 3600 * (dd_in - d_out) - 60 * m_out
return d_out, m_out, s_out
def arcs2rad(arcsec):
"""Converts ArcSec In to Radians"""
rads = radians(arcsec / 3600)
return rads
def e_2_fn(f):
"""Returns eccentricity of ellipsoid"""
e_2 = 2 * f - f ** 2
return e_2
def n_fn(a, e, lat):
"""
N = a / sqrt( 1 − e^2 sin^2 φ)
Where:
a = Semi-Major
e^2 = eccentricity
φ = latitude
"""
n_1 = a / sqrt(1 - e * (sin(lat)) ** 2)
return n_1
def h_fn(p_in, lat, n_in):
"""Calculates h for the iteration function"""
h = (p_in / cos(lat)) - n_in
return h
def p_fn(x, y):
"""Returns the value of p"""
p = (sqrt(pow(x, 2) + pow(y, 2)))
return p
def radius_at_lat(lat):
"""Calculates radius of Earth at specific latitude"""
ra = 6_378_137.0 # Equatorial Radius
rb = 6_356_752.3 # Polar Radius
lat = radians(lat)
rp = sqrt(
(pow(pow(ra, 2) * cos(lat), 2) +
pow(pow(rb, 2) * sin(lat), 2))
/
(pow(ra * cos(lat), 2) +
pow(rb * sin(lat), 2))
)
return rp
Мне нужно, чтобы это было точно для преобразования 0,01 м между любыми двумя датумами. Если вы поможете мне добраться туда, я буду исключительно благодарен.
Ура,