he2851/Masterprojekt
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Koordinatenumrechnungen funktionieren inkl. Randfälle, GHA2_num funktioniert mit Standard Ellipsoid

Browse Source
This commit is contained in:
Hendrik Möllersmann
2025-11-26 11:05:18 +01:00
parent 9031a12312
commit d76859d17b
5 changed files with 181 additions and 125 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
GHA_triaxial/approx_louville.py
View File

@@ -54,8 +54,8 @@ def num_update(ell: EllipsoidTriaxial, points, diffs):
def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int):
beta1, lamb1 = ell.cart2ell2(p1)
beta2, lamb2 = ell.cart2ell2(p2)
beta1, lamb1 = ell.cart2ell(p1)
beta2, lamb2 = ell.cart2ell(p2)
points = points_approx_gha2(ell.ax, beta1, lamb1, beta2, lamb2, 5)
for j in range(maxI):
@@ -82,7 +82,7 @@ def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int):
def show_points(ell: EllipsoidTriaxial, points):
points_cart = []
for point in points:
points_cart.append(ell.ell2cart2(point[0], point[1]))
points_cart.append(ell.ell2cart(point[0], point[1]))
points_cart = np.array(points_cart)
fig = plt.figure()

12
GHA_triaxial/panou.py
View File

@@ -152,16 +152,16 @@ if __name__ == "__main__":
lamb1 = wu.deg2rad(0)
beta2 = wu.deg2rad(60)
lamb2 = wu.deg2rad(175)
P1 = ell.ell2cart2(wu.deg2rad(60), wu.deg2rad(0))
P2 = ell.ell2cart2(wu.deg2rad(60), wu.deg2rad(175))
P1 = ell.ell2cart(wu.deg2rad(60), wu.deg2rad(0))
P2 = ell.ell2cart(wu.deg2rad(60), wu.deg2rad(175))
para1 = ell.cart2para(P1)
para2 = ell.cart2para(P2)
cart1 = ell.para2cart(para1[0], para1[1])
cart2 = ell.para2cart(para2[0], para2[1])
ell11 = ell.cart2ell2(P1)
ell21 = ell.cart2ell2(P2)
ell1 = ell.cart2ell2(cart1)
ell2 = ell.cart2ell2(cart2)
ell11 = ell.cart2ell(P1)
ell21 = ell.cart2ell(P2)
ell1 = ell.cart2ell(cart1)
ell2 = ell.cart2ell(cart2)
c = 0.06207487624
alpha0 = wu.gms2rad([2, 52, 26.2393])

67
GHA_triaxial/panou_2013_2GHA_num.py
View File

@@ -1,68 +1,9 @@
import numpy as np
#from ellipsoide import EllipsoidTriaxial
from ellipsoide import EllipsoidTriaxial
import Numerische_Integration.num_int_runge_kutta as rk
import ausgaben as aus
# Panou 2013
class EllipsoidTriaxial:
def __init__(self, ax: float, ay: float, b: float):
self.ax = ax
self.ay = ay
self.b = b
self.ex = np.sqrt((self.ax**2 - self.b**2) / self.ax**2)
self.ey = np.sqrt((self.ay**2 - self.b**2) / self.ay**2)
self.ee = np.sqrt((self.ax**2 - self.ay**2) / self.ax**2)
self.ex_ = np.sqrt((self.ax**2 - self.b**2) / self.b**2)
self.ey_ = np.sqrt((self.ay**2 - self.b**2) / self.b**2)
self.ee_ = np.sqrt((self.ax**2 - self.ay**2) / self.ay**2)
self.Ex = np.sqrt(self.ax**2 - self.b**2)
self.Ey = np.sqrt(self.ay**2 - self.b**2)
self.Ee = np.sqrt(self.ax**2 - self.ay**2)
@classmethod
def init_name(cls, name: str):
if name == "BursaFialova1993":
ax = 6378171.36
ay = 6378101.61
b = 6356751.84
return cls(ax, ay, b)
elif name == "BursaSima1980":
ax = 6378172
ay = 6378102.7
b = 6356752.6
return cls(ax, ay, b)
elif name == "BursaSima1980round":
# Panou 2013
ax = 6378172
ay = 6378103
b = 6356753
return cls(ax, ay, b)
elif name == "Eitschberger1978":
ax = 6378173.435
ay = 6378103.9
b = 6356754.4
return cls(ax, ay, b)
elif name == "Bursa1972":
ax = 6378173
ay = 6378104
b = 6356754
return cls(ax, ay, b)
elif name == "Bursa1970":
ax = 6378173
ay = 6378105
b = 6356754
return cls(ax, ay, b)
elif name == "Bessel-biaxial":
ax = 6377397.15509
ay = 6377397.15508
b = 6356078.96290
return cls(ax, ay, b)
def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30):
"""
@@ -383,8 +324,8 @@ if __name__ == "__main__":
beta2 = np.deg2rad(75)
lamb2 = np.deg2rad(66)
a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
print(np.rad2deg(a1))
print(np.rad2deg(a2))
print(aus.gms("a1", a1, 4))
print(aus.gms("a2", a2, 4))
print(s)

190
ellipsoide.py
View File

@@ -119,6 +119,12 @@ class EllipsoidTriaxial:
b = 6356751.84
return cls(ax, ay, b)
elif name == "BursaSima1980":
ax = 6378172
ay = 6378102.7
b = 6356752.6
return cls(ax, ay, b)
elif name == "BursaSima1980round":
# Panou 2013
ax = 6378172
ay = 6378103
b = 6356753
@@ -161,7 +167,7 @@ class EllipsoidTriaxial:
else:
return False
def ell2cart(self, beta: float, lamb: float, u: float) -> np.ndarray:
def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray:
"""
Panou 2014 12ff.
Ellipsoidische Breite+Länge sind nicht gleich der geodätischen
@@ -190,21 +196,34 @@ class EllipsoidTriaxial:
return np.array([x, y, z])
def ell2cart2(self, beta: float, lamb: float) -> np.ndarray:
def ell2cart(self, beta: float, lamb: float) -> np.ndarray:
"""
Panou, Korakitis 2019 2
:param beta: ellipsoidische Breite [rad]
:param lamb: ellipsoidische Länge [rad]
:return: Punkt in kartesischen Koordinaten
"""
B = self.Ex**2 * np.cos(beta)**2 + self.Ee**2 * np.sin(beta)**2
L = self.Ex**2 - self.Ee**2 * np.cos(lamb)**2
x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb)
y = self.ay * np.cos(beta) * np.sin(lamb)
z = self.b / self.Ex * np.sin(beta) * np.sqrt(L)
return np.array([x, y, z])
if beta == -np.pi/2:
return np.array([0, 0, -self.b])
elif beta == np.pi/2:
return np.array([0, 0, self.b])
elif beta == 0 and lamb == -np.pi/2:
return np.array([0, -self.ay, 0])
elif beta == 0 and lamb == np.pi/2:
return np.array([0, self.ay, 0])
elif beta == 0 and lamb == 0:
return np.array([self.ax, 0, 0])
elif beta == 0 and lamb == np.pi:
return np.array([-self.ax, 0, 0])
else:
B = self.Ex**2 * np.cos(beta)**2 + self.Ee**2 * np.sin(beta)**2
L = self.Ex**2 - self.Ee**2 * np.cos(lamb)**2
x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb)
y = self.ay * np.cos(beta) * np.sin(lamb)
z = self.b / self.Ex * np.sin(beta) * np.sqrt(L)
return np.array([x, y, z])
def cart2ell(self, point: np.ndarray) -> tuple[float, float, float]:
def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]:
"""
Panou 2014 15ff.
:param point: Punkt in kartesischen Koordinaten
@@ -232,21 +251,65 @@ class EllipsoidTriaxial:
return beta, lamb, u
def cart2ell2(self, point: np.ndarray) -> tuple[float, float]:
def cart2ell(self, point: np.ndarray) -> tuple[float, float]:
"""
Panou, Korakitis 2019 2f.
:param point: Punkt in kartesischen Koordinaten
:return: ellipsoidische Breite, ellipsoidische Länge
"""
x, y, z = point
eps = 1e-9
if abs(x) < eps and abs(y) < eps: # Punkt in der z-Achse
beta = np.pi / 2 if z > 0 else -np.pi / 2
lamb = 0.0
return beta, lamb
elif abs(x) < eps and abs(z) < eps: # Punkt in der y-Achse
beta = 0.0
lamb = np.pi / 2 if y > 0 else -np.pi / 2
return beta, lamb
elif abs(y) < eps and abs(z) < eps: # Punkt in der x-Achse
beta = 0.0
lamb = 0.0 if x > 0 else np.pi
return beta, lamb
# ---- Allgemeiner Fall -----
c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
(self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
self.ax ** 2 + self.ay ** 2) * z ** 2)
t2 = (-c1 + np.sqrt(c1 ** 2 - 4 * c0)) / 2
t1 = c0 / t2
beta = np.arctan(np.sqrt((t1 - self.b**2) / (self.ay**2 - t1)))
lamb = np.arctan(np.sqrt((t2 - self.ay**2) / (self.ax**2 - t2)))
num_beta = max(t1 - self.b ** 2, 0)
den_beta = max(self.ay ** 2 - t1, 0)
num_lamb = max(t2 - self.ay ** 2, 0)
den_lamb = max(self.ax ** 2 - t2, 0)
beta = np.arctan(np.sqrt(num_beta / den_beta))
lamb = np.arctan(np.sqrt(num_lamb / den_lamb))
if z < 0:
beta = -beta
if y < 0:
lamb = -lamb
if x < 0:
lamb = np.pi - lamb
if abs(x) < eps:
lamb = -np.pi/2 if y < 0 else np.pi/2
elif abs(y) < eps:
lamb = 0 if x > 0 else np.pi
elif abs(z) < eps:
beta = 0
return beta, lamb
def cart2geod(self, mode: str, point: np.ndarray, maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
@@ -259,6 +322,45 @@ class EllipsoidTriaxial:
:return: phi, lambda, h
"""
xG, yG, zG = point
eps = 1e-9
if abs(xG) < eps and abs(yG) < eps: # Punkt in der z-Achse
phi = np.pi / 2 if zG > 0 else -np.pi / 2
lamb = 0.0
h = abs(zG) - ell.b
return phi, lamb, h
elif abs(xG) < eps and abs(zG) < eps: # Punkt in der y-Achse
phi = 0.0
lamb = np.pi / 2 if yG > 0 else -np.pi / 2
h = abs(yG) - ell.ay
return phi, lamb, h
elif abs(yG) < eps and abs(zG) < eps: # Punkt in der x-Achse
phi = 0.0
lamb = 0.0 if xG > 0 else np.pi
h = abs(xG) - ell.ax
return phi, lamb, h
# elif abs(zG) < eps: # Punkt in der xy-Ebene
# phi = 0
# lamb = np.arctan2(yG / ell.ay**2, xG / ell.ax**2)
# rG = np.sqrt(xG ** 2 + yG ** 2)
# pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
# rE = np.sqrt(pE[0] ** 2 + pE[1] ** 2)
# h = rG - rE
# return phi, lamb, h
#
# elif abs(yG) < eps: # Punkt in der xz-Ebene
# phi = np.arctan2(zG / ell.b**2, xG / ell.ax**2)
# lamb = 0 if xG > 0 else np.pi
# rG = np.sqrt(xG ** 2 + zG ** 2)
# pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
# rE = np.sqrt(pE[0] ** 2 + pE[2] ** 2)
# h = rG - rE
# return phi, lamb, h
rG = np.sqrt(xG ** 2 + yG ** 2 + zG ** 2)
pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
@@ -286,6 +388,15 @@ class EllipsoidTriaxial:
lamb = np.arctan(1/(1-self.ee**2) * pE[1]/pE[0])
h = np.sign(zG - pE[2]) * np.sign(pE[2]) * np.sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2)
if xG < 0 and yG < 0:
lamb = -np.pi + lamb
elif xG < 0:
lamb = np.pi + lamb
if abs(zG) < eps:
phi = 0
return phi, lamb, h
def geod2cart(self, phi: float, lamb: float, h: float) -> np.ndarray:
@@ -342,12 +453,13 @@ class EllipsoidTriaxial:
:return: parametrische Koordinaten
"""
x, y, z = point
u_check1 = z*np.sqrt(1 - self.ee**2)
u_check2 = np.sqrt(x**2 * (1-self.ee**2) + y**2) * np.sqrt(1-self.ex**2)
if u_check1 <= u_check2:
u = np.arctan2(u_check1, u_check2)
else:
u = np.pi/2 * np.arctan2(u_check2, u_check1)
u = np.pi/2 - np.arctan2(u_check2, u_check1)
v_check1 = y
v_check2 = x*np.sqrt(1-self.ee**2)
@@ -356,6 +468,7 @@ class EllipsoidTriaxial:
v = 2 * np.arctan2(v_check1, v_check2 + v_factor)
else:
v = np.pi/2 - 2 * np.arctan2(v_check2, v_check1 + v_factor)
return u, v
def p_q(self, x, y, z) -> dict:
@@ -370,7 +483,7 @@ class EllipsoidTriaxial:
n = np.array([x / np.sqrt(H), y / ((1 - self.ee ** 2) * np.sqrt(H)), z / ((1 - self.ex ** 2) * np.sqrt(H))])
beta, lamb, u = self.cart2ell(np.array([x, y, z]))
beta, lamb, u = self.cart2ellu(np.array([x, y, z]))
B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2
L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2
@@ -403,37 +516,34 @@ class EllipsoidTriaxial:
if __name__ == "__main__":
ell = EllipsoidTriaxial.init_name("Eitschberger1978")
diff_list = []
for beta_deg in range(-150, 210, 30):
for lamb_deg in range(-150, 210, 30):
beta = wu.deg2rad(beta_deg)
lamb = wu.deg2rad(lamb_deg)
point = ell.ell2cart(beta, lamb)
point = np.array([5672455, 2698193, 1103177])
pointell, _, _, _ = ell.cartonell(point)
values1 = ell.cart2ell(point)
values2 = ell.cart2ell2(pointell)
carts1 = ell.ell2cart(values2[0], values2[1], ell.b)
carts2 = ell.ell2cart2(values2[0], values2[1])
elli = ell.cart2ell(point)
cart_elli = ell.ell2cart(elli[0], elli[1])
diff_ell = np.sum(np.abs(point-cart_elli))
para = ell.cart2para(point)
cart_para = ell.para2cart(para[0], para[1])
diff_para = np.sum(np.abs(point-cart_para))
pass
geod = ell.cart2geod("ligas1", point)
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
diff_geod1 = np.sum(np.abs(point-cart_geod))
# cart_x, cart_y, cart_z = np.array([5672455, 2698193, 1103177])
# cart_x, cart_y, cart_z = np.array([3415031.1337320395, 3414993.9029805865, 588647.548260936])
# cart_x, cart_y, cart_z = np.array([6378173.435, 0, 0])
# cart_x, cart_y, cart_z = np.array([0, 6378103.9, 0])
# cart_x, cart_y, cart_z = np.array([0, 0, 6356754.4])
cart_x, cart_y, cart_z = np.array([1000, 1000, 1000])
print(aus.xyz(cart_x, cart_y, cart_z, 10), " Startwerte")
geod = ell.cart2geod("ligas2", point)
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
diff_geod2 = np.sum(np.abs(point-cart_geod))
beta, lamb_ell, u = ell.cart2ell(np.array([cart_x, cart_y, cart_z]))
phi, lamb_geod, h = ell.cart2geod("ligas3", np.array([cart_x, cart_y, cart_z]))
u_para, v_para = ell.cart2para(np.array([cart_x, cart_y, cart_z]))
geod = ell.cart2geod("ligas3", point)
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
diff_geod3 = np.sum(np.abs(point-cart_geod))
e_x, e_y, e_z = ell.ell2cart(beta, lamb_ell, u)
diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod1, diff_geod2, diff_geod3])
print(aus.xyz(e_x, e_y, e_z, 10), f", Distanz ellipsoidisch: {np.linalg.norm(np.array([cart_x, cart_y, cart_z]) - np.array([e_x, e_y, e_z]))} m")
g_x, g_y, g_z = ell.geod2cart(phi, lamb_geod, h)
print(aus.xyz(g_x, g_y, g_z, 10), f", Distanz geodätisch: {np.linalg.norm(np.array([cart_x, cart_y, cart_z]) - np.array([g_x, g_y, g_z]))} m")
p_x, p_y, p_z = ell.para2cart(u_para, v_para)
print(aus.xyz(p_x, p_y, p_z, 10), f", Distanz parametrisch: {np.linalg.norm(np.array([cart_x, cart_y, cart_z]) - np.array([p_x, p_y, p_z]))} m")
pass
carts = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(20), 6356754.4)
diff_list = np.array(diff_list)
pass

31
jacobian_Ligas.py
View File

@@ -35,13 +35,16 @@ def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
j33 = (pE[1] - pG[1]) * G - F * pE[1]
detJ = j11 * j22 * j33 - j21 * j12 * j33 + j21 * j13 * j32
invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22],
[-j21*j33, j11*j33, j13*j21],
[j21*j32, -j11*j32, j11*j22-j12*j21]])
if detJ == 0:
invJ, fxE = case3(E, F, G, pG, pE)
else:
invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22],
[-j21*j33, j11*j33, j13*j21],
[j21*j32, -j11*j32, j11*j22-j12*j21]])
fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
(pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0],
(pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
(pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0],
(pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
return invJ, fxE
@@ -57,13 +60,15 @@ def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
j33 = (pE[1] - pG[1]) * G - F * pE[1]
detJ = -j11 * j23 * j32 - j21 * j12 * j33 + j21 * j13 * j32
if detJ == 0:
invJ, fxE = case2(E, F, G, pG, pE)
else:
invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23],
[-j21*j33, j11*j33, -(j11*j23-j13*j21)],
[j21*j32, -j11*j32, -j12*j21]])
invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23],
[-j21*j33, j11*j33, -(j11*j23-j13*j21)],
[j21*j32, -j11*j32, -j12*j21]])
fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
(pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0],
(pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
(pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0],
(pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
return invJ, fxE