|
| 1 | +# %% |
| 2 | +import json |
| 3 | +import multiprocessing as mp |
| 4 | +import os |
| 5 | +from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor, as_completed |
| 6 | +from datetime import datetime, timedelta, timezone |
| 7 | +from threading import Lock, Thread |
| 8 | + |
| 9 | +import fsspec |
| 10 | +import numpy as np |
| 11 | +import pandas as pd |
| 12 | +import pyproj |
| 13 | +from obspy import read_inventory |
| 14 | +from obspy.clients.fdsn import Client |
| 15 | +from sklearn.cluster import DBSCAN |
| 16 | +from tqdm import tqdm |
| 17 | +from args import parse_args |
| 18 | +from glob import glob |
| 19 | +import matplotlib.pyplot as plt |
| 20 | +from utils.plotting import plotting_ransac |
| 21 | + |
| 22 | +# %% |
| 23 | +if __name__ == "__main__": |
| 24 | + |
| 25 | + args = parse_args() |
| 26 | + root_path = args.root_path |
| 27 | + region = args.region |
| 28 | + iter = args.iter |
| 29 | + |
| 30 | + data_path = f"{region}/adloc" |
| 31 | + result_path = f"{region}/adloc" |
| 32 | + figure_path = f"{region}/adloc/figures" |
| 33 | + if not os.path.exists(figure_path): |
| 34 | + os.makedirs(figure_path) |
| 35 | + |
| 36 | + # %% |
| 37 | + # protocol = "gs" |
| 38 | + # token_json = f"{os.environ['HOME']}/.config/gcloud/application_default_credentials.json" |
| 39 | + # with open(token_json, "r") as fp: |
| 40 | + # token = json.load(fp) |
| 41 | + # fs = fsspec.filesystem(protocol, token=token) |
| 42 | + |
| 43 | + # %% |
| 44 | + event_csvs = sorted(glob(f"{root_path}/{data_path}/????/????.???.events_sst_{iter}.csv")) |
| 45 | + |
| 46 | + # %% |
| 47 | + events = [] |
| 48 | + picks = [] |
| 49 | + stations = [] |
| 50 | + for event_csv in tqdm(event_csvs, desc="Load event csvs"): |
| 51 | + pick_csv = event_csv.replace(f"events_sst_{iter}.csv", f"picks_sst_{iter}.csv") |
| 52 | + station_csv = event_csv.replace(f"events_sst_{iter}.csv", f"stations_sst_{iter}.csv") |
| 53 | + |
| 54 | + year, jday = event_csv.split("/")[-1].split(".")[:2] |
| 55 | + events_ = pd.read_csv(event_csv, dtype=str) |
| 56 | + picks_ = pd.read_csv(pick_csv, dtype=str) |
| 57 | + stations_ = pd.read_csv(station_csv) |
| 58 | + events_["year"] = year |
| 59 | + events_["jday"] = jday |
| 60 | + picks_["year"] = year |
| 61 | + picks_["jday"] = jday |
| 62 | + stations_["year"] = year |
| 63 | + stations_["jday"] = jday |
| 64 | + events.append(events_) |
| 65 | + picks.append(picks_) |
| 66 | + stations.append(stations_) |
| 67 | + |
| 68 | + events = pd.concat(events, ignore_index=True) |
| 69 | + picks = pd.concat(picks, ignore_index=True) |
| 70 | + stations = pd.concat(stations, ignore_index=True) |
| 71 | + |
| 72 | + station_terms = ( |
| 73 | + stations.groupby(["station_id"]) |
| 74 | + .apply( |
| 75 | + lambda x: pd.Series( |
| 76 | + { |
| 77 | + "station_term_time_p": ( |
| 78 | + (x.station_term_time_p * x.num_pick_p).sum() / x.num_pick_p.sum() |
| 79 | + if x.num_pick_p.sum() > 0 |
| 80 | + else 0 |
| 81 | + ), |
| 82 | + "station_term_time_s": ( |
| 83 | + (x.station_term_time_s * x.num_pick_s).sum() / x.num_pick_s.sum() |
| 84 | + if x.num_pick_s.sum() > 0 |
| 85 | + else 0 |
| 86 | + ), |
| 87 | + "station_term_amplitude": ( |
| 88 | + (x.station_term_amplitude * x.num_pick).sum() / x.num_pick.sum() if x.num_pick.sum() > 0 else 0 |
| 89 | + ), |
| 90 | + } |
| 91 | + ) |
| 92 | + ) |
| 93 | + .reset_index() |
| 94 | + ) |
| 95 | + if iter > 0: |
| 96 | + stations_prev = pd.read_csv(f"{root_path}/{result_path}/adloc_stations_sst_{iter-1}.csv") |
| 97 | + stations_prev.set_index("station_id", inplace=True) |
| 98 | + |
| 99 | + station_terms["station_term_time_p"] += ( |
| 100 | + station_terms["station_id"].map(stations_prev["station_term_time_p"]).fillna(0) |
| 101 | + ) |
| 102 | + station_terms["station_term_time_s"] += ( |
| 103 | + station_terms["station_id"].map(stations_prev["station_term_time_s"]).fillna(0) |
| 104 | + ) |
| 105 | + station_terms["station_term_amplitude"] += ( |
| 106 | + station_terms["station_id"].map(stations_prev["station_term_amplitude"]).fillna(0) |
| 107 | + ) |
| 108 | + |
| 109 | + stations = stations.groupby(["station_id"]).first().reset_index() |
| 110 | + stations.drop(["station_term_time_p", "station_term_time_s", "station_term_amplitude"], axis=1, inplace=True) |
| 111 | + stations = stations.merge(station_terms, on="station_id") |
| 112 | + |
| 113 | + events["dummy_id"] = events["year"] + "." + events["jday"] + "." + events["event_index"] |
| 114 | + picks["dummy_id"] = picks["year"] + "." + picks["jday"] + "." + picks["event_index"] |
| 115 | + |
| 116 | + events["event_index"] = np.arange(len(events)) |
| 117 | + picks = picks.drop("event_index", axis=1) |
| 118 | + picks = picks.merge(events[["dummy_id", "event_index"]], on="dummy_id") |
| 119 | + |
| 120 | + events.drop(["year", "jday", "dummy_id"], axis=1, inplace=True) |
| 121 | + picks.drop(["year", "jday", "dummy_id"], axis=1, inplace=True) |
| 122 | + stations.drop(["year", "jday"], axis=1, inplace=True) |
| 123 | + |
| 124 | + events.to_csv(f"{root_path}/{result_path}/adloc_events_sst_{iter}.csv", index=False) |
| 125 | + picks.to_csv(f"{root_path}/{result_path}/adloc_picks_sst_{iter}.csv", index=False) |
| 126 | + stations.to_csv(f"{root_path}/{result_path}/adloc_stations_sst_{iter}.csv", index=False) |
| 127 | + |
| 128 | + # %% |
| 129 | + |
| 130 | + events = pd.read_csv(f"{root_path}/{result_path}/adloc_events_sst_{iter}.csv") |
| 131 | + picks = pd.read_csv(f"{root_path}/{result_path}/adloc_picks_sst_{iter}.csv") |
| 132 | + stations = pd.read_csv(f"{root_path}/{result_path}/adloc_stations_sst_{iter}.csv") |
| 133 | + |
| 134 | + fig, ax = plt.subplots(3, 3, figsize=(12, 10)) |
| 135 | + ax[0, 0].scatter(events["longitude"], events["latitude"], c=events["depth_km"], s=1, cmap="viridis_r") |
| 136 | + ax[0, 0].set_title(f"Events {len(events)}") |
| 137 | + ax[0, 1].scatter(events["longitude"], events["depth_km"], c=events["depth_km"], s=1, cmap="viridis_r") |
| 138 | + ax[0, 1].invert_yaxis() |
| 139 | + ax[0, 1].set_title(f"Events depth") |
| 140 | + ax[0, 2].scatter(events["latitude"], events["depth_km"], c=events["depth_km"], s=1, cmap="viridis_r") |
| 141 | + ax[0, 2].invert_yaxis() |
| 142 | + ax[0, 2].set_title(f"Events latitude") |
| 143 | + ax[1, 0].scatter( |
| 144 | + stations["longitude"], stations["latitude"], c=stations["station_term_time_p"], marker="^", cmap="viridis_r" |
| 145 | + ) |
| 146 | + ax[1, 0].set_title(f"Station term time P {stations['station_term_time_p'].mean():.2f} s") |
| 147 | + ax[1, 1].scatter( |
| 148 | + stations["longitude"], stations["latitude"], c=stations["station_term_time_s"], marker="^", cmap="viridis_r" |
| 149 | + ) |
| 150 | + ax[1, 1].set_title(f"Station term time S {stations['station_term_time_s'].mean():.2f} s") |
| 151 | + ax[1, 2].scatter( |
| 152 | + stations["longitude"], stations["latitude"], c=stations["station_term_amplitude"], marker="^", cmap="viridis_r" |
| 153 | + ) |
| 154 | + ax[1, 2].set_title(f"Station term amplitude {stations['station_term_amplitude'].mean():.2f} m") |
| 155 | + ax[2, 0].hist(events["adloc_residual_time"], bins=30, edgecolor="white") |
| 156 | + ax[2, 0].set_title(f"Event residual time") |
| 157 | + ax[2, 1].hist(events["adloc_residual_amplitude"], bins=30, edgecolor="white") |
| 158 | + ax[2, 1].set_title(f"Event residual amplitude") |
| 159 | + idx = picks["adloc_mask"] == 1 |
| 160 | + ax[2, 2].hist(picks.loc[idx, "adloc_residual_time"], bins=30, edgecolor="white") |
| 161 | + ax[2, 2].set_title(f"Pick residual time") |
| 162 | + # ax[2, 2].hist(picks["adloc_residual_amplitude"], bins=30, edgecolor="white") |
| 163 | + # ax[2, 2].set_title(f"Pick residual amplitude") |
| 164 | + plt.tight_layout() |
| 165 | + plt.savefig(f"{root_path}/{figure_path}/adloc_summary_{iter}.png") |
| 166 | + plt.close() |
0 commit comments