zusammenfügen 02.2.

This commit is contained in:
2026-02-02 20:01:37 +01:00
parent 774acc3854
commit 4e7c55500b
11 changed files with 62593 additions and 2121 deletions

View File

@@ -13,11 +13,56 @@ from pyproj import CRS, Transformer
class Transformationen:
"""Koordinatentransformationen und Helmert-Transformation (Euler-Winkel) zwischen Referenzsystemen.
Die Klasse stellt Methoden zur Verfügung für:
- Aufbau einer Rotationsmatrix aus Eulerwinkeln,
- Schätzung von 7-Parameter-Transformationsparametern (dX, dY, dZ, Maßstab m, Eulerwinkel e1/e2/e3)
aus identischen Punkten zwischen lokalem Horizontsystem (LH) und geozentrisch-kartesischem System (ECEF),
- Anwendung der geschätzten Helmerttransformation auf Punkte, die nur im Ausgangssystem vorliegen,
- Transformation zwischen ETRS89 / UTM (+ DHHN2016 Normalhöhe) und ETRS89 geozentrisch-kartesisch (ECEF),
inkl. Nutzung einer BKG-Quasigeoidundulations-Datei (GeoTIFF) für PROJ.
Die grundlegende Funktionsweise der Transformationsschätzung lautet:
1) Identische Punkte aus Ausgangs- und Zielsystem ermitteln.
2) Näherung für Maßstab m0 aus mittleren Streckenverhältnissen bilden.
3) Näherungs-Rotation R0 aus lokalen Basen (u/v/w und U/V/W) bestimmen und daraus Euler-Näherungen ableiten.
4) Iterative Parameterschätzung (Gauss-Newton) auf Basis der Beobachtungsgleichung:
P = T + m * R(e1,e2,e3) * p
"""
def __init__(self, pfad_datenbank: str) -> None:
"""Initialisiert die Transformationsklasse.
Speichert den Pfad zur SQLite-Datenbank und initialisiert den Datenbankzugriff.
:param pfad_datenbank: Pfad zur SQLite-Datenbank.
:type pfad_datenbank: str
:return: None
:rtype: None
"""
self.pfad_datenbank = pfad_datenbank
self.db_zugriff = Datenbank.Datenbankzugriff(self.pfad_datenbank)
@staticmethod
def R_matrix_aus_euler(e1: float, e2: float, e3: float) -> sp.Matrix:
def R_matrix_aus_eulerwinkeln(e1: float, e2: float, e3: float) -> sp.Matrix:
"""Erstellt eine 3x3-Rotationsmatrix aus Eulerwinkeln.
Die Rotationsmatrix wird symbolisch (SymPy) aufgebaut. Die Eulerwinkel e1, e2, e3 werden
direkt in trigonometrische Ausdrücke eingesetzt und eine orthogonale Rotationsmatrix R(e1,e2,e3)
zur Verwendung in Helmert-Transformationen zurückgegeben.
:param e1: Eulerwinkel 1 (Radiant).
:type e1: float
:param e2: Eulerwinkel 2 (Radiant).
:type e2: float
:param e3: Eulerwinkel 3 (Radiant).
:type e3: float
:return: Rotationsmatrix R als SymPy-Matrix (3x3).
:rtype: sp.Matrix
"""
return sp.Matrix([
[
sp.cos(e2) * sp.cos(e3),
@@ -36,32 +81,48 @@ class Transformationen:
]
])
def Helmerttransformation_Euler_Transformationsparameter_berechne(self) -> dict[Any, float]:
db = Datenbank.Datenbankzugriff(self.pfad_datenbank)
dict_ausgangssystem = db.get_koordinaten("naeherung_lh", "Dict")
dict_zielsystem = db.get_koordinaten("naeherung_us", "Dict")
def Helmerttransformation_Euler_Transformationsparameter_berechnen(self) -> dict[Any, float]:
"""Schätzt die Helmert-Transformationsparameter aus identischen Punkten.
Aus der Datenbank werden Näherungskoordinaten des lokalen Horizontsystems (naeherung_lh)
und des geozentrisch-kartesischen Systems (naeherung_us) geladen. Für die Schnittmenge der
Punkte werden die 7 Helmertparameter geschätzt:
- Translation (dX, dY, dZ),
- Maßstab m,
- Eulerwinkel (e1, e2, e3).
Näherungen:
- m0: Mittelwert der Streckenverhältnisse aus allen Punktpaaren,
- R0: Anfangsrotationsmatrix aus lokalen Basisvektoren (u/v/w und U/V/W),
- Translation0: aus Schwerpunkten und m0/R0.
Die Parameterschätzung erfolgt iterativ mit P = I.
Abbruchkriterium: |dx_i| < schwellenwert in zwei aufeinanderfolgenden Iterationen oder max. 100 Iterationen.
:param: None
:return: Dictionary der finalen Parameter mit SymPy-Symbolen als Keys und float-Werten als Values
(Keys: dX, dY, dZ, m, e1, e2, e3).
:rtype: dict[Any, float]
"""
# Koordinaten des lokalen Horizontsystems des Tachymeters und der geozentrisch Kartesischen Näherungskoordinaten aus den statischen GNSS-Messungen aus der Tabelle Netzpunkte abfragen
dict_ausgangssystem = self.db_zugriff.get_koordinaten("naeherung_lh", "Dict")
dict_zielsystem = self.db_zugriff.get_koordinaten("naeherung_us", "Dict")
# Identische Punkte ermitteln
gemeinsame_punktnummern = sorted(set(dict_ausgangssystem.keys()) & set(dict_zielsystem.keys()))
anzahl_gemeinsame_punkte = len(gemeinsame_punktnummern)
liste_punkte_ausgangssystem = [dict_ausgangssystem[i] for i in gemeinsame_punktnummern]
liste_punkte_zielsystem = [dict_zielsystem[i] for i in gemeinsame_punktnummern]
print("Anzahl gemeinsame Punkte:", anzahl_gemeinsame_punkte)
print("Anzahl verwendete Punkte für die Helmerttransformation:", anzahl_gemeinsame_punkte)
print("\nErste Zielpunkte:")
for pn, P in list(zip(gemeinsame_punktnummern, liste_punkte_zielsystem))[:5]:
print(pn, [float(P[0]), float(P[1]), float(P[2])])
print("\nErste Ausgangspunkte:")
for pn, p in list(zip(gemeinsame_punktnummern, liste_punkte_ausgangssystem))[:5]:
print(pn, [float(p[0]), float(p[1]), float(p[2])])
# --- Näherungswerte (minimal erweitert) ---
p1, p2, p3 = liste_punkte_ausgangssystem[0], liste_punkte_ausgangssystem[1], liste_punkte_ausgangssystem[2]
P1, P2, P3 = liste_punkte_zielsystem[0], liste_punkte_zielsystem[1], liste_punkte_zielsystem[2]
# 1) Näherungswert Maßstab: Mittelwert aus allen Punktpaaren
# Näherungswert für dem Maßstab berechnen aus dem Mittelwert aller Punktpaare
ratios = []
for i, j in combinations(range(anzahl_gemeinsame_punkte), 2):
dp = (liste_punkte_ausgangssystem[j] - liste_punkte_ausgangssystem[i]).norm()
@@ -72,12 +133,7 @@ class Transformationen:
m0 = sum(ratios) / len(ratios)
if ratios:
print("min/mean/max:",
min(ratios),
sum(ratios) / len(ratios),
max(ratios))
# Näherungswert für die Translation berechnen
U = (P2 - P1) / (P2 - P1).norm()
W = (U.cross(P3 - P1)) / (U.cross(P3 - P1)).norm()
V = W.cross(U)
@@ -93,28 +149,18 @@ class Transformationen:
Translation0 = XS - m0 * R0 * xS
# 2) Test auf orthonormale Drehmatrix bei 3 Nachkommastellen!
if R0.T.applyfunc(lambda x: round(float(x), 3)) == R0.inv().applyfunc(lambda x: round(float(x), 3)) \
and (R0.T * R0).applyfunc(lambda x: round(float(x), 3)) == sp.eye(3).applyfunc(
lambda x: round(float(x), 3)) \
and ((round(R0.det(), 3) == 1.000 or round(R0.det(), 3) == -1.000)):
print("R ist Orthonormal!")
else:
print("R ist nicht Orthonormal!")
# 3) Euler-Näherungswerte aus R0
# Euler-Näherungswerte aus der Anfangsrotationsmatrix
e2_0 = sp.asin(R0[2, 0])
# Schutz gegen Division durch 0 wenn cos(e2) ~ 0:
cos_e2_0 = sp.cos(e2_0)
e1_0 = sp.acos(R0[2, 2] / cos_e2_0)
e3_0 = sp.acos(R0[0, 0] / cos_e2_0)
# --- Symbolische Unbekannte (klassische 7 Parameter) ---
# Symbolische Unbekannte
dX, dY, dZ, m, e1, e2, e3 = sp.symbols('dX dY dZ m e1 e2 e3')
R_symbolisch = self.R_matrix_aus_euler(e1, e2, e3)
R_symbolisch = self.R_matrix_aus_eulerwinkeln(e1, e2, e3)
# 4) Funktionales Modell
# Funktionales Modell
f_zeilen = []
for punkt in liste_punkte_ausgangssystem:
punkt_vektor = sp.Matrix([punkt[0], punkt[1], punkt[2]])
@@ -134,7 +180,7 @@ class Transformationen:
l_vektor = sp.Matrix([koord for P in liste_punkte_zielsystem for koord in P])
l = l_vektor
P_mat = sp.eye(3 * anzahl_gemeinsame_punkte)
P_matrix = sp.eye(3 * anzahl_gemeinsame_punkte)
l_berechnet_0 = None
while True:
@@ -156,8 +202,8 @@ class Transformationen:
dl_0 = l_vektor - l_berechnet_0
A_0 = A_ohne_zahlen.subs(zahlen_0).evalf(n=30)
N = A_0.T * P_mat * A_0
n_0 = A_0.T * P_mat * dl_0
N = A_0.T * P_matrix * A_0
n_0 = A_0.T * P_matrix * dl_0
Qxx_0 = N.inv()
dx = Qxx_0 * n_0
x = x0 + dx
@@ -165,7 +211,6 @@ class Transformationen:
anzahl_iterationen += 1
print(f"Iteration Nr.{anzahl_iterationen} abgeschlossen")
print(dx.evalf(n=3))
else:
zahlen_i = {
@@ -182,16 +227,17 @@ class Transformationen:
dl_i = l_vektor - l_berechnet_i
A_i = A_ohne_zahlen.subs(zahlen_i).evalf(n=30)
N_i = A_i.T * P_mat * A_i
N_i = A_i.T * P_matrix * A_i
Qxx_i = N_i.inv()
n_i = A_i.T * P_mat * dl_i
n_i = A_i.T * P_matrix * dl_i
dx = Qxx_i * n_i
x = sp.Matrix(x + dx)
anzahl_iterationen += 1
print(f"Iteration Nr.{anzahl_iterationen} abgeschlossen")
print(dx.evalf(n=3))
alle_kleiner = True
for i in range(dx.rows):
@@ -204,11 +250,7 @@ class Transformationen:
alle_kleiner_vorherige_iteration = alle_kleiner
print(l.evalf(n=3))
print(l_berechnet_0.evalf(n=3))
print(f"x = {x.evalf(n=3)}")
# --- Neuberechnung Zielsystem ---
# Neuberechnung Zielsystem
zahlen_final = {
dX: float(x[0]),
dY: float(x[1]),
@@ -228,12 +270,11 @@ class Transformationen:
Zi = l_berechnet_final[3 * i + 2]
liste_l_berechnet_final.append(sp.Matrix([Xi, Yi, Zi]))
print("")
print("l_berechnet_final:")
for punktnummer, l_fin in zip(gemeinsame_punktnummern, liste_l_berechnet_final):
print(f"{punktnummer}: {float(l_fin[0]):.3f}, {float(l_fin[1]):.3f}, {float(l_fin[2]):.3f}")
print("Streckendifferenzen:")
print("Streckendifferenzen zwischen Näherungskoordinate aus statischer GNSS-Messung und ergebnis der Helmerttransformation:")
streckendifferenzen = [
(punkt_zielsys - l_final).norm()
for punkt_zielsys, l_final in zip(liste_punkte_zielsystem, liste_l_berechnet_final)
@@ -245,30 +286,51 @@ class Transformationen:
Schwerpunktsdifferenz = Schwerpunkt_Zielsystem - Schwerpunkt_berechnet
print("\nDifferenz Schwerpunkt (Vektor):")
print("\nDifferenz Schwerpunkt zwischen Näherungskoordinate aus statischer GNSS-Messung und ergebnis der Helmerttransformation::")
print(Schwerpunktsdifferenz.evalf(3))
print("Betrag der Schwerpunkt-Differenz:")
print("Betrag der Schwerpunkt-Differenz zwischen Näherungskoordinate aus statischer GNSS-Messung und ergebnis der Helmerttransformation::")
print(f"{float(Schwerpunktsdifferenz.norm()):.3f}m")
return zahlen_final
def Helmerttransformation(self, transformationsparameter: dict) -> dict[Any, Any]:
db = Datenbank.Datenbankzugriff(self.pfad_datenbank)
dict_ausgangssystem = db.get_koordinaten("naeherung_lh", "Dict")
dict_zielsystem = db.get_koordinaten("naeherung_us", "Dict")
"""Wendet eine Helmerttransformation auf Punkte des Ausgangssystems an.
Aus der Datenbank werden Koordinaten des Ausgangssystems (naeherung_lh) und Zielsystems (naeherung_us) geladen.
Transformiert werden genau die Punkte, die:
- im Ausgangssystem vorhanden sind,
- im Zielsystem fehlen (symmetrische Differenz der Punktmengen, anschließend Filter auf Ausgangssystem).
Die Transformation erfolgt gemäß:
P = [dX, dY, dZ]^T + m * R(e1,e2,e3) * p
:param transformationsparameter: Transformationsparameter als Dictionary mit SymPy-Symbolen als Keys
(dX, dY, dZ, m, e1, e2, e3) und numerischen Werten als Values.
:type transformationsparameter: dict
:return: Dictionary {punktnummer: sp.Matrix([X, Y, Z])} der transformierten geozentrisch-kartesischen Koordinaten.
:rtype: dict[Any, Any]
"""
# Koordinaten des lokalen Horizontsystems des Tachymeters und der geozentrisch Kartesischen Näherungskoordinaten aus den statischen GNSS-Messungen aus der Tabelle Netzpunkte abfragen
dict_ausgangssystem = self.db_zugriff.get_koordinaten("naeherung_lh", "Dict")
dict_zielsystem = self.db_zugriff.get_koordinaten("naeherung_us", "Dict")
# Symbole definieren
dX, dY, dZ, m, e1, e2, e3 = sp.symbols('dX dY dZ m e1 e2 e3')
# Unterschiedliche Punkte zwischen Ausgangs- und Zielsystem ermitteln
unterschiedliche_punktnummern = sorted(set(dict_ausgangssystem.keys()) ^ set(dict_zielsystem.keys()))
punktnummern_transformieren = [
punktnummer for punktnummer in unterschiedliche_punktnummern if punktnummer in dict_ausgangssystem
]
liste_punkte_ausgangssystem = [dict_ausgangssystem[punktnummer] for punktnummer in punktnummern_transformieren]
R = self.R_matrix_aus_euler(transformationsparameter[e1], transformationsparameter[e2], transformationsparameter[e3])
# Rotationsmatrix aufstellen
R = self.R_matrix_aus_eulerwinkeln(transformationsparameter[e1], transformationsparameter[e2], transformationsparameter[e3])
f_zeilen = []
# Helmertransformation durchführen und Koordinaten speichern
for punkt in liste_punkte_ausgangssystem:
punkt_vektor = sp.Matrix([punkt[0], punkt[1], punkt[2]])
f_zeile_i = sp.Matrix([transformationsparameter[dX], transformationsparameter[dY], transformationsparameter[dZ]]) + transformationsparameter[m] * R * punkt_vektor
@@ -289,9 +351,28 @@ class Transformationen:
return dict_transformiert
def utm_to_XYZ(self, pfad_tif_quasigeoidundolation: str, liste_utm: list) -> dict[Any, Any]:
"""Rechnet UTM-Koordinaten (ETRS89 / UTM + DHHN2016) in ECEF-Koordinaten (ETRS89 geozentrisch-kartesisch) um.
Es wird ein PROJ-Transformer von:
- Quelle: EPSG:25832 + EPSG:7837 (ETRS89 / UTM Zone 32N + DHHN2016 Normalhöhe),
- Ziel: EPSG:4936 (ETRS89 geozentrisch-kartesisch)
initialisiert. Zusätzlich wird ein BKG-GeoTIFF (Quasigeoidunndulation) in den PROJ-Datenpfad eingebunden,
indem eine Kopie mit dem erwarteten Dateinamen "de_bkg_gcg2016.tif" im selben Ordner erzeugt wird.
:param pfad_tif_quasigeoidundolation: Pfad zur BKG-GeoTIFF-Datei (Quasigeoidundulation).
:type pfad_tif_quasigeoidundolation: str
:param liste_utm: Liste von UTM-Koordinaten in der Form [(punktnummer, E, N, Normalhoehe), ...].
:type liste_utm: list
:return: Dictionary {punktnummer: sp.Matrix([X, Y, Z])} mit ECEF-Koordinaten (Meter).
:rtype: dict[Any, Any]
"""
# tif vom BKG zur Quasigeoidundolation übergeben
pfad_gcg_tif = Path(pfad_tif_quasigeoidundolation)
pfad_gcg_tif_proj = pfad_gcg_tif.with_name("de_bkg_gcg2016.tif")
# Kopie des TIF anlegen (Dies ist voraussetzung für die Transformer-Bibliothek
if (not pfad_gcg_tif_proj.exists()) or (pfad_gcg_tif_proj.stat().st_size != pfad_gcg_tif.stat().st_size):
shutil.copy2(pfad_gcg_tif, pfad_gcg_tif_proj)
@@ -299,8 +380,9 @@ class Transformationen:
utm_epsg = 25832
crs_src = CRS.from_user_input(f"EPSG:{utm_epsg}+EPSG:7837") # ETRS89/DREF91 + DHHN2016
crs_dst = CRS.from_epsg(4936) # ETRS89 geozentrisch (ECEF)
crs_dst = CRS.from_epsg(4936) # ETRS89 geozentrisch kartesisch
# Umrechnungsvorgaben übergeben
tr_best = Transformer.from_crs(
crs_src,
crs_dst,
@@ -308,52 +390,54 @@ class Transformationen:
allow_ballpark=False,
)
# Koordinaten rechnen und in Dictionary speichern
dict_geozentrisch_kartesisch = {}
for Punktnummer, E, N, Normalhoehe in liste_utm:
X, Y, Z = tr_best.transform(E, N, Normalhoehe)
dict_geozentrisch_kartesisch[Punktnummer] = sp.Matrix([X, Y, Z])
# geographisch 3D + zeta
#crs_geog3d = CRS.from_epsg(4937) # ETRS89 (lon, lat, h)
#tr_h = Transformer.from_crs(
# crs_src,
# crs_geog3d,
# always_xy=True,
# allow_ballpark=False,
#)
#lon, lat, h = tr_h.transform(E, N, H)
#print("lon/lat/h:", lon, lat, h)
#print("zeta (h-H):", h - H)
return dict_geozentrisch_kartesisch
def ecef_to_utm(
self,
dict_koordinaten: dict,
pfad_gcg_tif: str | Path | None = None,
zone: int = 32,
):
pfad_gcg_tif: str | Path | None = None):
"""Rechnet ECEF-Koordinaten (ETRS89 geozentrisch-kartesisch) nach nach UTM (+ DHHN2016 Normalhöhe).
if pfad_gcg_tif is not None:
pfad_gcg_tif = Path(pfad_gcg_tif).resolve()
if not pfad_gcg_tif.exists():
raise FileNotFoundError(f"Quasigeoid-Datei nicht gefunden: {pfad_gcg_tif}")
Es wird ein PROJ-Transformer von:
pfad_proj_grid = pfad_gcg_tif.with_name("de_bkg_gcg2016.tif")
if (
not pfad_proj_grid.exists()
or pfad_proj_grid.stat().st_size != pfad_gcg_tif.stat().st_size
):
shutil.copy2(pfad_gcg_tif, pfad_proj_grid)
- Quelle: EPSG:4936 (ETRS89 geozentrisch-kartesisch),
- Ziel: EPSG:25832 + EPSG:7837 (ETRS89 / UTM Zone 32N + DHHN2016 Normalhöhe)
datadir.append_data_dir(str(pfad_proj_grid.parent))
initialisiert. Zusätzlich wird die BKG-GeoTIFF-Datei (Quasigeoidundulation) als PROJ-Grid eingebunden,
indem eine Kopie mit dem erwarteten Namen "de_bkg_gcg2016.tif" im selben Ordner erzeugt wird.
crs_src = CRS.from_epsg(4936) # ETRS89 geocentric (ECEF)
Die Methode akzeptiert Koordinatenwerte in verschiedenen Formen (SymPy-Matrix, numpy.ndarray,
Liste/Tuple, Skalar) und extrahiert daraus drei Werte (X, Y, Z). Die Ergebnisse (E, N, H) werden auf 8 Nachkommastellen gerundet.
:param dict_koordinaten: Dictionary {punktnummer: koordinate}, wobei koordinate X/Y/Z enthält.
:type dict_koordinaten: dict
:param pfad_gcg_tif: Pfad zur BKG-GeoTIFF-Datei (Quasigeoidundulation) als str.
:type pfad_gcg_tif: str | Path | None
:return: Dictionary {punktnummer: (E, N, H)} mit UTM-Koordinaten (Meter) und Normalhöhe.
:rtype: dict
"""
# Kopie des TIF vom BKG mit der Quasigeoidundolation erstellen
pfad_gcg_tif = Path(pfad_gcg_tif).resolve()
pfad_proj_grid = pfad_gcg_tif.with_name("de_bkg_gcg2016.tif")
if (
not pfad_proj_grid.exists()
or pfad_proj_grid.stat().st_size != pfad_gcg_tif.stat().st_size
):
shutil.copy2(pfad_gcg_tif, pfad_proj_grid)
datadir.append_data_dir(str(pfad_proj_grid.parent))
# EPSG-Codes feslegen
crs_src = CRS.from_epsg(4936) # ETRS89 geozentrisch-kartesisch
# Ziel-CRS: ETRS89 / UTM Zone 32/33 + DHHN2016 Normalhöhe
# EPSG:25832/25833 = ETRS89 / UTM; EPSG:7837 = DHHN2016 height
utm_epsg = 25800 + zone # 25832 oder 25833
utm_epsg = 25832
crs_dst = CRS.from_user_input(f"EPSG:{utm_epsg}+EPSG:7837")
tr = Transformer.from_crs(
@@ -365,6 +449,7 @@ class Transformationen:
tr_geo = Transformer.from_crs(CRS.from_epsg(4936), CRS.from_epsg(4979), always_xy=True)
# Koordinaten an Dictionary übergeben
dict_koordinaten_utm = {}
for punktnummer, koordinate in dict_koordinaten.items():
werte = []
@@ -399,20 +484,8 @@ class Transformationen:
# Skalar
werte.append(float(v))
if len(werte) < 3:
raise ValueError(f"Zu wenig skalare Werte gefunden: {werte}")
X, Y, Z = werte[0], werte[1], werte[2]
try:
E, N, H = tr.transform(X, Y, Z, errcheck=True)
except ProjError as e:
lon, lat, h_ell = tr_geo.transform(X, Y, Z, errcheck=True)
raise ProjError(
f"transform error (outside grid) | pn={punktnummer} | "
f"X,Y,Z={X},{Y},{Z} | lon/lat={lon},{lat} | h_ell={h_ell} | {e}"
)
E, N, H = tr.transform(X, Y, Z, errcheck=True)
# Runden, weil ansonsten aufgrund begrenzter Rechenkapazität falsche Werte Resultieren
dict_koordinaten_utm[punktnummer] = (round(E, 8), round(N, 8), round(H, 8))