seperate station term

Author: zhuwq0

scripts/cut_templates_cc.py

@@ -48,7 +48,12 @@ def fillin_missing_picks(picks, events, stations, config):
     picks_ps.drop(columns=["P_source", "S_source"], inplace=True)
 
     picks_ps = picks_ps.merge(events[["event_index", "event_timestamp"]], on="event_index")
-    picks_ps = picks_ps.merge(stations[["station_id", "station_term_time"]], on="station_id")
+    # picks_ps = picks_ps.merge(stations[["station_id", "station_term_time"]], on="station_id")
+    picks_ps = picks_ps.merge(stations[["station_id", "station_term_time_p", "station_term_time_s"]], on="station_id")
+    picks_ps["station_term_time"] = picks_ps.apply(
+        lambda row: row[f"station_term_time_{row['phase_type'].lower()}"], axis=1
+    )  ## Separate P and S station term
+    picks_ps.drop(columns=["station_term_time_p", "station_term_time_s"], inplace=True)
     picks_ps["phase_timestamp"] = picks_ps["event_timestamp"] + picks_ps["traveltime"] + picks_ps["station_term_time"]
     picks_ps["phase_time"] = reference_t0 + pd.to_timedelta(picks_ps["phase_timestamp"], unit="s")
     picks_ps["phase_score"] = min_phase_score
@@ -85,9 +90,16 @@ def predict_full_picks(picks, events, stations, config):
     s_picks = picks_full.assign(phase_type="S")
     picks_full = pd.concat([p_picks, s_picks])
     phase_type = picks_full["phase_type"].values
+    # picks_full = picks_full.merge(
+    #     stations[["station_id", "station_term_time", "x_km", "y_km", "z_km"]], on="station_id"
+    # )
     picks_full = picks_full.merge(
-        stations[["station_id", "station_term_time", "x_km", "y_km", "z_km"]], on="station_id"
+        stations[["station_id", "station_term_time_p", "station_term_time_s", "x_km", "y_km", "z_km"]], on="station_id"
     )
+    picks_full["station_term_time"] = picks_full.apply(
+        lambda row: row[f"station_term_time_{row['phase_type'].lower()}"], axis=1
+    )
+    picks_full.drop(columns=["station_term_time_p", "station_term_time_s"], inplace=True)
     station_dt = picks_full["station_term_time"].values
     station_locs = picks_full[["x_km", "y_km", "z_km"]].values
     picks_full.rename(columns={"x_km": "station_x_km", "y_km": "station_y_km", "z_km": "station_z_km"}, inplace=True)
@@ -379,6 +391,14 @@ def cut_templates(root_path, region, config):
     vp_vs_ratio = 1.73
     vs = [v / vp_vs_ratio for v in vp]
     h = 0.3
+
+    if os.path.exists(f"{root_path}/{region}/obspy/velocity.csv"):
+        velocity = pd.read_csv(f"{root_path}/{region}/obspy/velocity.csv")
+        zz = velocity["z_km"].values
+        vp = velocity["vp"].values
+        vs = velocity["vs"].values
+        h = 0.1
+
     vel = {"Z": zz, "P": vp, "S": vs}
     eikonal = {
         "vel": vel,
@@ -396,9 +416,16 @@ def cut_templates(root_path, region, config):
     min_phase_score = picks["phase_score"].min()
 
     picks = picks.merge(events[["event_index", "event_timestamp"]], on="event_index")
-    picks = picks.merge(stations[["station_id", "station_term_time"]], on="station_id")
+    # picks = picks.merge(stations[["station_id", "station_term_time"]], on="station_id")
+    picks = picks.merge(stations[["station_id", "station_term_time_p", "station_term_time_s"]], on="station_id")
+    picks["station_term_time"] = picks.apply(
+        lambda x: x[f"station_term_time_{x['phase_type'].lower()}"], axis=1
+    )  ## Separate P and S station term
+    picks.drop(columns=["station_term_time_p", "station_term_time_s"], inplace=True)
     picks["phase_timestamp"] = picks["phase_time"].apply(lambda x: (x - reference_t0).total_seconds())
-    picks["traveltime"] = picks["phase_timestamp"] - picks["event_timestamp"] - picks["station_term_time"]
+    picks["traveltime"] = (
+        picks["phase_timestamp"] - picks["event_timestamp"] - picks["station_term_time"]
+    )  ## Separate P and S station term
 
     picks = fillin_missing_picks(
         picks,

scripts/plot_catalog.py

@@ -9,16 +9,9 @@
 import numpy as np
 import pandas as pd
 import plotly.graph_objects as go
-
+from args import parse_args
 
 # %%
-def parse_args():
-    parser = argparse.ArgumentParser()
-    parser.add_argument("root_path", nargs="?", type=str, default="local", help="root path")
-    parser.add_argument("region", nargs="?", type=str, default="demo", help="region")
-    return parser.parse_args()
-
-
 args = parse_args()
 
 # %%
@@ -49,6 +42,7 @@ def parse_args():
 routine_catalog = f"{root_path}/{region}/obspy/catalog.csv"
 routine_exist = False
 if os.path.exists(routine_catalog):
+    print(f"Reading {routine_catalog}")
     routine_exist = True
     routine_catalog = pd.read_csv(routine_catalog, parse_dates=["time"])
 
@@ -57,6 +51,7 @@ def parse_args():
 gamma_file = f"{root_path}/{region}/gamma/gamma_events_000_001.csv"
 gamma_exist = False
 if os.path.exists(gamma_file):
+    print(f"Reading {gamma_file}")
     gamma_exist = True
     gamma_catalog = pd.read_csv(gamma_file, parse_dates=["time"])
     # gamma_catalog["depth_km"] = gamma_catalog["depth(m)"] / 1e3
@@ -66,6 +61,7 @@ def parse_args():
 adloc_file = f"{root_path}/{region}/adloc/ransac_events.csv"
 adloc_exist = False
 if os.path.exists(adloc_file):
+    print(f"Reading {adloc_file}")
     adloc_exist = True
     adloc_catalog = pd.read_csv(adloc_file, parse_dates=["time"])
     # adloc_catalog["magnitude"] = 0.0
@@ -75,20 +71,23 @@ def parse_args():
 adloc_dt_file = f"{root_path}/{region}/adloc_dd/adloc_dt_events.csv"
 adloc_dt_exist = False
 if os.path.exists(adloc_dt_file):
+    print(f"Reading {adloc_dt_file}")
     adloc_dt_exist = True
     adloc_dt_catalog = pd.read_csv(adloc_dt_file, parse_dates=["time"])
 
 # %%
 adloc_dtcc_file = f"{root_path}/{region}/adloc_dd/adloc_dtcc_events.csv"
 adloc_dtcc_exist = False
 if os.path.exists(adloc_dtcc_file):
+    print(f"Reading {adloc_dtcc_file}")
     adloc_dtcc_exist = True
     adloc_dtcc_catalog = pd.read_csv(adloc_dtcc_file, parse_dates=["time"])
 
 # %%
 hypodd_file = f"{root_path}/{region}/hypodd/hypodd_ct.reloc"
 hypodd_ct_exist = False
 if os.path.exists(hypodd_file):
+    print(f"Reading {hypodd_file}")
     hypodd_ct_exist = True
     columns = [
         "ID",
@@ -140,6 +139,7 @@ def parse_args():
 hypodd_file = f"{root_path}/{region}/hypodd/hypodd_cc.reloc"
 hypodd_cc_exist = False
 if os.path.exists(hypodd_file):
+    print(f"Reading {hypodd_file}")
     hypodd_cc_exist = True
     columns = [
         "ID",
@@ -191,6 +191,7 @@ def parse_args():
 growclust_file = f"{root_path}/{region}/growclust/growclust_ct_catalog.txt"
 growclust_ct_exist = False
 if os.path.exists(growclust_file):
+    print(f"Reading {growclust_file}")
     growclust_ct_exist = True
     columns = [
         "yr",
@@ -235,6 +236,7 @@ def parse_args():
 growclust_file = f"{root_path}/{region}/growclust/growclust_cc_catalog.txt"
 growclust_cc_exist = False
 if os.path.exists(growclust_file):
+    print(f"Reading {growclust_file}")
     growclust_cc_exist = True
     columns = [
         "yr",
@@ -746,6 +748,7 @@ def parse_args():
     ax[2, 0].set_title(f"HypoDD (CT): {len(catalog_ct_hypodd)}")
     ax[2, 0].set_xlim(xlim)
     ax[2, 0].set_ylim(ylim)
+
 if hypodd_cc_exist and (len(catalog_cc_hypodd) > 0):
     ax[2, 1].scatter(
         catalog_cc_hypodd["LAT"],
@@ -950,6 +953,25 @@ def plot3d(x, y, z, config, fig_name):
             config_plot3d,
             f"{root_path}/{figure_path}/earthquake_location_adloc.html",
         )
+
+    if adloc_dt_exist and len(adloc_dt_catalog) > 0:
+        plot3d(
+            adloc_dt_catalog["longitude"],
+            adloc_dt_catalog["latitude"],
+            adloc_dt_catalog["depth_km"],
+            config_plot3d,
+            f"{root_path}/{figure_path}/earthquake_location_adloc_dt.html",
+        )
+
+    if adloc_dtcc_exist and len(adloc_dtcc_catalog) > 0:
+        plot3d(
+            adloc_dtcc_catalog["longitude"],
+            adloc_dtcc_catalog["latitude"],
+            adloc_dtcc_catalog["depth_km"],
+            config_plot3d,
+            f"{root_path}/{figure_path}/earthquake_location_adloc_dtcc.html",
+        )
+
     if hypodd_ct_exist and len(catalog_ct_hypodd) > 0:
         plot3d(
             catalog_ct_hypodd["LON"],

scripts/run_adloc.py

@@ -86,8 +86,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(
@@ -117,6 +119,13 @@ def run_adloc(
     # %%
     config["eikonal"] = None
 
+    if os.path.exists(f"{root_path}/{region}/obspy/velocity.csv"):
+        velocity = pd.read_csv(f"{root_path}/{region}/obspy/velocity.csv")
+        zz = velocity["z_km"].values
+        vp = velocity["vp"].values
+        vs = velocity["vs"].values
+        h = 0.1
+
     vel = {"Z": zz, "P": vp, "S": vs}
     config["eikonal"] = {
         "vel": vel,
@@ -194,38 +203,37 @@ def run_adloc(
     for iter in range(MAX_SST_ITER):
         # picks, events = invert_location_iter(picks, stations, config, estimator, events_init=events_init, iter=iter)
         picks, events = invert_location(picks, stations, config, estimator, events_init=events_init, iter=iter)
-        # 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()
-        station_term_time = (
-            picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_time": "median"}).reset_index()
-        )
-        # station_term_amp = (
-        #     picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_amplitude": "mean"}).reset_index()
-        # )
         station_term_amp = (
             picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_amplitude": "median"}).reset_index()
         )
-        stations["station_term_time"] += (
-            stations["idx_sta"].map(station_term_time.set_index("idx_sta")["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)
+
+        ## Same P and S station term
+        # station_term_time = picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_time": "mean"}).reset_index()
+        # stations["station_term_time_p"] += (
+        #     stations["idx_sta"].map(station_term_time.set_index("idx_sta")["residual_time"]).fillna(0)
         # )
-        # stations["station_term_s"] = (
-        #     stations["idx_sta"]
-        #     .map(station_term[station_term["phase_type"] == 1].set_index("idx_sta")["residual"])
-        #     .fillna(0)
+        # stations["station_term_time_s"] += (
+        #     stations["idx_sta"].map(station_term_time.set_index("idx_sta")["residual_time"]).fillna(0)
         # )
 
+        ## Separate P and S station term
+        station_term_time = (
+            picks[picks["mask"] == 1.0].groupby(["idx_sta", "phase_type"]).agg({"residual_time": "mean"}).reset_index()
+        )
+        stations["station_term_time_p"] += (
+            stations["idx_sta"]
+            .map(station_term_time[station_term_time["phase_type"] == 0].set_index("idx_sta")["residual_time"])
+            .fillna(0)
+        )
+        stations["station_term_time_s"] += (
+            stations["idx_sta"]
+            .map(station_term_time[station_term_time["phase_type"] == 1].set_index("idx_sta")["residual_time"])
+            .fillna(0)
+        )
+
         plotting_ransac(stations, figure_path, config, picks, events_init, events, suffix=f"_ransac_sst_{iter}")
 
         if "event_index" not in events.columns: