Merge branch 'master' of github.com:kartverket/where

analysis/check_nyale13s/plot_nyale13s.py

@@ -210,23 +210,38 @@ def plot_statistics(dates, dof, variance_factor, colors):
     #fig.tight_layout()
     plt.savefig(f"img/{dset_id}/Statistics_{dset_id}_{start:%Y-%m-%d}_{end:%Y-%m-%d}.png", bbox_inches='tight')
 
-def plot_baseline(dates, baseline_length, baseline_length_ferr, num_obs_bs, sta_1, sta_2, local_tie):
+def plot_baseline(dates, baseline_length, baseline_length_ferr, num_obs_bs, vgos, sta_1, sta_2, local_tie):
     t = Time(dates, fmt="datetime", scale="utc")
     _wblr, wmean = wblr(baseline_length,t.mjd, baseline_length_ferr)
+
+    vgos = np.array(vgos)
+    dates = np.array(dates)
+    baseline_length = np.array(baseline_length)
+    baseline_length_ferr = np.array(baseline_length_ferr)
+    num_obs_bs = np.array(num_obs_bs)
+
+    num_vgos = np.sum(vgos)
+    num_sx = np.sum(~vgos)
+
     if local_tie:
         print(f"{dset_id}: Weighted mean: {wmean: 6.4f} [m], Local tie offset: {(wmean - local_tie)*Unit.m2mm: 6.2f} [mm]")
     else:
         print(f"{dset_id}: Weighted mean: {wmean: 16.4f} [m]")
     print(f"{dset_id}: Weighted blr (using weighted mean): {_wblr*Unit.m2mm: 6.2f} [mm]")
+    print(f"Number of S/X sessions: {num_sx}")
+    print(f"Number of VGOS sessions: {num_vgos}")
 
     x_padding = 1.5 # days
     y_padding = 0.01 # meter
 
     fig = plt.figure(figsize=(12, 8), dpi=150)
     xlim = (min(dates) - timedelta(days=x_padding), max(dates) + timedelta(days=x_padding))
     ylim = (wmean - y_padding, wmean + y_padding)
-    plt.errorbar(dates, baseline_length, yerr=baseline_length_ferr, fmt="o", marker=None, zorder=0, mew=0, ecolor="tab:gray")
-    im = plt.scatter(dates, baseline_length, c=num_obs_bs, zorder=100)
+    plt.errorbar(dates, baseline_length, yerr=baseline_length_ferr, fmt=",", marker=None, zorder=0, mew=0, ecolor="tab:gray")
+    if num_sx > 0:
+        im = plt.scatter(dates[~vgos], baseline_length[~vgos], c=num_obs_bs[~vgos], marker="o", zorder=100, label="S/X")
+    if num_vgos > 0:
+        im = plt.scatter(dates[vgos], baseline_length[vgos], c=num_obs_bs[vgos], marker="d", zorder=100, label="VGOS")
     plt.axhline(wmean, c="red",label="Weighted mean")
     if local_tie:
        plt.axhline(local_tie, c="blue", label="Local tie vector")
@@ -329,6 +344,7 @@ def get_state_values(dset, fieldname, fill_value=np.nan):
 baseline_length = []
 baseline_length_ferr = []
 baseline_dates = []
+vgos = []
 
 rms = dset_ts.rms_residual_estimate[idx]
 dof = dset_ts.degrees_of_freedom[idx]
@@ -385,6 +401,10 @@ def get_state_values(dset, fieldname, fill_value=np.nan):
         baseline_length.append(bl['baseline_length'])
         baseline_length_ferr.append(bl['baseline_length_ferr'])
         num_obs_bs.append(dset_session.num(baseline=baseline1) + dset_session.num(baseline=baseline2))
+        if session_code.startswith("v"):
+            vgos.append(True)
+        else:
+            vgos.append(False)
 
 
     if args.plot_residuals:
@@ -399,7 +419,7 @@ def get_state_values(dset, fieldname, fill_value=np.nan):
 local_tie = data.get(baseline1) or data.get(baseline2)
 
 plot_statistics(dates, dof, variance_factor, colors[0])
-plot_baseline(baseline_dates, baseline_length, baseline_length_ferr, num_obs_bs, station1, station2, local_tie)
+plot_baseline(baseline_dates, baseline_length, baseline_length_ferr, num_obs_bs, vgos, station1, station2, local_tie)
 
 #plot_residual_rms(dates, dates_sta_1, dates_sta_2, rms, rms_sta_1, rms_sta_2, station1, station2, colors)
 

analysis/vlbi_scale/vlbi_scale.py

@@ -174,7 +174,9 @@
         r_3 = dset_ts.neq_helmert_gamma[idx_2][sort_idx_2]
         volume = dset_ts.network_volume[idx_2][sort_idx_2]
         for rd, sc, s, tx, ty, tz, r1, r2, r3, v in zip(rundate, session_code, scale, t_x, t_y, t_z, r_1, r_2, r_3, volume):
-            session_idx = dset_ts.filter(rundate=rd)
+            session_idx = dset_ts.filter(rundate=rd, session_code=sc)
             station_idx = dset_ts.num_obs_estimate[session_idx] > 0
-            stations = "/".join(list(dset_ts.station[session_idx][station_idx])[1:]) # First station entry is 'all' for each rundate
+            stations = list(dset_ts.station[session_idx][station_idx])
+            stations.remove("all") 
+            stations = "/".join(stations)
             fid.write(f"{rd}, {sc}, {s}, {tx}, {ty}, {tz}, {r1}, {r2}, {r3}, {v}, {stations}\n")

config/trf.conf

@@ -11,18 +11,18 @@
 # VLBI
 #
 
-[-26-__vlbi]
-cdp           = -26-
+[1526__vlbi]
+cdp           = 1526
 pos_itrs      = -2354891.0740, -4647166.2620, 3668871.7118
 name          = DSS26
 
-[-36-__vlbi]
-cdp           = -36-
+[1536__vlbi]
+cdp           = 1536
 pos_itrs      = -4461169.0272, 2682814.7152, -3674083.1649
 name          = DSS36
 
-[-56-__vlbi]
-cdp           = -56-
+[1556__vlbi]
+cdp           = 1556
 pos_itrs      = 4849421.3311, -360549.3026, 4114647.0931
 name          = DSS56
 
@@ -39,7 +39,15 @@ vel_itrs       =     -0.01280        0.00563        0.01021
 ref_epoch     = 2018.0
 name          = NYALE13N
 
-
+[7399__vlbi]
+cdp           = 1399
+#pos_itrs      = 5085504.75069 ,  2668369.305995, -2768492.965349
+# from where estimation
+pos_itrs      = 5085503.170463,  2668368.531475, -2768492.810447
+#pos_itrs      = 5085505.433,     2668370.025,    -2768494.876
+#vel_itrs      = 
+#ref_epoch     =
+name          = HARTVGS
 # SLR
 #
 # example: 

config/where_pipeline_vlbi.conf

@@ -21,6 +21,7 @@ site:add_sections
 
 # 2023b: yearly ITRF updates without non-tidal atm loading
 [vlbi__2023b]
+filekey_suffix              = vlbi, 2023b
 site                        = ocean_tides, solid_tides, solid_ptides, ocean_ptides, atm_tides,
                               eccentricity
 site:add_sections
@@ -92,7 +93,7 @@ estimate_constant:add_sections
 estimate_constant:help      = Parameters to estimate as constant
 estimate_constant:type      = List[enum[...]]
 
-estimate_obs_rejectors      = rms, vlbi_obs_per_station, vlbi_obs_per_source
+estimate_obs_rejectors      = random, rms, vlbi_obs_per_station, vlbi_obs_per_source
 estimate_obs_rejectors:help = Which algorithms to use to discard observations between iterations
 estimate_obs_rejectors:add_sections
 estimate_obs_rejectors_independent = False
@@ -112,7 +113,7 @@ output:add_sections
 files_to_publish            = output_sinex
 files_to_publish:help       = List of file keys to publish
 
-[vlbi_nyales_test]
+[vlbi_random]
 estimate_obs_rejectors      = random, rms, vlbi_obs_per_station, vlbi_obs_per_source
 
 [vlbi__mastertest]