Merge pull request #15 from GeoscienceAustralia/orbitq-stat-fix

fixing stats + residual estimation

gnssanalysis/gn_diffaux.py

@@ -626,12 +626,13 @@ def rac_df_to_rms_df(rac_df):
     RMS of dX, dY, dZ and R (Radial), A (Along-track), C (Cross-track) per satellite (PRN).
     Additionally, a 1D mean difference, (dX + dY + dZ)/3, and 3D difference, sqrt(R^2 + A^2 + C^2), are computed for each PRN.
     """
-
     merged_data = rac_df.join(rac_df.attrs["diff_eci"])
     rms_df = merged_data.pow(2).groupby("PRN").mean().pow(0.5)
+    std_df = merged_data.groupby("PRN").std(ddof=0)
 
     rms_df["EST", "MEAN"] = rms_df.EST.mean(axis=1).groupby("PRN").mean()
     rms_df["EST", "AVG"] = rac_df.attrs["diff_eci"].EST.mean(axis=1).groupby("PRN").mean()
+    std_df["EST", "3D_RMS"] = rac_df.attrs["diff_eci"].EST.pow(2).sum(axis=1).pow(0.5).groupby("PRN").std(ddof=0)
     rms_df["EST", "3D_RMS"] = rac_df.attrs["diff_eci"].EST.pow(2).sum(axis=1).groupby("PRN").mean().pow(0.5)
     rms_df = rms_df.droplevel(0, axis=1).rename(
         columns={
@@ -643,10 +644,19 @@ def rac_df_to_rms_df(rac_df):
             "Z": "dZ_RMS",
         }
     )
-
+    std_df = std_df.droplevel(0, axis=1).rename(
+        columns={
+            "Radial": "R_RMS",
+            "Along-track": "A_RMS",
+            "Cross-track": "C_RMS",
+            "X": "dX_RMS",
+            "Y": "dY_RMS",
+            "Z": "dZ_RMS",
+        }
+    )
     # summarising over all SVs
     summary_df = _pd.DataFrame(
-        [rms_df.mean(axis=0), rms_df.std(axis=0), rms_df.pow(2).mean(axis=0).pow(0.5)], index=["AVG", "STD", "RMS"]
+        [rms_df.mean(axis=0), std_df.mean(axis=0), rms_df.pow(2).mean(axis=0).pow(0.5)], index=["AVG", "STD", "RMS"]
     )
 
     rms_df.attrs["summary"] = summary_df

gnssanalysis/gn_transform.py

@@ -36,7 +36,7 @@ def get_helmert7(pt1:_np.ndarray, pt2:_np.ndarray, scale_in_ppm:bool = True):
     """inversion of 7 Helmert parameters between 2 sets of points. pt1@hlm -> pt2"""
     A, rhs = gen_helm_aux(pt1, pt2)
     sol = list(_np.linalg.lstsq(A, rhs, rcond=-1))  # parameters
-    sol[0] = sol[0].flatten()  # flattening the HLM params arr to [Tx, Ty, Tz, Rx, Ry, Rz, Scale/mu]
+    sol[0] = sol[0].flatten()*-1.0  # flattening the HLM params arr to [Tx, Ty, Tz, Rx, Ry, Rz, Scale/mu]
     if scale_in_ppm:
         sol[0][-1] *= 1e6 # scale in ppm
     res = rhs - A @ sol[0]  # computing residuals for pt2