update station term for california

Author: zhuwq0

examples/california/plotting.py

@@ -33,7 +33,7 @@ def plotting(stations, figure_path, config, picks, events_old, locations, statio
     im = ax[0, 1].scatter(
         stations["x_km"],
         stations["y_km"],
-        c=stations["station_term"],
+        c=stations["station_term_time_p"],
         cmap="viridis_r",
         s=100,
         marker="^",
@@ -43,7 +43,7 @@ def plotting(stations, figure_path, config, picks, events_old, locations, statio
     ax[0, 1].set_ylim([ymin, ymax])
     cbar = fig.colorbar(im, ax=ax[0, 1])
     cbar.set_label("Residual (s)")
-    ax[0, 1].set_title(f"Station term: {np.mean(np.abs(stations['station_term'].values)):.4f} s")
+    ax[0, 1].set_title(f"Station term (P): {np.mean(np.abs(stations['station_term_time_p'].values)):.4f} s")
 
     im = ax[1, 0].scatter(
         locations["x_km"],
@@ -265,14 +265,21 @@ def plotting_ransac(stations, figure_path, config, picks, events_init, events, s
     plt.savefig(os.path.join(figure_path, f"error{suffix}.png"), bbox_inches="tight", dpi=300)
     plt.close(fig)
 
-    xmin, xmax = config["xlim_km"]
-    ymin, ymax = config["ylim_km"]
-    zmin, zmax = config["zlim_km"]
-    vmin, vmax = config["zlim_km"]
+    # xmin, xmax = config["xlim_km"]
+    # ymin, ymax = config["ylim_km"]
+    # zmin, zmax = config["zlim_km"]
+    # vmin, vmax = config["zlim_km"]
+    xmin = events["x_km"].min()
+    xmax = events["x_km"].max()
+    ymin = events["y_km"].min()
+    ymax = events["y_km"].max()
+    zmin = events["z_km"].min()
+    zmax = events["z_km"].max()
+    vmin, vmax = zmin, zmax
     events = events.sort_values("time", ascending=True)
     s = max(0.1, min(10, 5000 / len(events)))
     alpha = 0.8
-    fig, ax = plt.subplots(2, 3, figsize=(18, 8), gridspec_kw={"height_ratios": [2, 1]})
+    fig, ax = plt.subplots(2, 4, figsize=(25, 8), gridspec_kw={"height_ratios": [2, 1]})
     # fig, ax = plt.subplots(2, 3, figsize=(15, 8), gridspec_kw={"height_ratios": [2, 1]})
     im = ax[0, 0].scatter(
         events["x_km"],
@@ -299,7 +306,7 @@ def plotting_ransac(stations, figure_path, config, picks, events_init, events, s
     im = ax[0, 1].scatter(
         stations["x_km"],
         stations["y_km"],
-        c=stations["station_term_time"],
+        c=stations["station_term_time_p"],
         cmap="viridis_r",
         s=100,
         marker="^",
@@ -312,12 +319,12 @@ def plotting_ransac(stations, figure_path, config, picks, events_init, events, s
     ax[0, 1].set_ylabel("Y (km)")
     cbar = fig.colorbar(im, ax=ax[0, 1])
     cbar.set_label("Residual (s)")
-    ax[0, 1].set_title(f"Station term: {np.mean(np.abs(stations['station_term_time'].values)):.4f} s")
+    ax[0, 1].set_title(f"Station term (P): {np.mean(np.abs(stations['station_term_time_p'].values)):.4f} s")
 
     im = ax[0, 2].scatter(
         stations["x_km"],
         stations["y_km"],
-        c=stations["station_term_amplitude"],
+        c=stations["station_term_time_s"],
         cmap="viridis_r",
         s=100,
         marker="^",
@@ -329,8 +336,26 @@ def plotting_ransac(stations, figure_path, config, picks, events_init, events, s
     ax[0, 2].set_xlabel("X (km)")
     ax[0, 2].set_ylabel("Y (km)")
     cbar = fig.colorbar(im, ax=ax[0, 2])
+    cbar.set_label("Residual (s)")
+    ax[0, 2].set_title(f"Station term (S): {np.mean(np.abs(stations['station_term_time_s'].values)):.4f} s")
+
+    im = ax[0, 3].scatter(
+        stations["x_km"],
+        stations["y_km"],
+        c=stations["station_term_amplitude"],
+        cmap="viridis_r",
+        s=100,
+        marker="^",
+        alpha=alpha,
+    )
+    ax[0, 3].set_aspect("equal", "box")
+    ax[0, 3].set_xlim([xmin, xmax])
+    ax[0, 3].set_ylim([ymin, ymax])
+    ax[0, 3].set_xlabel("X (km)")
+    ax[0, 3].set_ylabel("Y (km)")
+    cbar = fig.colorbar(im, ax=ax[0, 3])
     cbar.set_label("Residual (log10 cm/s)")
-    ax[0, 2].set_title(f"Station term: {np.mean(np.abs(stations['station_term_amplitude'].values)):.4f} s")
+    ax[0, 3].set_title(f"Station term: {np.mean(np.abs(stations['station_term_amplitude'].values)):.4f} (log10 cm/s)")
 
     ## Separate P and S station term
     # im = ax[0, 1].scatter(
@@ -402,5 +427,146 @@ def plotting_ransac(stations, figure_path, config, picks, events_init, events, s
     # ax[1, 1].set_ylabel("Depth (km)")
     cbar = fig.colorbar(im, ax=ax[1, 1])
     cbar.set_label("Depth (km)")
+
+    if "magnitude" in events.columns:
+        im = ax[1, 2].hist(events["magnitude"], bins=30, edgecolor="white")
+        ax[1, 2].set_yscale("log")
+        ax[1, 2].set_xlabel("Magnitude")
+        ax[1, 2].set_ylabel("Count")
+
     plt.savefig(os.path.join(figure_path, f"location{suffix}.png"), bbox_inches="tight", dpi=300)
     plt.close(fig)
+
+
+# %%
+def plotting_dd(events, stations, config, figure_path, events_old, suffix=""):
+
+    xmin, xmax = config["xlim_km"]
+    ymin, ymax = config["ylim_km"]
+    zmin, zmax = config["zlim_km"]
+    vmin, vmax = zmin, zmax
+
+    s = max(0.1, min(10, 5000 / len(events)))
+    alpha = 0.8
+
+    fig, ax = plt.subplots(3, 2, figsize=(10, 10), gridspec_kw={"height_ratios": [2, 1, 1]})
+    im = ax[0, 0].scatter(
+        events_old["x_km"],
+        events_old["y_km"],
+        c=events_old["z_km"],
+        cmap="viridis_r",
+        s=s,
+        marker="o",
+        vmin=vmin,
+        vmax=vmax,
+        alpha=alpha,
+        linewidth=0.0,
+    )
+    ax[0, 0].set_xlim([xmin, xmax])
+    ax[0, 0].set_ylim([ymin, ymax])
+    cbar = fig.colorbar(im, ax=ax[0, 0])
+    cbar.set_label("Depth (km)")
+    ax[0, 0].set_title(f"ADLoc: {len(events_old)} events")
+
+    im = ax[0, 1].scatter(
+        events["x_km"],
+        events["y_km"],
+        c=events["z_km"],
+        cmap="viridis_r",
+        s=s,
+        marker="o",
+        vmin=vmin,
+        vmax=vmax,
+        alpha=alpha,
+        linewidth=0.0,
+    )
+    ax[0, 1].set_xlim([xmin, xmax])
+    ax[0, 1].set_ylim([ymin, ymax])
+    cbar = fig.colorbar(im, ax=ax[0, 1])
+    cbar.set_label("Depth (km)")
+    ax[0, 1].set_title(f"ADLoc DD: {len(events)} events")
+
+    # im = ax[1, 0].scatter(
+    #     events_new["x_km"],
+    #     events_new["z_km"],
+    #     c=events_new["z_km"],
+    #     cmap="viridis_r",
+    #     s=1,
+    #     marker="o",
+    #     vmin=vmin,
+    #     vmax=vmax,
+    # )
+    # ax[1, 0].set_xlim([xmin, xmax])
+    # ax[1, 0].set_ylim([zmax, zmin])
+    # cbar = fig.colorbar(im, ax=ax[1, 0])
+    # cbar.set_label("Depth (km)")
+
+    im = ax[1, 0].scatter(
+        events_old["x_km"],
+        events_old["z_km"],
+        c=events_old["z_km"],
+        cmap="viridis_r",
+        s=s,
+        marker="o",
+        vmin=vmin,
+        vmax=vmax,
+        alpha=alpha,
+        linewidth=0.0,
+    )
+    ax[1, 0].set_xlim([xmin, xmax])
+    ax[1, 0].set_ylim([zmax, zmin])
+    cbar = fig.colorbar(im, ax=ax[1, 0])
+    cbar.set_label("Depth (km)")
+
+    im = ax[1, 1].scatter(
+        events["x_km"],
+        events["z_km"],
+        c=events["z_km"],
+        cmap="viridis_r",
+        s=s,
+        marker="o",
+        vmin=vmin,
+        vmax=vmax,
+        alpha=alpha,
+        linewidth=0.0,
+    )
+    ax[1, 1].set_xlim([xmin, xmax])
+    ax[1, 1].set_ylim([zmax, zmin])
+    cbar = fig.colorbar(im, ax=ax[1, 1])
+    cbar.set_label("Depth (km)")
+
+    im = ax[2, 0].scatter(
+        events_old["y_km"],
+        events_old["z_km"],
+        c=events_old["z_km"],
+        cmap="viridis_r",
+        s=s,
+        marker="o",
+        vmin=vmin,
+        vmax=vmax,
+        alpha=alpha,
+        linewidth=0.0,
+    )
+    ax[2, 0].set_xlim([ymin, ymax])
+    ax[2, 0].set_ylim([zmax, zmin])
+    cbar = fig.colorbar(im, ax=ax[2, 0])
+    cbar.set_label("Depth (km)")
+
+    im = ax[2, 1].scatter(
+        events["y_km"],
+        events["z_km"],
+        c=events["z_km"],
+        cmap="viridis_r",
+        s=s,
+        marker="o",
+        vmin=vmin,
+        vmax=vmax,
+        alpha=alpha,
+        linewidth=0.0,
+    )
+    ax[2, 1].set_xlim([ymin, ymax])
+    ax[2, 1].set_ylim([zmax, zmin])
+    cbar = fig.colorbar(im, ax=ax[2, 1])
+    cbar.set_label("Depth (km)")
+
+    plt.savefig(os.path.join(figure_path, f"location{suffix}.png"), bbox_inches="tight", dpi=300)

examples/california/run_adloc.py

@@ -40,7 +40,7 @@ def run_adloc(
         jday = int(jday.split(".")[1])
 
         # %%
-        data_path = f"{region}/gamma_bo_debug/{year:04d}"
+        data_path = f"{region}/gamma_bo/{year:04d}"
         result_path = f"{region}/adloc/{year:04d}"
         if not os.path.exists(f"{root_path}/{result_path}"):
             os.makedirs(f"{root_path}/{result_path}")
@@ -110,8 +110,10 @@ def run_adloc(
 
         # %%
         stations["depth_km"] = -stations["elevation_m"] / 1000
-        if "station_term_time" not in stations.columns:
-            stations["station_term_time"] = 0.0
+        if "station_term_time_p" not in stations.columns:
+            stations["station_term_time_p"] = 0.0
+        if "station_term_time_s" not in stations.columns:
+            stations["station_term_time_s"] = 0.0
         if "station_term_amplitude" not in stations.columns:
             stations["station_term_amplitude"] = 0.0
         stations[["x_km", "y_km"]] = stations.apply(
@@ -229,16 +231,15 @@ def run_adloc(
         #     suffix=f"_ransac_sst_{iter}",
         # )
 
-        stations["station_term"] = 0.0
-        plotting(
-            stations,
-            f"{root_path}/{figure_path}",
-            config,
-            picks,
-            events_init,
-            events_init,
-            suffix=f"_ransac_sst_init",
-        )
+        # plotting(
+        #     stations,
+        #     f"{root_path}/{figure_path}",
+        #     config,
+        #     picks,
+        #     events_init,
+        #     events_init,
+        #     suffix=f"_ransac_sst_init",
+        # )
 
         for iter in range(MAX_SST_ITER):
             if iter == 0:
@@ -257,7 +258,8 @@ def run_adloc(
                 config["min_score"] = 0.5
                 config["min_p_picks"] = 0  # for filtering
                 config["min_s_picks"] = 0  # for filtering
-            stations["station_term_time"] = 0.0  # no station term
+            stations["station_term_time_p"] = 0.0  # no station term
+            stations["station_term_time_s"] = 0.0  # no station term
             # picks, events = invert_location_iter(picks, stations, config, estimator, events_init=events_init, iter=iter)
             if iter == 0:
                 picks, events = invert_location(picks, stations, config, estimator, events_init=events_init, iter=iter)
@@ -268,31 +270,28 @@ def run_adloc(
                 )
             # station_term = picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_time": "mean"}).reset_index()
             station_term_time = (
-                picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_time": "mean"}).reset_index()
+                picks[picks["mask"] == 1.0]
+                .groupby(["idx_sta", "phase_type"])
+                .agg({"residual_time": "mean"})
+                .reset_index()
             )
+            station_term_time.set_index("idx_sta", inplace=True)
             station_term_amp = (
                 picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_amplitude": "mean"}).reset_index()
             )
-            stations["station_term_time"] += (
-                stations["idx_sta"].map(station_term_time.set_index("idx_sta")["residual_time"]).fillna(0)
+            stations["station_term_time_p"] += (
+                stations["idx_sta"]
+                .map(station_term_time[station_term_time["phase_type"] == 0]["residual_time"])
+                .fillna(0)
+            )
+            stations["station_term_time_s"] += (
+                stations["idx_sta"]
+                .map(station_term_time[station_term_time["phase_type"] == 1]["residual_time"])
+                .fillna(0)
             )
             stations["station_term_amplitude"] += (
                 stations["idx_sta"].map(station_term_amp.set_index("idx_sta")["residual_amplitude"]).fillna(0)
             )
-            ## Separate P and S station term
-            # station_term = (
-            #     picks[picks["mask"] == 1.0].groupby(["idx_sta", "phase_type"]).agg({"residual": "mean"}).reset_index()
-            # )
-            # stations["station_term_p"] = (
-            #     stations["idx_sta"]
-            #     .map(station_term[station_term["phase_type"] == 0].set_index("idx_sta")["residual"])
-            #     .fillna(0)
-            # )
-            # stations["station_term_s"] = (
-            #     stations["idx_sta"]
-            #     .map(station_term[station_term["phase_type"] == 1].set_index("idx_sta")["residual"])
-            #     .fillna(0)
-            # )
 
             # plotting_ransac(
             #     stations,