Push

Netzqualität_Zuverlässigkeit.py

@@ -49,9 +49,8 @@ class Zuverlaessigkeit:
             interpretation = (
                 "Nullhypothese H₀ verworfen!\n"
                 "Dies kann folgende Gründe haben:\n"
-                "→ Es befinden sich grobe Fehler im Datenmaterial.\n"
-                "→ Das funktionale Modell ist fehlerhaft.\n"
-                "→ Das stochastische Modell ist zu optimistisch."
+                "→ 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 {
@@ -111,3 +110,107 @@ class Zuverlaessigkeit:
             "δ0": nzp,
         })
         return Lokaltest_innere_Zuv
+
+
+    def aeussere_zuverlaessigkeit_EF_EP(Lokaltest, labels, Qxx, A, P, s0_apost, unbekannten_liste, x):
+        df = Lokaltest.copy()
+        labels = list(labels)
+        Qxx = np.asarray(Qxx, float)
+        A = np.asarray(A, float)
+        P = np.asarray(P, float)
+        x = np.asarray(x, float).reshape(-1)
+        ri = df["r_i"].astype(float).to_numpy()
+        GF = df["GF_i"].astype(float).to_numpy()
+        s_vi = df["s_vi"].astype(float).to_numpy()
+        GRZW = df["GRZW_i"].astype(float).to_numpy()
+        nzp = df["δ0"].astype(float).to_numpy()
+        n = A.shape[0]  # Anzahl Beobachtungen
+        u = A.shape[1]  # Anzahl Unbekannte
+
+        # Einflussfaktor EF berechnen
+        EF = np.zeros(n, dtype=float)
+
+        for i in range(n):
+            # 1) ∇l_i aufstellen
+            nabla_l = np.zeros((n, 1))
+            nabla_l[i, 0] = GRZW[i]
+
+            # 2) ∇x_i = Qxx * A^T * P * ∇l_i
+            nabla_x = Qxx @ (A.T @ (P @ nabla_l))
+
+            # 3) EF_i^2 = (∇x_i^T * Qxx^{-1} * ∇x_i) / s0^2
+            Qxx_inv_nabla_x = np.linalg.solve(Qxx, nabla_x)  # = Qxx^{-1} ∇x_i
+            #EF2 = float((nabla_x.T @ Qxx_inv_nabla_x) / (float(s0_apost) ** 2)).item()
+            EF2 = ((nabla_x.T @ Qxx_inv_nabla_x) / (float(s0_apost) ** 2)).item()
+
+            EF[i] = np.sqrt(EF2)
+
+        # Koordinaten-Dict aus x
+        coords = {}
+        j = 0
+        while j < len(unbekannten_liste):
+            name = str(unbekannten_liste[j])
+            if name.startswith("X"):
+                pn = name[1:]
+                coords[pn] = (x[j], x[j + 1], x[j + 2])
+                j += 3
+            else:
+                j += 1
+
+        # EP + Standpunkte
+        EP_m = np.full(len(labels), np.nan, dtype=float)
+        standpunkte = [""] * len(labels)
+
+        for i, lbl in enumerate(labels):
+            parts = lbl.split("_")
+            sp = None
+            zp = None
+
+            # Tachymeter: ID_SD_GRP_SP_ZP / ID_R_GRP_SP_ZP / ID_ZW_GRP_SP_ZP
+            if ("_SD_" in lbl) or ("_R_" in lbl) or ("_ZW_" in lbl):
+                if len(parts) >= 5:
+                    sp = parts[3].strip()
+                    zp = parts[4].strip()
+
+            # GNSS: *_gnssbx_SP_ZP etc.
+            if ("gnss" in lbl) and (len(parts) >= 4):
+                sp = parts[-2].strip()
+                zp = parts[-1].strip()
+
+            standpunkte[i] = sp if sp is not None else ""
+            one_minus_r = (1.0 - ri[i])
+
+            # SD + GNSS: direkt in m
+            if ("_SD_" in lbl) or ("gnss" in lbl):
+                EP_m[i] = one_minus_r * GF[i]
+
+            # R / ZW: Winkel -> Streckenäquivalent über s
+            elif ("_R_" in lbl) or ("_ZW_" in lbl):
+                if sp and zp and (sp in coords) and (zp in coords):
+                    X1, Y1, Z1 = coords[sp]
+                    X2, Y2, Z2 = coords[zp]
+                    s = float(np.sqrt((X2 - X1) ** 2 + (Y2 - Y1) ** 2 + (Z2 - Z1) ** 2))
+                    EP_m[i] = one_minus_r * ((GF[i]) * s)
+
+        # SP am Standpunkt (2D)
+        diagQ = np.diag(Qxx)
+        SP_cache_mm = {}
+
+        for sp in set([s for s in standpunkte if s]):
+            idx_x = [k for k, sym in enumerate(unbekannten_liste) if str(sym) == f"X{sp}"][0]
+            qx = diagQ[idx_x]
+            qy = diagQ[idx_x + 1]
+            SP_cache_mm[sp] = float(s0_apost) * np.sqrt(qx + qy) * 1000.0
+
+        SP_mm = np.array([SP_cache_mm.get(sp, np.nan) for sp in standpunkte], dtype=float)
+
+        out = pd.DataFrame({
+            "Beobachtung": labels,
+            "Stand-Pkt": standpunkte,
+            "EF": EF,
+            "EP [mm]": EP_m * 1000.0,
+            "SP [mm]": SP_mm,
+            "EF*SP [mm]": EF * SP_mm,
+        })
+        return out
+