Push 05.02.2026
This commit is contained in:
157
Berechnungen.py
157
Berechnungen.py
@@ -1,7 +1,4 @@
|
||||
from decimal import Decimal
|
||||
import sympy as sp
|
||||
from typing import Any
|
||||
import math
|
||||
import numpy as np
|
||||
|
||||
import Datenbank
|
||||
@@ -383,92 +380,30 @@ class Berechnungen:
|
||||
return schraegdistanz_boden, zw_boden
|
||||
|
||||
|
||||
class Einheitenumrechnung:
|
||||
"""Einheitenumrechnungen für Winkel- und Längeneinheiten.
|
||||
class ENU:
|
||||
"""Transformationen zwischen XYZ und lokalem ENU-System.
|
||||
|
||||
Die Klasse stellt Methoden zur Verfügung für:
|
||||
|
||||
- Umrechnung von Millibogensekunden (mas) in Radiant,
|
||||
- Umrechnung von Millimetern (mm) in Meter,
|
||||
- Umrechnung von Gon und Milligon (mgon) in Radiant (Decimal-basiert).
|
||||
- Bestimmung eines ENU-Referenzpunkts über den Schwerpunkt gegebener XYZ-Koordinaten,
|
||||
- Aufbau der lokalen Rotationsmatrix R0 (XYZ nach ENU) aus geodätischer Breite B und Länge L,
|
||||
- Aufbau einer Rotationsmatrix R_ENU für einen gesamten Unbekanntenvektor,
|
||||
- Transformation der Kovarianz-Matrix Qxx in das ENU-System,
|
||||
- Transformation von Punktkoordinaten (XYZ) in lokale ENU-Koordinaten relativ zum Schwerpunkt.
|
||||
"""
|
||||
def mas_to_rad(mas: float) -> float:
|
||||
"""Rechnet Millibogensekunden (mas) in Radiant um.
|
||||
|
||||
Es gilt: rad = mas * (pi / (180 * 3600 * 1000)).
|
||||
|
||||
:param mas: Winkel in Millibogensekunden (mas).
|
||||
:type mas: float
|
||||
:return: Winkel in Radiant.
|
||||
:rtype: float
|
||||
"""
|
||||
umrechnungsfaktor = 1 / 1000 * 1 / 3600 * sp.pi / 180
|
||||
grad = mas * umrechnungsfaktor
|
||||
return grad
|
||||
|
||||
def mm_to_m(mm: float) -> float:
|
||||
"""Rechnet Millimeter in Meter um.
|
||||
|
||||
Es gilt: m = mm / 1000.
|
||||
|
||||
:param mm: Länge in Millimeter.
|
||||
:type mm: float
|
||||
:return: Länge in Meter.
|
||||
:rtype: float
|
||||
"""
|
||||
m = mm / 1000
|
||||
return m
|
||||
|
||||
def gon_to_rad_Decimal(gon: float) -> Decimal:
|
||||
"""Rechnet Gon in Radiant um (Decimal-basiert).
|
||||
|
||||
Es gilt: 400 gon = 2*pi und damit rad = (gon / 200) * pi.
|
||||
|
||||
:param gon: Winkel in Gon.
|
||||
:type gon: float
|
||||
:return: Winkel in Radiant als Decimal.
|
||||
:rtype: Decimal
|
||||
"""
|
||||
gon = Decimal(gon)
|
||||
pi = Decimal(str(math.pi))
|
||||
rad = (gon / Decimal(200)) * pi
|
||||
return rad
|
||||
|
||||
def mgon_to_rad_Decimal(gon: float) -> Decimal:
|
||||
"""Rechnet Milligon (mgon) in Radiant um (Decimal-basiert).
|
||||
|
||||
Es gilt: 1 mgon = 0.001 gon und damit rad = (mgon / 200000) * pi.
|
||||
|
||||
:param gon: Winkel in Milligon (mgon).
|
||||
:type gon: float
|
||||
:return: Winkel in Radiant als Decimal.
|
||||
:rtype: Decimal
|
||||
"""
|
||||
gon = Decimal(gon)
|
||||
pi = Decimal(str(math.pi))
|
||||
rad = (gon / Decimal(200000)) * pi
|
||||
return rad
|
||||
|
||||
def rad_to_gon_Decimal(rad: float) -> Decimal:
|
||||
"""Rechnet Radiant in Gon um (Decimal-basiert).
|
||||
|
||||
Es gilt: 400 gon = 2*pi und damit rad = (gon / 200) * pi.
|
||||
|
||||
:param rad: Winkel in Rad.
|
||||
:type rad: float
|
||||
:return: Winkel in Gon als Decimal.
|
||||
:rtype: Decimal
|
||||
"""
|
||||
rad = Decimal(rad)
|
||||
pi = Decimal(str(math.pi))
|
||||
gon = (rad / pi) * Decimal(200)
|
||||
return gon
|
||||
|
||||
|
||||
class ENU:
|
||||
|
||||
@staticmethod
|
||||
def berechne_schwerpunkt_fuer_enu(berechnungen, dict_xyz):
|
||||
"""
|
||||
Berechnet die ENU-Referenz (B0, L0) aus dem Schwerpunkt gegebener XYZ-Koordinaten.
|
||||
|
||||
:param berechnungen: Hilfsobjekt mit Methoden zur Berechnung von Breite B und Länge L.
|
||||
:type berechnungen: Any
|
||||
:param dict_xyz: Dictionary der Koordinaten.
|
||||
:type dict_xyz: dict
|
||||
:return: Tuple aus geodätischer Breite B0 und Länge L0 des Schwerpunktes.
|
||||
:rtype: tuple[float, float]
|
||||
"""
|
||||
XYZ = np.array(list(dict_xyz.values()), dtype=float)
|
||||
X0, Y0, Z0 = XYZ.mean(axis=0)
|
||||
B0 = float(berechnungen.B(X0, Y0, Z0))
|
||||
@@ -478,6 +413,18 @@ class ENU:
|
||||
|
||||
@staticmethod
|
||||
def berechne_R0_ENU(berechnungen, B, L):
|
||||
"""
|
||||
Erzeugt die 3×3-Rotationsmatrix R0 für die Transformation von XYZ nach ENU.
|
||||
|
||||
:param berechnungen: Hilfsobjekt mit Methoden zur Berechnung der ENU-Komponenten.
|
||||
:type berechnungen: Any
|
||||
:param B: Geodätische Breite des ENU-Referenzpunkts.
|
||||
:type B: float
|
||||
:param L: Geodätische Länge des ENU-Referenzpunkts.
|
||||
:type L: float
|
||||
:return: Rotationsmatrix R0 (3×3) für XYZ nach ENU.
|
||||
:rtype: numpy.ndarray
|
||||
"""
|
||||
# East
|
||||
r11 = berechnungen.E(L, 1, 0)
|
||||
r12 = berechnungen.E(L, 0, 1)
|
||||
@@ -504,6 +451,20 @@ class ENU:
|
||||
|
||||
@staticmethod
|
||||
def berechne_R_ENU(unbekannten_liste, R0):
|
||||
"""
|
||||
Erstellt eine Rotationsmatrix R für die Umrechnung eines Unbekanntenvektors ins ENU-System.
|
||||
|
||||
Für jede Punkt-ID, die über Symbolnamen (z.B. X1, Y1, Z1) in unbekannten_liste erkannt wird,
|
||||
wird die 3×3-Rotation R0 in die passende Position der Gesamtmatrix eingesetzt.
|
||||
|
||||
:param unbekannten_liste: Liste der Unbekannten
|
||||
:type unbekannten_liste: list
|
||||
:param R0: Rotationsmatrix (3×3) für XYZ nach ENU.
|
||||
:type R0: numpy.ndarray
|
||||
:return: Rotationsmarix R (n×n) für den gesamten Unbekanntenvektor.
|
||||
:rtype: numpy.ndarray
|
||||
"""
|
||||
|
||||
names = [str(s) for s in unbekannten_liste]
|
||||
n = len(names)
|
||||
R = np.eye(n, dtype=float)
|
||||
@@ -526,6 +487,27 @@ class ENU:
|
||||
|
||||
@staticmethod
|
||||
def transform_Qxx_zu_QxxENU(Qxx, unbekannten_liste, berechnungen, dict_xyz):
|
||||
"""
|
||||
Transformiert die Kofaktor-Matrix Qxx in das ENU-System.
|
||||
|
||||
Es wird zunächst eine ENU-Referenz über den Schwerpunkt der Koordinaten bestimmt.
|
||||
Anschließend wird R0 (XYZ nach ENU) und die Matrix R_ENU aufgebaut. Die Transformation
|
||||
erfolgt über:
|
||||
|
||||
Qxx_ENU = R_ENU · Qxx · R_ENUᵀ
|
||||
|
||||
:param Qxx: Kofaktor-Matrix der Unbekannten im XYZ-System.
|
||||
:type Qxx: numpy.ndarray
|
||||
:param unbekannten_liste: Liste der Unbekannten.
|
||||
:type unbekannten_liste: list
|
||||
:param berechnungen: Hilfsobjekt für Breite/Länge und ENU-Komponenten.
|
||||
:type berechnungen: Any
|
||||
:param dict_xyz: Dictionary der Koordinaten.
|
||||
:type dict_xyz: dict
|
||||
:return: Qxx_ENU, (B0, L0), R0 mit Schwerpunkt-Referenz und lokaler Rotationsmatrix.
|
||||
:rtype: tuple[numpy.ndarray, tuple[float, float], numpy.ndarray]
|
||||
"""
|
||||
|
||||
B0, L0 = ENU.berechne_schwerpunkt_fuer_enu(berechnungen, dict_xyz)
|
||||
R0 = ENU.berechne_R0_ENU(berechnungen, B0, L0)
|
||||
R_ENU = ENU.berechne_R_ENU(unbekannten_liste, R0)
|
||||
@@ -538,6 +520,17 @@ class ENU:
|
||||
|
||||
@staticmethod
|
||||
def transform_Koord_zu_KoordENU(dict_xyz, R0):
|
||||
"""
|
||||
Transformiert Punktkoordinaten (XYZ) in ENU-Koordinaten relativ zum Schwerpunkt.
|
||||
|
||||
:param dict_xyz: Dictionary der Koordinaten.
|
||||
:type dict_xyz: dict
|
||||
:param R0: Rotationsmatrix (3×3) für XYZ nach ENU.
|
||||
:type R0: numpy.ndarray
|
||||
:return: Dictionary der ENU-Koordinaten.
|
||||
:rtype: dict
|
||||
"""
|
||||
|
||||
XYZ = np.asarray(list(dict_xyz.values()), dtype=float).reshape(-1, 3)
|
||||
XYZ0 = XYZ.mean(axis=0).reshape(3, )
|
||||
|
||||
|
||||
37410
Campusnetz.ipynb
37410
Campusnetz.ipynb
File diff suppressed because it is too large
Load Diff
104
Datenbank.py
104
Datenbank.py
@@ -1,11 +1,12 @@
|
||||
from decimal import Decimal
|
||||
import os
|
||||
import pandas as pd
|
||||
import sqlite3
|
||||
import sympy as sp
|
||||
from sympy import MutableDenseMatrix
|
||||
from typing import Any
|
||||
|
||||
from Berechnungen import Einheitenumrechnung
|
||||
from Einheitenumrechnung import Einheitenumrechnung
|
||||
|
||||
class Datenbank_anlegen:
|
||||
"""Legt die SQLite-Datenbank für die Ausgleichungsrechnung an.
|
||||
@@ -147,6 +148,7 @@ class Datenbankzugriff:
|
||||
|
||||
- Schreiben in die Datenbank durch alle set_* Methoden
|
||||
- Lesen aus der Datenbank durch alle get_* Methoden
|
||||
- Darstellen der aus der Datenbenak gelesenen Daten in tabellarischer Form durch alle tabelle_* Methoden
|
||||
|
||||
:ivar pfad_datenbank: Pfad zur SQLite-Datenbankdatei.
|
||||
:vartype pfad_datenbank: str
|
||||
@@ -241,7 +243,7 @@ class Datenbankzugriff:
|
||||
else:
|
||||
id_instrument = cursor.execute(
|
||||
"SELECT instrumenteID FROM Instrumente WHERE typ = ? AND name =?", (typ, name))
|
||||
print(f"Das Instrument {name} ist bereits in der Datenbank vorhanden.\nEs hat die ID {id_instrument.fetchone()[0]}")
|
||||
print(f"Das Instrument {name} ist bereits in der Datenbank vorhanden.")
|
||||
con.commit()
|
||||
cursor.close()
|
||||
con.close()
|
||||
@@ -252,7 +254,8 @@ class Datenbankzugriff:
|
||||
|
||||
Prüft, ob instrumenteID existiert und ob mindestens eine Genauigkeitsangabe übergeben wurde.
|
||||
Je nach Beobachtungsart werden Einheitenumrechnungen durchgeführt (z. B. mgon → rad bzw. mm → m).
|
||||
Der Eintrag wird nur ergänzt, wenn in Genauigkeiten kein identischer Datensatz vorhanden ist.
|
||||
Der Eintrag wird neu hinzugefügt, wenn in Genauigkeiten kein identischer Datensatz vorhanden ist.
|
||||
Ansonsten wird der bestehende Datensatz in der Datenbank aktualisiert.
|
||||
|
||||
:param instrumenteID: ID des Instruments in der Tabelle Instrumente.
|
||||
:type instrumenteID: int
|
||||
@@ -300,51 +303,36 @@ class Datenbankzugriff:
|
||||
if isinstance(stabw_apriori_streckenprop, Decimal):
|
||||
stabw_apriori_streckenprop = float(stabw_apriori_streckenprop)
|
||||
|
||||
# Überprüfen, ob die Genauigkeitsinformation für die jeweilige Beobachtungsart des Instruments bereits vorhanden ist. Wenn nein, wird die Benutzereingabe in die Datenbank gespeichert.
|
||||
sql = "SELECT 1 FROM Genauigkeiten WHERE instrumenteID = ? AND beobachtungsart = ?"
|
||||
params = [instrumenteID, beobachtungsart]
|
||||
# Prüfen, ob für dieses Instrument und diese Art bereits ein Eintrag existiert
|
||||
vorhanden = cursor.execute("SELECT 1 FROM Genauigkeiten WHERE instrumenteID = ? AND beobachtungsart = ?",
|
||||
(instrumenteID, beobachtungsart)).fetchone()
|
||||
|
||||
if stabw_apriori_konstant is None:
|
||||
sql += " AND stabw_apriori_konstant IS NULL"
|
||||
if vorhanden:
|
||||
# Bestehenden Datensatz aktualisieren
|
||||
cursor.execute(
|
||||
"UPDATE Genauigkeiten SET stabw_apriori_konstant = ?, stabw_apriori_streckenprop = ? WHERE instrumenteID = ? AND beobachtungsart = ?",
|
||||
(stabw_apriori_konstant, stabw_apriori_streckenprop, instrumenteID, beobachtungsart)
|
||||
)
|
||||
print(f"Die Genauigkeitsangabe für {beobachtungsart} (Instrument: {instrumentenname}) wurde aktualisiert.")
|
||||
else:
|
||||
sql += " AND stabw_apriori_konstant = ?"
|
||||
params.append(stabw_apriori_konstant)
|
||||
|
||||
if stabw_apriori_streckenprop is None:
|
||||
sql += " AND stabw_apriori_streckenprop IS NULL"
|
||||
else:
|
||||
sql += " AND stabw_apriori_streckenprop = ?"
|
||||
params.append(stabw_apriori_streckenprop)
|
||||
|
||||
liste_genauigkeiten = cursor.execute(sql, tuple(params)).fetchall()
|
||||
|
||||
if liste_genauigkeiten == []:
|
||||
# Neuen Datensatz anlegen
|
||||
if stabw_apriori_konstant is not None and stabw_apriori_streckenprop is not None:
|
||||
cursor.execute(
|
||||
"INSERT INTO Genauigkeiten (instrumenteID, beobachtungsart, stabw_apriori_konstant, stabw_apriori_streckenprop) VALUES (?, ?, ?, ?)",
|
||||
(instrumenteID, beobachtungsart, stabw_apriori_konstant, stabw_apriori_streckenprop)
|
||||
)
|
||||
print(
|
||||
f"Die Genauigkeitsangabe für die Beobachtungsart {beobachtungsart} des Instrumentes {instrumentenname} wurde erfolgreich hinzugefügt.")
|
||||
|
||||
elif stabw_apriori_konstant is None and stabw_apriori_streckenprop is not None:
|
||||
elif stabw_apriori_konstant is None:
|
||||
cursor.execute(
|
||||
"INSERT INTO Genauigkeiten (instrumenteID, beobachtungsart, stabw_apriori_streckenprop) VALUES (?, ?, ?)",
|
||||
(instrumenteID, beobachtungsart, stabw_apriori_streckenprop)
|
||||
)
|
||||
print(
|
||||
f"Die Genauigkeitsangabe für die Beobachtungsart {beobachtungsart} des Instrumentes {instrumentenname} wurde erfolgreich hinzugefügt.")
|
||||
|
||||
elif stabw_apriori_streckenprop is None and stabw_apriori_konstant is not None:
|
||||
else:
|
||||
cursor.execute(
|
||||
"INSERT INTO Genauigkeiten (instrumenteID, beobachtungsart, stabw_apriori_konstant) VALUES (?, ?, ?)",
|
||||
(instrumenteID, beobachtungsart, stabw_apriori_konstant)
|
||||
)
|
||||
print(
|
||||
f"Die Genauigkeitsangabe für die Beobachtungsart {beobachtungsart} des Instrumentes {instrumentenname} wurde erfolgreich hinzugefügt.")
|
||||
|
||||
else:
|
||||
print("Die Genauigkeitsangabe ist bereits in der Datenbank vorhanden.")
|
||||
print(
|
||||
f"Die Genauigkeitsangabe für {beobachtungsart} (Instrument: {instrumentenname}) wurde erfolgreich hinzugefügt.")
|
||||
|
||||
con.commit()
|
||||
cursor.close()
|
||||
@@ -711,7 +699,7 @@ class Datenbankzugriff:
|
||||
return liste_instrumente
|
||||
|
||||
def get_alle_instrumente_liste(self: str) -> list:
|
||||
"""Liest alle Instrumente aus der Tabelle Instrumente.
|
||||
"""Liest alle Instrumente aus der Tabelle Instrumente auf.
|
||||
|
||||
Gibt eine Liste der gefundenen Instrumente zurück. Falls keine Instrumente vorhanden sind,
|
||||
wird eine Textausgabe mit verfügbaren Typen zurückgegeben.
|
||||
@@ -977,4 +965,52 @@ class Datenbankzugriff:
|
||||
con.close()
|
||||
return liste_varianzkomponenten
|
||||
|
||||
def tabelle_instrumente_aus_db(self) -> None:
|
||||
"""Stellt die in der Datenbank gespeicherten Instrumente tabellarisch dar.
|
||||
|
||||
Die Methode liest alle Datensätze aus der Tabelle Instrumente über die
|
||||
zugehörige Datenbankzugriffsmethode und erzeugt daraus einen pandas-DataFrame
|
||||
zur Anzeige im Notebook. Sind keine Instrumente vorhanden, wird ein Hinweistext
|
||||
ausgegeben, dass die Instrumente zuerst anzulegen sind.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
liste_instrumente_in_db = self.get_alle_instrumente_liste()
|
||||
|
||||
# Prüfen, ob Datensätze in der Tabelle Instrumente enthalten sind
|
||||
if isinstance(liste_instrumente_in_db, list) and len(liste_instrumente_in_db) > 0:
|
||||
df_instrumente = pd.DataFrame(liste_instrumente_in_db, columns=['InstrumenteID', 'Typ', 'Bezeichnung'])
|
||||
display(df_instrumente.style.hide(axis='index'))
|
||||
|
||||
else:
|
||||
print(
|
||||
"Es wurden noch keine Instrumente angelegt. Bitte in der folgenden Zelle nachholen und diese Zelle erneut ausführen!")
|
||||
|
||||
def tabelle_genauigkeiten_aus_db(self) -> None:
|
||||
"""Stellt die a-priori Genauigkeiten der Beobachtungsgruppen tabellarisch dar.
|
||||
|
||||
Die Methode liest alle Einträge aus der Tabelle Genauigkeiten über die
|
||||
zugehörige Datenbankzugriffsmethode und erstellt eine
|
||||
übersichtliche Anzeige im Notebook. Liegen keine Genauigkeitsangaben vor,
|
||||
wird ein Hinweistext ausgegeben, dass die a-priori Genauigkeiten zunächst zu
|
||||
erfassen sind.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
genauigkeiten_dict = self.get_genauigkeiten_dict()
|
||||
# Prüfen, ob Datensätze in der Tabelle Instrumente enthalten sind
|
||||
if genauigkeiten_dict == {}:
|
||||
print(
|
||||
"Es wurden noch keine apriori Genauigkeiten zu den Beobachtungsgruppen erfasst. Bitte in der folgenden Zelle nachholen und diese Zelle erneut ausführen.")
|
||||
else:
|
||||
formatierte_daten = list(genauigkeiten_dict.values())
|
||||
spalten = [
|
||||
'instrumenteID',
|
||||
'beobachtungsart',
|
||||
'stabw_apriori_konstant',
|
||||
'stabw_apriori_streckenprop'
|
||||
]
|
||||
df = pd.DataFrame(formatierte_daten, columns=spalten)
|
||||
display(df.style.hide(axis='index'))
|
||||
|
||||
@@ -9,6 +9,32 @@ class Datumsfestlegung:
|
||||
@staticmethod
|
||||
def weiches_datum(Q_ll: np.ndarray, Q_AA: np.ndarray, varianzkompontenschaetzung_erfolgt: bool,
|
||||
dict_indizes_beobachtungsgruppen: dict) -> np.ndarray:
|
||||
"""
|
||||
Erstellt die erweiterte Kofaktor- und Gewichtsmatrix für eine weiche Datumsfestlegung.
|
||||
|
||||
Aus den Kofaktormatrizen der Beobachtungen Q_ll und der Kofaktormatrix der Anschlusspunkte Q_AA
|
||||
wird eine erweiterte Kofaktormatrix Q_ext aufgebaut. Zusätzlich wird die zugehörige Gewichtsmatrix P erzeugt.
|
||||
|
||||
Falls keine Varianzkomponentenschätzung durchgeführt wurde, wird P als Blockmatrix aus den
|
||||
Inversen (Gewichten) von Q_ll und Q_AA aufgebaut.
|
||||
|
||||
Falls eine Varianzkomponentenschätzung durchgeführt wurde, wird Q_ext entsprechend den in definierten
|
||||
Beobachtungsgruppen in Teilblöcke zerlegt (z.B. SD, R, ZW, gnss, niv, lA). Für jeden Block wird die Gewichtsmatrix
|
||||
berechnet und anschließend zu einer Gesamtgewichtsmatrix zusammengesetzt.
|
||||
|
||||
:param Q_ll: a-priori Kofaktormatrix der Beobachtungen.
|
||||
:type Q_ll: numpy.ndarray
|
||||
:param Q_AA: a-priori Kofaktormatrix der Anschlusspunkte.
|
||||
:type Q_AA: numpy.ndarray
|
||||
:param varianzkompontenschaetzung_erfolgt: Kennzeichen, ob eine Varianzkomponentenschätzung berücksichtigt werden soll.
|
||||
:type varianzkompontenschaetzung_erfolgt: bool
|
||||
:param dict_indizes_beobachtungsgruppen: Dictionary mit Indexbereichen je Beobachtungsgruppe zur Blockzerlegung.
|
||||
:type dict_indizes_beobachtungsgruppen: dict
|
||||
:return: Tuple aus erweiterter Kofaktormatrix Q_ext und zugehöriger Gewichtsmatrix P.
|
||||
:rtype: tuple[numpy.ndarray, numpy.ndarray]
|
||||
:raises ValueError: Wenn Q_ll oder Q_AA keine quadratische Matrix ist.
|
||||
"""
|
||||
|
||||
if Stochastisches_Modell.StochastischesModell.berechne_P(Q_ll).ndim != 2 or \
|
||||
Stochastisches_Modell.StochastischesModell.berechne_P(Q_ll).shape[0] != \
|
||||
Stochastisches_Modell.StochastischesModell.berechne_P(Q_ll).shape[1]:
|
||||
@@ -95,21 +121,37 @@ class Datumsfestlegung:
|
||||
Z(aufgeteilt_lA_invertiert, aufgeteilt_gnss_invertiert),
|
||||
Z(aufgeteilt_lA_invertiert, aufgeteilt_niv_invertiert), aufgeteilt_lA_invertiert]
|
||||
])
|
||||
|
||||
# print(aufgeteilt)
|
||||
# print(beobachtungsgruppe, indizes)
|
||||
# Export.matrix_to_csv(
|
||||
# fr"Zwischenergebnisse\_{beobachtungsgruppe}.csv",
|
||||
# [""],
|
||||
# labels,
|
||||
# aufgeteilt,
|
||||
# f"{beobachtungsgruppe}"
|
||||
# )
|
||||
|
||||
return Q_ext, P
|
||||
|
||||
@staticmethod
|
||||
def indizes_beobachtungsvektor_nach_beobachtungsgruppe(Jacobimatrix_symbolisch_liste_beobachtungsvektor):
|
||||
"""
|
||||
Ermittelt Indexbereiche des Beobachtungsvektors getrennt nach Beobachtungsgruppen.
|
||||
|
||||
Die Funktion analysiert die Bezeichner des symbolischen Beobachtungsvektors (z.B.
|
||||
aus der Jacobi-Matrix) und ordnet jede Beobachtung anhand ihres Kennzeichens
|
||||
(Beobachtungsart) einer Gruppe zu. Für jede Beobachtungsgruppe wird anschließend
|
||||
der minimale und maximale Index im Beobachtungsvektor bestimmt.
|
||||
|
||||
Unterstützte Beobachtungsgruppen sind u.a.:
|
||||
|
||||
- SD : Tachymeter-Strecken,
|
||||
- R : Tachymeter-Richtungen,
|
||||
- ZW : Tachymeter-Zenitwinkel,
|
||||
- gnss: GNSS-Basislinienkomponenten (bx/by/bz),
|
||||
- niv : Geometrisches Nivellement,
|
||||
- lA : Pseudobeobachtungen.
|
||||
|
||||
Die zurückgegebenen Indexbereiche dienen insbesondere zur Blockzerlegung von
|
||||
Kofaktor- oder Gewichtsmatrizen (z.B. bei Varianzkomponentenschätzung oder
|
||||
weicher Datumsfestlegung).
|
||||
|
||||
:param Jacobimatrix_symbolisch_liste_beobachtungsvektor: Liste der Beobachtungen.
|
||||
:type Jacobimatrix_symbolisch_liste_beobachtungsvektor: list
|
||||
:return: Dictionary mit Indexbereichen je Beobachtungsgruppe.
|
||||
:rtype: dict
|
||||
:raises ValueError: Wenn für eine Beobachtungsgruppe keine Indizes ermittelt werden können.
|
||||
"""
|
||||
liste_strecken_indizes = []
|
||||
liste_richtungen_indizes = []
|
||||
liste_zenitwinkel_indizes = []
|
||||
@@ -156,7 +198,7 @@ class Datumsfestlegung:
|
||||
return names, {n: i for i, n in enumerate(names)}
|
||||
|
||||
@staticmethod
|
||||
def build_G_from_names(x0, unbekannten_liste, liste_punktnummern=None, mit_massstab=True):
|
||||
def erstelle_G(x0, unbekannten_liste, liste_punktnummern=None, mit_massstab=True):
|
||||
"""
|
||||
Baut G (u x d) in den vollen Unbekanntenraum.
|
||||
Wenn liste_punktnummern=None, werden alle Punkt-IDs aus unbekannten_liste
|
||||
@@ -235,11 +277,8 @@ class Datumsfestlegung:
|
||||
return E
|
||||
|
||||
|
||||
def loese_geraendert_mit_Qxx(N, n, G):
|
||||
"""
|
||||
löst [N G; G^T 0] [dx;k] = [n;0]
|
||||
und liefert zusätzlich Q_xx als oberen linken Block von inv(K).
|
||||
"""
|
||||
def berechne_dx(N, n, G):
|
||||
|
||||
N = np.asarray(N, float)
|
||||
n = np.asarray(n, float).reshape(-1, 1)
|
||||
G = np.asarray(G, float)
|
||||
|
||||
85
Einheitenumrechnung.py
Normal file
85
Einheitenumrechnung.py
Normal file
@@ -0,0 +1,85 @@
|
||||
from decimal import Decimal
|
||||
import sympy as sp
|
||||
import math
|
||||
|
||||
class Einheitenumrechnung:
|
||||
"""Einheitenumrechnungen für Winkel- und Längeneinheiten.
|
||||
|
||||
Die Klasse stellt Methoden zur Verfügung für:
|
||||
|
||||
- Umrechnung von Millibogensekunden (mas) in Radiant,
|
||||
- Umrechnung von Millimetern (mm) in Meter,
|
||||
- Umrechnung von Gon und Milligon (mgon) in Radiant (Decimal-basiert).
|
||||
"""
|
||||
|
||||
def mas_to_rad(mas: float) -> float:
|
||||
"""Rechnet Millibogensekunden (mas) in Radiant um.
|
||||
|
||||
Es gilt: rad = mas * (pi / (180 * 3600 * 1000)).
|
||||
|
||||
:param mas: Winkel in Millibogensekunden (mas).
|
||||
:type mas: float
|
||||
:return: Winkel in Radiant.
|
||||
:rtype: float
|
||||
"""
|
||||
umrechnungsfaktor = 1 / 1000 * 1 / 3600 * sp.pi / 180
|
||||
grad = mas * umrechnungsfaktor
|
||||
return grad
|
||||
|
||||
def mm_to_m(mm: float) -> float:
|
||||
"""Rechnet Millimeter in Meter um.
|
||||
|
||||
Es gilt: m = mm / 1000.
|
||||
|
||||
:param mm: Länge in Millimeter.
|
||||
:type mm: float
|
||||
:return: Länge in Meter.
|
||||
:rtype: float
|
||||
"""
|
||||
m = mm / 1000
|
||||
return m
|
||||
|
||||
def gon_to_rad_Decimal(gon: float) -> Decimal:
|
||||
"""Rechnet Gon in Radiant um (Decimal-basiert).
|
||||
|
||||
Es gilt: 400 gon = 2*pi und damit rad = (gon / 200) * pi.
|
||||
|
||||
:param gon: Winkel in Gon.
|
||||
:type gon: float
|
||||
:return: Winkel in Radiant als Decimal.
|
||||
:rtype: Decimal
|
||||
"""
|
||||
gon = Decimal(gon)
|
||||
pi = Decimal(str(math.pi))
|
||||
rad = (gon / Decimal(200)) * pi
|
||||
return rad
|
||||
|
||||
def mgon_to_rad_Decimal(gon: float) -> Decimal:
|
||||
"""Rechnet Milligon (mgon) in Radiant um (Decimal-basiert).
|
||||
|
||||
Es gilt: 1 mgon = 0.001 gon und damit rad = (mgon / 200000) * pi.
|
||||
|
||||
:param gon: Winkel in Milligon (mgon).
|
||||
:type gon: float
|
||||
:return: Winkel in Radiant als Decimal.
|
||||
:rtype: Decimal
|
||||
"""
|
||||
gon = Decimal(gon)
|
||||
pi = Decimal(str(math.pi))
|
||||
rad = (gon / Decimal(200000)) * pi
|
||||
return rad
|
||||
|
||||
def rad_to_gon_Decimal(rad: float) -> Decimal:
|
||||
"""Rechnet Radiant in Gon um (Decimal-basiert).
|
||||
|
||||
Es gilt: 400 gon = 2*pi und damit rad = (gon / 200) * pi.
|
||||
|
||||
:param rad: Winkel in Rad.
|
||||
:type rad: float
|
||||
:return: Winkel in Gon als Decimal.
|
||||
:rtype: Decimal
|
||||
"""
|
||||
rad = Decimal(rad)
|
||||
pi = Decimal(str(math.pi))
|
||||
gon = (rad / pi) * Decimal(200)
|
||||
return gon
|
||||
11
Import.py
11
Import.py
@@ -6,6 +6,7 @@ import xml.etree.ElementTree as ET
|
||||
|
||||
from Berechnungen import Berechnungen
|
||||
import Berechnungen
|
||||
from Einheitenumrechnung import Einheitenumrechnung
|
||||
|
||||
|
||||
class Import:
|
||||
@@ -625,7 +626,7 @@ class Import:
|
||||
|
||||
if Import_fortsetzen:
|
||||
nummer_zielpunkt = 0
|
||||
# Abfragen der aktuell höschten Nummer im Attribut beobachtungsgruppeID der Tabelle Beobachtungen
|
||||
# Abfragen der aktuell höchsten Nummer im Attribut beobachtungsgruppeID der Tabelle Beobachtungen
|
||||
try:
|
||||
nummer_beobachtungsgruppeID = max(liste_beobachtungsgruppeID)
|
||||
except:
|
||||
@@ -741,7 +742,7 @@ class Import:
|
||||
richtung2 = self.string_to_decimal(liste[5]) - Decimal(200)
|
||||
zenitwinkel_vollsatz_gon = (self.string_to_decimal(liste_aktueller_zielpunkt[6]) - self.string_to_decimal(
|
||||
liste[6]) + 400) / 2
|
||||
zenitwinkel_vollsatz_rad = Berechnungen.Einheitenumrechnung.gon_to_rad_Decimal(zenitwinkel_vollsatz_gon)
|
||||
zenitwinkel_vollsatz_rad = Einheitenumrechnung.gon_to_rad_Decimal(zenitwinkel_vollsatz_gon)
|
||||
distanz_vollsatz = (self.string_to_decimal(liste_aktueller_zielpunkt[7]) + self.string_to_decimal(
|
||||
liste[7])) / 2
|
||||
if richtung2 < 0 and richtung1 != Decimal(0):
|
||||
@@ -749,7 +750,7 @@ class Import:
|
||||
elif richtung2 > 400:
|
||||
richtung2 -= Decimal(400)
|
||||
richtung_vollsatz_gon = (richtung1 + richtung2) / 2
|
||||
richtung_vollsatz_rad = Berechnungen.Einheitenumrechnung.gon_to_rad_Decimal(richtung_vollsatz_gon)
|
||||
richtung_vollsatz_rad = Einheitenumrechnung.gon_to_rad_Decimal(richtung_vollsatz_gon)
|
||||
|
||||
if liste_aktueller_zielpunkt[8] == liste[8]:
|
||||
prismenhoehe = liste_aktueller_zielpunkt[8]
|
||||
@@ -891,7 +892,7 @@ class Import:
|
||||
|
||||
# Es werden nur Höhendifferenzen für Punkte berechnet, für die Näherungskoordinaten in der Datenbank vorliegen.
|
||||
if Import_fortsetzen:
|
||||
print(f"Für folgende Nivellementpunkte werden die Höhen in der Ausgleichung berechnet: {liste_punktnummern_in_db}\nFür folgende Punkte wird aktuell keine Höhe in der Ausgleichung berechnet: {liste_punktnummern_nicht_in_db}. Bei Bedarf im folgenden Schritt ändern!")
|
||||
print(f"Folgende Stand- und Zielpunkte des geometrischen Nivellements werden für die Beobachtungsgruppe ausgeglichen: {liste_punktnummern_in_db}\nFür folgende Punkte wird aktuell keine Höhe in der Ausgleichung berechnet: {liste_punktnummern_nicht_in_db}. Bei Bedarf im folgenden Schritt ändern!")
|
||||
return dict_punkt_mittelwert_punkthoehen, liste_punktnummern_in_db
|
||||
|
||||
def import_beobachtungen_nivellement_naeherung_punkthoehen(self, dict_punkt_mittelwert_punkthoehen: dict,
|
||||
@@ -1173,7 +1174,7 @@ class Import:
|
||||
return f"Die Beobachtungen aus der Datei {pfad_datei} wurden erfolgreich importiert."
|
||||
|
||||
else:
|
||||
print(f"Anzahl nicht RVVR durch 4 teilbar. Bitte die Datei {pfad_datei} überprüfen! Der Import wurde abgebrochen.")
|
||||
print(f"Anzahl RVVR durch 4 teilbar. Bitte die Datei {pfad_datei} überprüfen! Der Import wurde abgebrochen.")
|
||||
Import_fortsetzen = False
|
||||
|
||||
def import_koordinaten_gnss(self, pfad_datei: str, liste_sapos_stationen_genauigkeiten: list) -> str:
|
||||
|
||||
@@ -270,7 +270,7 @@ class Transformationen:
|
||||
Zi = l_berechnet_final[3 * i + 2]
|
||||
liste_l_berechnet_final.append(sp.Matrix([Xi, Yi, Zi]))
|
||||
|
||||
print("l_berechnet_final:")
|
||||
print("Koordinaten berechnet aus Helmerttransformation:")
|
||||
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}")
|
||||
|
||||
|
||||
564
Netzqualitaet_Genauigkeit.py
Normal file
564
Netzqualitaet_Genauigkeit.py
Normal file
@@ -0,0 +1,564 @@
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from scipy.stats import f
|
||||
import pandas as pd
|
||||
from IPython.display import HTML
|
||||
from IPython.display import display, clear_output
|
||||
import Berechnungen
|
||||
import Einheitenumrechnung
|
||||
|
||||
|
||||
class Genauigkeitsmaße:
|
||||
"""Berechnung von Genauigkeitsmaße zur Bewertung der erreichten Netzqualität.
|
||||
|
||||
Die Klasse stellt Methoden zur Verfügung für:
|
||||
|
||||
- Berechnung der a-posteriori Standardabweichung der Gewichtseinheit s₀
|
||||
- Berechnung des Helmertschen Punktfehlers (2D/3D),
|
||||
- Berechnung der Standardellipse (Helmertschen Fehlerellipse),
|
||||
- Berechnung der Konfidenzellipse auf Basis eines Konfidenzniveaus (alpha) mit Skalierung über die F-Verteilung,
|
||||
- Berechnung von Konfidenzellipsen im lokalen ENU-System durch Transformation von Qxx → Qxx_ENU,
|
||||
inkl. Ausgabe/Export tabellarischer Ergebnisse.
|
||||
"""
|
||||
|
||||
|
||||
@staticmethod
|
||||
def berechne_s0apost(v: np.ndarray, P: np.ndarray, r: int) -> float:
|
||||
"""
|
||||
Berechnet die a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
|
||||
Die a-posteriori Standardabweichung dient als Qualitätsmaß für die
|
||||
Ausgleichung nach der Methode der kleinsten Quadrate. Dabei beschreibt
|
||||
r die Redundanz (Freiheitsgrade).
|
||||
|
||||
:param v: Residuenvektor der Beobachtungen.
|
||||
:type v: numpy.ndarray
|
||||
:param P: Gewichtsmatrix der Beobachtungen.
|
||||
:type P: numpy.ndarray
|
||||
:param r: Redundanz bzw. Anzahl der Freiheitsgrade der Ausgleichung.
|
||||
:type r: int
|
||||
:return: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:rtype: float
|
||||
"""
|
||||
|
||||
vTPv_matrix = v.T @ P @ v
|
||||
vTPv = vTPv_matrix.item()
|
||||
s0apost = np.sqrt(vTPv / r)
|
||||
print(f"s0 a posteriori beträgt: {s0apost:.4f}")
|
||||
return s0apost
|
||||
|
||||
|
||||
@staticmethod
|
||||
def helmert_punktfehler(Qxx, s0_apost, unbekannten_liste):
|
||||
"""
|
||||
Berechnet den Helmertschen Punktfehler (2D/3D) anhand der Standardabweichungen der Koordinaten der Punkte.
|
||||
|
||||
Aus der Kofaktor-Matrix der Unbekannten Qxx werden die Kofaktoren punktweise ausgelesen. Durch Multiplikation
|
||||
mit der a-posteriori Standardabweichung der Gewichtseinheit s₀ werden die Standardabweichungen je Koordinate
|
||||
(σx, σy, σz) sowie der Helmertsche Punktfehler σP berechnet:
|
||||
|
||||
Die Punktzuordnung erfolgt über die Symbolnamen der Unbekanntenliste (z.B. X1, Y1, Z1).
|
||||
Die Dimension (2D/3D) wird interaktiv per Eingabe abgefragt. Zusätzlich werden die
|
||||
Ergebnisse als Tabelle ausgegeben und in eine Excel-Datei exportiert.
|
||||
|
||||
:param Qxx: Kofaktor-Matrix der Unbekannten.
|
||||
:type Qxx: numpy.ndarray
|
||||
:param s0_apost: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:type s0_apost: float
|
||||
:param unbekannten_liste: Liste der Unbekannten.
|
||||
:type unbekannten_liste: list
|
||||
:return: Tabelle mit Standardabweichungen und Helmertschem Punktfehler je Punkt.
|
||||
:rtype: pandas.DataFrame
|
||||
:raises ValueError: Wenn eine ungültige Dimension (nicht 2 oder 3) eingegeben wird.
|
||||
"""
|
||||
dim = int(input("Helmertscher Punktfehler (2 = 2D, 3 = 3D): "))
|
||||
diagQ = np.diag(Qxx)
|
||||
daten = []
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
punkt_ids = []
|
||||
for n in namen_str:
|
||||
if n.upper().startswith('X'):
|
||||
punkt_ids.append(n[1:])
|
||||
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
|
||||
idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
|
||||
|
||||
qx = diagQ[idx_x]
|
||||
qy = diagQ[idx_y]
|
||||
qz = 0.0
|
||||
|
||||
if dim == 3:
|
||||
try:
|
||||
idx_z = next(i for i, n in enumerate(namen_str) if n.upper() == f"Z{pid}".upper())
|
||||
qz = diagQ[idx_z]
|
||||
except StopIteration:
|
||||
qz = 0.0
|
||||
|
||||
sx = s0_apost * np.sqrt(qx)
|
||||
sy = s0_apost * np.sqrt(qy)
|
||||
sz = s0_apost * np.sqrt(qz) if dim == 3 else 0
|
||||
|
||||
sP = s0_apost * np.sqrt(qx + qy + qz)
|
||||
|
||||
daten.append([pid, float(sx), float(sy), float(sz), float(sP)])
|
||||
except:
|
||||
continue
|
||||
|
||||
helmert_punktfehler = pd.DataFrame(daten, columns=["Punkt", "σx", "σy", "σz", f"σP_{dim}D"])
|
||||
display(HTML(helmert_punktfehler.to_html(index=False)))
|
||||
helmert_punktfehler.to_excel(r"Zwischenergebnisse\Standardabweichungen_Helmertscher_Punktfehler.xlsx",index=False)
|
||||
return helmert_punktfehler
|
||||
|
||||
|
||||
|
||||
@staticmethod
|
||||
def standardellipse(Qxx, s0_apost, unbekannten_liste):
|
||||
"""
|
||||
Berechnet die Standardellipse (Helmertsche Fehlerellipse) für die Punkte aus Qxx und s₀ a posteriori.
|
||||
|
||||
Für jeden Punkt werden aus der Kofaktor-Matrix der Unbekannten Qxx die
|
||||
Kofaktoren von X und Y ausgelesen (qxx, qyy, qyx).
|
||||
Daraus werden Standardabweichungen σx, σy sowie die Kovarianz σxy bestimmt und
|
||||
anschließend die Parameter der Standardellipse berechnet:
|
||||
|
||||
- Große und kleine Halbachse der Standardellipse,
|
||||
- Richtungswinkel θ der großen Halbachse in gon.
|
||||
|
||||
Die Punktzuordnung erfolgt über die Symbolnamen der Unbekanntenliste (z.B. X1, Y1).
|
||||
Zusätzlich werden die Ergebnisse tabellarisch ausgegeben und in eine Excel-Datei expoertiert.
|
||||
|
||||
:param Qxx: Kofaktor-Matrix der Unbekannten.
|
||||
:type Qxx: numpy.ndarray
|
||||
:param s0_apost: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:type s0_apost: float
|
||||
:param unbekannten_liste: Liste der Unbekannten.
|
||||
:type unbekannten_liste: list
|
||||
:return: Tabelle mit Standardabweichungen und Parametern der Standardellipse je Punkt.
|
||||
:rtype: pandas.DataFrame
|
||||
"""
|
||||
Qxx = np.asarray(Qxx, float)
|
||||
daten = []
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
punkt_ids = []
|
||||
for n in namen_str:
|
||||
if n.upper().startswith('X'):
|
||||
punkt_ids.append(n[1:])
|
||||
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
|
||||
idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
|
||||
|
||||
qxx = Qxx[idx_x, idx_x]
|
||||
qyy = Qxx[idx_y, idx_y]
|
||||
qyx = Qxx[idx_y, idx_x]
|
||||
|
||||
# Standardabweichungen
|
||||
sx = s0_apost * np.sqrt(qxx)
|
||||
sy = s0_apost * np.sqrt(qyy)
|
||||
sxy = (s0_apost ** 2) * qyx
|
||||
|
||||
k = np.sqrt((qxx - qyy) ** 2 + 4 * (qyx ** 2))
|
||||
|
||||
# Q_dmax/min = 0.5 * (Qyy + Qxx +/- k)
|
||||
q_dmax = 0.5 * (qyy + qxx + k)
|
||||
q_dmin = 0.5 * (qyy + qxx - k)
|
||||
|
||||
# Halbachsen
|
||||
s_max = s0_apost * np.sqrt(q_dmax)
|
||||
s_min = s0_apost * np.sqrt(q_dmin)
|
||||
|
||||
# Richtungswinkel theta in gon:
|
||||
zaehler = 2 * qyx
|
||||
nenner = qxx - qyy
|
||||
t_grund = 0.5 * np.arctan(abs(zaehler) / abs(nenner)) * (200 / np.pi)
|
||||
|
||||
# Quadrantenabfrage
|
||||
if nenner > 0 and qyx > 0: # Qxx - Qyy > 0 und Qyx > 0
|
||||
t_gon = t_grund # 0 - 50 gon
|
||||
elif nenner < 0 and qyx > 0: # Qxx - Qyy < 0 und Qyx > 0
|
||||
t_gon = 100 - t_grund # 50 - 100 gon
|
||||
elif nenner < 0 and qyx < 0: # Qxx - Qyy < 0 und Qyx < 0
|
||||
t_gon = 100 + t_grund # 100 - 150 gon
|
||||
elif nenner > 0 and qyx < 0: # Qxx - Qyy > 0 und Qyx < 0
|
||||
t_gon = 200 - t_grund # 150 - 200 gon
|
||||
else:
|
||||
t_gon = 0.0
|
||||
|
||||
daten.append([
|
||||
pid,
|
||||
float(sx), float(sy), float(sxy),
|
||||
float(s_max), float(s_min),
|
||||
float(t_gon)
|
||||
])
|
||||
except:
|
||||
continue
|
||||
standardellipse = pd.DataFrame(daten, columns=["Punkt", "σx [m]", "σy [m]", "σxy [m]", "Große Halbachse [m]", "Kleine Halbachse [m]", "θ [gon]"])
|
||||
standardellipse["σx [m]"] = standardellipse["σx [m]"].astype(float).round(4)
|
||||
standardellipse["σy [m]"] = standardellipse["σy [m]"].astype(float).round(4)
|
||||
standardellipse["Große Halbachse [m]"] = standardellipse["Große Halbachse [m]"].astype(float).round(4)
|
||||
standardellipse["Kleine Halbachse [m]"] = standardellipse["Kleine Halbachse [m]"].astype(float).round(4)
|
||||
standardellipse["θ [gon]"] = standardellipse["θ [gon]"].astype(float).round(3)
|
||||
display(HTML(standardellipse.to_html(index=False)))
|
||||
standardellipse.to_excel(r"Zwischenergebnisse\Standardellipse.xlsx", index=False)
|
||||
return standardellipse
|
||||
|
||||
|
||||
|
||||
@staticmethod
|
||||
def konfidenzellipse(Qxx, s0_apost, unbekannten_liste, R, ausgabe_erfolgt):
|
||||
"""
|
||||
Berechnet die Konfidenzellipse für Punkte aus Qxx und einem Konfidenzniveau.
|
||||
|
||||
Auf Basis der Kovarianz-Matrix der Unbekannten Qxx und der a-posteriori
|
||||
Standardabweichung der Gewichtseinheit s₀ werden für jeden Punkt die Parameter
|
||||
der Konfidenzellipse berechnet. Das Konfidenzniveau wird mittels einer Eingabe
|
||||
über alpha gewählt (Standard: 0.05 für 95%).
|
||||
|
||||
Die Punktzuordnung erfolgt über die Symbolnamen der Unbekanntenliste (z.B. X1, Y1).
|
||||
Optional wird die Tabelle ausgegeben und als Excel-Datei exportiert, abhängig von
|
||||
ausgabe_erfolgt.
|
||||
|
||||
:param Qxx: Kofaktor-Matrix der geschätzten Unbekannten.
|
||||
:type Qxx: numpy.ndarray
|
||||
:param s0_apost: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:type s0_apost: float
|
||||
:param unbekannten_liste: Liste der Unbekannten.
|
||||
:type unbekannten_liste: list
|
||||
:param R: Redundanz (Freiheitsgrade) für die F-Verteilung.
|
||||
:type R: int
|
||||
:param ausgabe_erfolgt: Steuert, ob eine Ausgabe/Dateischreibung erfolgen soll (False = Ausgabe).
|
||||
:type ausgabe_erfolgt: bool
|
||||
:return: Tabelle der Konfidenzellipse je Punkt, verwendetes alpha.
|
||||
:rtype: tuple[pandas.DataFrame, float]
|
||||
:raises ValueError: Wenn alpha nicht in (0, 1) liegt oder nicht in float umgewandelt werden kann.
|
||||
"""
|
||||
alpha_input = input("Konfidenzniveau wählen (z.B. 0.05 für 95%, 0.01 für 99%) [Standard=0.05]: ")
|
||||
if alpha_input.strip() == "":
|
||||
alpha = 0.05
|
||||
else:
|
||||
alpha = float(alpha_input)
|
||||
print(f"→ Verwende alpha = {alpha} (Konfidenz = {(1 - alpha) * 100:.1f}%)")
|
||||
Qxx = np.asarray(Qxx, float)
|
||||
daten = []
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
punkt_ids = [n[1:] for n in namen_str if n.upper().startswith('X')]
|
||||
|
||||
# Faktor für Konfidenzellipse (F-Verteilung)
|
||||
kk = float(np.sqrt(2.0 * f.ppf(1.0 - alpha, 2, R)))
|
||||
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
|
||||
idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
|
||||
|
||||
qxx = Qxx[idx_x, idx_x]
|
||||
qyy = Qxx[idx_y, idx_y]
|
||||
qyx = Qxx[idx_y, idx_x]
|
||||
|
||||
# Standardabweichungen
|
||||
sx = s0_apost * np.sqrt(qxx)
|
||||
sy = s0_apost * np.sqrt(qyy)
|
||||
sxy = (s0_apost ** 2) * qyx
|
||||
|
||||
k = np.sqrt((qxx - qyy) ** 2 + 4 * (qyx ** 2))
|
||||
|
||||
# Q_dmax/min = 0.5 * (Qyy + Qxx +/- k)
|
||||
q_dmax = 0.5 * (qyy + qxx + k)
|
||||
q_dmin = 0.5 * (qyy + qxx - k)
|
||||
|
||||
# Halbachsen der Standardellipse
|
||||
s_max = s0_apost * np.sqrt(q_dmax)
|
||||
s_min = s0_apost * np.sqrt(q_dmin)
|
||||
|
||||
# Halbachsen der Konfidenzellipse
|
||||
A_K = kk * s_max
|
||||
B_K = kk * s_min
|
||||
|
||||
# Richtungswinkel theta in gon:
|
||||
zaehler = 2 * qyx
|
||||
nenner = qxx - qyy
|
||||
t_grund = 0.5 * np.arctan(abs(zaehler) / abs(nenner)) * (200 / np.pi)
|
||||
|
||||
# Quadrantenabfrage
|
||||
if nenner > 0 and qyx > 0:
|
||||
t_gon = t_grund # 0 - 50 gon
|
||||
elif nenner < 0 and qyx > 0:
|
||||
t_gon = 100 - t_grund # 50 - 100 gon
|
||||
elif nenner < 0 and qyx < 0:
|
||||
t_gon = 100 + t_grund # 100 - 150 gon
|
||||
elif nenner > 0 and qyx < 0:
|
||||
t_gon = 200 - t_grund # 150 - 200 gon
|
||||
else:
|
||||
t_gon = 0.0
|
||||
|
||||
daten.append([
|
||||
pid,
|
||||
float(sx), float(sy), float(sxy),
|
||||
float(A_K), float(B_K),
|
||||
float(t_gon)
|
||||
])
|
||||
|
||||
except:
|
||||
continue
|
||||
|
||||
konfidenzellipse = pd.DataFrame(daten, columns=["Punkt", "σx [m]", "σy [m]", "σxy [m]", "Große Halbachse [m]", "Kleine Halbachse [m]", "θ [gon]"])
|
||||
konfidenzellipse["Große Halbachse [m]"] = konfidenzellipse["Große Halbachse [m]"].round(4)
|
||||
konfidenzellipse["Kleine Halbachse [m]"] = konfidenzellipse["Kleine Halbachse [m]"].round(4)
|
||||
konfidenzellipse["θ [gon]"] = konfidenzellipse["θ [gon]"].round(3)
|
||||
if ausgabe_erfolgt == False:
|
||||
display(HTML(konfidenzellipse.to_html(index=False)))
|
||||
konfidenzellipse.to_excel(r"Zwischenergebnisse\Konfidenzellipse.xlsx", index=False)
|
||||
return konfidenzellipse, alpha
|
||||
|
||||
|
||||
@staticmethod
|
||||
def konfidenzellipsen_enu(a, b, ausgabe_parameterschaetzung, liste_unbekannte, ausgleichungsergebnis, s0apost, r_gesamt):
|
||||
"""
|
||||
Berechnet Konfidenzellipsen im lokalen ENU-System aus einer ins ENU-System transformierten Qxx-Matrix.
|
||||
|
||||
Die Funktion transformiert zunächst die Kofaktor-Matrix der Unbekannten Qxx
|
||||
in ein East-North-Up-System (ENU) bezogen auf den Schwerpunkt der verwendeten
|
||||
Punkte (B0, L0). Anschließend wird auf Basis der transformierten Matrix die
|
||||
Konfidenzellipse über die Funktion "konfidenzellipse" bestimmt.
|
||||
Zum Schluss werden Spaltennamen an die ENU-Notation angepasst, Werte gerundet,
|
||||
tabellarisch ausgegeben und als Excel-Datei exportiert.
|
||||
|
||||
:param a: Große Halbachse a des Referenzellipsoids (z.B. WGS84/GRS80) in Metern.
|
||||
:type a: float
|
||||
:param b: Große Halbachse b des Referenzellipsoids (z.B. WGS84/GRS80) in Metern.
|
||||
:type b: float
|
||||
:param ausgabe_parameterschaetzung: Dictonary der Ergebnisse der Parameterschätzung, muss "Q_xx" enthalten.
|
||||
:type ausgabe_parameterschaetzung: dict
|
||||
:param liste_unbekannte: Liste der Unbekannten.
|
||||
:type liste_unbekannte: list
|
||||
:param ausgleichungsergebnis: Dictionary der geschätzten Punktkoordinaten (XYZ) zur ENU-Referenzbildung.
|
||||
:type ausgleichungsergebnis: dict
|
||||
:param s0apost: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:type s0apost: float
|
||||
:param r_gesamt: Redundanz (Freiheitsgrade) für die Konfidenzberechnung.
|
||||
:type r_gesamt: int
|
||||
:return: Tabelle der Konfidenzellipse im ENU-System, Rotationsmatrix R0 der ENU-Transformation.
|
||||
:rtype: tuple[pandas.DataFrame, numpy.ndarray]
|
||||
:raises KeyError: Wenn ``ausgabe_parameterschaetzung`` keinen Eintrag ``"Q_xx"`` enthält.
|
||||
"""
|
||||
|
||||
berechnungen = Berechnungen.Berechnungen(a, b)
|
||||
|
||||
# 1) Qxx ins ENU-System transformieren
|
||||
Qxx_enu, (B0, L0), R0 = Berechnungen.ENU.transform_Qxx_zu_QxxENU(
|
||||
Qxx=ausgabe_parameterschaetzung["Q_xx"],
|
||||
unbekannten_liste= liste_unbekannte,
|
||||
berechnungen=berechnungen,
|
||||
dict_xyz= ausgleichungsergebnis,
|
||||
)
|
||||
|
||||
print(
|
||||
f"ENU-Referenz (Schwerpunkt): B0={Einheitenumrechnung.rad_to_gon_Decimal(B0):.8f} rad, L0={Einheitenumrechnung.rad_to_gon_Decimal(L0):.8f} rad")
|
||||
|
||||
# 2) Konfidenzellipse im ENU-System
|
||||
Konfidenzellipse_ENU, alpha = Genauigkeitsmaße.konfidenzellipse(
|
||||
Qxx_enu,
|
||||
s0apost,
|
||||
liste_unbekannte,
|
||||
r_gesamt,
|
||||
ausgabe_erfolgt = True
|
||||
)
|
||||
|
||||
# 3) Spaltennamen anpassen
|
||||
Konfidenzellipse_ENU = Konfidenzellipse_ENU.rename(columns={
|
||||
"σx [m]": "σE [m]",
|
||||
"σy [m]": "σN [m]",
|
||||
"σxy [m]": "σEN [m]",
|
||||
"θ [gon]": "θ_EN [gon]"
|
||||
})
|
||||
|
||||
# 4) Runden und Anzeigen
|
||||
Konfidenzellipse_ENU["σE [m]"] = Konfidenzellipse_ENU["σE [m]"].round(4)
|
||||
Konfidenzellipse_ENU["σN [m]"] = Konfidenzellipse_ENU["σN [m]"].round(4)
|
||||
Konfidenzellipse_ENU["Große Halbachse [m]"] = Konfidenzellipse_ENU["Große Halbachse [m]"].round(4)
|
||||
Konfidenzellipse_ENU["Kleine Halbachse [m]"] = Konfidenzellipse_ENU["Kleine Halbachse [m]"].round(4)
|
||||
Konfidenzellipse_ENU["θ_EN [gon]"] = Konfidenzellipse_ENU["θ_EN [gon]"].round(4)
|
||||
|
||||
display(HTML(Konfidenzellipse_ENU.to_html(index=False)))
|
||||
|
||||
# 5) Export
|
||||
Konfidenzellipse_ENU.to_excel(r"Zwischenergebnisse\Konfidenzellipse_ENU.xlsx", index=False)
|
||||
return Konfidenzellipse_ENU, R0
|
||||
|
||||
|
||||
|
||||
class Plot:
|
||||
"""Visualisierung geodätischer Netze und Genauigkeitsmaße.
|
||||
|
||||
Die Klasse stellt Methoden zur Verfügung für:
|
||||
|
||||
- grafische Darstellung von geodätischen Netzen im lokalen ENU-System,
|
||||
- Visualisierung von Beobachtungen als Verbindungslinien,
|
||||
- Darstellung von Konfidenzellipsen,
|
||||
- interaktive Netzdarstellung mit Plotly inklusive Hover-Informationen,
|
||||
- Skalierung und Layout-Anpassung zur anschaulichen Präsentation von
|
||||
Lagegenauigkeiten.
|
||||
|
||||
Die Klasse dient ausschließlich der Ergebnisvisualisierung und nimmt keine
|
||||
numerischen Berechnungen vor.
|
||||
"""
|
||||
|
||||
@staticmethod
|
||||
def netzplot_ellipsen(
|
||||
Koord_ENU,
|
||||
unbekannten_labels,
|
||||
beobachtungs_labels,
|
||||
df_konf_ellipsen_enu,
|
||||
skalierung=1000,
|
||||
n_ellipse_pts=60,
|
||||
titel="Netzplot im ENU-System mit Konfidenzellipsen"
|
||||
):
|
||||
"""
|
||||
Erstellt einen Netzplot im ENU-System inklusive Konfidenzellipsen, Netzpunkten und Beobachtungslinien.
|
||||
|
||||
Die Funktion visualisiert das geodätische Netz im East-North-Up-System (ENU)
|
||||
mit Plotly. Dabei werden:
|
||||
|
||||
- Beobachtungen als Verbindungslinien zwischen Punkten dargestellt, deren Ansicht aus- und eingeschaltet werden kann,
|
||||
- Konfidenzellipsen je Punkt (Halbachsen und Richtungswinkel),
|
||||
- Netzpunkte mit Punkt-ID und Koordinaten im Hover-Text angezeigt.
|
||||
|
||||
Die Ellipsen werden zur besseren Sichtbarkeit mit einem Faktor "skalierung" vergrößert. Dieser kann angepasst werden.
|
||||
Der Richtungswinkel wird in gon erwartet und intern nach Radiant umgerechnet.
|
||||
|
||||
:param Koord_ENU: Dictionary der Punktkoordinaten im ENU-System.
|
||||
:type Koord_ENU: dict
|
||||
:param unbekannten_labels: Liste der Unbekannten zur Ableitung der Punkt-IDs (z.B. X1, Y1, Z1).
|
||||
:type unbekannten_labels: list
|
||||
:param beobachtungs_labels: Liste der Beobachtungen zur Ableitung von Verbindungslinien.
|
||||
:type beobachtungs_labels: list
|
||||
:param df_konf_ellipsen_enu: DataFrame mit Konfidenzellipsenparametern je Punkt.
|
||||
:type df_konf_ellipsen_enu: pandas.DataFrame
|
||||
:param skalierung: Faktor zur visuellen Vergrößerung der Ellipsen im Plot.
|
||||
:type skalierung: float
|
||||
:param n_ellipse_pts: Anzahl der Stützpunkte zur Approximation der Ellipse.
|
||||
:type n_ellipse_pts: int
|
||||
:param titel: Titel des Plots.
|
||||
:type titel: str
|
||||
:return: None
|
||||
:rtype: None
|
||||
:raises ValueError: Wenn weder "θ_EN [gon]" noch "θ [gon]" im DataFrame vorhanden ist.
|
||||
"""
|
||||
|
||||
names = [str(s).strip() for s in unbekannten_labels]
|
||||
|
||||
if "θ_EN [gon]" in df_konf_ellipsen_enu.columns:
|
||||
theta_col = "θ_EN [gon]"
|
||||
elif "θ [gon]" in df_konf_ellipsen_enu.columns:
|
||||
theta_col = "θ [gon]"
|
||||
else:
|
||||
raise ValueError("Spalte 'θ_EN [gon]' oder 'θ [gon]' fehlt im DataFrame.")
|
||||
|
||||
punkt_ids = sorted({nm[1:] for nm in names if nm and nm[0].upper() in ("X", "Y", "Z")})
|
||||
|
||||
fig = go.Figure()
|
||||
|
||||
# 1) Darstellungen der Beobachtungen
|
||||
beob_typen = {
|
||||
'GNSS-Basislinien': {'pattern': 'gnss', 'color': 'rgba(255, 100, 0, 0.4)'},
|
||||
'Tachymeter-Beob': {'pattern': '', 'color': 'rgba(100, 100, 100, 0.3)'}
|
||||
}
|
||||
|
||||
for typ, info in beob_typen.items():
|
||||
x_l, y_l = [], []
|
||||
for bl in beobachtungs_labels:
|
||||
bl_str = str(bl).lower()
|
||||
is_typ = ((info['pattern'] in bl_str and info['pattern'] != '') or
|
||||
(info['pattern'] == '' and 'gnss' not in bl_str and 'niv' not in bl_str))
|
||||
if not is_typ:
|
||||
continue
|
||||
|
||||
bl_raw = str(bl)
|
||||
pts = []
|
||||
for pid in punkt_ids:
|
||||
if (f"_{pid}" in bl_raw) or bl_raw.startswith(f"{pid}_"):
|
||||
if pid in Koord_ENU:
|
||||
pts.append(pid)
|
||||
|
||||
if len(pts) >= 2:
|
||||
p1, p2 = pts[0], pts[1]
|
||||
x_l.extend([Koord_ENU[p1][0], Koord_ENU[p2][0], None]) # E
|
||||
y_l.extend([Koord_ENU[p1][1], Koord_ENU[p2][1], None]) # N
|
||||
|
||||
if x_l:
|
||||
fig.add_trace(go.Scatter(x=x_l, y=y_l, mode='lines', name=typ,
|
||||
line=dict(color=info['color'], width=1)))
|
||||
|
||||
# 2) Darstellung der Konfidenzellipsen
|
||||
t = np.linspace(0, 2 * np.pi, n_ellipse_pts)
|
||||
first = True
|
||||
for _, row in df_konf_ellipsen_enu.iterrows():
|
||||
pid = str(row["Punkt"])
|
||||
if pid not in Koord_ENU:
|
||||
continue
|
||||
|
||||
a = float(row["Große Halbachse [m]"]) * skalierung
|
||||
b = float(row["Kleine Halbachse [m]"]) * skalierung
|
||||
theta = float(row[theta_col]) * np.pi / 200.0 # gon->rad
|
||||
|
||||
ex = a * np.cos(t)
|
||||
ey = b * np.sin(t)
|
||||
|
||||
c, s = np.cos(theta), np.sin(theta)
|
||||
xr = c * ex - s * ey
|
||||
yr = s * ex + c * ey
|
||||
|
||||
E0, N0, _ = Koord_ENU[pid]
|
||||
|
||||
fig.add_trace(go.Scatter(
|
||||
x=E0 + xr, y=N0 + yr,
|
||||
mode="lines",
|
||||
line=dict(color="red", width=1.5),
|
||||
name=f"Ellipsen (×{skalierung})",
|
||||
legendgroup="Ellipsen",
|
||||
showlegend=first,
|
||||
hoverinfo="skip"
|
||||
))
|
||||
first = False
|
||||
|
||||
# 3) Darstellung der Punkte
|
||||
xs, ys, texts, hovers = [], [], [], []
|
||||
for pid in punkt_ids:
|
||||
if pid not in Koord_ENU:
|
||||
continue
|
||||
E, N, U = Koord_ENU[pid]
|
||||
xs.append(E);
|
||||
ys.append(N);
|
||||
texts.append(pid)
|
||||
hovers.append(f"Punkt {pid}<br>E={E:.4f} m<br>N={N:.4f} m<br>U={U:.4f} m")
|
||||
|
||||
fig.add_trace(go.Scatter(
|
||||
x=xs, y=ys, mode="markers+text",
|
||||
text=texts, textposition="top center",
|
||||
marker=dict(size=8, color="black"),
|
||||
name="Netzpunkte",
|
||||
hovertext=hovers, hoverinfo="text"
|
||||
))
|
||||
|
||||
fig.update_layout(
|
||||
title=f"{titel} (Ellipsen ×{skalierung})",
|
||||
xaxis=dict(title="E [m]", scaleanchor="y", scaleratio=1, showgrid=True, gridcolor="lightgrey"),
|
||||
yaxis=dict(title="N [m]", showgrid=True, gridcolor="lightgrey"),
|
||||
width=1100, height=900,
|
||||
template="plotly_white",
|
||||
plot_bgcolor="white"
|
||||
)
|
||||
|
||||
fig.add_annotation(
|
||||
text=f"<b>Maßstab Ellipsen:</b><br>Dargestellte Größe = Konfidenzellipse × {skalierung}",
|
||||
align='left', showarrow=False, xref='paper', yref='paper', x=0.02, y=0.05,
|
||||
bgcolor="white", bordercolor="black", borderwidth=1
|
||||
)
|
||||
|
||||
fig.show(config={'scrollZoom': True})
|
||||
828
Netzqualitaet_Zuverlaessigkeit.py
Normal file
828
Netzqualitaet_Zuverlaessigkeit.py
Normal file
@@ -0,0 +1,828 @@
|
||||
from dataclasses import dataclass
|
||||
import numpy as np
|
||||
from scipy import stats
|
||||
from scipy.stats import norm
|
||||
import pandas as pd
|
||||
from IPython.display import HTML
|
||||
from IPython.display import display, clear_output
|
||||
import ipywidgets as widgets
|
||||
import itables
|
||||
from itables.widget import ITable
|
||||
|
||||
|
||||
@dataclass
|
||||
class Zuverlaessigkeit:
|
||||
|
||||
@staticmethod
|
||||
def gesamtredundanz(n, u):
|
||||
"""
|
||||
Berechnet die Gesamtredundanz des Netzes.
|
||||
|
||||
Die Gesamtredundanz ergibt sich aus der Differenz zwischen der Anzahl der
|
||||
Beobachtungen n und der Anzahl der Unbekannten u. Sie entspricht der Anzahl
|
||||
der Freiheitsgrade.
|
||||
|
||||
:param n: Anzahl der Beobachtungen.
|
||||
:type n: int
|
||||
:param u: Anzahl der Unbekannten.
|
||||
:type u: int
|
||||
:return: Gesamtredundanz des Netzes.
|
||||
:rtype: int
|
||||
"""
|
||||
r_gesamt = n - u
|
||||
print(f"Die Gesamtredundanz des Netzes beträgt: {r_gesamt}")
|
||||
return r_gesamt
|
||||
|
||||
|
||||
@staticmethod
|
||||
def berechne_R(Q_vv, P):
|
||||
"""
|
||||
Berechnet die Redundanzmatrix R aus Qvv und der Gewichtsmatrix P.
|
||||
|
||||
Die Redundanzmatrix wird definiert als:
|
||||
R = Qvv · P
|
||||
|
||||
:param Q_vv: Kofaktor-Matrix der Residuen.
|
||||
:type Q_vv: numpy.ndarray
|
||||
:param P: Gewichtsmatrix der Beobachtungen.
|
||||
:type P: numpy.ndarray
|
||||
:return: Redundanzmatrix R.
|
||||
:rtype: numpy.ndarray
|
||||
"""
|
||||
R = Q_vv @ P
|
||||
return R
|
||||
|
||||
|
||||
@staticmethod
|
||||
def berechne_ri(R):
|
||||
"""
|
||||
Berechnet die Redundanzanteile einzelner Beobachtungen.
|
||||
|
||||
Die Redundanzanteile rᵢ ergeben sich aus den Diagonalelementen der Redundanzmatrix R.
|
||||
Zusätzlich werden die effektiven Redundanzanteile EVi in Prozent berechnet:
|
||||
|
||||
EVi = 100 · rᵢ
|
||||
|
||||
:param R: Redundanzmatrix R.
|
||||
:type R: numpy.ndarray
|
||||
:return: Tuple aus Redundanzanteilen rᵢ und effektiven Redundanzanteilen EVi in Prozent.
|
||||
:rtype: tuple[numpy.ndarray, numpy.ndarray]
|
||||
"""
|
||||
ri = np.diag(R)
|
||||
EVi = 100.0 * ri
|
||||
return ri, EVi
|
||||
|
||||
|
||||
@staticmethod
|
||||
def klassifiziere_ri(ri):
|
||||
"""
|
||||
Klassifiziert einen Redundanzanteil rᵢ nach seiner Kontrollierbarkeit.
|
||||
|
||||
Der Redundanzanteil wird anhand üblicher geodätischer Schwellenwerte
|
||||
qualitativ bewertet.
|
||||
|
||||
:param ri: Redundanzanteil einer einzelnen Beobachtung.
|
||||
:type ri: float
|
||||
:return: Qualitative Bewertung der Kontrollierbarkeit.
|
||||
:rtype: str
|
||||
"""
|
||||
if ri < 0.01:
|
||||
return "nicht kontrollierbar"
|
||||
elif ri < 0.10:
|
||||
return "schlecht kontrollierbar"
|
||||
elif ri < 0.30:
|
||||
return "ausreichend kontrollierbar"
|
||||
elif ri < 0.70:
|
||||
return "gut kontrollierbar"
|
||||
else:
|
||||
return "nahezu vollständig redundant"
|
||||
|
||||
|
||||
@staticmethod
|
||||
def redundanzanteile_ri(Qvv, P, liste_beob):
|
||||
"""
|
||||
Berechnet und dokumentiert Redundanzanteile rᵢ und EVᵢ für alle Beobachtungen.
|
||||
|
||||
Die Ergebnisse werden als DataFrame ausgegeben, als HTML-Tabelle angezeigt und
|
||||
als Excel-Datei exportiert.
|
||||
|
||||
:param Qvv: Kofaktor-Matrix der Residuen.
|
||||
:type Qvv: numpy.ndarray
|
||||
:param P: Gewichtsmatrix der Beobachtungen.
|
||||
:type P: numpy.ndarray
|
||||
:param liste_beob: Liste der Beobachtungslabels (Zeilenbeschriftungen) zur Zuordnung der Ergebnisse.
|
||||
:type liste_beob: list
|
||||
:return: Redundanzmatrix R, Redundanzanteile rᵢ, effektive Redundanzanteile EVᵢ, Ergebnistabelle als DataFrame.
|
||||
:rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, pandas.DataFrame]
|
||||
"""
|
||||
R = Zuverlaessigkeit.berechne_R(Qvv, P)
|
||||
ri, EVi = Zuverlaessigkeit.berechne_ri(R)
|
||||
ri = np.asarray(ri).reshape(-1)
|
||||
EVi = np.asarray(EVi).reshape(-1)
|
||||
|
||||
labels = [str(s) for s in liste_beob]
|
||||
klassen = [Zuverlaessigkeit.klassifiziere_ri(r) for r in ri]
|
||||
Redundanzanteile = pd.DataFrame({"Beobachtung": labels, "r_i": ri, "EV_i [%]": EVi, "Klassifikation": klassen, })
|
||||
display(HTML(Redundanzanteile.to_html(index=False)))
|
||||
Redundanzanteile.to_excel(r"Zwischenergebnisse\Redundanzanteile.xlsx", index=False)
|
||||
return R, ri, EVi, Redundanzanteile
|
||||
|
||||
|
||||
@staticmethod
|
||||
def globaltest(r_gesamt, sigma0_apost, sigma0_apriori=1):
|
||||
"""
|
||||
Führt den Globaltest zur Prüfung des Ausgleichungsmodells durch.
|
||||
|
||||
Der Globaltest überprüft, ob die a-posteriori Standardabweichung der Gewichtseinheit σ̂₀
|
||||
mit der a-priori Annahme σ₀ vereinbar ist. Als Testgröße wird verwendet:
|
||||
|
||||
T_G = (σ̂₀²) / (σ₀²)
|
||||
|
||||
Die Entscheidung erfolgt über die F-Verteilung. Das Signifikanzniveau alpha wird interaktiv abgefragt
|
||||
(Standard: 0.001). Zusätzlich wird eine Ergebnis-Tabelle und eine Interpretation ausgegeben.
|
||||
|
||||
:param r_gesamt: Gesamtredundanz bzw. Freiheitsgrade.
|
||||
:type r_gesamt: int
|
||||
:param sigma0_apost: a-posteriori Standardabweichung der Gewichtseinheit σ̂₀.
|
||||
:type sigma0_apost: float
|
||||
:param sigma0_apriori: a-priori Standardabweichung der Gewichtseinheit σ₀ (Standard=1).
|
||||
:type sigma0_apriori: float
|
||||
:return: Dictionary mit Testparametern, Testergebnis (H₀ angenommen/verworfen) und Interpretation.
|
||||
:rtype: dict[str, Any]
|
||||
:raises ValueError: Wenn alpha nicht in (0, 1) liegt oder nicht in float umgewandelt werden kann.
|
||||
"""
|
||||
|
||||
alpha_input = input("Irrtumswahrscheinlichkeit α wählen (z.B. 0.05, 0.01) [Standard=0.001]: ").strip()
|
||||
alpha = 0.001 if alpha_input == "" else float(alpha_input)
|
||||
T_G = (sigma0_apost ** 2) / (sigma0_apriori ** 2)
|
||||
F_krit = stats.f.ppf(1 - alpha, r_gesamt, 10 ** 9)
|
||||
H0 = T_G < F_krit
|
||||
|
||||
if H0:
|
||||
interpretation = (
|
||||
"Nullhypothese H₀ angenommen.\n"
|
||||
)
|
||||
else:
|
||||
interpretation = (
|
||||
"Nullhypothese H₀ verworfen!\n"
|
||||
"Dies kann folgende Gründe haben:\n"
|
||||
"→ Es befinden sich grobe Fehler im Datenmaterial. Bitte Lokaltest durchführen und ggf. grobe Fehler im Datenmaterial entfernen.\n"
|
||||
"→ Das stochastische Modell ist zu optimistisch. Bitte Gewichte überprüfen und ggf. anpassen."
|
||||
)
|
||||
globaltest = pd.DataFrame([
|
||||
["Freiheitsgrad", r_gesamt],
|
||||
["σ̂₀ a posteriori", sigma0_apost],
|
||||
["σ₀ a priori", sigma0_apriori],
|
||||
["Signifikanzniveau α", alpha, ],
|
||||
["Testgröße T_G", T_G, ],
|
||||
["Kritischer Wert Fₖ", F_krit],
|
||||
["Nullhypothese H₀", "angenommen" if H0 else "verworfen"], ], columns=["Größe", "Wert"])
|
||||
|
||||
display(HTML(globaltest.to_html(index=False)))
|
||||
print(interpretation)
|
||||
|
||||
return {
|
||||
"r_gesamt": r_gesamt,
|
||||
"sigma0_apost": sigma0_apost,
|
||||
"sigma0_apriori": sigma0_apriori,
|
||||
"alpha": alpha,
|
||||
"T_G": T_G,
|
||||
"F_krit": F_krit,
|
||||
"H0_angenommen": H0,
|
||||
"Interpretation": interpretation,
|
||||
}
|
||||
|
||||
|
||||
def lokaltest_innere_Zuverlaessigkeit(v, Q_vv, ri, labels, s0_apost, alpha, beta):
|
||||
"""
|
||||
Führt den Lokaltest zur Grobfehlerdetektion je Beobachtung durch.
|
||||
|
||||
Auf Basis der Residuen v, der Kofaktor-Matrix der Residuen Qvv und der Redundanzanteile rᵢ
|
||||
werden für jede Beobachtung statistische Kennwerte zur Detektion grober Fehler berechnet. Dazu zählen:
|
||||
|
||||
- Grobfehlerabschätzung: GFᵢ = − vᵢ / rᵢ
|
||||
- Standardabweichung der Residuen: s_vᵢ = s₀ · √q_vᵢ (mit q_vᵢ = diag(Qvv))
|
||||
- Normierte Verbesserung: NVᵢ = |vᵢ| / s_vᵢ
|
||||
- Nichtzentralitätsparameter: δ₀ = k + k_A
|
||||
mit k aus dem zweiseitigen Normalquantil (α) und k_A aus der Testmacht (1−β)
|
||||
- Grenzwert der Aufdeckbarkeit (Minimal detektierbarer Grobfehler): GRZWᵢ = (s_vᵢ / rᵢ) · δ₀
|
||||
|
||||
Beobachtungen werden als auffällig markiert, wenn NVᵢ > δ₀. Für rᵢ = 0 wird die Grobfehlerabschätzung
|
||||
und der Grenzwert als NaN gesetzt.
|
||||
|
||||
:param v: Residuenvektor der Beobachtungen.
|
||||
:type v: array_like
|
||||
:param Q_vv: Kofaktor-Matrix der Residuen.
|
||||
:type Q_vv: array_like
|
||||
:param ri: Redundanzanteile der Beobachtungen.
|
||||
:type ri: array_like
|
||||
:param labels: Liste der Beobachtungen zur Zuordnung in der Ergebnistabelle.
|
||||
:type labels: list
|
||||
:param s0_apost: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:type s0_apost: float
|
||||
:param alpha: Irrtumswahrscheinlichkeit α (Signifikanzniveau, zweiseitiger Test).
|
||||
:type alpha: float
|
||||
:param beta: Wahrscheinlichkeit β für einen Fehler 2. Art (Testmacht = 1−β).
|
||||
:type beta: float
|
||||
:return: DataFrame mit NVᵢ, Auffälligkeit, Grobfehlerabschätzung GFᵢ und Grenzwert GRZWᵢ je Beobachtung.
|
||||
:rtype: pandas.DataFrame
|
||||
:raises ValueError: Wenn alpha oder beta nicht im Intervall (0, 1) liegen.
|
||||
"""
|
||||
v = np.asarray(v, float).reshape(-1)
|
||||
Q_vv = np.asarray(Q_vv, float)
|
||||
ri = np.asarray(ri, float).reshape(-1)
|
||||
labels = list(labels)
|
||||
|
||||
# Grobfehlerabschätzung:
|
||||
ri_ = np.where(ri == 0, np.nan, ri)
|
||||
GF = -v / ri_
|
||||
|
||||
# Standardabweichungen der Residuen
|
||||
qv = np.diag(Q_vv).astype(float)
|
||||
s_vi = float(s0_apost) * np.sqrt(qv)
|
||||
|
||||
# Normierte Verbesserung NV
|
||||
NV = np.abs(v) / s_vi
|
||||
|
||||
# Quantile k und kA (zweiseitig),
|
||||
k = float(norm.ppf(1 - alpha / 2))
|
||||
kA = float(norm.ppf(1 - beta)) # (Testmacht 1-β)
|
||||
|
||||
# Nichtzentralitätsparameter δ0
|
||||
nzp = k + kA
|
||||
|
||||
# Grenzwert für die Aufdeckbarkeit eines GF (GRZW)
|
||||
GRZW_i = (s_vi / ri_) * nzp
|
||||
|
||||
auffaellig = NV > nzp
|
||||
|
||||
Lokaltest_innere_Zuv = pd.DataFrame({
|
||||
"Beobachtung": labels,
|
||||
"v_i": v,
|
||||
"r_i": ri,
|
||||
"s_vi": s_vi,
|
||||
"k": k,
|
||||
"NV_i": NV,
|
||||
"auffaellig": auffaellig,
|
||||
"GF_i": GF,
|
||||
"GRZW_i": GRZW_i,
|
||||
"alpha": alpha,
|
||||
"beta": beta,
|
||||
"kA": kA,
|
||||
"δ0": nzp,
|
||||
})
|
||||
return Lokaltest_innere_Zuv
|
||||
|
||||
|
||||
def aufruf_lokaltest(liste_beob, alpha, ausgabe_parameterschaetzung, ri, s0_aposteriori):
|
||||
"""Startet den Lokaltest und erzeugt die interaktive Tabelle.
|
||||
|
||||
:param liste_beob: Liste der Beobachtungslabels.
|
||||
:type liste_beob: list
|
||||
:param alpha: Signifikanzniveau.
|
||||
:type alpha: float
|
||||
:param ausgabe_parameterschaetzung: Dictionary mit den Ergebnissen der letzten Iteration der Parameterschätzung.
|
||||
:type ausgabe_parameterschaetzung: dict
|
||||
:param ri: Redundanz.
|
||||
:type ri: Any
|
||||
:param s0_aposteriori: a-posteriori Standardabweichung.
|
||||
:type s0_aposteriori: float
|
||||
:return: ausschalten_dict
|
||||
:rtype: dict
|
||||
"""
|
||||
|
||||
# Initialisieren einer interaktiven Tabelle für die Benutzereingaben
|
||||
itables.init_notebook_mode()
|
||||
labels = [str(s) for s in liste_beob]
|
||||
|
||||
# Benutzereingabe von β
|
||||
beta_input = input("Macht des Tests (1-β) wählen [Standard: 80 % -> 0.80]: ").strip()
|
||||
beta = 0.80 if beta_input == "" else float(beta_input)
|
||||
|
||||
# Berechnungen für den Lokaltest
|
||||
Lokaltest = Zuverlaessigkeit.lokaltest_innere_Zuverlaessigkeit(
|
||||
v=ausgabe_parameterschaetzung["v"],
|
||||
Q_vv=ausgabe_parameterschaetzung["Q_vv"],
|
||||
ri=ri,
|
||||
labels=labels,
|
||||
s0_apost=s0_aposteriori,
|
||||
alpha=alpha,
|
||||
beta=beta
|
||||
)
|
||||
|
||||
if "v_i" in Lokaltest.columns:
|
||||
Lokaltest["v_i"] = Lokaltest["v_i"].round(6)
|
||||
if "r_i" in Lokaltest.columns:
|
||||
Lokaltest["r_i"] = Lokaltest["r_i"].round(4)
|
||||
if "s_vi" in Lokaltest.columns:
|
||||
Lokaltest["s_vi"] = Lokaltest["s_vi"].round(6)
|
||||
if "GF_i" in Lokaltest.columns:
|
||||
Lokaltest["GF_i"] = Lokaltest["GF_i"].round(6)
|
||||
if "GRZW_i" in Lokaltest.columns:
|
||||
Lokaltest["GRZW_i"] = Lokaltest["GRZW_i"].round(6)
|
||||
|
||||
# Anlegen des Dataframes
|
||||
df = Lokaltest.copy()
|
||||
|
||||
if "Beobachtung" not in df.columns:
|
||||
if df.index.name == "Beobachtung":
|
||||
df = df.reset_index()
|
||||
else:
|
||||
df = df.reset_index().rename(columns={"index": "Beobachtung"})
|
||||
|
||||
if "Beobachtung_ausschalten" not in df.columns:
|
||||
df.insert(0, "Beobachtung_ausschalten", "")
|
||||
else:
|
||||
zeile = df.pop("Beobachtung_ausschalten")
|
||||
df.insert(0, "Beobachtung_ausschalten", zeile)
|
||||
|
||||
gui = LokaltestInnereZuverlaessigkeitGUI(df)
|
||||
gui.ausgabe_erstellen()
|
||||
gui.zeige_tabelle()
|
||||
|
||||
Lokaltest.to_excel(r"Zwischenergebnisse\Lokaltest_innere_Zuverlaessugkeit.xlsx", index=False)
|
||||
return gui.ausschalten_dict, beta
|
||||
|
||||
|
||||
|
||||
def aeussere_zuverlaessigkeit(
|
||||
Lokaltest, bezeichnung, Qxx, A, P, s0_apost, unbekannten_liste, x,
|
||||
ausschliessen=("lA_",),
|
||||
):
|
||||
"""
|
||||
Berechnet Parameter der äußeren Zuverlässigkeit (EP/EF) je Beobachtung.
|
||||
|
||||
Auf Basis der Ergebnisse des Lokaltests werden für jede Beobachtung Maße der äußeren
|
||||
Zuverlässigkeit bestimmt. Dazu zählen:
|
||||
|
||||
- Einfluss auf die (relative) Punktlage EP:
|
||||
- aus geschätzter Modellstörung: EP_GF,i = |(1 - r_i) · GF_i|
|
||||
- aus Grenzwert der nicht mehr aufdeckbaren Modellstörung: EP_GRZW,i = |(1 - r_i) · GRZW_i|
|
||||
Für Winkelbeobachtungen (R/ZW) wird EP in eine äquivalente Querabweichung (in m) umgerechnet: q = EP · s
|
||||
wobei EP als Winkelstörung im Bogenmaß (rad) und s als räumliche Strecke zwischen Stand- und
|
||||
Zielpunkt verwendet wird.
|
||||
- Einflussfaktor / Netzverzerrung EF (Worst-Case-Einfluss einer nicht detektierten Störung):
|
||||
Es wird eine Einzelstörung Δl_i = GRZW_i angesetzt (alle anderen Δl_j = 0) und in den
|
||||
Unbekanntenraum übertragen: Δx = Q_xx · A^T · P · Δl
|
||||
Der Einflussfaktor wird lokal (nur für die von der Beobachtung berührten Punktkoordinaten,
|
||||
i.d.R. Stand- und Zielpunkt) über die gewichtete Norm berechnet: EF_i^2 = (Δx_loc^T · Q_loc^{-1} · Δx_loc) / s0^2
|
||||
mit s0 = a posteriori Standardabweichung der Gewichtseinheit.
|
||||
- Punktstreuungsmaß SP_3D und maximale Verfälschung EF·SP:
|
||||
Für die berührten Punkte wird je Punkt der 3×3-Block aus Q_xx betrachtet: als Maß wird die maximale Spur
|
||||
verwendet: SP_3D,loc = s0 · sqrt( max( tr(Q_P) ) )
|
||||
und daraus: (EF·SP)_i = EF_i · SP_3D,loc
|
||||
|
||||
Pseudobeobachtungen (z.B. Lagerungs-/Anschlussgleichungen) können über Präfixe in
|
||||
"ausschliessen" aus der Auswertung entfernt werden. Es wird geprüft, ob die Anzahl
|
||||
der Bezeichnungen und die Zeilenanzahl des Lokaltests zur Beobachtungsanzahl von A passen.
|
||||
|
||||
:param Lokaltest: DataFrame des Lokaltests`.
|
||||
:type Lokaltest: pandas.DataFrame
|
||||
:param bezeichnung: Bezeichnungen der Beobachtungen.
|
||||
:type bezeichnung: list
|
||||
:param Qxx: Kofaktor-Matrix der Unbekannten.
|
||||
:type Qxx: numpy.ndarray
|
||||
:param A: Jacobi-Matrix (A-Matrix).
|
||||
:type A: numpy.ndarray
|
||||
:param P: Gewichtsmatrix der Beobachtungen.
|
||||
:type P: numpy.ndarray
|
||||
:param s0_apost: a-posteriori Standardabweichung der Gewichtseinheit s₀.
|
||||
:type s0_apost: float
|
||||
:param unbekannten_liste: Liste der Unbekannten.
|
||||
:type unbekannten_liste: list
|
||||
:param x: Unbekanntenvektor.
|
||||
:type x: array_like
|
||||
:param ausschliessen: Präfixe von Beobachtungsbezeichnungen, die aus der Auswertung entfernt werden sollen
|
||||
(Standard: ("lA_",) für Lagerungs-/Pseudobeobachtungen).
|
||||
:type ausschliessen: tuple
|
||||
:return: DataFrame mit Stand/Zielpunkt, Redundanzanteil rᵢ, EP (aus GF und GRZW), EF sowie SP_3D und EF·SP_3D.
|
||||
:rtype: pandas.DataFrame
|
||||
:raises ValueError: Wenn die Anzahl der Bezeichnungen oder die Zeilenanzahl des Lokaltests nicht zu A passt.
|
||||
"""
|
||||
|
||||
lokaltest_daten = Lokaltest.copy()
|
||||
bezeichnung = [str(l) for l in list(bezeichnung)]
|
||||
|
||||
Qxx = np.asarray(Qxx, float)
|
||||
A = np.asarray(A, float)
|
||||
P = np.asarray(P, float)
|
||||
x = np.asarray(x, float).reshape(-1)
|
||||
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
# Konsistenzprüfung
|
||||
n = A.shape[0]
|
||||
if len(bezeichnung) != n:
|
||||
raise ValueError(f"len(labels)={len(bezeichnung)} passt nicht zu A.shape[0]={n}.")
|
||||
if len(lokaltest_daten) != n:
|
||||
raise ValueError(f"Lokaltest hat {len(lokaltest_daten)} Zeilen, A hat {n} Beobachtungen.")
|
||||
|
||||
# Pseudobeobachtungen lA rausfiltern
|
||||
beobachtungen = np.ones(n, dtype=bool)
|
||||
if ausschliessen:
|
||||
for i, bez in enumerate(bezeichnung):
|
||||
if any(bez.startswith(pref) for pref in ausschliessen):
|
||||
beobachtungen[i] = False
|
||||
|
||||
lokaltest_daten = lokaltest_daten.loc[beobachtungen].reset_index(drop=True)
|
||||
bezeichnung = [bez for (bez, k) in zip(bezeichnung, beobachtungen) if k]
|
||||
A = A[beobachtungen, :]
|
||||
P = P[np.ix_(beobachtungen, beobachtungen)]
|
||||
n = A.shape[0]
|
||||
|
||||
# Daten aus dem Lokaltest
|
||||
ri = lokaltest_daten["r_i"].astype(float).to_numpy()
|
||||
GF = lokaltest_daten["GF_i"].astype(float).to_numpy()
|
||||
GRZW = lokaltest_daten["GRZW_i"].astype(float).to_numpy()
|
||||
s0 = float(s0_apost)
|
||||
|
||||
# Punktkoordinaten
|
||||
koordinaten = {}
|
||||
punkt_ids = sorted({name[1:] for name in namen_str
|
||||
if name[:1].upper() in ("X", "Y", "Z") and len(name) > 1})
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
ix = namen_str.index(f"X{pid}")
|
||||
iy = namen_str.index(f"Y{pid}")
|
||||
iz = namen_str.index(f"Z{pid}")
|
||||
koordinaten[pid] = (x[ix], x[iy], x[iz])
|
||||
except ValueError:
|
||||
continue
|
||||
|
||||
# Standpunkt/Zielpunkt
|
||||
standpunkte = [""] * n
|
||||
zielpunkte = [""] * n
|
||||
for i, bez in enumerate(bezeichnung):
|
||||
parts = bez.split("_")
|
||||
sp, zp = None, None
|
||||
|
||||
if any(k in bez for k in ["_SD_", "_R_", "_ZW_"]):
|
||||
if len(parts) >= 5:
|
||||
sp, zp = parts[3].strip(), parts[4].strip()
|
||||
elif "gnss" in bez.lower():
|
||||
if len(parts) >= 2:
|
||||
sp, zp = parts[-2].strip(), parts[-1].strip()
|
||||
elif "niv" in bez.lower():
|
||||
if len(parts) >= 4:
|
||||
sp = parts[3].strip()
|
||||
if len(parts) >= 5:
|
||||
zp = parts[4].strip()
|
||||
else:
|
||||
sp = parts[-1].strip()
|
||||
standpunkte[i] = sp or ""
|
||||
zielpunkte[i] = zp or ""
|
||||
|
||||
# Berechnung des EP
|
||||
EP_GF = np.abs((1.0 - ri) * GF)
|
||||
EP_grzw = np.abs((1.0 - ri) * GRZW)
|
||||
|
||||
EP_hat_m = np.full(n, np.nan, float)
|
||||
EP_grzw_m = np.full(n, np.nan, float)
|
||||
|
||||
for i, bez in enumerate(bezeichnung):
|
||||
sp = standpunkte[i]
|
||||
zp = zielpunkte[i]
|
||||
|
||||
wenn_winkel = ("_R_" in bez) or ("_ZW_" in bez)
|
||||
if not wenn_winkel:
|
||||
EP_hat_m[i] = EP_GF[i]
|
||||
EP_grzw_m[i] = EP_grzw[i]
|
||||
continue
|
||||
|
||||
# Wenn Winkel: Querabweichung = Winkel * Strecke (3D)
|
||||
if sp in koordinaten and zp in koordinaten:
|
||||
X1, Y1, Z1 = koordinaten[sp]
|
||||
X2, Y2, Z2 = koordinaten[zp]
|
||||
s = np.sqrt((X2 - X1) ** 2 + (Y2 - Y1) ** 2 + (Z2 - Z1) ** 2)
|
||||
EP_hat_m[i] = EP_GF[i] * s
|
||||
EP_grzw_m[i] = EP_grzw[i] * s
|
||||
|
||||
# Berechnung von EF
|
||||
EF = np.full(n, np.nan, float)
|
||||
SP_m = np.full(n, np.nan, float)
|
||||
EF_SP_m = np.full(n, np.nan, float)
|
||||
|
||||
for i in range(n):
|
||||
sp = standpunkte[i]
|
||||
zp = zielpunkte[i]
|
||||
bloecke = []
|
||||
idx = []
|
||||
try:
|
||||
if sp:
|
||||
b = [
|
||||
namen_str.index(f"X{sp}"),
|
||||
namen_str.index(f"Y{sp}"),
|
||||
namen_str.index(f"Z{sp}")
|
||||
]
|
||||
bloecke.append(b)
|
||||
idx += b
|
||||
if zp:
|
||||
b = [
|
||||
namen_str.index(f"X{zp}"),
|
||||
namen_str.index(f"Y{zp}"),
|
||||
namen_str.index(f"Z{zp}")
|
||||
]
|
||||
bloecke.append(b)
|
||||
idx += b
|
||||
except ValueError:
|
||||
continue
|
||||
if not bloecke:
|
||||
continue
|
||||
idx = list(dict.fromkeys(idx))
|
||||
|
||||
dl = np.zeros((n, 1)) # dl ungestört
|
||||
dl[i, 0] = GRZW[i] # dl gestört durch GRZW
|
||||
dx = Qxx @ (A.T @ (P @ dl)) # dx mit Störung
|
||||
|
||||
dx_pkt = dx[idx, :]
|
||||
Q_pkt = Qxx[np.ix_(idx, idx)]
|
||||
|
||||
# EF
|
||||
EF2 = (dx_pkt.T @ np.linalg.solve(Q_pkt, dx_pkt)).item() / (s0 ** 2)
|
||||
EF[i] = np.sqrt(max(0.0, EF2))
|
||||
|
||||
# SP 3D: Spur der 3x3 Matrizen in einer Liste
|
||||
spur_matrix_liste = [np.trace(Qxx[np.ix_(b, b)]) for b in bloecke]
|
||||
if not spur_matrix_liste:
|
||||
continue
|
||||
|
||||
# SP des schlechtesten Punktes bestimmen
|
||||
sigma_max = s0 * np.sqrt(max(spur_matrix_liste))
|
||||
SP_m[i] = sigma_max
|
||||
|
||||
EF_SP_m[i] = EF[i] * sigma_max
|
||||
|
||||
aeussere_zuverlaessigkeit = pd.DataFrame({
|
||||
"Beobachtung": bezeichnung,
|
||||
"Stand-Pkt": standpunkte,
|
||||
"Ziel-Pkt": zielpunkte,
|
||||
"r_i": ri,
|
||||
"EP_GF [mm]": EP_hat_m * 1000.0,
|
||||
"EP_grzw [mm]": EP_grzw_m * 1000.0,
|
||||
"EF": EF,
|
||||
"SP_3D [mm]": SP_m * 1000.0,
|
||||
"EF*SP_3D [mm]": EF_SP_m * 1000.0,
|
||||
})
|
||||
return aeussere_zuverlaessigkeit
|
||||
|
||||
class LokaltestInnereZuverlaessigkeitGUI:
|
||||
"""Interaktive Auswahloberfläche für den Lokaltest (innere Zuverlässigkeit).
|
||||
|
||||
Die Klasse erzeugt eine ITable-Tabelle auf Basis des Lokaltest-DataFrames und stellt
|
||||
eine Mehrfachauswahl bereit. Für GNSS-Basislinien wird sichergestellt, dass bei Auswahl
|
||||
einer Komponente (bx/by/bz) automatisch das gesamte Trio gewählt bzw. abgewählt wird.
|
||||
"""
|
||||
|
||||
def __init__(self, df):
|
||||
"""Initialisiert die GUI-Objekte.
|
||||
|
||||
:param df: DataFrame des Lokaltests (inkl. Spalte "Beobachtung").
|
||||
:type df: pandas.DataFrame
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
self.df = df
|
||||
try:
|
||||
if not (self.df.index.equals(pd.RangeIndex(start=0, stop=len(self.df), step=1))):
|
||||
self.df = self.df.reset_index(drop=True)
|
||||
except:
|
||||
self.df = self.df.reset_index(drop=True)
|
||||
|
||||
self.tabelle = None
|
||||
|
||||
self.dict_gnss = {}
|
||||
self.dict_gnss_erweitert = {}
|
||||
|
||||
self.auswahl_zeilen_vorher = set()
|
||||
self.update_durch_code = False
|
||||
|
||||
self.ausschalten_dict = {}
|
||||
|
||||
self.output = widgets.Output()
|
||||
self.btn_auswahl_speichern = widgets.Button(description="Auswahl speichern", icon="download")
|
||||
self.btn_auswahl_zuruecksetzen = widgets.Button(description="Rückgängig", icon="refresh")
|
||||
|
||||
@staticmethod
|
||||
def gnss_komponenten_extrahieren(beobachtung: str):
|
||||
"""Extrahiert GNSS-Komponente und einen eindeutigen Key für bx/by/bz-Trio.
|
||||
|
||||
:param beobachtung: Text aus Spalte "Beobachtung".
|
||||
:type beobachtung: str
|
||||
:return: (komponente, key) oder (None, None)
|
||||
:rtype: tuple[str | None, str | None]
|
||||
"""
|
||||
beobachtung = str(beobachtung).strip()
|
||||
for gnss_komponente in ["bx", "by", "bz"]:
|
||||
bezeichnung = f"_gnss{gnss_komponente}_"
|
||||
if bezeichnung in beobachtung:
|
||||
key = beobachtung.replace(bezeichnung, "_gnss_")
|
||||
return gnss_komponente, key
|
||||
return None, None
|
||||
|
||||
def gnss_dictionary_erstellen(self) -> None:
|
||||
"""Exportiert die Tabelleneinträge in ein Dictionary auf Basis der Tabellenzeilen.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
liste_beobachtungen = self.tabelle.df["Beobachtung"].astype(str).tolist()
|
||||
# Als Instanzvariable speichern
|
||||
self.dict_gnss = {}
|
||||
|
||||
for i, beobachtung in enumerate(liste_beobachtungen):
|
||||
beobachtung = str(beobachtung).strip()
|
||||
|
||||
if "_gnssbx_" in beobachtung:
|
||||
key = beobachtung.split("_gnssbx_", 1)[1].strip()
|
||||
if key not in self.dict_gnss:
|
||||
self.dict_gnss[key] = {}
|
||||
self.dict_gnss[key]["bx"] = i
|
||||
|
||||
if "_gnssby_" in beobachtung:
|
||||
key = beobachtung.split("_gnssby_", 1)[1].strip()
|
||||
if key not in self.dict_gnss:
|
||||
self.dict_gnss[key] = {}
|
||||
self.dict_gnss[key]["by"] = i
|
||||
|
||||
if "_gnssbz_" in beobachtung:
|
||||
key = beobachtung.split("_gnssbz_", 1)[1].strip()
|
||||
if key not in self.dict_gnss:
|
||||
self.dict_gnss[key] = {}
|
||||
self.dict_gnss[key]["bz"] = i
|
||||
|
||||
def gnss_dictionary_erweitert_erstellen(self) -> None:
|
||||
"""Baut ein Dictionary mit alles GNSS-Komponten auf.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
self.dict_gnss_erweitert = {}
|
||||
|
||||
for idx, row in self.df.iterrows():
|
||||
value, key = self.gnss_komponenten_extrahieren(row["Beobachtung"])
|
||||
if key:
|
||||
if key not in self.dict_gnss_erweitert:
|
||||
self.dict_gnss_erweitert[key] = []
|
||||
self.dict_gnss_erweitert[key].append(idx)
|
||||
|
||||
def export_ausschalten_dict(self, Eintrag_Auswahl: str = "beobachtung_ausschalten", Wert_nicht_ausgewaehlt: str = "") -> dict:
|
||||
"""Exportiert die aktuelle Auswahl in ein Dictionary.
|
||||
|
||||
:param Eintrag_Auswahl: Wert für ausgewählte Beobachtungen (beobachtung_ausschalten).
|
||||
:type Eintrag_Auswahl: str
|
||||
:param Wert_nicht_ausgewaehlt: Wert für nicht ausgewählte Beobachtungen ("").
|
||||
:type Wert_nicht_ausgewaehlt: str
|
||||
:return: Dict {Beobachtung: "beobachtung_ausschalten" oder ""}
|
||||
:rtype: dict
|
||||
"""
|
||||
auswahl = set(self.tabelle.selected_rows or [])
|
||||
liste_beobachtungen = self.tabelle.df["Beobachtung"].astype(str).tolist()
|
||||
|
||||
# Zurückgabe des Ergebnisdictionarys für die Weiterverarbeitung
|
||||
return {
|
||||
liste_beobachtungen[i]: (Eintrag_Auswahl if i in auswahl else Wert_nicht_ausgewaehlt)
|
||||
for i in range(len(liste_beobachtungen))
|
||||
}
|
||||
|
||||
def aktualisiere_ausschalten_dict(self) -> None:
|
||||
"""Aktualisiert das ausschalten_dict in der Instanzvariablen.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
neu = self.export_ausschalten_dict()
|
||||
self.ausschalten_dict.clear()
|
||||
self.ausschalten_dict.update(neu)
|
||||
|
||||
def ausgabe_aktualisieren(self) -> None:
|
||||
"""Aktualisiert Ausgabe und schreibt self.ausschalten_dict neu.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
self.aktualisiere_ausschalten_dict()
|
||||
|
||||
with self.output:
|
||||
clear_output(wait=True)
|
||||
auswahl = self.tabelle.selected_rows or []
|
||||
print(f"AUSGESCHALTET: {len(auswahl)}")
|
||||
|
||||
if len(auswahl) > 0:
|
||||
zeilen = [c for c in ["Beobachtung", "v_i", "r_i", "auffaellig", "GF_i", "GRZW_i"] if c in self.tabelle.df.columns]
|
||||
display(self.tabelle.df.iloc[auswahl][zeilen].head(30))
|
||||
|
||||
def auswahl_exportieren(self, _=None) -> None:
|
||||
"""Button-Callback: Ausgabe der ausgeschalteten Beobachtungen.
|
||||
|
||||
:param _: Button-Event
|
||||
:type _: Any
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
self.aktualisiere_ausschalten_dict()
|
||||
|
||||
with self.output:
|
||||
clear_output(wait=True)
|
||||
print("ausschalten_dict ist aktualisiert.")
|
||||
ausgeschaltet = [k for k, v in self.ausschalten_dict.items() if v == "X"]
|
||||
print(f"Nur ausgeschaltete Beobachtungen ({len(ausgeschaltet)}):")
|
||||
display(ausgeschaltet[:300])
|
||||
|
||||
def auswahl_zuruecksetzen(self, _=None) -> None:
|
||||
"""Button-aktion: setzt Auswahl zurück.
|
||||
|
||||
:param _: Button-Event
|
||||
:type _: Any
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
self.tabelle.selected_rows = []
|
||||
self.ausgabe_aktualisieren()
|
||||
|
||||
def gnss_auswahl_synchronisieren(self, aenderungen: dict) -> None:
|
||||
"""Synchronisiert die GNSS-bx/by/bz Auswahl.
|
||||
|
||||
:param aenderungen: Dictionary mit Änderungen.
|
||||
:type aenderungen: dict
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
if self.update_durch_code:
|
||||
return
|
||||
|
||||
auswahl_aktuell = set(aenderungen["new"] or [])
|
||||
|
||||
hinzufuegen = auswahl_aktuell - self.auswahl_zeilen_vorher
|
||||
entfernt = self.auswahl_zeilen_vorher - auswahl_aktuell
|
||||
auswahl_final = set(auswahl_aktuell)
|
||||
|
||||
# Hinzufügen -> Alle Komponten auswählen
|
||||
for index in hinzufuegen:
|
||||
if index not in self.df.index:
|
||||
continue
|
||||
|
||||
beob_name = str(self.df.loc[index, "Beobachtung"])
|
||||
value, key = self.gnss_komponenten_extrahieren(beob_name)
|
||||
|
||||
if key in self.dict_gnss_erweitert:
|
||||
for p_idx in self.dict_gnss_erweitert[key]:
|
||||
auswahl_final.add(p_idx)
|
||||
|
||||
# Entfernen -> alle Komponenten abwählen
|
||||
for index in entfernt:
|
||||
if index not in self.df.index:
|
||||
continue
|
||||
|
||||
beob_name = str(self.df.loc[index, "Beobachtung"])
|
||||
value, key = self.gnss_komponenten_extrahieren(beob_name)
|
||||
|
||||
if key in self.dict_gnss_erweitert:
|
||||
for p_idx in self.dict_gnss_erweitert[key]:
|
||||
if p_idx in auswahl_final:
|
||||
auswahl_final.remove(p_idx)
|
||||
|
||||
# Nur bei Änderungen zurück in die Ausgabe schreiben
|
||||
if auswahl_final != auswahl_aktuell:
|
||||
self.update_durch_code = True
|
||||
self.tabelle.selected_rows = sorted(list(auswahl_final))
|
||||
self.update_durch_code = False
|
||||
|
||||
self.auswahl_zeilen_vorher = set(self.tabelle.selected_rows or [])
|
||||
self.ausgabe_aktualisieren()
|
||||
|
||||
|
||||
def ausgabe_erstellen(self) -> None:
|
||||
"""Erstellt die Tabelle und verbindet die Buttons.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
self.tabelle = ITable(
|
||||
self.df,
|
||||
maxBytes=5 * 1024 * 1024, # 5 MB
|
||||
columnDefs=[
|
||||
{"targets": 0, "orderable": False, "className": "select-checkbox", "width": "26px"},
|
||||
],
|
||||
select={"style": "multi", "selector": "td:first-child"},
|
||||
order=[[1, "asc"]],
|
||||
)
|
||||
|
||||
self.gnss_dictionary_erstellen()
|
||||
self.gnss_dictionary_erweitert_erstellen()
|
||||
|
||||
self.auswahl_zeilen_vorher = set(self.tabelle.selected_rows or [])
|
||||
self.aktualisiere_ausschalten_dict()
|
||||
|
||||
self.btn_auswahl_speichern.on_click(self.auswahl_exportieren)
|
||||
self.btn_auswahl_zuruecksetzen.on_click(self.auswahl_zuruecksetzen)
|
||||
|
||||
self.tabelle.observe(self.gnss_auswahl_synchronisieren, names="selected_rows")
|
||||
|
||||
def zeige_tabelle(self) -> None:
|
||||
"""Zeigt die Tabelle an.
|
||||
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
display(widgets.VBox([self.tabelle, widgets.HBox([self.btn_auswahl_speichern, self.btn_auswahl_zuruecksetzen]), self.output]))
|
||||
self.ausgabe_aktualisieren()
|
||||
@@ -1,321 +0,0 @@
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from scipy.stats import f
|
||||
import pandas as pd
|
||||
import Berechnungen
|
||||
|
||||
|
||||
class Genauigkeitsmaße:
|
||||
|
||||
|
||||
@staticmethod
|
||||
def berechne_s0apost(v: np.ndarray, P: np.ndarray, r: int) -> float:
|
||||
vTPv_matrix = v.T @ P @ v
|
||||
vTPv = float(vTPv_matrix.item())
|
||||
s0apost = np.sqrt(vTPv / r)
|
||||
return float(s0apost)
|
||||
|
||||
|
||||
@staticmethod
|
||||
def helmert_punktfehler(Qxx, s0_apost, unbekannten_liste, dim=3):
|
||||
diagQ = np.diag(Qxx)
|
||||
daten = []
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
punkt_ids = []
|
||||
for n in namen_str:
|
||||
if n.upper().startswith('X'):
|
||||
punkt_ids.append(n[1:])
|
||||
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
|
||||
idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
|
||||
|
||||
qx = diagQ[idx_x]
|
||||
qy = diagQ[idx_y]
|
||||
qz = 0.0
|
||||
|
||||
if dim == 3:
|
||||
try:
|
||||
idx_z = next(i for i, n in enumerate(namen_str) if n.upper() == f"Z{pid}".upper())
|
||||
qz = diagQ[idx_z]
|
||||
except StopIteration:
|
||||
qz = 0.0
|
||||
|
||||
sx = s0_apost * np.sqrt(qx)
|
||||
sy = s0_apost * np.sqrt(qy)
|
||||
sz = s0_apost * np.sqrt(qz) if dim == 3 else 0
|
||||
|
||||
sP = s0_apost * np.sqrt(qx + qy + qz)
|
||||
|
||||
daten.append([pid, float(sx), float(sy), float(sz), float(sP)])
|
||||
except:
|
||||
continue
|
||||
|
||||
helmert_punktfehler = pd.DataFrame(daten, columns=["Punkt", "σx", "σy", "σz", f"σP_{dim}D"])
|
||||
return helmert_punktfehler
|
||||
|
||||
|
||||
|
||||
@staticmethod
|
||||
def standardellipse(Qxx, s0_apost, unbekannten_liste):
|
||||
Qxx = np.asarray(Qxx, float)
|
||||
daten = []
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
punkt_ids = []
|
||||
for n in namen_str:
|
||||
if n.upper().startswith('X'):
|
||||
punkt_ids.append(n[1:])
|
||||
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
|
||||
idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
|
||||
|
||||
qxx = Qxx[idx_x, idx_x]
|
||||
qyy = Qxx[idx_y, idx_y]
|
||||
qyx = Qxx[idx_y, idx_x]
|
||||
|
||||
# Standardabweichungen
|
||||
sx = s0_apost * np.sqrt(qxx)
|
||||
sy = s0_apost * np.sqrt(qyy)
|
||||
sxy = (s0_apost ** 2) * qyx
|
||||
|
||||
k = np.sqrt((qxx - qyy) ** 2 + 4 * (qyx ** 2))
|
||||
|
||||
# Q_dmax/min = 0.5 * (Qyy + Qxx +/- k)
|
||||
q_dmax = 0.5 * (qyy + qxx + k)
|
||||
q_dmin = 0.5 * (qyy + qxx - k)
|
||||
|
||||
# Halbachsen
|
||||
s_max = s0_apost * np.sqrt(q_dmax)
|
||||
s_min = s0_apost * np.sqrt(q_dmin)
|
||||
|
||||
# Richtungswinkel theta in gon:
|
||||
zaehler = 2 * qyx
|
||||
nenner = qxx - qyy
|
||||
t_grund = 0.5 * np.arctan(abs(zaehler) / abs(nenner)) * (200 / np.pi)
|
||||
|
||||
# Quadrantenabfrage
|
||||
if nenner > 0 and qyx > 0: # Qxx - Qyy > 0 und Qyx > 0
|
||||
t_gon = t_grund # 0 - 50 gon
|
||||
elif nenner < 0 and qyx > 0: # Qxx - Qyy < 0 und Qyx > 0
|
||||
t_gon = 100 - t_grund # 50 - 100 gon
|
||||
elif nenner < 0 and qyx < 0: # Qxx - Qyy < 0 und Qyx < 0
|
||||
t_gon = 100 + t_grund # 100 - 150 gon
|
||||
elif nenner > 0 and qyx < 0: # Qxx - Qyy > 0 und Qyx < 0
|
||||
t_gon = 200 - t_grund # 150 - 200 gon
|
||||
else:
|
||||
t_gon = 0.0
|
||||
|
||||
daten.append([
|
||||
pid,
|
||||
float(sx), float(sy), float(sxy),
|
||||
float(s_max), float(s_min),
|
||||
float(t_gon)
|
||||
])
|
||||
except:
|
||||
continue
|
||||
standardellipse = pd.DataFrame(daten, columns=["Punkt", "σx [m]", "σy [m]", "σxy [m]", "Große Halbachse [m]", "Kleine Halbachse [m]", "θ [gon]"])
|
||||
return standardellipse
|
||||
|
||||
|
||||
|
||||
@staticmethod
|
||||
def konfidenzellipse(Qxx, s0_apost, unbekannten_liste, R, alpha):
|
||||
Qxx = np.asarray(Qxx, float)
|
||||
daten = []
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
punkt_ids = [n[1:] for n in namen_str if n.upper().startswith('X')]
|
||||
|
||||
# Faktor für Konfidenzellipse (F-Verteilung)
|
||||
kk = float(np.sqrt(2.0 * f.ppf(1.0 - alpha, 2, R)))
|
||||
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
|
||||
idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
|
||||
|
||||
qxx = Qxx[idx_x, idx_x]
|
||||
qyy = Qxx[idx_y, idx_y]
|
||||
qyx = Qxx[idx_y, idx_x]
|
||||
|
||||
# Standardabweichungen
|
||||
sx = s0_apost * np.sqrt(qxx)
|
||||
sy = s0_apost * np.sqrt(qyy)
|
||||
sxy = (s0_apost ** 2) * qyx
|
||||
|
||||
k = np.sqrt((qxx - qyy) ** 2 + 4 * (qyx ** 2))
|
||||
|
||||
# Q_dmax/min = 0.5 * (Qyy + Qxx +/- k)
|
||||
q_dmax = 0.5 * (qyy + qxx + k)
|
||||
q_dmin = 0.5 * (qyy + qxx - k)
|
||||
|
||||
# Halbachsen der Standardellipse
|
||||
s_max = s0_apost * np.sqrt(q_dmax)
|
||||
s_min = s0_apost * np.sqrt(q_dmin)
|
||||
|
||||
# Halbachsen der Konfidenzellipse
|
||||
A_K = kk * s_max
|
||||
B_K = kk * s_min
|
||||
|
||||
# Richtungswinkel theta in gon:
|
||||
zaehler = 2 * qyx
|
||||
nenner = qxx - qyy
|
||||
t_grund = 0.5 * np.arctan(abs(zaehler) / abs(nenner)) * (200 / np.pi)
|
||||
|
||||
# Quadrantenabfrage
|
||||
if nenner > 0 and qyx > 0:
|
||||
t_gon = t_grund # 0 - 50 gon
|
||||
elif nenner < 0 and qyx > 0:
|
||||
t_gon = 100 - t_grund # 50 - 100 gon
|
||||
elif nenner < 0 and qyx < 0:
|
||||
t_gon = 100 + t_grund # 100 - 150 gon
|
||||
elif nenner > 0 and qyx < 0:
|
||||
t_gon = 200 - t_grund # 150 - 200 gon
|
||||
else:
|
||||
t_gon = 0.0
|
||||
|
||||
daten.append([
|
||||
pid,
|
||||
float(sx), float(sy), float(sxy),
|
||||
float(A_K), float(B_K),
|
||||
float(t_gon)
|
||||
])
|
||||
|
||||
except:
|
||||
continue
|
||||
|
||||
konfidenzellipse = pd.DataFrame(daten, columns=["Punkt", "σx [m]", "σy [m]", "σxy [m]", "Große Halbachse [m]",
|
||||
"Kleine Halbachse [m]", "θ [gon]"])
|
||||
return konfidenzellipse
|
||||
|
||||
|
||||
|
||||
class Plot:
|
||||
|
||||
|
||||
@staticmethod
|
||||
def netzplot_ellipsen(
|
||||
Koord_ENU,
|
||||
unbekannten_labels,
|
||||
beobachtungs_labels,
|
||||
df_konf_ellipsen_enu,
|
||||
v_faktor=1000,
|
||||
n_ellipse_pts=60,
|
||||
title="Netzplot im ENU-System mit Konfidenzellipsen"
|
||||
):
|
||||
names = [str(s).strip() for s in unbekannten_labels]
|
||||
|
||||
if "θ_EN [gon]" in df_konf_ellipsen_enu.columns:
|
||||
theta_col = "θ_EN [gon]"
|
||||
elif "θ [gon]" in df_konf_ellipsen_enu.columns:
|
||||
theta_col = "θ [gon]"
|
||||
else:
|
||||
raise ValueError("Spalte 'θ_EN [gon]' oder 'θ [gon]' fehlt im DataFrame.")
|
||||
|
||||
punkt_ids = sorted({nm[1:] for nm in names if nm and nm[0].upper() in ("X", "Y", "Z")})
|
||||
|
||||
fig = go.Figure()
|
||||
|
||||
# 1) Darstellungen der Beobachtungen
|
||||
beob_typen = {
|
||||
'GNSS-Basislinien': {'pattern': 'gnss', 'color': 'rgba(255, 100, 0, 0.4)'},
|
||||
'Tachymeter-Beob': {'pattern': '', 'color': 'rgba(100, 100, 100, 0.3)'}
|
||||
}
|
||||
|
||||
for typ, info in beob_typen.items():
|
||||
x_l, y_l = [], []
|
||||
for bl in beobachtungs_labels:
|
||||
bl_str = str(bl).lower()
|
||||
is_typ = ((info['pattern'] in bl_str and info['pattern'] != '') or
|
||||
(info['pattern'] == '' and 'gnss' not in bl_str and 'niv' not in bl_str))
|
||||
if not is_typ:
|
||||
continue
|
||||
|
||||
bl_raw = str(bl)
|
||||
pts = []
|
||||
for pid in punkt_ids:
|
||||
if (f"_{pid}" in bl_raw) or bl_raw.startswith(f"{pid}_"):
|
||||
if pid in Koord_ENU:
|
||||
pts.append(pid)
|
||||
|
||||
if len(pts) >= 2:
|
||||
p1, p2 = pts[0], pts[1]
|
||||
x_l.extend([Koord_ENU[p1][0], Koord_ENU[p2][0], None]) # E
|
||||
y_l.extend([Koord_ENU[p1][1], Koord_ENU[p2][1], None]) # N
|
||||
|
||||
if x_l:
|
||||
fig.add_trace(go.Scatter(x=x_l, y=y_l, mode='lines', name=typ,
|
||||
line=dict(color=info['color'], width=1)))
|
||||
|
||||
# 2) Darstellung der Konfidenzellipsen
|
||||
t = np.linspace(0, 2 * np.pi, n_ellipse_pts)
|
||||
first = True
|
||||
for _, row in df_konf_ellipsen_enu.iterrows():
|
||||
pid = str(row["Punkt"])
|
||||
if pid not in Koord_ENU:
|
||||
continue
|
||||
|
||||
a = float(row["a_K"]) * v_faktor
|
||||
b = float(row["b_K"]) * v_faktor
|
||||
theta = float(row[theta_col]) * np.pi / 200.0 # gon->rad
|
||||
|
||||
ex = a * np.cos(t)
|
||||
ey = b * np.sin(t)
|
||||
|
||||
c, s = np.cos(theta), np.sin(theta)
|
||||
xr = c * ex - s * ey
|
||||
yr = s * ex + c * ey
|
||||
|
||||
E0, N0, _ = Koord_ENU[pid]
|
||||
|
||||
fig.add_trace(go.Scatter(
|
||||
x=E0 + xr, y=N0 + yr,
|
||||
mode="lines",
|
||||
line=dict(color="red", width=1.5),
|
||||
name=f"Ellipsen (×{v_faktor})",
|
||||
legendgroup="Ellipsen",
|
||||
showlegend=first,
|
||||
hoverinfo="skip"
|
||||
))
|
||||
first = False
|
||||
|
||||
# 3) Darstellung der Punkte
|
||||
xs, ys, texts, hovers = [], [], [], []
|
||||
for pid in punkt_ids:
|
||||
if pid not in Koord_ENU:
|
||||
continue
|
||||
E, N, U = Koord_ENU[pid]
|
||||
xs.append(E);
|
||||
ys.append(N);
|
||||
texts.append(pid)
|
||||
hovers.append(f"Punkt {pid}<br>E={E:.4f} m<br>N={N:.4f} m<br>U={U:.4f} m")
|
||||
|
||||
fig.add_trace(go.Scatter(
|
||||
x=xs, y=ys, mode="markers+text",
|
||||
text=texts, textposition="top center",
|
||||
marker=dict(size=8, color="black"),
|
||||
name="Netzpunkte",
|
||||
hovertext=hovers, hoverinfo="text"
|
||||
))
|
||||
|
||||
fig.update_layout(
|
||||
title=f"{title} (Ellipsen ×{v_faktor})",
|
||||
xaxis=dict(title="E [m]", scaleanchor="y", scaleratio=1, showgrid=True, gridcolor="lightgrey"),
|
||||
yaxis=dict(title="N [m]", showgrid=True, gridcolor="lightgrey"),
|
||||
width=1100, height=900,
|
||||
template="plotly_white",
|
||||
plot_bgcolor="white"
|
||||
)
|
||||
|
||||
fig.add_annotation(
|
||||
text=f"<b>Maßstab Ellipsen:</b><br>Dargestellte Größe = Konfidenzellipse × {v_faktor}",
|
||||
align='left', showarrow=False, xref='paper', yref='paper', x=0.02, y=0.05,
|
||||
bgcolor="white", bordercolor="black", borderwidth=1
|
||||
)
|
||||
|
||||
fig.show(config={'scrollZoom': True})
|
||||
@@ -1,309 +0,0 @@
|
||||
from dataclasses import dataclass
|
||||
import numpy as np
|
||||
from scipy import stats
|
||||
from scipy.stats import norm
|
||||
import pandas as pd
|
||||
|
||||
@dataclass
|
||||
class Zuverlaessigkeit:
|
||||
|
||||
@staticmethod
|
||||
def gesamtredundanz(n, u):
|
||||
r = n - u
|
||||
return r
|
||||
|
||||
|
||||
@staticmethod
|
||||
def berechne_R(Q_vv, P):
|
||||
R = Q_vv @ P
|
||||
return R #Redundanzmatrix
|
||||
|
||||
|
||||
@staticmethod
|
||||
def berechne_ri(R):
|
||||
ri = np.diag(R)
|
||||
EVi = 100.0 * ri
|
||||
return ri, EVi #Redundanzanteile
|
||||
|
||||
|
||||
@staticmethod
|
||||
def klassifiziere_ri(ri): #Klassifizierung der Redundanzanteile
|
||||
if ri < 0.01:
|
||||
return "nicht kontrollierbar"
|
||||
elif ri < 0.10:
|
||||
return "schlecht kontrollierbar"
|
||||
elif ri < 0.30:
|
||||
return "ausreichend kontrollierbar"
|
||||
elif ri < 0.70:
|
||||
return "gut kontrollierbar"
|
||||
else:
|
||||
return "nahezu vollständig redundant"
|
||||
|
||||
@staticmethod
|
||||
def globaltest(r_gesamt, sigma0_apost, sigma0_apriori, alpha):
|
||||
T_G = (sigma0_apost ** 2) / (sigma0_apriori ** 2)
|
||||
F_krit = stats.f.ppf(1 - alpha, r_gesamt, 10 ** 9)
|
||||
H0 = T_G < F_krit
|
||||
|
||||
if H0:
|
||||
interpretation = (
|
||||
"Nullhypothese H₀ angenommen.\n"
|
||||
)
|
||||
else:
|
||||
interpretation = (
|
||||
"Nullhypothese H₀ verworfen!\n"
|
||||
"Dies kann folgende Gründe haben:\n"
|
||||
"→ Es befinden sich grobe Fehler im Datenmaterial. Bitte Lokaltest durchführen und ggf. grobe Fehler im Datenmaterial entfernen.\n"
|
||||
"→ Das stochastische Modell ist zu optimistisch. Bitte Gewichte überprüfen und ggf. anpassen."
|
||||
)
|
||||
|
||||
return {
|
||||
"r_gesamt": r_gesamt,
|
||||
"sigma0_apost": sigma0_apost,
|
||||
"sigma0_apriori": sigma0_apriori,
|
||||
"alpha": alpha,
|
||||
"T_G": T_G,
|
||||
"F_krit": F_krit,
|
||||
"H0_angenommen": H0,
|
||||
"Interpretation": interpretation,
|
||||
}
|
||||
|
||||
|
||||
def lokaltest_innere_Zuverlaessigkeit(v, Q_vv, ri, labels, s0_apost, alpha, beta):
|
||||
v = np.asarray(v, float).reshape(-1)
|
||||
Q_vv = np.asarray(Q_vv, float)
|
||||
ri = np.asarray(ri, float).reshape(-1)
|
||||
labels = list(labels)
|
||||
|
||||
# Grobfehlerabschätzung:
|
||||
ri_ = np.where(ri == 0, np.nan, ri)
|
||||
GF = -v / ri_
|
||||
|
||||
# Standardabweichungen der Residuen
|
||||
qv = np.diag(Q_vv).astype(float)
|
||||
s_vi = float(s0_apost) * np.sqrt(qv)
|
||||
|
||||
# Normierte Verbesserung NV
|
||||
NV = np.abs(v) / s_vi
|
||||
|
||||
# Quantile k und kA (zweiseitig),
|
||||
k = float(norm.ppf(1 - alpha / 2))
|
||||
kA = float(norm.ppf(1 - beta)) # (Testmacht 1-β)
|
||||
|
||||
# Nichtzentralitätsparameter δ0
|
||||
nzp = k + kA
|
||||
|
||||
# Grenzwert für die Aufdeckbarkeit eines GF (GRZW)
|
||||
GRZW_i = (s_vi / ri_) * nzp
|
||||
|
||||
auffaellig = NV > nzp
|
||||
|
||||
Lokaltest_innere_Zuv = pd.DataFrame({
|
||||
"Beobachtung": labels,
|
||||
"v_i": v,
|
||||
"r_i": ri,
|
||||
"s_vi": s_vi,
|
||||
"k": k,
|
||||
"NV_i": NV,
|
||||
"auffaellig": auffaellig,
|
||||
"GF_i": GF,
|
||||
"GRZW_i": GRZW_i,
|
||||
"alpha": alpha,
|
||||
"beta": beta,
|
||||
"kA": kA,
|
||||
"δ0": nzp,
|
||||
})
|
||||
return Lokaltest_innere_Zuv
|
||||
|
||||
|
||||
|
||||
def aeussere_zuverlaessigkeit(
|
||||
Lokaltest, labels, Qxx, A, P, s0_apost, unbekannten_liste, x,
|
||||
angle_units="rad",
|
||||
ep_use_abs=True,
|
||||
exclude_prefixes=("lA_",),
|
||||
):
|
||||
df = Lokaltest.copy()
|
||||
labels = [str(l) for l in list(labels)]
|
||||
|
||||
Qxx = np.asarray(Qxx, float)
|
||||
A = np.asarray(A, float)
|
||||
P = np.asarray(P, float)
|
||||
x = np.asarray(x, float).reshape(-1)
|
||||
|
||||
namen_str = [str(sym) for sym in unbekannten_liste]
|
||||
|
||||
n = A.shape[0]
|
||||
if len(labels) != n:
|
||||
raise ValueError(f"len(labels)={len(labels)} passt nicht zu A.shape[0]={n}.")
|
||||
if len(df) != n:
|
||||
raise ValueError(f"Lokaltest hat {len(df)} Zeilen, A hat {n} Beobachtungen.")
|
||||
|
||||
# Pseudobeobachtungen rausfiltern
|
||||
keep = np.ones(n, dtype=bool)
|
||||
if exclude_prefixes:
|
||||
for i, lbl in enumerate(labels):
|
||||
if any(lbl.startswith(pref) for pref in exclude_prefixes):
|
||||
keep[i] = False
|
||||
|
||||
# alles konsistent kürzen (wichtig: auch A & P!)
|
||||
df = df.loc[keep].reset_index(drop=True)
|
||||
labels = [lbl for (lbl, k) in zip(labels, keep) if k]
|
||||
A = A[keep, :]
|
||||
P = P[np.ix_(keep, keep)]
|
||||
|
||||
# neue n
|
||||
n = A.shape[0]
|
||||
|
||||
# Daten aus dem Lokaltest
|
||||
ri = df["r_i"].astype(float).to_numpy()
|
||||
GF = df["GF_i"].astype(float).to_numpy()
|
||||
GRZW = df["GRZW_i"].astype(float).to_numpy()
|
||||
|
||||
s0 = float(s0_apost)
|
||||
|
||||
def to_rad(val):
|
||||
if angle_units == "rad":
|
||||
return val
|
||||
if angle_units == "gon":
|
||||
return val * (np.pi / 200.0)
|
||||
if angle_units == "deg":
|
||||
return val * (np.pi / 180.0)
|
||||
raise ValueError("angle_units muss 'rad', 'gon' oder 'deg' sein.")
|
||||
|
||||
# Punktkoordinaten aus x (für Streckenäquivalent bei Winkel-EP)
|
||||
coords = {}
|
||||
punkt_ids = sorted({name[1:] for name in namen_str
|
||||
if name[:1].upper() in ("X", "Y", "Z") and len(name) > 1})
|
||||
for pid in punkt_ids:
|
||||
try:
|
||||
ix = namen_str.index(f"X{pid}")
|
||||
iy = namen_str.index(f"Y{pid}")
|
||||
iz = namen_str.index(f"Z{pid}")
|
||||
coords[pid] = (x[ix], x[iy], x[iz])
|
||||
except ValueError:
|
||||
continue
|
||||
|
||||
# Standpunkt/Zielpunkt
|
||||
standpunkte = [""] * n
|
||||
zielpunkte = [""] * n
|
||||
for i, lbl in enumerate(labels):
|
||||
parts = lbl.split("_")
|
||||
sp, zp = None, None
|
||||
|
||||
if any(k in lbl for k in ["_SD_", "_R_", "_ZW_"]):
|
||||
if len(parts) >= 5:
|
||||
sp, zp = parts[3].strip(), parts[4].strip()
|
||||
elif "gnss" in lbl.lower():
|
||||
if len(parts) >= 2:
|
||||
sp, zp = parts[-2].strip(), parts[-1].strip()
|
||||
elif "niv" in lbl.lower():
|
||||
if len(parts) >= 4:
|
||||
sp = parts[3].strip()
|
||||
if len(parts) >= 5:
|
||||
zp = parts[4].strip()
|
||||
else:
|
||||
sp = parts[-1].strip()
|
||||
|
||||
standpunkte[i] = sp or ""
|
||||
zielpunkte[i] = zp or ""
|
||||
|
||||
# Berechnung des EPs
|
||||
EP_GF = (1.0 - ri) * GF
|
||||
EP_grzw = (1.0 - ri) * GRZW
|
||||
if ep_use_abs:
|
||||
EP_GF = np.abs(EP_GF)
|
||||
EP_grzw = np.abs(EP_grzw)
|
||||
|
||||
EP_hat_m = np.full(n, np.nan, float)
|
||||
EP_grzw_m = np.full(n, np.nan, float)
|
||||
|
||||
for i, lbl in enumerate(labels):
|
||||
sp = standpunkte[i]
|
||||
zp = zielpunkte[i]
|
||||
|
||||
is_angle = ("_R_" in lbl) or ("_ZW_" in lbl)
|
||||
if not is_angle:
|
||||
EP_hat_m[i] = EP_GF[i]
|
||||
EP_grzw_m[i] = EP_grzw[i]
|
||||
continue
|
||||
|
||||
# Winkel -> Querabweichung = Winkel(rad) * Strecke (3D)
|
||||
if sp in coords and zp in coords:
|
||||
X1, Y1, Z1 = coords[sp]
|
||||
X2, Y2, Z2 = coords[zp]
|
||||
s = np.sqrt((X2 - X1) ** 2 + (Y2 - Y1) ** 2 + (Z2 - Z1) ** 2)
|
||||
|
||||
EP_hat_m[i] = to_rad(EP_GF[i]) * s
|
||||
EP_grzw_m[i] = to_rad(EP_grzw[i]) * s
|
||||
|
||||
# 3x3 Blöcke
|
||||
def idx_xyz(pid):
|
||||
return [
|
||||
namen_str.index(f"X{pid}"),
|
||||
namen_str.index(f"Y{pid}"),
|
||||
namen_str.index(f"Z{pid}")
|
||||
]
|
||||
|
||||
# EF lokal + SP lokal (3D)
|
||||
EF = np.full(n, np.nan, float)
|
||||
SP_loc_m = np.full(n, np.nan, float)
|
||||
EFSP_loc_m = np.full(n, np.nan, float)
|
||||
|
||||
for i in range(n):
|
||||
sp = standpunkte[i]
|
||||
zp = zielpunkte[i]
|
||||
|
||||
blocks = []
|
||||
idx = []
|
||||
|
||||
try:
|
||||
if sp:
|
||||
b = idx_xyz(sp)
|
||||
blocks.append(b)
|
||||
idx += b
|
||||
if zp:
|
||||
b = idx_xyz(zp)
|
||||
blocks.append(b)
|
||||
idx += b
|
||||
except ValueError:
|
||||
continue
|
||||
|
||||
if not blocks:
|
||||
continue
|
||||
|
||||
idx = list(dict.fromkeys(idx)) # unique
|
||||
|
||||
# Δx_i aus Grenzstörung
|
||||
dl = np.zeros((n, 1))
|
||||
dl[i, 0] = GRZW[i]
|
||||
dx = Qxx @ (A.T @ (P @ dl))
|
||||
|
||||
dx_loc = dx[idx, :]
|
||||
Q_loc = Qxx[np.ix_(idx, idx)]
|
||||
|
||||
# EF lokal
|
||||
EF2 = (dx_loc.T @ np.linalg.solve(Q_loc, dx_loc)).item() / (s0 ** 2)
|
||||
EF[i] = np.sqrt(max(0.0, EF2))
|
||||
|
||||
# SP lokal 3D: max trace der 3x3 Punktblöcke
|
||||
tr_list = [np.trace(Qxx[np.ix_(b, b)]) for b in blocks]
|
||||
if not tr_list:
|
||||
continue
|
||||
|
||||
sigmaPmax_loc = s0 * np.sqrt(max(tr_list))
|
||||
SP_loc_m[i] = sigmaPmax_loc
|
||||
EFSP_loc_m[i] = EF[i] * sigmaPmax_loc
|
||||
|
||||
ausgabe_zuv = pd.DataFrame({
|
||||
"Beobachtung": labels,
|
||||
"Stand-Pkt": standpunkte,
|
||||
"Ziel-Pkt": zielpunkte,
|
||||
"r_i": ri,
|
||||
"EP_GF [mm]": EP_hat_m * 1000.0,
|
||||
"EP_grzw [mm]": EP_grzw_m * 1000.0,
|
||||
"EF": EF,
|
||||
"SP_loc_3D [mm]": SP_loc_m * 1000.0,
|
||||
"EF*SP_loc_3D [mm]": EFSP_loc_m * 1000.0,
|
||||
})
|
||||
return ausgabe_zuv
|
||||
@@ -1,25 +1,43 @@
|
||||
from Stochastisches_Modell import StochastischesModell
|
||||
from Netzqualität_Genauigkeit import Genauigkeitsmaße
|
||||
|
||||
from Datumsfestlegung import Datumsfestlegung
|
||||
import numpy as np
|
||||
import Export
|
||||
import importlib
|
||||
import Datenbank
|
||||
import importlib
|
||||
import Export
|
||||
from Export import Export as Exporter
|
||||
import Stochastisches_Modell
|
||||
import Funktionales_Modell
|
||||
import Koordinatentransformationen
|
||||
import Parameterschaetzung
|
||||
import Netzqualität_Genauigkeit
|
||||
import Datumsfestlegung
|
||||
import numpy as np
|
||||
import sympy as sp
|
||||
import pandas as pd
|
||||
|
||||
|
||||
def ausgleichung_global(A, dl, Q_ext, P):
|
||||
"""
|
||||
Führt eine Ausgleichung nach kleinsten Quadraten durch.
|
||||
|
||||
Aus der Designmatrix A, dem Verbesserungsvektor dl und der Gewichtsmatrix P wird das Normalgleichungssystem
|
||||
aufgestellt und gelöst. Anschließend werden Residuen sowie Kofaktor-Matrizen der Unbekannten und Beobachtungen berechnet.
|
||||
|
||||
Es werden folgende Berechnungsschitte durchgeführt:
|
||||
|
||||
1) Normalgleichungsmatrix N = Aᵀ · P · A und Absolutglied n = Aᵀ · P · dl
|
||||
2) Lösung dx = N⁻¹ · n
|
||||
3) Residuen v = A · dx − dl
|
||||
4) Kofaktormatrix der Unbekannten Q_xx
|
||||
5) Kofaktormatrix der ausgeglichenen Beobachtungen Q_ll_dach
|
||||
6) Kofaktormatrix der Verbesserungen Q_vv
|
||||
|
||||
:param A: Jacobi-Matrix (A-Matrix).
|
||||
:type A: array_like
|
||||
:param dl: Verbesserungsvektor bzw. Beobachtungsabweichungen.
|
||||
:type dl: array_like
|
||||
:param Q_ext: a-priori Kofaktormatrix der Beobachtungen.
|
||||
:type Q_ext: array_like
|
||||
:param P: Gewichtsmatrix der Beobachtungen.
|
||||
:type P: array_like
|
||||
:return: Dictionary mit Ausgleichungsergebnissen, Zuschlagsvektor dx.
|
||||
:rtype: tuple[dict[str, Any], numpy.ndarray]
|
||||
:raises numpy.linalg.LinAlgError: Wenn das Normalgleichungssystem singulär ist und nicht gelöst werden kann.
|
||||
"""
|
||||
A=np.asarray(A, float)
|
||||
dl = np.asarray(dl, float).reshape(-1, 1)
|
||||
Q_ext = np.asarray(Q_ext, float)
|
||||
@@ -62,7 +80,7 @@ def ausgleichung_global(A, dl, Q_ext, P):
|
||||
|
||||
|
||||
|
||||
def ausgleichung_lokal1(A, dl, Q_ll):
|
||||
def ausgleichung_lokal(A, dl, Q_ll):
|
||||
A = np.asarray(A, dtype=float)
|
||||
dl = np.asarray(dl, dtype=float).reshape(-1, 1)
|
||||
Q_ll = np.asarray(Q_ll, dtype=float)
|
||||
@@ -110,175 +128,6 @@ def ausgleichung_lokal1(A, dl, Q_ll):
|
||||
return dict_ausgleichung, dx
|
||||
|
||||
|
||||
|
||||
def ausgleichung_lokal(
|
||||
A, dl, Q_ll,
|
||||
*,
|
||||
datumfestlegung="hart",
|
||||
x0=None,
|
||||
unbekannten_liste=None,
|
||||
liste_punktnummern=None,
|
||||
datenbank=None,
|
||||
mit_massstab=True
|
||||
):
|
||||
A = np.asarray(A, float)
|
||||
dl = np.asarray(dl, float).reshape(-1, 1)
|
||||
Q_ll = np.asarray(Q_ll, float)
|
||||
|
||||
# 1) Gewichtsmatrix
|
||||
P = np.linalg.inv(Q_ll)
|
||||
|
||||
# 2) Normalgleichungen
|
||||
N = A.T @ P @ A
|
||||
n = A.T @ P @ dl
|
||||
u = N.shape[0]
|
||||
|
||||
# 3) Datumsfestlegung
|
||||
if datumfestlegung == "weiche Lagerung":
|
||||
# hier wurde Q_ll bereits extern erweitert → ganz normal lösen
|
||||
dx = np.linalg.solve(N, n)
|
||||
|
||||
elif datumfestlegung == "gesamtspur":
|
||||
if x0 is None or unbekannten_liste is None or liste_punktnummern is None:
|
||||
raise ValueError("gesamtspur benötigt x0, unbekannten_liste, liste_punktnummern")
|
||||
|
||||
G = Datumsfestlegung.build_G_from_names(
|
||||
x0,
|
||||
unbekannten_liste,
|
||||
liste_punktnummern,
|
||||
mit_massstab=mit_massstab
|
||||
)
|
||||
dx, k = Datumsfestlegung.berechne_dx_geraendert(N, n, G)
|
||||
|
||||
elif datumfestlegung == "teilspur":
|
||||
if x0 is None or unbekannten_liste is None or liste_punktnummern is None or datenbank is None:
|
||||
raise ValueError("teilspur benötigt x0, unbekannten_liste, liste_punktnummern, datenbank")
|
||||
|
||||
# G über alle Punkte
|
||||
G = Datumsfestlegung.build_G_from_names(
|
||||
x0,
|
||||
unbekannten_liste,
|
||||
liste_punktnummern,
|
||||
mit_massstab=mit_massstab
|
||||
)
|
||||
|
||||
# Auswahl aus DB
|
||||
liste_datumskoordinaten = datenbank.get_datumskoordinate()
|
||||
if not liste_datumskoordinaten:
|
||||
raise ValueError("Teilspur gewählt, aber keine Datumskoordinaten in der DB gesetzt.")
|
||||
|
||||
aktive = Datumsfestlegung.aktive_indices_from_selection(
|
||||
[(s[1:], s[0]) for s in liste_datumskoordinaten],
|
||||
unbekannten_liste
|
||||
)
|
||||
|
||||
E = Datumsfestlegung.auswahlmatrix_E(u, aktive)
|
||||
Gi = E @ G
|
||||
|
||||
dx, k = Datumsfestlegung.berechne_dx_geraendert(N, n, Gi)
|
||||
|
||||
else:
|
||||
raise ValueError(f"Unbekannte Datumsfestlegung: {datumfestlegung}")
|
||||
|
||||
# 4) Residuen
|
||||
v = dl - A @ dx
|
||||
|
||||
# 5) Kofaktormatrix der Unbekannten
|
||||
Q_xx = StochastischesModell.berechne_Q_xx(N)
|
||||
|
||||
# 6) Kofaktormatrix der Beobachtungen
|
||||
Q_ll_dach = StochastischesModell.berechne_Q_ll_dach(A, Q_xx)
|
||||
|
||||
# 7) Kofaktormatrix der Verbesserungen
|
||||
Q_vv = StochastischesModell.berechne_Qvv(Q_ll, Q_ll_dach)
|
||||
|
||||
dict_ausgleichung = {
|
||||
"dx": dx,
|
||||
"v": v,
|
||||
"P": P,
|
||||
"N": N,
|
||||
"Q_xx": Q_xx,
|
||||
"Q_ll_dach": Q_ll_dach,
|
||||
"Q_vv": Q_vv,
|
||||
"Q_ll": Q_ll,
|
||||
}
|
||||
|
||||
return dict_ausgleichung, dx
|
||||
|
||||
|
||||
def ausgleichung_spurminimierung(A, dl, Q_ll, *, datumfestlegung, x0, unbekannten_liste, datenbank=None, mit_massstab=True):
|
||||
"""
|
||||
Freie Ausgleichung mit Gesamtspur- oder Teilspurminimierung.
|
||||
Rückgabe-Layout wie ausgleichung_global.
|
||||
"""
|
||||
A = np.asarray(A, float)
|
||||
dl = np.asarray(dl, float).reshape(-1, 1)
|
||||
Q_ll = np.asarray(Q_ll, float)
|
||||
|
||||
# 1) Gewichtsmatrix P
|
||||
P = StochastischesModell.berechne_P(Q_ll)
|
||||
|
||||
# 2) Normalgleichungen
|
||||
N = A.T @ P @ A
|
||||
n = A.T @ P @ dl
|
||||
|
||||
# 3) G über ALLE Punkte aus unbekannten_liste (automatisch)
|
||||
G = Datumsfestlegung.build_G_from_names(
|
||||
x0=x0,
|
||||
unbekannten_liste=unbekannten_liste,
|
||||
liste_punktnummern=None, # <- wichtig: auto
|
||||
mit_massstab=mit_massstab
|
||||
)
|
||||
|
||||
# 4) Gesamtspur / Teilspur
|
||||
if datumfestlegung == "gesamtspur":
|
||||
Gi = G
|
||||
k = None # wird unten überschrieben
|
||||
|
||||
elif datumfestlegung == "teilspur":
|
||||
if datenbank is None:
|
||||
raise ValueError("teilspur benötigt datenbank mit get_datumskoordinate()")
|
||||
liste_datumskoordinaten = datenbank.get_datumskoordinate()
|
||||
if not liste_datumskoordinaten:
|
||||
raise ValueError("Teilspur gewählt, aber keine Datumskoordinaten in der DB gesetzt.")
|
||||
|
||||
# ["X10034","Y10034"] -> [("10034","X"),("10034","Y")]
|
||||
auswahl = [(s[1:], s[0]) for s in liste_datumskoordinaten]
|
||||
aktive = Datumsfestlegung.aktive_indices_from_selection(auswahl, unbekannten_liste)
|
||||
|
||||
E = Datumsfestlegung.auswahlmatrix_E(N.shape[0], aktive)
|
||||
Gi = E @ G
|
||||
|
||||
else:
|
||||
raise ValueError("datumfestlegung muss 'gesamtspur' oder 'teilspur' sein")
|
||||
|
||||
# 5) Lösung per Ränderung + korrektes Q_xx
|
||||
dx, k, Q_xx = Datumsfestlegung.loese_geraendert_mit_Qxx(N, n, Gi)
|
||||
|
||||
# 6) Residuen (wie bei dir)
|
||||
v = A @ dx - dl
|
||||
|
||||
# 7) Kofaktormatrix der Beobachtungen
|
||||
Q_ll_dach = StochastischesModell.berechne_Q_ll_dach(A, Q_xx)
|
||||
|
||||
# 8) Kofaktormatrix der Verbesserungen
|
||||
Q_vv = StochastischesModell.berechne_Qvv(Q_ll, Q_ll_dach)
|
||||
|
||||
dict_ausgleichung = {
|
||||
"dx": dx,
|
||||
"v": v,
|
||||
"P": P,
|
||||
"N": N,
|
||||
"n": n,
|
||||
"k": k,
|
||||
"Q_xx": Q_xx,
|
||||
"Q_ll_dach": Q_ll_dach,
|
||||
"Q_vv": Q_vv,
|
||||
"Q_ll": Q_ll,
|
||||
"datumfestlegung": datumfestlegung,
|
||||
}
|
||||
return dict_ausgleichung, dx
|
||||
|
||||
class Iterationen:
|
||||
"""Iterative Ausgleichung auf Basis des funktionalen und stochastischen Modells.
|
||||
|
||||
|
||||
40
Proben.py
40
Proben.py
@@ -1,6 +1,23 @@
|
||||
import numpy as np
|
||||
# d
|
||||
|
||||
def atpv_probe(A, P, v, tol=1e-7):
|
||||
"""
|
||||
Führt die ATPv-Probe zur Kontrolle der Lösung des Normalgleichungssystems durch.
|
||||
|
||||
Die Funktion überprüft, ob der Ausdruck Aᵀ · P · v näherungsweise Null ist.
|
||||
Die Prüfung erfolgt unter Verwendung einer vorgegebenen Toleranz.
|
||||
|
||||
:param A: Jacobi-Matrix (A-Matrix).
|
||||
:type A: array_like
|
||||
:param P: Gewichtsmatrix der Beobachtungen.
|
||||
:type P: array_like
|
||||
:param v: Residuenvektor der Beobachtungen.
|
||||
:type v: array_like
|
||||
:param tol: Absolute Toleranz für den Vergleich mit Null.
|
||||
:type tol: float
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
A = np.asarray(A, float)
|
||||
P = np.asarray(P, float)
|
||||
v = np.asarray(v, float).reshape(-1, 1)
|
||||
@@ -14,6 +31,27 @@ def atpv_probe(A, P, v, tol=1e-7):
|
||||
|
||||
|
||||
def hauptprobe(A, x, l, v, tol=1e-7):
|
||||
"""
|
||||
Führt die Hauptprobe zur Überprüfung der berechneten Residuen durch.
|
||||
|
||||
Die Hauptprobe kontrolliert, ob die Residuen v mit der Beziehung
|
||||
v = A · x − l übereinstimmen. Stimmen der berechnete Residuenvektor
|
||||
und der über das funktionale Modell rekonstruierte Residuenvektor
|
||||
innerhalb der Toleranz überein, gilt die Ausgleichung als konsistent.
|
||||
|
||||
:param A: Jacobi-Matrix (A-Matrix).
|
||||
:type A: array_like
|
||||
:param x: Lösungsvektor der Unbekannten.
|
||||
:type x: array_like
|
||||
:param l: Beobachtungsvektor.
|
||||
:type l: array_like
|
||||
:param v: Residuenvektor aus der Ausgleichung.
|
||||
:type v: array_like
|
||||
:param tol: Absolute Toleranz für den Vergleich der Residuen.
|
||||
:type tol: float
|
||||
:return: None
|
||||
:rtype: None
|
||||
"""
|
||||
A = np.asarray(A, float)
|
||||
x = np.asarray(x, float).reshape(-1, 1)
|
||||
l = np.asarray(l, float).reshape(-1, 1)
|
||||
|
||||
@@ -44,7 +44,7 @@ class StochastischesModell:
|
||||
self.db_zugriff = Datenbankzugriff(self.pfad_datenbank)
|
||||
|
||||
def Qll_symbolisch(self, liste_beobachtungen_symbolisch: list) -> sp.Matrix:
|
||||
"""Erstellt die symbolische Varianz-Kovarianz-Matrix Qll der Beobachtungen.
|
||||
"""Erstellt die symbolische Kofaktormatrix der Beobachtungen.
|
||||
|
||||
Aus den symbolischen Beobachtungskennungen wird die Beobachtungsart abgeleitet (Tachymeter: SD/R/ZW,
|
||||
GNSS: gnssbx/gnssby/gnssbz, Geometrisches Nivellement: niv). Für jede Beobachtung wird eine symbolische Varianzgleichung
|
||||
@@ -61,7 +61,7 @@ class StochastischesModell:
|
||||
|
||||
:param liste_beobachtungen_symbolisch: Liste der symbolischen Beobachtungskennungen.
|
||||
:type liste_beobachtungen_symbolisch: list
|
||||
:return: Symbolische Varianz-Kovarianz-Matrix Qll.
|
||||
:return: Symbolische Kofaktormatrix Qll.
|
||||
:rtype: sp.Matrix
|
||||
"""
|
||||
liste_standardabweichungen_symbole = []
|
||||
@@ -204,7 +204,7 @@ class StochastischesModell:
|
||||
return Qll
|
||||
|
||||
def Qll_numerisch(self, Qll_Matrix_Symbolisch: sp.Matrix, liste_beobachtungen_symbolisch: list) -> np.Matrix:
|
||||
"""Erstellt eine numerische Varianz-Kovarianz-Matrix aus einer symbolischen Qll-Matrix.
|
||||
"""Erstellt eine numerische Kofaktormatrix der Beobachtungen aus einer symbolischen Qll-Matrix.
|
||||
|
||||
Es werden die zur Substitution benötigten Werte aus der Datenbank abgefragt und den in Qll vorkommenden Symbolen zugeordnet,
|
||||
u. a.:
|
||||
@@ -221,11 +221,11 @@ class StochastischesModell:
|
||||
|
||||
Die numerische Matrix wird als CSV-Datei in Zwischenergebnisse\\Qll_Numerisch.csv exportiert.
|
||||
|
||||
:param Qll_Matrix_Symbolisch: Symbolische Varianz-Kovarianz-Matrix Qll.
|
||||
:param Qll_Matrix_Symbolisch: Symbolische Kofaktormatrix der Beobachtungen Qll.
|
||||
:type Qll_Matrix_Symbolisch: sp.Matrix
|
||||
:param liste_beobachtungen_symbolisch: Liste der symbolischen Beobachtungskennungen.
|
||||
:type liste_beobachtungen_symbolisch: list
|
||||
:return: Numerische Varianz-Kovarianz-Matrix Qll als Numpy-Array.
|
||||
:return: Numerische Kofaktormatrix Qll als Numpy-Array.
|
||||
:rtype: np.Matrix
|
||||
:raises ValueError: Falls Symbole in Qll_Matrix_Symbolisch enthalten sind, für die keine Substitutionen vorhanden sind.
|
||||
"""
|
||||
@@ -380,7 +380,7 @@ class StochastischesModell:
|
||||
return Qll_numerisch
|
||||
|
||||
def QAA_symbolisch(self, liste_beobachtungen_symbolisch: list) -> np.Matrix:
|
||||
"""Erstellt die symbolische Varianz-Kovarianz-Matrix QAA der Anschlusspunkte (weiche Lagerung).
|
||||
"""Erstellt die symbolische Kofaktormatrix der Anschlusspunkte (weiche Lagerung).
|
||||
|
||||
Es werden ausschließlich Beobachtungen berücksichtigt, deren Kennung mit "lA_" beginnt. Für jede Anschlussbedingung
|
||||
wird eine (symbolische) Standardabweichung StabwAA_* angesetzt und mit der Varianzkomponente der Beobachtungsgruppe
|
||||
@@ -390,7 +390,7 @@ class StochastischesModell:
|
||||
|
||||
:param liste_beobachtungen_symbolisch: Liste der symbolischen Beobachtungskennungen.
|
||||
:type liste_beobachtungen_symbolisch: list
|
||||
:return: Symbolische Varianz-Kovarianz-Matrix QAA.
|
||||
:return: Symbolische Kofaktormatrix QAA.
|
||||
:rtype: sp.Matrix
|
||||
"""
|
||||
liste_standardabweichungen_symbole = []
|
||||
@@ -428,11 +428,11 @@ class StochastischesModell:
|
||||
|
||||
Die numerische Matrix wird als CSV-Datei in Zwischenergebnisse\\QAA_Numerisch.csv exportiert.
|
||||
|
||||
:param QAA_Matrix_Symbolisch: Symbolische Varianz-Kovarianz-Matrix QAA.
|
||||
:param QAA_Matrix_Symbolisch: Symbolische Kofaktormatrix QAA.
|
||||
:type QAA_Matrix_Symbolisch: sp.Matrix
|
||||
:param liste_beobachtungen_symbolisch: Liste der symbolischen Beobachtungskennungen.
|
||||
:type liste_beobachtungen_symbolisch: list
|
||||
:return: Numerische Varianz-Kovarianz-Matrix QAA als Numpy-Array.
|
||||
:return: Numerische VKofaktormatrix QAA als Numpy-Array.
|
||||
:rtype: np.Matrix
|
||||
"""
|
||||
# Symbolische Listen
|
||||
@@ -483,11 +483,11 @@ class StochastischesModell:
|
||||
|
||||
@staticmethod
|
||||
def berechne_P(Q_ll: np.ndarray) -> np.ndarray:
|
||||
"""Berechnet die Gewichtsmatrix P aus einer Varianz-Kovarianz-Matrix Qll.
|
||||
"""Berechnet die Gewichtsmatrix P aus der Kofaktormatrix.
|
||||
|
||||
Die Gewichtsmatrix wird als Inverse von Qll gebildet: P = inv(Qll).
|
||||
|
||||
:param Q_ll: Varianz-Kovarianz-Matrix der Beobachtungen.
|
||||
:param Q_ll: Kofaktormatrix der Beobachtungen.
|
||||
:type Q_ll: np.ndarray
|
||||
:return: Gewichtsmatrix P.
|
||||
:rtype: np.ndarray
|
||||
@@ -515,7 +515,7 @@ class StochastischesModell:
|
||||
|
||||
@staticmethod
|
||||
def berechne_Q_ll_dach(A: np.ndarray, Q_xx: np.ndarray) -> np.ndarray:
|
||||
"""Berechnet die (geschätzte) Varianz-Kovarianz-Matrix der Beobachtungen Qll_dach.
|
||||
"""Berechnet die (geschätzte) Kofaktormatrix der Beobachtungen Qll_dach.
|
||||
|
||||
Die Matrix wird gemäß Qll_dach = A @ Qxx @ A.T gebildet.
|
||||
|
||||
@@ -523,7 +523,7 @@ class StochastischesModell:
|
||||
:type A: np.ndarray
|
||||
:param Q_xx: Kofaktormatrix der Unbekannten.
|
||||
:type Q_xx: np.ndarray
|
||||
:return: Geschätzte Varianz-Kovarianz-Matrix der Beobachtungen Qll_dach.
|
||||
:return: Geschätzte VKofaktormatrix der Beobachtungen Qll_dach.
|
||||
:rtype: np.ndarray
|
||||
"""
|
||||
Q_ll_dach = A @ Q_xx @ A.T
|
||||
@@ -531,16 +531,16 @@ class StochastischesModell:
|
||||
|
||||
@staticmethod
|
||||
def berechne_Qvv(Q_ll: np.ndarray, Q_ll_dach: np.ndarray) -> np.ndarray:
|
||||
"""Berechnet die Varianz-Kovarianz-Matrix der Residuen Qvv.
|
||||
"""Berechnet die Kofaktormatrix der Residuen Qvv.
|
||||
|
||||
Die Residuenkovarianz wird als Differenz aus Beobachtungs-Kovarianz und dem durch das Modell erklärten Anteil gebildet:
|
||||
Qvv = Qll - Qll_dach.
|
||||
|
||||
:param Q_ll: Varianz-Kovarianz-Matrix der Beobachtungen.
|
||||
:param Q_ll: Kofaktormatrix der Beobachtungen.
|
||||
:type Q_ll: np.ndarray
|
||||
:param Q_ll_dach: Geschätzte Varianz-Kovarianz-Matrix der Beobachtungen.
|
||||
:param Q_ll_dach: Geschätzte Kofaktormatrix der Beobachtungen.
|
||||
:type Q_ll_dach: np.ndarray
|
||||
:return: Varianz-Kovarianz-Matrix der Residuen Qvv.
|
||||
:return: Kofaktormatrix der Residuen Qvv.
|
||||
:rtype: np.ndarray
|
||||
"""
|
||||
Q_vv = Q_ll - Q_ll_dach
|
||||
|
||||
@@ -5,7 +5,7 @@ import pandas as pd
|
||||
|
||||
import Datenbank
|
||||
import Datumsfestlegung
|
||||
from Netzqualität_Genauigkeit import Genauigkeitsmaße
|
||||
from Netzqualitaet_Genauigkeit import Genauigkeitsmaße
|
||||
|
||||
class VKS:
|
||||
"""Varianzkomponentenschätzung (VKS) mit Anpassung durch Benutzereingaben.
|
||||
|
||||
Reference in New Issue
Block a user