Earthquake Catalog Workshop¶
Contributors: Eric Beauce, Gabrielle Tepp, Clara Yoon, Ellen Yu, Weiqiang Zhu (alphabetical order)
Machine Learning Part¶
Download data using Obpsy and NCEDC/SCEDC AWS Public Dataset
PhaseNet for P/S phase picking
GaMMA for phase association
ADLoc for earthquake location
HypoDD for earthquake relocation
QuakeFlow
Setup Environment¶
In [ ]:
Copied!
# Run once per session (Colab or a fresh environment)
import os
os.environ["TF_USE_LEGACY_KERAS"] = "1" # PhaseNet uses tf.compat.v1 layers, which Keras 3 removed
!pip install -q obspy cartopy pyproj fsspec s3fs tf-keras git+https://github.com/AI4EPS/GaMMA.git git+https://github.com/AI4EPS/ADLoc.git
!git clone -q https://github.com/AI4EPS/PhaseNet.git QuakeFlow/PhaseNet
!wget -q -nc https://raw.githubusercontent.com/AI4EPS/Earthquake_Catalog_Workshop/main/notebooks/run_hypodd_ct.sh
!which gfortran || apt-get -qq install -y gfortran
# Run once per session (Colab or a fresh environment)
import os
os.environ["TF_USE_LEGACY_KERAS"] = "1" # PhaseNet uses tf.compat.v1 layers, which Keras 3 removed
!pip install -q obspy cartopy pyproj fsspec s3fs tf-keras git+https://github.com/AI4EPS/GaMMA.git git+https://github.com/AI4EPS/ADLoc.git
!git clone -q https://github.com/AI4EPS/PhaseNet.git QuakeFlow/PhaseNet
!wget -q -nc https://raw.githubusercontent.com/AI4EPS/Earthquake_Catalog_Workshop/main/notebooks/run_hypodd_ct.sh
!which gfortran || apt-get -qq install -y gfortran
In [2]:
Copied!
# %%
import json
import os
import pickle
from collections import defaultdict
from datetime import datetime
from typing import Dict
from glob import glob
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import fsspec
import matplotlib.pyplot as plt
import numpy as np
import obspy
import pandas as pd
from obspy.clients.fdsn.mass_downloader import (
CircularDomain,
MassDownloader,
Restrictions,
)
from pyproj import Proj
from tqdm import tqdm
from adloc.eikonal2d import init_eikonal2d
from adloc.sacloc2d import ADLoc
from adloc.utils import invert_location
from gamma.utils import association, estimate_eps
# %%
import json
import os
import pickle
from collections import defaultdict
from datetime import datetime
from typing import Dict
from glob import glob
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import fsspec
import matplotlib.pyplot as plt
import numpy as np
import obspy
import pandas as pd
from obspy.clients.fdsn.mass_downloader import (
CircularDomain,
MassDownloader,
Restrictions,
)
from pyproj import Proj
from tqdm import tqdm
from adloc.eikonal2d import init_eikonal2d
from adloc.sacloc2d import ADLoc
from adloc.utils import invert_location
from gamma.utils import association, estimate_eps
Setup configurations¶
In [3]:
Copied!
def set_config(root_path: str = "local", region: str = "demo") -> Dict:
if not os.path.exists(f"{root_path}/{region}"):
os.makedirs(f"{root_path}/{region}", exist_ok=True)
regions = {
"demo": {
"longitude0": -117.504,
"latitude0": 35.705,
"maxradius_degree": 0.5,
"mindepth": 0,
"maxdepth": 30,
"starttime": "2019-07-04T00:00:00",
"endtime": "2019-07-05T00:00:00",
"network": "CI",
"channel": "HH*,BH*,EH*,HN*",
"provider": [
"SCEDC"
],
},
"ridgecrest": {
"longitude0": -117.504,
"latitude0": 35.705,
"maxradius_degree": 0.5,
"mindepth": 0,
"maxdepth": 30,
"starttime": "2019-07-04T00:00:00",
"endtime": "2019-07-10T00:00:00",
"network": "CI",
"channel": "HH*,BH*,EH*,HN*",
"provider": [
"SCEDC"
],
},
}
## Set config
config = regions[region.lower()]
## PhaseNet
config["phasenet"] = {}
## GaMMA
config["gamma"] = {}
## ADLoc
config["adloc"] = {}
## HypoDD
config["hypodd"] = {}
with open(f"{root_path}/{region}/config.json", "w") as fp:
json.dump(config, fp, indent=2)
print(json.dumps(config, indent=4))
return config
def set_config(root_path: str = "local", region: str = "demo") -> Dict:
if not os.path.exists(f"{root_path}/{region}"):
os.makedirs(f"{root_path}/{region}", exist_ok=True)
regions = {
"demo": {
"longitude0": -117.504,
"latitude0": 35.705,
"maxradius_degree": 0.5,
"mindepth": 0,
"maxdepth": 30,
"starttime": "2019-07-04T00:00:00",
"endtime": "2019-07-05T00:00:00",
"network": "CI",
"channel": "HH*,BH*,EH*,HN*",
"provider": [
"SCEDC"
],
},
"ridgecrest": {
"longitude0": -117.504,
"latitude0": 35.705,
"maxradius_degree": 0.5,
"mindepth": 0,
"maxdepth": 30,
"starttime": "2019-07-04T00:00:00",
"endtime": "2019-07-10T00:00:00",
"network": "CI",
"channel": "HH*,BH*,EH*,HN*",
"provider": [
"SCEDC"
],
},
}
## Set config
config = regions[region.lower()]
## PhaseNet
config["phasenet"] = {}
## GaMMA
config["gamma"] = {}
## ADLoc
config["adloc"] = {}
## HypoDD
config["hypodd"] = {}
with open(f"{root_path}/{region}/config.json", "w") as fp:
json.dump(config, fp, indent=2)
print(json.dumps(config, indent=4))
return config
In [4]:
Copied!
region = "demo"
config = set_config(region = region)
region = "demo"
config = set_config(region = region)
{
"longitude0": -117.504,
"latitude0": 35.705,
"maxradius_degree": 0.5,
"mindepth": 0,
"maxdepth": 30,
"starttime": "2019-07-04T00:00:00",
"endtime": "2019-07-05T00:00:00",
"network": "CI",
"channel": "HH*,BH*,EH*,HN*",
"provider": [
"SCEDC"
],
"phasenet": {},
"gamma": {},
"adloc": {},
"hypodd": {}
}
Download the standard catalog for comparison¶
In [5]:
Copied!
def plot_catalog(catalog: pd.DataFrame, method="Standard", region: str = "demo", config: Dict = {}):
## Plot earthquake locations
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=config["longitude0"], central_latitude=config["latitude0"]))
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE, linewidth=0.2)
scatter = ax.scatter(
catalog['longitude'],
catalog['latitude'],
c=catalog['depth_km'],
cmap='viridis_r',
s=5,
alpha=0.6,
vmin = config["mindepth"],
vmax = config["maxdepth"]/2,
transform=ccrs.PlateCarree()
)
plt.colorbar(scatter, label='Depth (km)')
ax.set_extent([config["longitude0"] - config["maxradius_degree"]/2, config["longitude0"] + config["maxradius_degree"]/2, config["latitude0"] - config["maxradius_degree"]/2, config["latitude0"] + config["maxradius_degree"]/2])
# Add gridlines
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
plt.title(f'{method} Catalog ({len(catalog)})')
fig, axes = plt.subplots(2, 1, figsize=(8, 6))
# Histogram of magnitudes
axes[0].hist(catalog['magnitude'], bins=20, color='C0', edgecolor='white')
axes[0].set_xlabel('Magnitude')
axes[0].set_ylabel('Frequency')
axes[0].set_yscale('log')
# axes[0].set_title(f'{method} Catalog ({len(catalog)})')
axes[0].grid(linestyle='--')
# Scatter plot of magnitudes over time
axes[1].scatter(
pd.to_datetime(catalog['time']),
catalog['magnitude'],
s=2**catalog['magnitude'],
c="C0"
)
axes[1].set_xlim(pd.to_datetime(config["starttime"]), pd.to_datetime(config["endtime"]))
axes[1].set_xlabel('Time')
axes[1].set_ylabel('Magnitude')
# axes[1].set_title(f'{method} Catalog ({len(catalog)})')
axes[1].tick_params(axis='x', rotation=45)
axes[1].grid(linestyle='--')
plt.tight_layout()
plt.show()
def plot_stations(stations: pd.DataFrame, catalog: pd.DataFrame=None, region: str = "demo", config: Dict = {}):
## Plot earthquake locations
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=config["longitude0"], central_latitude=config["latitude0"]))
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE, linewidth=0.2)
scatter = ax.scatter(
stations['longitude'],
stations['latitude'],
c="C0",
s=40,
marker="^",
alpha=0.6,
transform=ccrs.PlateCarree()
)
if catalog is not None:
scatter = ax.scatter(
catalog['longitude'],
catalog['latitude'],
c=catalog['depth_km'],
cmap='viridis_r',
s=0.5,
alpha=0.6,
vmin = config["mindepth"],
vmax = config["maxdepth"]/2,
transform=ccrs.PlateCarree()
)
plt.colorbar(scatter, label='Depth (km)')
ax.set_extent([config["longitude0"] - config["maxradius_degree"], config["longitude0"] + config["maxradius_degree"], config["latitude0"] - config["maxradius_degree"], config["latitude0"] + config["maxradius_degree"]])
# Add gridlines
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
plt.title(f'Stations ({len(stations["station"].unique())})')
plt.show()
def download_catalog(root_path: str = "local", region: str = "demo", config: Dict = {}):
result_path = f"{region}/obspy"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
# print(json.dumps(config, indent=4))
## Download catalog
client = obspy.clients.fdsn.Client("usgs")
events = client.get_events(
starttime=config["starttime"],
endtime=config["endtime"],
latitude=config["latitude0"],
longitude=config["longitude0"],
maxradius=config["maxradius_degree"],
)
print(f"Number of events: {len(events)}")
## Save catalog
catalog = defaultdict(list)
for event in events:
if len(event.magnitudes) > 0:
catalog["time"].append(event.origins[0].time.datetime.strftime("%Y-%m-%dT%H:%M:%S.%f"))
catalog["magnitude"].append(event.magnitudes[0].mag)
catalog["longitude"].append(event.origins[0].longitude)
catalog["latitude"].append(event.origins[0].latitude)
catalog["depth_km"].append(event.origins[0].depth/1e3)
catalog = pd.DataFrame.from_dict(catalog).sort_values(["time"])
catalog.to_csv(f"{root_path}/{result_path}/catalog.csv", index=False)
return catalog
def plot_catalog(catalog: pd.DataFrame, method="Standard", region: str = "demo", config: Dict = {}):
## Plot earthquake locations
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=config["longitude0"], central_latitude=config["latitude0"]))
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE, linewidth=0.2)
scatter = ax.scatter(
catalog['longitude'],
catalog['latitude'],
c=catalog['depth_km'],
cmap='viridis_r',
s=5,
alpha=0.6,
vmin = config["mindepth"],
vmax = config["maxdepth"]/2,
transform=ccrs.PlateCarree()
)
plt.colorbar(scatter, label='Depth (km)')
ax.set_extent([config["longitude0"] - config["maxradius_degree"]/2, config["longitude0"] + config["maxradius_degree"]/2, config["latitude0"] - config["maxradius_degree"]/2, config["latitude0"] + config["maxradius_degree"]/2])
# Add gridlines
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
plt.title(f'{method} Catalog ({len(catalog)})')
fig, axes = plt.subplots(2, 1, figsize=(8, 6))
# Histogram of magnitudes
axes[0].hist(catalog['magnitude'], bins=20, color='C0', edgecolor='white')
axes[0].set_xlabel('Magnitude')
axes[0].set_ylabel('Frequency')
axes[0].set_yscale('log')
# axes[0].set_title(f'{method} Catalog ({len(catalog)})')
axes[0].grid(linestyle='--')
# Scatter plot of magnitudes over time
axes[1].scatter(
pd.to_datetime(catalog['time']),
catalog['magnitude'],
s=2**catalog['magnitude'],
c="C0"
)
axes[1].set_xlim(pd.to_datetime(config["starttime"]), pd.to_datetime(config["endtime"]))
axes[1].set_xlabel('Time')
axes[1].set_ylabel('Magnitude')
# axes[1].set_title(f'{method} Catalog ({len(catalog)})')
axes[1].tick_params(axis='x', rotation=45)
axes[1].grid(linestyle='--')
plt.tight_layout()
plt.show()
def plot_stations(stations: pd.DataFrame, catalog: pd.DataFrame=None, region: str = "demo", config: Dict = {}):
## Plot earthquake locations
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=config["longitude0"], central_latitude=config["latitude0"]))
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE, linewidth=0.2)
scatter = ax.scatter(
stations['longitude'],
stations['latitude'],
c="C0",
s=40,
marker="^",
alpha=0.6,
transform=ccrs.PlateCarree()
)
if catalog is not None:
scatter = ax.scatter(
catalog['longitude'],
catalog['latitude'],
c=catalog['depth_km'],
cmap='viridis_r',
s=0.5,
alpha=0.6,
vmin = config["mindepth"],
vmax = config["maxdepth"]/2,
transform=ccrs.PlateCarree()
)
plt.colorbar(scatter, label='Depth (km)')
ax.set_extent([config["longitude0"] - config["maxradius_degree"], config["longitude0"] + config["maxradius_degree"], config["latitude0"] - config["maxradius_degree"], config["latitude0"] + config["maxradius_degree"]])
# Add gridlines
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
plt.title(f'Stations ({len(stations["station"].unique())})')
plt.show()
def download_catalog(root_path: str = "local", region: str = "demo", config: Dict = {}):
result_path = f"{region}/obspy"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
# print(json.dumps(config, indent=4))
## Download catalog
client = obspy.clients.fdsn.Client("usgs")
events = client.get_events(
starttime=config["starttime"],
endtime=config["endtime"],
latitude=config["latitude0"],
longitude=config["longitude0"],
maxradius=config["maxradius_degree"],
)
print(f"Number of events: {len(events)}")
## Save catalog
catalog = defaultdict(list)
for event in events:
if len(event.magnitudes) > 0:
catalog["time"].append(event.origins[0].time.datetime.strftime("%Y-%m-%dT%H:%M:%S.%f"))
catalog["magnitude"].append(event.magnitudes[0].mag)
catalog["longitude"].append(event.origins[0].longitude)
catalog["latitude"].append(event.origins[0].latitude)
catalog["depth_km"].append(event.origins[0].depth/1e3)
catalog = pd.DataFrame.from_dict(catalog).sort_values(["time"])
catalog.to_csv(f"{root_path}/{result_path}/catalog.csv", index=False)
return catalog
In [6]:
Copied!
standard_catalog = download_catalog(region=region, config=config)
plot_catalog(standard_catalog, method="Standard", config=config)
standard_catalog = download_catalog(region=region, config=config)
plot_catalog(standard_catalog, method="Standard", config=config)
Number of events: 709
Download stations¶
In [7]:
Copied!
def download_station(root_path: str = "local", region: str = "demo", config: Dict = {}):
result_dir = f"{region}/obspy"
if not os.path.exists(f"{root_path}/{result_dir}"):
os.makedirs(f"{root_path}/{result_dir}")
if not os.path.exists(f"{root_path}/{result_dir}/inventory/"):
os.makedirs(f"{root_path}/{result_dir}/inventory/")
## Download stations
stations = obspy.core.inventory.Inventory()
for provider in config["provider"]:
client = obspy.clients.fdsn.Client(provider)
stations += client.get_stations(
network=config["network"],
station="*",
starttime=config["starttime"],
endtime=config["endtime"],
latitude=config["latitude0"],
longitude=config["longitude0"],
maxradius=config["maxradius_degree"],
channel=config["channel"],
level="response",
)
stations.write(f"{root_path}/{result_dir}/inventory.xml", format="STATIONXML")
print("Number of stations: {}".format(sum([len(x) for x in stations])))
## Save stations
station_dict = defaultdict(dict)
for network in stations:
for station in network:
inv = stations.select(network=network.code, station=station.code)
inv.write(f"{root_path}/{result_dir}/inventory/{network.code}.{station.code}.xml", format="STATIONXML")
for channel in station:
sid = f"{network.code}.{station.code}.{channel.location_code}.{channel.code}"
station_dict[sid] = {
"network": network.code,
"station": station.code,
"location": channel.location_code,
"channel": channel.code,
"longitude": channel.longitude,
"latitude": channel.latitude,
"elevation_m": channel.elevation,
"response": round(channel.response.instrument_sensitivity.value, 2),
}
# with open(f"{root_path}/{result_dir}/stations.json", "w") as fp:
# json.dump(station_dict, fp, indent=2)
with open(f"{root_path}/{result_dir}/stations.pkl", "wb") as fp:
pickle.dump(stations, fp)
stations = pd.DataFrame.from_dict(station_dict, orient="index")
stations.to_csv(f"{root_path}/{result_dir}/stations.csv", index=False)
return stations
def download_station(root_path: str = "local", region: str = "demo", config: Dict = {}):
result_dir = f"{region}/obspy"
if not os.path.exists(f"{root_path}/{result_dir}"):
os.makedirs(f"{root_path}/{result_dir}")
if not os.path.exists(f"{root_path}/{result_dir}/inventory/"):
os.makedirs(f"{root_path}/{result_dir}/inventory/")
## Download stations
stations = obspy.core.inventory.Inventory()
for provider in config["provider"]:
client = obspy.clients.fdsn.Client(provider)
stations += client.get_stations(
network=config["network"],
station="*",
starttime=config["starttime"],
endtime=config["endtime"],
latitude=config["latitude0"],
longitude=config["longitude0"],
maxradius=config["maxradius_degree"],
channel=config["channel"],
level="response",
)
stations.write(f"{root_path}/{result_dir}/inventory.xml", format="STATIONXML")
print("Number of stations: {}".format(sum([len(x) for x in stations])))
## Save stations
station_dict = defaultdict(dict)
for network in stations:
for station in network:
inv = stations.select(network=network.code, station=station.code)
inv.write(f"{root_path}/{result_dir}/inventory/{network.code}.{station.code}.xml", format="STATIONXML")
for channel in station:
sid = f"{network.code}.{station.code}.{channel.location_code}.{channel.code}"
station_dict[sid] = {
"network": network.code,
"station": station.code,
"location": channel.location_code,
"channel": channel.code,
"longitude": channel.longitude,
"latitude": channel.latitude,
"elevation_m": channel.elevation,
"response": round(channel.response.instrument_sensitivity.value, 2),
}
# with open(f"{root_path}/{result_dir}/stations.json", "w") as fp:
# json.dump(station_dict, fp, indent=2)
with open(f"{root_path}/{result_dir}/stations.pkl", "wb") as fp:
pickle.dump(stations, fp)
stations = pd.DataFrame.from_dict(station_dict, orient="index")
stations.to_csv(f"{root_path}/{result_dir}/stations.csv", index=False)
return stations
In [8]:
Copied!
stations = download_station(region=region, config=config)
plot_stations(stations, catalog = standard_catalog, region=region, config=config)
stations = download_station(region=region, config=config)
plot_stations(stations, catalog = standard_catalog, region=region, config=config)
Number of stations: 15
Download waveform data¶
In [9]:
Copied!
def map_remote_path(provider, bucket, starttime, network, station, location, channel):
starttime = pd.Timestamp(starttime).round("h").to_pydatetime() # in case of 2021-01-01T23:59:xxx
if provider.lower() == "scedc":
year = starttime.strftime("%Y")
dayofyear = starttime.strftime("%j")
if location == "":
location = "__"
path = f"s3://{bucket}/{year}/{year}_{dayofyear}/{network}{station:_<5}{channel}{location:_<2}_{year}{dayofyear}.ms"
elif provider.lower() == "ncedc":
year = starttime.strftime("%Y")
dayofyear = starttime.strftime("%j")
path = f"s3://{bucket}/{network}/{year}/{year}.{dayofyear}/{station}.{network}.{channel}.{location}.D.{year}.{dayofyear}"
else:
raise ValueError(f"Unknown provider: {provider}")
return path
def download_waveform(root_path: str = "local", region: str = "demo", config: Dict = {}):
waveform_dir = f"{region}/waveforms"
if not os.path.exists(f"{root_path}/{waveform_dir}"):
os.makedirs(f"{root_path}/{waveform_dir}")
# print(json.dumps(config, indent=4))
## Download from cloud
for provider in config["provider"]:
if provider.lower() in ["scedc", "ncedc"]:
cloud = {
"provider": provider.lower(),
"bucket": f"{provider.lower()}-pds/continuous_waveforms",
}
else:
continue
DELTATIME = "1D"
starttime = datetime.fromisoformat(config["starttime"]).strftime("%Y-%m-%d")
starttimes = pd.date_range(starttime, config["endtime"], freq=DELTATIME, tz="UTC", inclusive="left").to_list()
# with open(f'{root_path}/{region}/obspy/stations.json', 'r') as f:
# stations = json.load(f)
stations = pd.read_csv(f"{root_path}/{region}/obspy/stations.csv", na_filter=False)
stations["instrument"] = stations['channel'].apply(lambda x: x[:-1])
stations["component"] = stations['channel'].apply(lambda x: x[-1])
location_priorities = ('', '00', '10', '01', '20', '02', '30', '03', '40', '04',
'50', '05', '60', '06', '70', '07', '80', '08', '90', '09')
location_priority_map = {loc: i for i, loc in enumerate(location_priorities)}
instrument_priorities = ('HH', 'BH', 'MH', 'EH', 'LH', 'HL', 'BL', 'ML', 'EL', 'LL', 'SH')
instrument_priority_map = {ch: i for i, ch in enumerate(instrument_priorities)}
component_priorities = ('E', 'N', 'Z', '1', '2', '3')
component_priority_map = {ch: i for i, ch in enumerate(component_priorities)}
stations['location_priority'] = stations['location'].map(location_priority_map)
stations['instrument_priority'] = stations['instrument'].apply(lambda x: instrument_priority_map.get(x, len(instrument_priorities)))
stations['component_priority'] = stations['component'].apply(lambda x: component_priority_map.get(x, len(component_priorities)))
stations.sort_values(['network', 'station', 'location_priority', 'instrument_priority', "component_priority"], inplace=True)
stations.drop(['location_priority', 'instrument_priority', 'component_priority'], axis=1, inplace=True)
for starttime in starttimes:
for _, station in tqdm(stations.groupby(["network", "station"]), desc=f"Downloading {starttime}"):
prev = ""
nch = set()
for _, row in station.iterrows():
if len(nch) >= 3:
break
network, station, location, channel, instrument = row["network"], row["station"], row["location"], row["channel"], row["instrument"]
if instrument != prev:
prev = instrument
nch = set()
mseed_path = map_remote_path(
cloud["provider"],
cloud["bucket"],
starttime,
network,
station,
location,
channel,
)
try:
if os.path.exists(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed"):
# print(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed downloaded.")
nch.add(channel[-1])
continue
if not os.path.exists(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}"):
os.makedirs(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}")
with fsspec.open(f"{mseed_path}", "rb", s3={"anon": True}) as f:
data = f.read()
with open(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed", "wb") as f:
f.write(data)
nch.add(channel[-1])
except Exception as e:
# print(f"Failed to download {e}")
pass
# %% Download from FDSN
domain = CircularDomain(
longitude=config["longitude0"],
latitude=config["latitude0"],
minradius=0,
maxradius=config["maxradius_degree"],
)
restrictions = Restrictions(
starttime=obspy.UTCDateTime(config["starttime"]),
endtime=obspy.UTCDateTime(config["endtime"]),
chunklength_in_sec=3600 * 24, # 1 day
network=config["network"] if "network" in config else None,
station=config["station"] if "station" in config else None,
minimum_interstation_distance_in_m=0,
minimum_length=0.1,
reject_channels_with_gaps=False,
)
def get_mseed_storage(network, station, location, channel, starttime, endtime):
mseed_name = f"{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed"
if os.path.exists(f"{root_path}/{waveform_dir}/{mseed_name}"):
# print(f"{root_path}/{waveform_dir}/{mseed_name} downloaded.")
return True
return f"{root_path}/{waveform_dir}/{mseed_name}"
mdl = MassDownloader(
providers=config["provider"],
# providers=["IRIS"],
)
mdl.download(
domain,
restrictions,
mseed_storage=get_mseed_storage,
stationxml_storage=f"{root_path}/{waveform_dir}/stations",
download_chunk_size_in_mb=20,
threads_per_client=3,
print_report=False,
)
return
def map_remote_path(provider, bucket, starttime, network, station, location, channel):
starttime = pd.Timestamp(starttime).round("h").to_pydatetime() # in case of 2021-01-01T23:59:xxx
if provider.lower() == "scedc":
year = starttime.strftime("%Y")
dayofyear = starttime.strftime("%j")
if location == "":
location = "__"
path = f"s3://{bucket}/{year}/{year}_{dayofyear}/{network}{station:_<5}{channel}{location:_<2}_{year}{dayofyear}.ms"
elif provider.lower() == "ncedc":
year = starttime.strftime("%Y")
dayofyear = starttime.strftime("%j")
path = f"s3://{bucket}/{network}/{year}/{year}.{dayofyear}/{station}.{network}.{channel}.{location}.D.{year}.{dayofyear}"
else:
raise ValueError(f"Unknown provider: {provider}")
return path
def download_waveform(root_path: str = "local", region: str = "demo", config: Dict = {}):
waveform_dir = f"{region}/waveforms"
if not os.path.exists(f"{root_path}/{waveform_dir}"):
os.makedirs(f"{root_path}/{waveform_dir}")
# print(json.dumps(config, indent=4))
## Download from cloud
for provider in config["provider"]:
if provider.lower() in ["scedc", "ncedc"]:
cloud = {
"provider": provider.lower(),
"bucket": f"{provider.lower()}-pds/continuous_waveforms",
}
else:
continue
DELTATIME = "1D"
starttime = datetime.fromisoformat(config["starttime"]).strftime("%Y-%m-%d")
starttimes = pd.date_range(starttime, config["endtime"], freq=DELTATIME, tz="UTC", inclusive="left").to_list()
# with open(f'{root_path}/{region}/obspy/stations.json', 'r') as f:
# stations = json.load(f)
stations = pd.read_csv(f"{root_path}/{region}/obspy/stations.csv", na_filter=False)
stations["instrument"] = stations['channel'].apply(lambda x: x[:-1])
stations["component"] = stations['channel'].apply(lambda x: x[-1])
location_priorities = ('', '00', '10', '01', '20', '02', '30', '03', '40', '04',
'50', '05', '60', '06', '70', '07', '80', '08', '90', '09')
location_priority_map = {loc: i for i, loc in enumerate(location_priorities)}
instrument_priorities = ('HH', 'BH', 'MH', 'EH', 'LH', 'HL', 'BL', 'ML', 'EL', 'LL', 'SH')
instrument_priority_map = {ch: i for i, ch in enumerate(instrument_priorities)}
component_priorities = ('E', 'N', 'Z', '1', '2', '3')
component_priority_map = {ch: i for i, ch in enumerate(component_priorities)}
stations['location_priority'] = stations['location'].map(location_priority_map)
stations['instrument_priority'] = stations['instrument'].apply(lambda x: instrument_priority_map.get(x, len(instrument_priorities)))
stations['component_priority'] = stations['component'].apply(lambda x: component_priority_map.get(x, len(component_priorities)))
stations.sort_values(['network', 'station', 'location_priority', 'instrument_priority', "component_priority"], inplace=True)
stations.drop(['location_priority', 'instrument_priority', 'component_priority'], axis=1, inplace=True)
for starttime in starttimes:
for _, station in tqdm(stations.groupby(["network", "station"]), desc=f"Downloading {starttime}"):
prev = ""
nch = set()
for _, row in station.iterrows():
if len(nch) >= 3:
break
network, station, location, channel, instrument = row["network"], row["station"], row["location"], row["channel"], row["instrument"]
if instrument != prev:
prev = instrument
nch = set()
mseed_path = map_remote_path(
cloud["provider"],
cloud["bucket"],
starttime,
network,
station,
location,
channel,
)
try:
if os.path.exists(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed"):
# print(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed downloaded.")
nch.add(channel[-1])
continue
if not os.path.exists(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}"):
os.makedirs(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}")
with fsspec.open(f"{mseed_path}", "rb", s3={"anon": True}) as f:
data = f.read()
with open(f"{root_path}/{waveform_dir}/{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed", "wb") as f:
f.write(data)
nch.add(channel[-1])
except Exception as e:
# print(f"Failed to download {e}")
pass
# %% Download from FDSN
domain = CircularDomain(
longitude=config["longitude0"],
latitude=config["latitude0"],
minradius=0,
maxradius=config["maxradius_degree"],
)
restrictions = Restrictions(
starttime=obspy.UTCDateTime(config["starttime"]),
endtime=obspy.UTCDateTime(config["endtime"]),
chunklength_in_sec=3600 * 24, # 1 day
network=config["network"] if "network" in config else None,
station=config["station"] if "station" in config else None,
minimum_interstation_distance_in_m=0,
minimum_length=0.1,
reject_channels_with_gaps=False,
)
def get_mseed_storage(network, station, location, channel, starttime, endtime):
mseed_name = f"{starttime.strftime('%Y/%j')}/{network}.{station}.{location}.{channel}.mseed"
if os.path.exists(f"{root_path}/{waveform_dir}/{mseed_name}"):
# print(f"{root_path}/{waveform_dir}/{mseed_name} downloaded.")
return True
return f"{root_path}/{waveform_dir}/{mseed_name}"
mdl = MassDownloader(
providers=config["provider"],
# providers=["IRIS"],
)
mdl.download(
domain,
restrictions,
mseed_storage=get_mseed_storage,
stationxml_storage=f"{root_path}/{waveform_dir}/stations",
download_chunk_size_in_mb=20,
threads_per_client=3,
print_report=False,
)
return
In [10]:
Copied!
download_waveform(region=region, config=config)
download_waveform(region=region, config=config)
Downloading 2019-07-04 00:00:00+00:00: 0%| | 0/15 [00:00<?, ?it/s]
Downloading 2019-07-04 00:00:00+00:00: 40%|████ | 6/15 [00:00<00:00, 12.27it/s]
Downloading 2019-07-04 00:00:00+00:00: 100%|██████████| 15/15 [00:00<00:00, 30.29it/s]
[2026-10-19 22:18:07,151] - obspy.clients.fdsn.mass_downloader - INFO: Initializing FDSN client(s) for SCEDC.
[2026-10-19 22:18:07,155] - obspy.clients.fdsn.mass_downloader - INFO: Successfully initialized 1 client(s): SCEDC.
[2026-10-19 22:18:07,156] - obspy.clients.fdsn.mass_downloader - INFO: Total acquired or preexisting stations: 0
[2026-10-19 22:18:07,156] - obspy.clients.fdsn.mass_downloader - INFO: Client 'SCEDC' - Requesting unreliable availability.
[2026-10-19 22:18:12,149] - obspy.clients.fdsn.mass_downloader - INFO: Client 'SCEDC' - Successfully requested availability (4.99 seconds)
[2026-10-19 22:18:12,160] - obspy.clients.fdsn.mass_downloader - INFO: Client 'SCEDC' - Found 14 stations (36 channels).
[2026-10-19 22:18:12,160] - obspy.clients.fdsn.mass_downloader - INFO: Client 'SCEDC' - Will attempt to download data from 14 stations.
[2026-10-19 22:18:12,161] - obspy.clients.fdsn.mass_downloader - INFO: Client 'SCEDC' - Status for 36 time intervals/channels before downloading: IGNORE
[2026-10-19 22:18:12,161] - obspy.clients.fdsn.mass_downloader - INFO: Client 'SCEDC' - No station information to download.
Run PhaseNet to pick P/S picks¶
In [11]:
Copied!
def run_phasenet(root_path: str = "local", region: str = "demo", config: Dict = {} ) -> str:
result_path = f"{region}/phasenet"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
# %%
waveform_dir = f"{region}/waveforms"
mseed_list = sorted(glob(f"{root_path}/{waveform_dir}/????/???/*.mseed"))
# %% group 3C channels
mseed_list = sorted(list(set([x.split(".mseed")[0][:-1] + "*.mseed" for x in mseed_list])))
# %%
with open(f"{root_path}/{result_path}/mseed_list.csv", "w") as fp:
fp.write("fname\n")
fp.write("\n".join(mseed_list))
# %%
model_path = "QuakeFlow/PhaseNet/"
cmd = f"python {model_path}/phasenet/predict.py --model={model_path}/model/190703-214543 --data_dir=./ --data_list={root_path}/{result_path}/mseed_list.csv --response_xml={root_path}/{region}/obspy/inventory.xml --format=mseed --amplitude --highpass_filter=1.0 --result_dir={root_path}/{result_path} --result_fname=phasenet_picks --batch_size=1"
# cmd += " --sampling_rate 100"
os.system(cmd)
return f"{root_path}/{result_path}/phasenet_picks.csv"
def run_phasenet(root_path: str = "local", region: str = "demo", config: Dict = {} ) -> str:
result_path = f"{region}/phasenet"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
# %%
waveform_dir = f"{region}/waveforms"
mseed_list = sorted(glob(f"{root_path}/{waveform_dir}/????/???/*.mseed"))
# %% group 3C channels
mseed_list = sorted(list(set([x.split(".mseed")[0][:-1] + "*.mseed" for x in mseed_list])))
# %%
with open(f"{root_path}/{result_path}/mseed_list.csv", "w") as fp:
fp.write("fname\n")
fp.write("\n".join(mseed_list))
# %%
model_path = "QuakeFlow/PhaseNet/"
cmd = f"python {model_path}/phasenet/predict.py --model={model_path}/model/190703-214543 --data_dir=./ --data_list={root_path}/{result_path}/mseed_list.csv --response_xml={root_path}/{region}/obspy/inventory.xml --format=mseed --amplitude --highpass_filter=1.0 --result_dir={root_path}/{result_path} --result_fname=phasenet_picks --batch_size=1"
# cmd += " --sampling_rate 100"
os.system(cmd)
return f"{root_path}/{result_path}/phasenet_picks.csv"
In [12]:
Copied!
phasenet_picks = run_phasenet(region=region, config=config)
phasenet_picks = run_phasenet(region=region, config=config)
2026-10-19 22:18:14,375 Pred log: local/demo/phasenet 2026-10-19 22:18:14,375 Dataset size: 17 2026-10-19 22:18:14,404 Model: depths 5, filters 8, filter size 7x1, pool size: 4x1, dilation rate: 1x1
WARNING: All log messages before absl::InitializeLog() is called are written to STDERR I0000 00:00:1792473494.753865 3636693 mlir_graph_optimization_pass.cc:437] MLIR V1 optimization pass is not enabled 2026-10-19 22:18:14,807 restoring model QuakeFlow/PhaseNet//model/190703-214543/model_95.ckpt Pred: 0%| | 0/17 [00:00<?, ?it/s]
Pred: 6%|▌ | 1/17 [00:02<00:45, 2.84s/it]
Pred: 12%|█▏ | 2/17 [00:04<00:32, 2.15s/it]
Pred: 18%|█▊ | 3/17 [00:06<00:27, 1.98s/it]
Pred: 24%|██▎ | 4/17 [00:08<00:24, 1.92s/it]
Pred: 29%|██▉ | 5/17 [00:09<00:22, 1.85s/it]
Pred: 35%|███▌ | 6/17 [00:11<00:19, 1.80s/it]
Pred: 41%|████ | 7/17 [00:13<00:18, 1.83s/it]
Pred: 47%|████▋ | 8/17 [00:15<00:16, 1.83s/it]
Pred: 53%|█████▎ | 9/17 [00:16<00:14, 1.78s/it]
Pred: 59%|█████▉ | 10/17 [00:18<00:12, 1.75s/it]
Pred: 65%|██████▍ | 11/17 [00:20<00:10, 1.68s/it]
Pred: 71%|███████ | 12/17 [00:21<00:08, 1.70s/it]
Pred: 76%|███████▋ | 13/17 [00:23<00:06, 1.70s/it]
Pred: 82%|████████▏ | 14/17 [00:25<00:04, 1.65s/it]
Pred: 88%|████████▊ | 15/17 [00:26<00:03, 1.67s/it]
Pred: 94%|█████████▍| 16/17 [00:28<00:01, 1.63s/it]
Pred: 100%|██████████| 17/17 [00:30<00:00, 1.69s/it] Pred: 100%|██████████| 17/17 [00:30<00:00, 1.78s/it]
Done with 14088 P-picks and 14637 S-picks
Run GaMMA to associate P/S picks¶
In [13]:
Copied!
def run_gamma(root_path: str = "local", region: str = "demo", config: Dict = {}):
data_path = f"{region}/phasenet"
result_path = f"{region}/gamma"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
picks_csv = f"{data_path}/phasenet_picks.csv"
gamma_events_csv = f"{result_path}/gamma_events.csv"
gamma_picks_csv = f"{result_path}/gamma_picks.csv"
station_json = f"{region}/obspy/stations.json"
## read picks
picks = pd.read_csv(f"{root_path}/{picks_csv}")
picks.drop(columns=["event_index"], inplace=True, errors="ignore")
picks["id"] = picks["station_id"]
picks["timestamp"] = picks["phase_time"]
picks["amp"] = picks["phase_amplitude"]
picks["type"] = picks["phase_type"]
picks["prob"] = picks["phase_score"]
## read stations
# stations = pd.read_json(f"{root_path}/{station_json}", orient="index")
stations = pd.read_csv(f"{root_path}/{region}/obspy/stations.csv", na_filter=False)
stations["id"] = stations.apply(lambda x: f"{x['network']}.{x['station']}.{x['location']}.{x['channel'][:-1]}", axis=1)
stations = stations.groupby("id").agg(lambda x: x.iloc[0] if len(set(x)) == 1 else sorted(list(x))).reset_index()
proj = Proj(f"+proj=aeqd +lon_0={config['longitude0']} +lat_0={config['latitude0']} +units=km")
stations[["x(km)", "y(km)"]] = stations.apply(
lambda x: pd.Series(proj(longitude=x.longitude, latitude=x.latitude)), axis=1
)
stations["z(km)"] = stations["elevation_m"].apply(lambda x: -x / 1e3)
# print(stations.to_string())
## setting GaMMA configs
config["use_dbscan"] = True
config["use_amplitude"] = True
config["method"] = "BGMM"
if config["method"] == "BGMM": ## BayesianGaussianMixture
config["oversample_factor"] = 5
if config["method"] == "GMM": ## GaussianMixture
config["oversample_factor"] = 1
# earthquake location
config["vel"] = {"p": 6.0, "s": 6.0 / 1.75}
config["dims"] = ["x(km)", "y(km)", "z(km)"]
minlat, maxlat = config["latitude0"] - config["maxradius_degree"], config["latitude0"] + config["maxradius_degree"]
minlon, maxlon = config["longitude0"] - config["maxradius_degree"], config["longitude0"] + config["maxradius_degree"]
xmin, ymin = proj(minlon, minlat)
xmax, ymax = proj(maxlon, maxlat)
# zmin, zmax = config["mindepth"], config["maxdepth"]
zmin = config["mindepth"] if "mindepth" in config else 0
zmax = config["maxdepth"] if "maxdepth" in config else 30
config["x(km)"] = (xmin, xmax)
config["y(km)"] = (ymin, ymax)
config["z(km)"] = (zmin, zmax)
config["bfgs_bounds"] = (
(config["x(km)"][0] - 1, config["x(km)"][1] + 1), # x
(config["y(km)"][0] - 1, config["y(km)"][1] + 1), # y
(0, config["z(km)"][1] + 1), # z
(None, None), # t
)
# DBSCAN
config["dbscan_eps"] = estimate_eps(stations, config["vel"]["p"]) # s
config["dbscan_min_samples"] = 3
## Eikonal for 1D velocity model
zz = [0.0, 5.5, 5.5, 16.0, 16.0, 32.0, 32.0]
vp = [5.5, 5.5, 6.3, 6.3, 6.7, 6.7, 7.8]
vp_vs_ratio = 1.73
vs = [v / vp_vs_ratio for v in vp]
h = 0.3
vel = {"z": zz, "p": vp, "s": vs}
config["eikonal"] = {"vel": vel, "h": h, "xlim": config["x(km)"], "ylim": config["y(km)"], "zlim": config["z(km)"]}
# filtering
config["min_picks_per_eq"] = 5
config["min_p_picks_per_eq"] = 0
config["min_s_picks_per_eq"] = 0
config["max_sigma11"] = 2.0 # s
config["max_sigma22"] = 1.0 # log10(m/s)
config["max_sigma12"] = 1.0 # covariance
## filter picks without amplitude measurements
if config["use_amplitude"]:
picks = picks[picks["amp"] != -1]
# for k, v in config.items():
# print(f"{k}: {v}")
print(f"Number of picks: {len(picks)}")
# %%
event_idx0 = 0 ## current earthquake index
assignments = []
events, assignments = association(picks, stations, config, event_idx0, config["method"])
if len(events) == 0:
return
## create catalog
events = pd.DataFrame(events)
events[["longitude", "latitude"]] = events.apply(
lambda x: pd.Series(proj(longitude=x["x(km)"], latitude=x["y(km)"], inverse=True)), axis=1
)
events["depth_km"] = events["z(km)"]
events.sort_values("time", inplace=True)
with open(f"{root_path}/{gamma_events_csv}", "w") as fp:
events.to_csv(fp, index=False, float_format="%.3f", date_format="%Y-%m-%dT%H:%M:%S.%f")
## add assignment to picks
assignments = pd.DataFrame(assignments, columns=["pick_index", "event_index", "gamma_score"])
picks = picks.join(assignments.set_index("pick_index")).fillna(-1).astype({"event_index": int})
picks.sort_values(["phase_time"], inplace=True)
with open(f"{root_path}/{gamma_picks_csv}", "w") as fp:
picks.to_csv(fp, index=False, date_format="%Y-%m-%dT%H:%M:%S.%f")
# return f"{root_path}/{result_path}/gamma_picks.csv", f"{root_path}/{result_path}/gamma_events.csv"
return events
def run_gamma(root_path: str = "local", region: str = "demo", config: Dict = {}):
data_path = f"{region}/phasenet"
result_path = f"{region}/gamma"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
picks_csv = f"{data_path}/phasenet_picks.csv"
gamma_events_csv = f"{result_path}/gamma_events.csv"
gamma_picks_csv = f"{result_path}/gamma_picks.csv"
station_json = f"{region}/obspy/stations.json"
## read picks
picks = pd.read_csv(f"{root_path}/{picks_csv}")
picks.drop(columns=["event_index"], inplace=True, errors="ignore")
picks["id"] = picks["station_id"]
picks["timestamp"] = picks["phase_time"]
picks["amp"] = picks["phase_amplitude"]
picks["type"] = picks["phase_type"]
picks["prob"] = picks["phase_score"]
## read stations
# stations = pd.read_json(f"{root_path}/{station_json}", orient="index")
stations = pd.read_csv(f"{root_path}/{region}/obspy/stations.csv", na_filter=False)
stations["id"] = stations.apply(lambda x: f"{x['network']}.{x['station']}.{x['location']}.{x['channel'][:-1]}", axis=1)
stations = stations.groupby("id").agg(lambda x: x.iloc[0] if len(set(x)) == 1 else sorted(list(x))).reset_index()
proj = Proj(f"+proj=aeqd +lon_0={config['longitude0']} +lat_0={config['latitude0']} +units=km")
stations[["x(km)", "y(km)"]] = stations.apply(
lambda x: pd.Series(proj(longitude=x.longitude, latitude=x.latitude)), axis=1
)
stations["z(km)"] = stations["elevation_m"].apply(lambda x: -x / 1e3)
# print(stations.to_string())
## setting GaMMA configs
config["use_dbscan"] = True
config["use_amplitude"] = True
config["method"] = "BGMM"
if config["method"] == "BGMM": ## BayesianGaussianMixture
config["oversample_factor"] = 5
if config["method"] == "GMM": ## GaussianMixture
config["oversample_factor"] = 1
# earthquake location
config["vel"] = {"p": 6.0, "s": 6.0 / 1.75}
config["dims"] = ["x(km)", "y(km)", "z(km)"]
minlat, maxlat = config["latitude0"] - config["maxradius_degree"], config["latitude0"] + config["maxradius_degree"]
minlon, maxlon = config["longitude0"] - config["maxradius_degree"], config["longitude0"] + config["maxradius_degree"]
xmin, ymin = proj(minlon, minlat)
xmax, ymax = proj(maxlon, maxlat)
# zmin, zmax = config["mindepth"], config["maxdepth"]
zmin = config["mindepth"] if "mindepth" in config else 0
zmax = config["maxdepth"] if "maxdepth" in config else 30
config["x(km)"] = (xmin, xmax)
config["y(km)"] = (ymin, ymax)
config["z(km)"] = (zmin, zmax)
config["bfgs_bounds"] = (
(config["x(km)"][0] - 1, config["x(km)"][1] + 1), # x
(config["y(km)"][0] - 1, config["y(km)"][1] + 1), # y
(0, config["z(km)"][1] + 1), # z
(None, None), # t
)
# DBSCAN
config["dbscan_eps"] = estimate_eps(stations, config["vel"]["p"]) # s
config["dbscan_min_samples"] = 3
## Eikonal for 1D velocity model
zz = [0.0, 5.5, 5.5, 16.0, 16.0, 32.0, 32.0]
vp = [5.5, 5.5, 6.3, 6.3, 6.7, 6.7, 7.8]
vp_vs_ratio = 1.73
vs = [v / vp_vs_ratio for v in vp]
h = 0.3
vel = {"z": zz, "p": vp, "s": vs}
config["eikonal"] = {"vel": vel, "h": h, "xlim": config["x(km)"], "ylim": config["y(km)"], "zlim": config["z(km)"]}
# filtering
config["min_picks_per_eq"] = 5
config["min_p_picks_per_eq"] = 0
config["min_s_picks_per_eq"] = 0
config["max_sigma11"] = 2.0 # s
config["max_sigma22"] = 1.0 # log10(m/s)
config["max_sigma12"] = 1.0 # covariance
## filter picks without amplitude measurements
if config["use_amplitude"]:
picks = picks[picks["amp"] != -1]
# for k, v in config.items():
# print(f"{k}: {v}")
print(f"Number of picks: {len(picks)}")
# %%
event_idx0 = 0 ## current earthquake index
assignments = []
events, assignments = association(picks, stations, config, event_idx0, config["method"])
if len(events) == 0:
return
## create catalog
events = pd.DataFrame(events)
events[["longitude", "latitude"]] = events.apply(
lambda x: pd.Series(proj(longitude=x["x(km)"], latitude=x["y(km)"], inverse=True)), axis=1
)
events["depth_km"] = events["z(km)"]
events.sort_values("time", inplace=True)
with open(f"{root_path}/{gamma_events_csv}", "w") as fp:
events.to_csv(fp, index=False, float_format="%.3f", date_format="%Y-%m-%dT%H:%M:%S.%f")
## add assignment to picks
assignments = pd.DataFrame(assignments, columns=["pick_index", "event_index", "gamma_score"])
picks = picks.join(assignments.set_index("pick_index")).fillna(-1).astype({"event_index": int})
picks.sort_values(["phase_time"], inplace=True)
with open(f"{root_path}/{gamma_picks_csv}", "w") as fp:
picks.to_csv(fp, index=False, date_format="%Y-%m-%dT%H:%M:%S.%f")
# return f"{root_path}/{result_path}/gamma_picks.csv", f"{root_path}/{result_path}/gamma_events.csv"
return events
In [14]:
Copied!
gamma_catalog = run_gamma(region=region, config=config)
gamma_catalog = run_gamma(region=region, config=config)
Number of picks: 28725 Eikonal Solver:
Iter 0, error = 999.945 Iter 1, error = 0.000 Time: 0.903 Eikonal Solver: Iter 0, error = 999.906 Iter 1, error = 0.000 Time: 0.022
Clustering (eps=6.19): 0%| | 0/638 [00:00<?, ?it/s]
Clustering (eps=6.19): 100%|██████████| 638/638 [00:00<00:00, 113600.18it/s]
Associating 637 clusters with 17 CPUs
/var/folders/0n/_fgybwj14vxfg5mgnsyxy3f40000gn/T/ipykernel_63161/3463100311.py:112: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
picks = picks.join(assignments.set_index("pick_index")).fillna(-1).astype({"event_index": int})
In [15]:
Copied!
plot_catalog(gamma_catalog, "GaMMA", region=region, config=config)
plot_catalog(gamma_catalog, "GaMMA", region=region, config=config)
Run ADLoc to locate absolute earthquake locations¶
In [16]:
Copied!
def run_adloc(root_path: str = "local", region: str = "demo", config: Dict = {}):
data_path = f"{root_path}/{region}/gamma"
result_path = f"{root_path}/{region}/adloc"
figure_path = f"{root_path}/{region}/adloc/figures"
if not os.path.exists(result_path):
os.makedirs(result_path)
if not os.path.exists(figure_path):
os.makedirs(figure_path)
picks_file = f"{data_path}/gamma_picks.csv"
events_file = f"{data_path}/gamma_events.csv"
stations_file = f"{root_path}/{region}/obspy/stations.csv"
proj = Proj(f"+proj=aeqd +lon_0={config['longitude0']} +lat_0={config['latitude0']} +units=km")
config["proj"] = proj # ADLoc >= 1.0.5 locates in this frame, from the GaMMA locations; without it the catalogue shifts
## read picks and associated events
picks = pd.read_csv(picks_file)
picks["phase_time"] = pd.to_datetime(picks["phase_time"])
# drop unnecessary columns
picks.drop(["id", "timestamp", "type", "amp", "prob", "event_idx"], axis=1, inplace=True, errors="ignore")
if os.path.exists(events_file):
events = pd.read_csv(events_file)
events["time"] = pd.to_datetime(events["time"])
events[["x_km", "y_km"]] = events.apply(
lambda x: pd.Series(proj(longitude=x.longitude, latitude=x.latitude)), axis=1
)
events["z_km"] = events["depth_km"] if "depth_km" in events.columns else 10.0
else:
events = None
## read stations
# stations = pd.read_json(stations_file, orient="index")
stations = pd.read_csv(stations_file, na_filter=False)
stations["station_id"] = stations.apply(lambda x: f"{x['network']}.{x['station']}.{x['location']}.{x['channel'][:-1]}", axis=1)
stations = stations.groupby("station_id").agg(lambda x: x.iloc[0] if len(set(x)) == 1 else sorted(list(x))).reset_index()
stations["depth_km"] = -stations["elevation_m"] / 1000
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(
lambda x: pd.Series(proj(longitude=x.longitude, latitude=x.latitude)), axis=1
)
stations["z_km"] = stations["elevation_m"].apply(lambda x: -x / 1e3)
## setting ADLoc configs
config["use_amplitude"] = True
minlat, maxlat = config["latitude0"] - config["maxradius_degree"], config["latitude0"] + config["maxradius_degree"]
minlon, maxlon = config["longitude0"] - config["maxradius_degree"], config["longitude0"] + config["maxradius_degree"]
xmin, ymin = proj(minlon, minlat)
xmax, ymax = proj(maxlon, maxlat)
zmin, zmax = config["mindepth"], config["maxdepth"]
config["xlim_km"] = (xmin, xmax)
config["ylim_km"] = (ymin, ymax)
config["zlim_km"] = (zmin, zmax)
## Eikonal for 1D velocity model
zz = [0.0, 5.5, 5.5, 16.0, 16.0, 32.0, 32.0]
vp = [5.5, 5.5, 6.3, 6.3, 6.7, 6.7, 7.8]
vp_vs_ratio = 1.73
vs = [v / vp_vs_ratio for v in vp]
# Northern California (Gil7)
# zz = [0.0, 1.0, 3.0, 4.0, 5.0, 17.0, 25.0, 62.0]
# vp = [3.2, 3.2, 4.5, 4.8, 5.51, 6.21, 6.89, 7.83]
# vs = [1.5, 1.5, 2.4, 2.78, 3.18, 3.40, 3.98, 4.52]
h = 0.3
vel = {"Z": zz, "P": vp, "S": vs}
config["eikonal"] = {
"vel": vel,
"h": h,
"xlim_km": config["xlim_km"],
"ylim_km": config["ylim_km"],
"zlim_km": config["zlim_km"],
}
config["eikonal"] = init_eikonal2d(config["eikonal"])
# RASAC
config["min_picks"] = 6
config["min_picks_ratio"] = 0.5
config["max_residual_time"] = 1.0
config["max_residual_amplitude"] = 1.0
config["min_score"] = 0.5
config["min_s_picks"] = 1.5
config["min_p_picks"] = 1.5
config["bfgs_bounds"] = (
(config["xlim_km"][0] - 1, config["xlim_km"][1] + 1), # x
(config["ylim_km"][0] - 1, config["ylim_km"][1] + 1), # y
(0, config["zlim_km"][1] + 1),
(None, None), # t
)
# %%
mapping_phase_type_int = {"P": 0, "S": 1}
picks["phase_type"] = picks["phase_type"].map(mapping_phase_type_int)
if "phase_amplitude" in picks.columns:
picks["phase_amplitude"] = picks["phase_amplitude"].apply(lambda x: np.log10(x) + 2.0) # convert to log10(cm/s)
# %%
stations["idx_sta"] = np.arange(len(stations))
if events is None:
picks = picks.merge(stations[["station_id", "x_km", "y_km", "z_km"]], on="station_id")
events = picks.groupby("event_index").agg({"x_km": "mean", "y_km": "mean", "z_km": "mean", "phase_time": "min"})
events["z_km"] = 10.0 # km default depth
events.rename({"phase_time": "time"}, axis=1, inplace=True)
events["event_index"] = events.index
events.reset_index(drop=True, inplace=True)
events["idx_eve"] = np.arange(len(events))
picks.drop(["x_km", "y_km", "z_km"], axis=1, inplace=True)
else:
events["idx_eve"] = np.arange(len(events))
picks = picks.merge(events[["event_index", "idx_eve"]], on="event_index")
picks = picks.merge(stations[["station_id", "idx_sta"]], on="station_id")
# for key, value in config.items():
# print(f"{key}: {value}")
print(f"Number of picks: {len(picks)}")
print(f"Number of events: {len(events)}")
# %%
estimator = ADLoc(config, eikonal=config["eikonal"])
# %%
MAX_SST_ITER = 8
events_init = events.copy()
for iter in range(MAX_SST_ITER):
picks, events = invert_location(picks, stations, config, estimator, events_init=events_init, iter=iter)
station_term_amp = (
picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_amplitude": "median"}).reset_index()
)
station_term_amp.set_index("idx_sta", inplace=True)
stations["station_term_amplitude"] += stations["idx_sta"].map(station_term_amp["residual_amplitude"]).fillna(0)
station_term_time = (
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)
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)
)
if "event_index" not in events.columns:
events["event_index"] = events.merge(picks[["idx_eve", "event_index"]], on="idx_eve")["event_index"]
events[["longitude", "latitude"]] = events.apply(
lambda x: pd.Series(proj(x["x_km"], x["y_km"], inverse=True)), axis=1
)
events["depth_km"] = events["z_km"]
picks["adloc_mask"] = picks["mask"]
picks["adloc_residual_time"] = picks["residual_time"]
picks["adloc_residual_amplitude"] = picks["residual_amplitude"]
picks.to_csv(os.path.join(result_path, f"adloc_picks_sst_{iter}.csv"), index=False)
events.to_csv(os.path.join(result_path, f"adloc_events_sst_{iter}.csv"), index=False)
stations.to_csv(os.path.join(result_path, f"adloc_stations_sst_{iter}.csv"), index=False)
# %%
if "event_index" not in events.columns:
events["event_index"] = events.merge(picks[["idx_eve", "event_index"]], on="idx_eve")["event_index"]
events[["longitude", "latitude"]] = events.apply(
lambda x: pd.Series(proj(x["x_km"], x["y_km"], inverse=True)), axis=1
)
events["depth_km"] = events["z_km"]
events.drop(["idx_eve", "x_km", "y_km", "z_km"], axis=1, inplace=True, errors="ignore")
events.sort_values(["time"], inplace=True)
picks["phase_type"] = picks["phase_type"].map({0: "P", 1: "S"})
picks.drop(
["idx_eve", "idx_sta", "mask", "residual_time", "residual_amplitude"], axis=1, inplace=True, errors="ignore"
)
picks.sort_values(["phase_time"], inplace=True)
stations.drop(["idx_sta", "x_km", "y_km", "z_km"], axis=1, inplace=True, errors="ignore")
picks.to_csv(os.path.join(result_path, "adloc_picks.csv"), index=False)
events.to_csv(os.path.join(result_path, "adloc_events.csv"), index=False)
stations.to_csv(os.path.join(result_path, "adloc_stations.csv"), index=False)
return events
def run_adloc(root_path: str = "local", region: str = "demo", config: Dict = {}):
data_path = f"{root_path}/{region}/gamma"
result_path = f"{root_path}/{region}/adloc"
figure_path = f"{root_path}/{region}/adloc/figures"
if not os.path.exists(result_path):
os.makedirs(result_path)
if not os.path.exists(figure_path):
os.makedirs(figure_path)
picks_file = f"{data_path}/gamma_picks.csv"
events_file = f"{data_path}/gamma_events.csv"
stations_file = f"{root_path}/{region}/obspy/stations.csv"
proj = Proj(f"+proj=aeqd +lon_0={config['longitude0']} +lat_0={config['latitude0']} +units=km")
config["proj"] = proj # ADLoc >= 1.0.5 locates in this frame, from the GaMMA locations; without it the catalogue shifts
## read picks and associated events
picks = pd.read_csv(picks_file)
picks["phase_time"] = pd.to_datetime(picks["phase_time"])
# drop unnecessary columns
picks.drop(["id", "timestamp", "type", "amp", "prob", "event_idx"], axis=1, inplace=True, errors="ignore")
if os.path.exists(events_file):
events = pd.read_csv(events_file)
events["time"] = pd.to_datetime(events["time"])
events[["x_km", "y_km"]] = events.apply(
lambda x: pd.Series(proj(longitude=x.longitude, latitude=x.latitude)), axis=1
)
events["z_km"] = events["depth_km"] if "depth_km" in events.columns else 10.0
else:
events = None
## read stations
# stations = pd.read_json(stations_file, orient="index")
stations = pd.read_csv(stations_file, na_filter=False)
stations["station_id"] = stations.apply(lambda x: f"{x['network']}.{x['station']}.{x['location']}.{x['channel'][:-1]}", axis=1)
stations = stations.groupby("station_id").agg(lambda x: x.iloc[0] if len(set(x)) == 1 else sorted(list(x))).reset_index()
stations["depth_km"] = -stations["elevation_m"] / 1000
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(
lambda x: pd.Series(proj(longitude=x.longitude, latitude=x.latitude)), axis=1
)
stations["z_km"] = stations["elevation_m"].apply(lambda x: -x / 1e3)
## setting ADLoc configs
config["use_amplitude"] = True
minlat, maxlat = config["latitude0"] - config["maxradius_degree"], config["latitude0"] + config["maxradius_degree"]
minlon, maxlon = config["longitude0"] - config["maxradius_degree"], config["longitude0"] + config["maxradius_degree"]
xmin, ymin = proj(minlon, minlat)
xmax, ymax = proj(maxlon, maxlat)
zmin, zmax = config["mindepth"], config["maxdepth"]
config["xlim_km"] = (xmin, xmax)
config["ylim_km"] = (ymin, ymax)
config["zlim_km"] = (zmin, zmax)
## Eikonal for 1D velocity model
zz = [0.0, 5.5, 5.5, 16.0, 16.0, 32.0, 32.0]
vp = [5.5, 5.5, 6.3, 6.3, 6.7, 6.7, 7.8]
vp_vs_ratio = 1.73
vs = [v / vp_vs_ratio for v in vp]
# Northern California (Gil7)
# zz = [0.0, 1.0, 3.0, 4.0, 5.0, 17.0, 25.0, 62.0]
# vp = [3.2, 3.2, 4.5, 4.8, 5.51, 6.21, 6.89, 7.83]
# vs = [1.5, 1.5, 2.4, 2.78, 3.18, 3.40, 3.98, 4.52]
h = 0.3
vel = {"Z": zz, "P": vp, "S": vs}
config["eikonal"] = {
"vel": vel,
"h": h,
"xlim_km": config["xlim_km"],
"ylim_km": config["ylim_km"],
"zlim_km": config["zlim_km"],
}
config["eikonal"] = init_eikonal2d(config["eikonal"])
# RASAC
config["min_picks"] = 6
config["min_picks_ratio"] = 0.5
config["max_residual_time"] = 1.0
config["max_residual_amplitude"] = 1.0
config["min_score"] = 0.5
config["min_s_picks"] = 1.5
config["min_p_picks"] = 1.5
config["bfgs_bounds"] = (
(config["xlim_km"][0] - 1, config["xlim_km"][1] + 1), # x
(config["ylim_km"][0] - 1, config["ylim_km"][1] + 1), # y
(0, config["zlim_km"][1] + 1),
(None, None), # t
)
# %%
mapping_phase_type_int = {"P": 0, "S": 1}
picks["phase_type"] = picks["phase_type"].map(mapping_phase_type_int)
if "phase_amplitude" in picks.columns:
picks["phase_amplitude"] = picks["phase_amplitude"].apply(lambda x: np.log10(x) + 2.0) # convert to log10(cm/s)
# %%
stations["idx_sta"] = np.arange(len(stations))
if events is None:
picks = picks.merge(stations[["station_id", "x_km", "y_km", "z_km"]], on="station_id")
events = picks.groupby("event_index").agg({"x_km": "mean", "y_km": "mean", "z_km": "mean", "phase_time": "min"})
events["z_km"] = 10.0 # km default depth
events.rename({"phase_time": "time"}, axis=1, inplace=True)
events["event_index"] = events.index
events.reset_index(drop=True, inplace=True)
events["idx_eve"] = np.arange(len(events))
picks.drop(["x_km", "y_km", "z_km"], axis=1, inplace=True)
else:
events["idx_eve"] = np.arange(len(events))
picks = picks.merge(events[["event_index", "idx_eve"]], on="event_index")
picks = picks.merge(stations[["station_id", "idx_sta"]], on="station_id")
# for key, value in config.items():
# print(f"{key}: {value}")
print(f"Number of picks: {len(picks)}")
print(f"Number of events: {len(events)}")
# %%
estimator = ADLoc(config, eikonal=config["eikonal"])
# %%
MAX_SST_ITER = 8
events_init = events.copy()
for iter in range(MAX_SST_ITER):
picks, events = invert_location(picks, stations, config, estimator, events_init=events_init, iter=iter)
station_term_amp = (
picks[picks["mask"] == 1.0].groupby("idx_sta").agg({"residual_amplitude": "median"}).reset_index()
)
station_term_amp.set_index("idx_sta", inplace=True)
stations["station_term_amplitude"] += stations["idx_sta"].map(station_term_amp["residual_amplitude"]).fillna(0)
station_term_time = (
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)
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)
)
if "event_index" not in events.columns:
events["event_index"] = events.merge(picks[["idx_eve", "event_index"]], on="idx_eve")["event_index"]
events[["longitude", "latitude"]] = events.apply(
lambda x: pd.Series(proj(x["x_km"], x["y_km"], inverse=True)), axis=1
)
events["depth_km"] = events["z_km"]
picks["adloc_mask"] = picks["mask"]
picks["adloc_residual_time"] = picks["residual_time"]
picks["adloc_residual_amplitude"] = picks["residual_amplitude"]
picks.to_csv(os.path.join(result_path, f"adloc_picks_sst_{iter}.csv"), index=False)
events.to_csv(os.path.join(result_path, f"adloc_events_sst_{iter}.csv"), index=False)
stations.to_csv(os.path.join(result_path, f"adloc_stations_sst_{iter}.csv"), index=False)
# %%
if "event_index" not in events.columns:
events["event_index"] = events.merge(picks[["idx_eve", "event_index"]], on="idx_eve")["event_index"]
events[["longitude", "latitude"]] = events.apply(
lambda x: pd.Series(proj(x["x_km"], x["y_km"], inverse=True)), axis=1
)
events["depth_km"] = events["z_km"]
events.drop(["idx_eve", "x_km", "y_km", "z_km"], axis=1, inplace=True, errors="ignore")
events.sort_values(["time"], inplace=True)
picks["phase_type"] = picks["phase_type"].map({0: "P", 1: "S"})
picks.drop(
["idx_eve", "idx_sta", "mask", "residual_time", "residual_amplitude"], axis=1, inplace=True, errors="ignore"
)
picks.sort_values(["phase_time"], inplace=True)
stations.drop(["idx_sta", "x_km", "y_km", "z_km"], axis=1, inplace=True, errors="ignore")
picks.to_csv(os.path.join(result_path, "adloc_picks.csv"), index=False)
events.to_csv(os.path.join(result_path, "adloc_events.csv"), index=False)
stations.to_csv(os.path.join(result_path, "adloc_stations.csv"), index=False)
return events
In [17]:
Copied!
adloc_catalog = run_adloc(region=region, config=config)
adloc_catalog = run_adloc(region=region, config=config)
Eikonal Solver:
Iter 0, error = 999.945 Iter 1, error = 0.000 Time: 0.641 Eikonal Solver: Iter 0, error = 999.906 Iter 1, error = 0.000 Time: 0.027 Number of picks: 24681 Number of events: 1274 Using 17 CPUs for ADLoc inversion.
Iter 0: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 0: 0%| | 1/1274 [00:01<37:40, 1.78s/it]
Iter 0: 0%| | 3/1274 [00:02<14:26, 1.47it/s]
Iter 0: 0%| | 6/1274 [00:03<08:26, 2.50it/s]
Iter 0: 1%|▏ | 17/1274 [00:03<02:08, 9.79it/s]
Iter 0: 2%|▏ | 29/1274 [00:03<01:03, 19.51it/s]
Iter 0: 4%|▎ | 45/1274 [00:03<00:34, 35.23it/s]
Iter 0: 5%|▌ | 65/1274 [00:03<00:20, 57.71it/s]
Iter 0: 7%|▋ | 87/1274 [00:03<00:14, 83.95it/s]
Iter 0: 9%|▉ | 112/1274 [00:03<00:10, 114.40it/s]
Iter 0: 11%|█ | 140/1274 [00:03<00:07, 148.85it/s]
Iter 0: 13%|█▎ | 164/1274 [00:03<00:06, 166.17it/s]
Iter 0: 15%|█▍ | 186/1274 [00:03<00:06, 172.14it/s]
Iter 0: 16%|█▋ | 210/1274 [00:04<00:05, 187.84it/s]
Iter 0: 18%|█▊ | 232/1274 [00:04<00:05, 185.16it/s]
Iter 0: 20%|█▉ | 253/1274 [00:04<00:05, 182.08it/s]
Iter 0: 22%|██▏ | 281/1274 [00:04<00:04, 204.95it/s]
Iter 0: 24%|██▍ | 303/1274 [00:04<00:04, 206.03it/s]
Iter 0: 26%|██▌ | 325/1274 [00:04<00:04, 204.01it/s]
Iter 0: 28%|██▊ | 351/1274 [00:04<00:04, 216.83it/s]
Iter 0: 30%|██▉ | 381/1274 [00:04<00:03, 236.44it/s]
Iter 0: 32%|███▏ | 407/1274 [00:04<00:03, 242.53it/s]
Iter 0: 34%|███▍ | 432/1274 [00:05<00:03, 239.94it/s]
Iter 0: 36%|███▌ | 457/1274 [00:05<00:03, 218.10it/s]
Iter 0: 38%|███▊ | 483/1274 [00:05<00:03, 229.35it/s]
Iter 0: 40%|███▉ | 507/1274 [00:05<00:03, 230.66it/s]
Iter 0: 42%|████▏ | 535/1274 [00:05<00:03, 243.54it/s]
Iter 0: 45%|████▍ | 569/1274 [00:05<00:02, 269.95it/s]
Iter 0: 47%|████▋ | 597/1274 [00:05<00:02, 265.64it/s]
Iter 0: 49%|████▉ | 628/1274 [00:05<00:02, 277.85it/s]
Iter 0: 52%|█████▏ | 659/1274 [00:05<00:02, 283.34it/s]
Iter 0: 54%|█████▍ | 688/1274 [00:06<00:02, 282.31it/s]
Iter 0: 57%|█████▋ | 720/1274 [00:06<00:01, 293.13it/s]
Iter 0: 59%|█████▉ | 750/1274 [00:06<00:01, 285.01it/s]
Iter 0: 61%|██████▏ | 783/1274 [00:06<00:01, 295.26it/s]
Iter 0: 64%|██████▍ | 813/1274 [00:06<00:01, 282.00it/s]
Iter 0: 66%|██████▌ | 842/1274 [00:06<00:01, 267.84it/s]
Iter 0: 68%|██████▊ | 870/1274 [00:06<00:01, 267.28it/s]
Iter 0: 70%|███████ | 897/1274 [00:06<00:01, 260.71it/s]
Iter 0: 73%|███████▎ | 925/1274 [00:06<00:01, 264.54it/s]
Iter 0: 75%|███████▍ | 952/1274 [00:07<00:01, 253.87it/s]
Iter 0: 77%|███████▋ | 983/1274 [00:07<00:01, 268.68it/s]
Iter 0: 79%|███████▉ | 1011/1274 [00:07<00:01, 243.64it/s]
Iter 0: 81%|████████▏ | 1037/1274 [00:07<00:00, 240.82it/s]
Iter 0: 83%|████████▎ | 1062/1274 [00:07<00:01, 195.56it/s]
Iter 0: 85%|████████▌ | 1084/1274 [00:07<00:00, 194.01it/s]
Iter 0: 87%|████████▋ | 1111/1274 [00:07<00:00, 208.94it/s]
Iter 0: 89%|████████▉ | 1133/1274 [00:07<00:00, 208.38it/s]
Iter 0: 91%|█████████ | 1155/1274 [00:08<00:00, 205.37it/s]
Iter 0: 93%|█████████▎| 1180/1274 [00:08<00:00, 217.02it/s]
Iter 0: 95%|█████████▍| 1204/1274 [00:08<00:00, 222.65it/s]
Iter 0: 96%|█████████▋| 1227/1274 [00:08<00:00, 221.04it/s]
Iter 0: 98%|█████████▊| 1250/1274 [00:08<00:00, 194.25it/s]
Iter 0: 100%|█████████▉| 1271/1274 [00:08<00:00, 180.64it/s]
Iter 0: 100%|██████████| 1274/1274 [00:09<00:00, 133.70it/s]
ADLoc locates 1049 events outof 1274 events using 22642 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 1: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 1: 0%| | 1/1274 [00:01<41:57, 1.98s/it]
Iter 1: 0%| | 3/1274 [00:02<12:43, 1.67it/s]
Iter 1: 0%| | 4/1274 [00:02<09:18, 2.28it/s]
Iter 1: 0%| | 5/1274 [00:02<07:34, 2.79it/s]
Iter 1: 0%| | 6/1274 [00:03<09:26, 2.24it/s]
Iter 1: 1%| | 12/1274 [00:03<02:46, 7.60it/s]
Iter 1: 2%|▏ | 29/1274 [00:03<00:46, 26.58it/s]
Iter 1: 4%|▍ | 52/1274 [00:03<00:22, 55.47it/s]
Iter 1: 6%|▌ | 77/1274 [00:03<00:13, 88.19it/s]
Iter 1: 9%|▊ | 111/1274 [00:03<00:08, 136.36it/s]
Iter 1: 11%|█ | 139/1274 [00:03<00:07, 161.77it/s]
Iter 1: 13%|█▎ | 169/1274 [00:03<00:05, 191.45it/s]
Iter 1: 16%|█▌ | 204/1274 [00:04<00:04, 226.93it/s]
Iter 1: 18%|█▊ | 234/1274 [00:04<00:04, 244.95it/s]
Iter 1: 21%|██ | 263/1274 [00:04<00:04, 250.72it/s]
Iter 1: 23%|██▎ | 292/1274 [00:04<00:03, 259.16it/s]
Iter 1: 25%|██▌ | 320/1274 [00:04<00:03, 263.51it/s]
Iter 1: 27%|██▋ | 348/1274 [00:04<00:03, 247.76it/s]
Iter 1: 29%|██▉ | 375/1274 [00:04<00:03, 252.12it/s]
Iter 1: 32%|███▏ | 402/1274 [00:04<00:03, 255.42it/s]
Iter 1: 34%|███▎ | 429/1274 [00:04<00:03, 216.48it/s]
Iter 1: 36%|███▌ | 459/1274 [00:05<00:03, 235.23it/s]
Iter 1: 38%|███▊ | 486/1274 [00:05<00:03, 239.76it/s]
Iter 1: 40%|████ | 515/1274 [00:05<00:03, 247.91it/s]
Iter 1: 43%|████▎ | 542/1274 [00:05<00:02, 253.44it/s]
Iter 1: 45%|████▍ | 568/1274 [00:05<00:02, 239.90it/s]
Iter 1: 47%|████▋ | 593/1274 [00:05<00:03, 208.44it/s]
Iter 1: 48%|████▊ | 617/1274 [00:05<00:03, 214.31it/s]
Iter 1: 51%|█████ | 651/1274 [00:05<00:02, 246.24it/s]
Iter 1: 54%|█████▎ | 682/1274 [00:05<00:02, 263.47it/s]
Iter 1: 56%|█████▌ | 711/1274 [00:06<00:02, 269.46it/s]
Iter 1: 58%|█████▊ | 744/1274 [00:06<00:01, 282.86it/s]
Iter 1: 61%|██████ | 778/1274 [00:06<00:01, 294.33it/s]
Iter 1: 64%|██████▎ | 809/1274 [00:06<00:01, 298.10it/s]
Iter 1: 66%|██████▌ | 840/1274 [00:06<00:01, 275.71it/s]
Iter 1: 68%|██████▊ | 869/1274 [00:06<00:01, 264.03it/s]
Iter 1: 70%|███████ | 898/1274 [00:06<00:01, 264.67it/s]
Iter 1: 73%|███████▎ | 928/1274 [00:06<00:01, 273.76it/s]
Iter 1: 75%|███████▌ | 956/1274 [00:06<00:01, 263.32it/s]
Iter 1: 77%|███████▋ | 983/1274 [00:07<00:01, 261.77it/s]
Iter 1: 79%|███████▉ | 1011/1274 [00:07<00:00, 266.64it/s]
Iter 1: 82%|████████▏ | 1040/1274 [00:07<00:00, 271.80it/s]
Iter 1: 84%|████████▍ | 1068/1274 [00:07<00:00, 250.68it/s]
Iter 1: 86%|████████▌ | 1094/1274 [00:07<00:00, 245.92it/s]
Iter 1: 88%|████████▊ | 1122/1274 [00:07<00:00, 254.57it/s]
Iter 1: 90%|█████████ | 1148/1274 [00:07<00:00, 250.20it/s]
Iter 1: 92%|█████████▏| 1174/1274 [00:07<00:00, 252.86it/s]
Iter 1: 95%|█████████▍| 1204/1274 [00:07<00:00, 265.00it/s]
Iter 1: 97%|█████████▋| 1231/1274 [00:08<00:00, 240.55it/s]
Iter 1: 99%|█████████▉| 1264/1274 [00:08<00:00, 263.85it/s]
Iter 1: 100%|██████████| 1274/1274 [00:08<00:00, 144.17it/s]
ADLoc locates 1031 events outof 1274 events using 22437 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 2: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 2: 0%| | 1/1274 [00:01<36:46, 1.73s/it]
Iter 2: 0%| | 2/1274 [00:01<17:11, 1.23it/s]
Iter 2: 0%| | 3/1274 [00:02<11:37, 1.82it/s]
Iter 2: 0%| | 5/1274 [00:02<06:04, 3.48it/s]
Iter 2: 0%| | 6/1274 [00:03<08:41, 2.43it/s]
Iter 2: 1%| | 11/1274 [00:03<03:07, 6.73it/s]
Iter 2: 2%|▏ | 21/1274 [00:03<01:13, 17.13it/s]
Iter 2: 3%|▎ | 41/1274 [00:03<00:29, 42.28it/s]
Iter 2: 5%|▌ | 67/1274 [00:03<00:15, 77.87it/s]
Iter 2: 7%|▋ | 94/1274 [00:03<00:10, 114.18it/s]
Iter 2: 10%|▉ | 125/1274 [00:03<00:07, 152.54it/s]
Iter 2: 12%|█▏ | 158/1274 [00:03<00:05, 193.07it/s]
Iter 2: 15%|█▍ | 188/1274 [00:03<00:04, 218.59it/s]
Iter 2: 17%|█▋ | 215/1274 [00:03<00:04, 229.48it/s]
Iter 2: 19%|█▉ | 243/1274 [00:04<00:04, 243.03it/s]
Iter 2: 21%|██▏ | 272/1274 [00:04<00:03, 252.74it/s]
Iter 2: 23%|██▎ | 299/1274 [00:04<00:03, 249.06it/s]
Iter 2: 26%|██▌ | 326/1274 [00:04<00:03, 253.63it/s]
Iter 2: 28%|██▊ | 353/1274 [00:04<00:03, 233.24it/s]
Iter 2: 30%|███ | 383/1274 [00:04<00:03, 246.43it/s]
Iter 2: 32%|███▏ | 412/1274 [00:04<00:03, 256.01it/s]
Iter 2: 34%|███▍ | 439/1274 [00:04<00:03, 242.84it/s]
Iter 2: 37%|███▋ | 467/1274 [00:04<00:03, 252.86it/s]
Iter 2: 39%|███▉ | 495/1274 [00:05<00:02, 260.02it/s]
Iter 2: 41%|████ | 524/1274 [00:05<00:02, 266.64it/s]
Iter 2: 43%|████▎ | 551/1274 [00:05<00:02, 265.09it/s]
Iter 2: 45%|████▌ | 578/1274 [00:05<00:02, 255.85it/s]
Iter 2: 47%|████▋ | 604/1274 [00:05<00:02, 251.16it/s]
Iter 2: 50%|████▉ | 634/1274 [00:05<00:02, 260.90it/s]
Iter 2: 52%|█████▏ | 665/1274 [00:05<00:02, 274.06it/s]
Iter 2: 54%|█████▍ | 693/1274 [00:05<00:02, 213.88it/s]
Iter 2: 56%|█████▋ | 717/1274 [00:06<00:02, 214.74it/s]
Iter 2: 58%|█████▊ | 745/1274 [00:06<00:02, 226.84it/s]
Iter 2: 61%|██████ | 773/1274 [00:06<00:02, 238.62it/s]
Iter 2: 63%|██████▎ | 799/1274 [00:06<00:01, 244.24it/s]
Iter 2: 65%|██████▍ | 826/1274 [00:06<00:01, 247.56it/s]
Iter 2: 67%|██████▋ | 852/1274 [00:06<00:01, 241.92it/s]
Iter 2: 69%|██████▉ | 877/1274 [00:06<00:01, 225.37it/s]
Iter 2: 71%|███████ | 906/1274 [00:06<00:01, 242.48it/s]
Iter 2: 73%|███████▎ | 931/1274 [00:06<00:01, 241.37it/s]
Iter 2: 75%|███████▌ | 960/1274 [00:06<00:01, 253.86it/s]
Iter 2: 77%|███████▋ | 987/1274 [00:07<00:01, 258.46it/s]
Iter 2: 80%|███████▉ | 1014/1274 [00:07<00:00, 260.29it/s]
Iter 2: 82%|████████▏ | 1044/1274 [00:07<00:00, 262.93it/s]
Iter 2: 84%|████████▍ | 1071/1274 [00:07<00:00, 260.45it/s]
Iter 2: 86%|████████▋ | 1102/1274 [00:07<00:00, 270.71it/s]
Iter 2: 89%|████████▉ | 1133/1274 [00:07<00:00, 280.97it/s]
Iter 2: 92%|█████████▏| 1166/1274 [00:07<00:00, 294.21it/s]
Iter 2: 94%|█████████▍| 1198/1274 [00:07<00:00, 290.94it/s]
Iter 2: 96%|█████████▋| 1228/1274 [00:07<00:00, 280.63it/s]
Iter 2: 99%|█████████▊| 1257/1274 [00:08<00:00, 279.81it/s]
Iter 2: 100%|██████████| 1274/1274 [00:08<00:00, 144.68it/s]
ADLoc locates 1026 events outof 1274 events using 22386 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 3: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 3: 0%| | 1/1274 [00:01<37:14, 1.76s/it]
Iter 3: 0%| | 3/1274 [00:02<11:57, 1.77it/s]
Iter 3: 0%| | 5/1274 [00:02<06:49, 3.10it/s]
Iter 3: 0%| | 6/1274 [00:02<09:10, 2.30it/s]
Iter 3: 2%|▏ | 20/1274 [00:03<01:33, 13.45it/s]
Iter 3: 4%|▎ | 45/1274 [00:03<00:32, 37.92it/s]
Iter 3: 5%|▌ | 68/1274 [00:03<00:19, 62.50it/s]
Iter 3: 8%|▊ | 100/1274 [00:03<00:11, 101.91it/s]
Iter 3: 10%|█ | 129/1274 [00:03<00:08, 135.48it/s]
Iter 3: 13%|█▎ | 163/1274 [00:03<00:06, 176.86it/s]
Iter 3: 15%|█▍ | 190/1274 [00:03<00:05, 195.46it/s]
Iter 3: 17%|█▋ | 221/1274 [00:03<00:04, 222.63it/s]
Iter 3: 20%|█▉ | 249/1274 [00:03<00:04, 236.13it/s]
Iter 3: 22%|██▏ | 277/1274 [00:04<00:04, 243.50it/s]
Iter 3: 24%|██▍ | 305/1274 [00:04<00:04, 233.03it/s]
Iter 3: 26%|██▌ | 331/1274 [00:04<00:03, 239.29it/s]
Iter 3: 28%|██▊ | 359/1274 [00:04<00:04, 197.34it/s]
Iter 3: 31%|███ | 397/1274 [00:04<00:03, 239.52it/s]
Iter 3: 33%|███▎ | 426/1274 [00:04<00:03, 246.98it/s]
Iter 3: 36%|███▌ | 453/1274 [00:04<00:03, 249.57it/s]
Iter 3: 38%|███▊ | 480/1274 [00:04<00:03, 254.54it/s]
Iter 3: 40%|███▉ | 507/1274 [00:04<00:02, 255.84it/s]
Iter 3: 42%|████▏ | 538/1274 [00:05<00:02, 269.59it/s]
Iter 3: 44%|████▍ | 566/1274 [00:05<00:02, 259.50it/s]
Iter 3: 47%|████▋ | 593/1274 [00:05<00:02, 253.27it/s]
Iter 3: 49%|████▊ | 619/1274 [00:05<00:02, 250.46it/s]
Iter 3: 51%|█████ | 645/1274 [00:05<00:02, 252.51it/s]
Iter 3: 53%|█████▎ | 676/1274 [00:05<00:02, 265.34it/s]
Iter 3: 55%|█████▌ | 704/1274 [00:05<00:02, 267.39it/s]
Iter 3: 58%|█████▊ | 734/1274 [00:05<00:02, 267.92it/s]
Iter 3: 60%|██████ | 767/1274 [00:05<00:01, 283.76it/s]
Iter 3: 63%|██████▎ | 799/1274 [00:06<00:01, 288.36it/s]
Iter 3: 65%|██████▍ | 828/1274 [00:06<00:01, 274.27it/s]
Iter 3: 67%|██████▋ | 856/1274 [00:06<00:01, 243.46it/s]
Iter 3: 69%|██████▉ | 882/1274 [00:06<00:01, 245.86it/s]
Iter 3: 72%|███████▏ | 912/1274 [00:06<00:01, 252.35it/s]
Iter 3: 74%|███████▎ | 939/1274 [00:06<00:01, 257.06it/s]
Iter 3: 76%|███████▌ | 967/1274 [00:06<00:01, 261.73it/s]
Iter 3: 78%|███████▊ | 995/1274 [00:06<00:01, 264.39it/s]
Iter 3: 80%|████████ | 1023/1274 [00:06<00:00, 265.50it/s]
Iter 3: 83%|████████▎ | 1052/1274 [00:07<00:00, 269.99it/s]
Iter 3: 85%|████████▍ | 1080/1274 [00:07<00:00, 256.80it/s]
Iter 3: 87%|████████▋ | 1111/1274 [00:07<00:00, 267.65it/s]
Iter 3: 89%|████████▉ | 1138/1274 [00:07<00:00, 267.15it/s]
Iter 3: 92%|█████████▏| 1167/1274 [00:07<00:00, 271.03it/s]
Iter 3: 94%|█████████▍| 1199/1274 [00:07<00:00, 281.64it/s]
Iter 3: 96%|█████████▋| 1228/1274 [00:07<00:00, 211.04it/s]
Iter 3: 100%|█████████▉| 1270/1274 [00:07<00:00, 258.29it/s]
Iter 3: 100%|██████████| 1274/1274 [00:08<00:00, 150.49it/s]
ADLoc locates 1028 events outof 1274 events using 22391 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 4: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 4: 0%| | 1/1274 [00:01<35:13, 1.66s/it]
Iter 4: 0%| | 2/1274 [00:01<15:48, 1.34it/s]
Iter 4: 0%| | 3/1274 [00:01<09:57, 2.13it/s]
Iter 4: 0%| | 5/1274 [00:02<05:11, 4.07it/s]
Iter 4: 0%| | 6/1274 [00:02<08:44, 2.42it/s]
Iter 4: 1%| | 13/1274 [00:02<02:29, 8.46it/s]
Iter 4: 3%|▎ | 41/1274 [00:03<00:31, 38.97it/s]
Iter 4: 5%|▌ | 68/1274 [00:03<00:17, 70.87it/s]
Iter 4: 8%|▊ | 105/1274 [00:03<00:09, 120.17it/s]
Iter 4: 11%|█ | 134/1274 [00:03<00:07, 151.75it/s]
Iter 4: 13%|█▎ | 165/1274 [00:03<00:05, 184.94it/s]
Iter 4: 15%|█▌ | 192/1274 [00:03<00:05, 203.28it/s]
Iter 4: 17%|█▋ | 222/1274 [00:03<00:04, 226.74it/s]
Iter 4: 20%|█▉ | 250/1274 [00:03<00:04, 233.14it/s]
Iter 4: 22%|██▏ | 280/1274 [00:03<00:04, 247.85it/s]
Iter 4: 24%|██▍ | 308/1274 [00:04<00:03, 256.61it/s]
Iter 4: 26%|██▋ | 336/1274 [00:04<00:03, 253.76it/s]
Iter 4: 29%|██▊ | 365/1274 [00:04<00:03, 262.58it/s]
Iter 4: 31%|███ | 396/1274 [00:04<00:03, 269.18it/s]
Iter 4: 34%|███▎ | 427/1274 [00:04<00:03, 279.46it/s]
Iter 4: 36%|███▌ | 456/1274 [00:04<00:03, 266.92it/s]
Iter 4: 38%|███▊ | 484/1274 [00:04<00:02, 270.52it/s]
Iter 4: 40%|████ | 512/1274 [00:04<00:02, 259.19it/s]
Iter 4: 42%|████▏ | 539/1274 [00:04<00:02, 260.96it/s]
Iter 4: 44%|████▍ | 566/1274 [00:04<00:02, 262.35it/s]
Iter 4: 47%|████▋ | 593/1274 [00:05<00:02, 251.81it/s]
Iter 4: 49%|████▊ | 619/1274 [00:05<00:03, 188.17it/s]
Iter 4: 52%|█████▏ | 661/1274 [00:05<00:02, 238.55it/s]
Iter 4: 54%|█████▍ | 689/1274 [00:05<00:02, 245.83it/s]
Iter 4: 57%|█████▋ | 724/1274 [00:05<00:02, 269.81it/s]
Iter 4: 59%|█████▉ | 753/1274 [00:05<00:01, 270.86it/s]
Iter 4: 62%|██████▏ | 784/1274 [00:05<00:01, 279.65it/s]
Iter 4: 64%|██████▍ | 813/1274 [00:05<00:01, 274.08it/s]
Iter 4: 66%|██████▌ | 842/1274 [00:06<00:01, 262.88it/s]
Iter 4: 68%|██████▊ | 869/1274 [00:06<00:01, 254.19it/s]
Iter 4: 70%|███████ | 895/1274 [00:06<00:01, 250.96it/s]
Iter 4: 72%|███████▏ | 921/1274 [00:06<00:01, 251.67it/s]
Iter 4: 74%|███████▍ | 948/1274 [00:06<00:01, 254.87it/s]
Iter 4: 77%|███████▋ | 976/1274 [00:06<00:01, 257.75it/s]
Iter 4: 79%|███████▉ | 1009/1274 [00:06<00:00, 275.93it/s]
Iter 4: 81%|████████▏ | 1037/1274 [00:06<00:00, 256.47it/s]
Iter 4: 84%|████████▎ | 1064/1274 [00:06<00:00, 255.48it/s]
Iter 4: 86%|████████▌ | 1094/1274 [00:07<00:00, 264.10it/s]
Iter 4: 88%|████████▊ | 1126/1274 [00:07<00:00, 277.58it/s]
Iter 4: 91%|█████████ | 1155/1274 [00:07<00:00, 280.40it/s]
Iter 4: 93%|█████████▎| 1184/1274 [00:07<00:00, 277.90it/s]
Iter 4: 95%|█████████▌| 1216/1274 [00:07<00:00, 288.86it/s]
Iter 4: 98%|█████████▊| 1245/1274 [00:07<00:00, 284.77it/s]
Iter 4: 100%|██████████| 1274/1274 [00:08<00:00, 152.67it/s]
ADLoc locates 1027 events outof 1274 events using 22390 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 5: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 5: 0%| | 1/1274 [00:01<37:15, 1.76s/it]
Iter 5: 0%| | 3/1274 [00:01<10:56, 1.94it/s]
Iter 5: 0%| | 4/1274 [00:02<08:00, 2.64it/s]
Iter 5: 0%| | 5/1274 [00:02<06:45, 3.13it/s]
Iter 5: 0%| | 6/1274 [00:03<10:55, 1.93it/s]
Iter 5: 1%| | 12/1274 [00:03<03:09, 6.66it/s]
Iter 5: 2%|▏ | 29/1274 [00:03<00:53, 23.44it/s]
Iter 5: 4%|▍ | 53/1274 [00:03<00:24, 50.32it/s]
Iter 5: 6%|▌ | 77/1274 [00:03<00:15, 79.26it/s]
Iter 5: 8%|▊ | 101/1274 [00:03<00:10, 107.42it/s]
Iter 5: 9%|▉ | 119/1274 [00:03<00:09, 121.27it/s]
Iter 5: 11%|█ | 141/1274 [00:03<00:08, 137.76it/s]
Iter 5: 13%|█▎ | 163/1274 [00:04<00:07, 156.13it/s]
Iter 5: 15%|█▍ | 189/1274 [00:04<00:06, 175.39it/s]
Iter 5: 17%|█▋ | 213/1274 [00:04<00:05, 187.18it/s]
Iter 5: 18%|█▊ | 234/1274 [00:04<00:08, 128.85it/s]
Iter 5: 20%|██ | 261/1274 [00:04<00:06, 155.46it/s]
Iter 5: 22%|██▏ | 281/1274 [00:04<00:06, 165.05it/s]
Iter 5: 24%|██▍ | 305/1274 [00:04<00:05, 180.80it/s]
Iter 5: 26%|██▌ | 326/1274 [00:04<00:05, 172.96it/s]
Iter 5: 27%|██▋ | 348/1274 [00:05<00:05, 184.53it/s]
Iter 5: 29%|██▉ | 369/1274 [00:05<00:04, 190.66it/s]
Iter 5: 31%|███ | 391/1274 [00:05<00:04, 193.94it/s]
Iter 5: 32%|███▏ | 414/1274 [00:05<00:04, 200.35it/s]
Iter 5: 34%|███▍ | 435/1274 [00:05<00:04, 195.01it/s]
Iter 5: 36%|███▌ | 455/1274 [00:05<00:04, 181.40it/s]
Iter 5: 37%|███▋ | 474/1274 [00:05<00:04, 181.43it/s]
Iter 5: 39%|███▊ | 493/1274 [00:05<00:04, 181.33it/s]
Iter 5: 40%|████ | 515/1274 [00:05<00:04, 182.69it/s]
Iter 5: 42%|████▏ | 534/1274 [00:06<00:04, 183.11it/s]
Iter 5: 43%|████▎ | 554/1274 [00:06<00:04, 170.91it/s]
Iter 5: 45%|████▍ | 573/1274 [00:06<00:04, 173.56it/s]
Iter 5: 46%|████▋ | 591/1274 [00:06<00:04, 158.14it/s]
Iter 5: 48%|████▊ | 608/1274 [00:06<00:04, 160.29it/s]
Iter 5: 49%|████▉ | 625/1274 [00:06<00:04, 161.75it/s]
Iter 5: 51%|█████ | 646/1274 [00:06<00:03, 172.66it/s]
Iter 5: 52%|█████▏ | 667/1274 [00:06<00:03, 174.29it/s]
Iter 5: 54%|█████▍ | 687/1274 [00:06<00:03, 181.39it/s]
Iter 5: 55%|█████▌ | 707/1274 [00:07<00:03, 186.31it/s]
Iter 5: 58%|█████▊ | 733/1274 [00:07<00:02, 198.50it/s]
Iter 5: 59%|█████▉ | 753/1274 [00:07<00:02, 197.66it/s]
Iter 5: 61%|██████ | 778/1274 [00:07<00:02, 205.85it/s]
Iter 5: 63%|██████▎ | 800/1274 [00:07<00:02, 204.89it/s]
Iter 5: 64%|██████▍ | 821/1274 [00:07<00:02, 196.71it/s]
Iter 5: 66%|██████▌ | 841/1274 [00:07<00:02, 186.88it/s]
Iter 5: 68%|██████▊ | 860/1274 [00:07<00:02, 178.07it/s]
Iter 5: 69%|██████▉ | 880/1274 [00:07<00:02, 180.54it/s]
Iter 5: 71%|███████ | 899/1274 [00:08<00:02, 182.19it/s]
Iter 5: 72%|███████▏ | 918/1274 [00:08<00:02, 171.07it/s]
Iter 5: 74%|███████▎ | 937/1274 [00:08<00:01, 168.92it/s]
Iter 5: 75%|███████▍ | 955/1274 [00:08<00:01, 171.54it/s]
Iter 5: 76%|███████▋ | 973/1274 [00:08<00:01, 173.37it/s]
Iter 5: 78%|███████▊ | 993/1274 [00:08<00:01, 179.85it/s]
Iter 5: 80%|███████▉ | 1014/1274 [00:08<00:01, 177.58it/s]
Iter 5: 81%|████████ | 1034/1274 [00:08<00:01, 176.91it/s]
Iter 5: 83%|████████▎ | 1052/1274 [00:08<00:01, 175.65it/s]
Iter 5: 84%|████████▍ | 1072/1274 [00:09<00:01, 179.49it/s]
Iter 5: 86%|████████▌ | 1090/1274 [00:09<00:01, 176.36it/s]
Iter 5: 87%|████████▋ | 1108/1274 [00:09<00:01, 127.89it/s]
Iter 5: 88%|████████▊ | 1125/1274 [00:09<00:01, 134.46it/s]
Iter 5: 90%|█████████ | 1149/1274 [00:09<00:00, 157.85it/s]
Iter 5: 92%|█████████▏| 1171/1274 [00:09<00:00, 172.54it/s]
Iter 5: 94%|█████████▍| 1197/1274 [00:09<00:00, 194.15it/s]
Iter 5: 96%|█████████▌| 1218/1274 [00:09<00:00, 184.92it/s]
Iter 5: 97%|█████████▋| 1238/1274 [00:10<00:00, 180.60it/s]
Iter 5: 99%|█████████▊| 1257/1274 [00:10<00:00, 180.34it/s]
Iter 5: 100%|██████████| 1274/1274 [00:11<00:00, 112.54it/s]
ADLoc locates 1026 events outof 1274 events using 22366 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 6: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 6: 0%| | 1/1274 [00:03<1:05:12, 3.07s/it]
Iter 6: 0%| | 3/1274 [00:03<18:47, 1.13it/s]
Iter 6: 0%| | 4/1274 [00:03<14:55, 1.42it/s]
Iter 6: 0%| | 5/1274 [00:03<11:58, 1.77it/s]
Iter 6: 0%| | 6/1274 [00:05<16:46, 1.26it/s]
Iter 6: 1%|▏ | 18/1274 [00:05<02:42, 7.72it/s]
Iter 6: 3%|▎ | 32/1274 [00:05<01:12, 17.22it/s]
Iter 6: 4%|▍ | 49/1274 [00:05<00:39, 31.04it/s]
Iter 6: 5%|▌ | 65/1274 [00:05<00:26, 45.41it/s]
Iter 6: 7%|▋ | 83/1274 [00:05<00:18, 63.58it/s]
Iter 6: 8%|▊ | 104/1274 [00:05<00:13, 87.82it/s]
Iter 6: 9%|▉ | 121/1274 [00:05<00:11, 103.25it/s]
Iter 6: 11%|█ | 139/1274 [00:06<00:09, 118.55it/s]
Iter 6: 12%|█▏ | 158/1274 [00:06<00:08, 134.29it/s]
Iter 6: 14%|█▍ | 178/1274 [00:06<00:07, 146.94it/s]
Iter 6: 16%|█▌ | 202/1274 [00:06<00:06, 163.37it/s]
Iter 6: 17%|█▋ | 221/1274 [00:06<00:06, 168.01it/s]
Iter 6: 19%|█▉ | 242/1274 [00:06<00:05, 178.31it/s]
Iter 6: 21%|██ | 263/1274 [00:06<00:05, 184.04it/s]
Iter 6: 22%|██▏ | 283/1274 [00:06<00:05, 174.86it/s]
Iter 6: 24%|██▍ | 307/1274 [00:06<00:05, 177.87it/s]
Iter 6: 26%|██▌ | 327/1274 [00:07<00:05, 182.88it/s]
Iter 6: 27%|██▋ | 346/1274 [00:07<00:05, 157.94it/s]
Iter 6: 28%|██▊ | 363/1274 [00:07<00:05, 154.32it/s]
Iter 6: 30%|███ | 385/1274 [00:07<00:05, 170.43it/s]
Iter 6: 32%|███▏ | 407/1274 [00:07<00:04, 180.50it/s]
Iter 6: 33%|███▎ | 426/1274 [00:07<00:06, 132.63it/s]
Iter 6: 35%|███▍ | 442/1274 [00:07<00:06, 138.23it/s]
Iter 6: 36%|███▌ | 459/1274 [00:08<00:05, 141.03it/s]
Iter 6: 38%|███▊ | 478/1274 [00:08<00:05, 152.18it/s]
Iter 6: 39%|███▉ | 499/1274 [00:08<00:04, 161.60it/s]
Iter 6: 41%|████ | 516/1274 [00:08<00:04, 159.94it/s]
Iter 6: 42%|████▏ | 533/1274 [00:08<00:04, 156.64it/s]
Iter 6: 43%|████▎ | 550/1274 [00:08<00:04, 145.38it/s]
Iter 6: 44%|████▍ | 565/1274 [00:08<00:05, 133.66it/s]
Iter 6: 45%|████▌ | 579/1274 [00:08<00:05, 135.09it/s]
Iter 6: 47%|████▋ | 593/1274 [00:08<00:05, 127.65it/s]
Iter 6: 48%|████▊ | 609/1274 [00:09<00:05, 128.20it/s]
Iter 6: 49%|████▉ | 622/1274 [00:09<00:05, 120.52it/s]
Iter 6: 50%|████▉ | 635/1274 [00:09<00:05, 117.14it/s]
Iter 6: 51%|█████ | 648/1274 [00:09<00:05, 118.95it/s]
Iter 6: 52%|█████▏ | 662/1274 [00:09<00:04, 124.02it/s]
Iter 6: 53%|█████▎ | 675/1274 [00:09<00:04, 124.07it/s]
Iter 6: 54%|█████▍ | 688/1274 [00:09<00:04, 125.04it/s]
Iter 6: 55%|█████▌ | 703/1274 [00:09<00:04, 130.98it/s]
Iter 6: 56%|█████▋ | 717/1274 [00:09<00:04, 130.02it/s]
Iter 6: 57%|█████▋ | 732/1274 [00:10<00:04, 134.54it/s]
Iter 6: 59%|█████▊ | 746/1274 [00:10<00:04, 122.86it/s]
Iter 6: 60%|█████▉ | 759/1274 [00:10<00:04, 121.24it/s]
Iter 6: 61%|██████ | 772/1274 [00:10<00:04, 123.63it/s]
Iter 6: 62%|██████▏ | 785/1274 [00:10<00:04, 117.97it/s]
Iter 6: 63%|██████▎ | 797/1274 [00:10<00:04, 109.67it/s]
Iter 6: 64%|██████▎ | 809/1274 [00:10<00:04, 106.17it/s]
Iter 6: 64%|██████▍ | 820/1274 [00:10<00:04, 100.06it/s]
Iter 6: 65%|██████▌ | 831/1274 [00:11<00:04, 95.77it/s]
Iter 6: 66%|██████▌ | 843/1274 [00:11<00:04, 101.95it/s]
Iter 6: 67%|██████▋ | 854/1274 [00:11<00:04, 97.30it/s]
Iter 6: 68%|██████▊ | 868/1274 [00:11<00:03, 104.09it/s]
Iter 6: 69%|██████▉ | 880/1274 [00:11<00:03, 107.84it/s]
Iter 6: 70%|███████ | 896/1274 [00:11<00:03, 120.99it/s]
Iter 6: 72%|███████▏ | 911/1274 [00:11<00:02, 125.99it/s]
Iter 6: 73%|███████▎ | 924/1274 [00:11<00:02, 123.50it/s]
Iter 6: 74%|███████▎ | 939/1274 [00:11<00:02, 126.95it/s]
Iter 6: 75%|███████▌ | 956/1274 [00:12<00:02, 138.36it/s]
Iter 6: 76%|███████▌ | 970/1274 [00:12<00:02, 133.38it/s]
Iter 6: 77%|███████▋ | 984/1274 [00:12<00:02, 129.12it/s]
Iter 6: 78%|███████▊ | 1000/1274 [00:12<00:02, 135.67it/s]
Iter 6: 80%|███████▉ | 1014/1274 [00:12<00:01, 133.30it/s]
Iter 6: 81%|████████ | 1031/1274 [00:12<00:01, 140.89it/s]
Iter 6: 82%|████████▏ | 1049/1274 [00:12<00:01, 150.41it/s]
Iter 6: 84%|████████▎ | 1065/1274 [00:12<00:01, 140.26it/s]
Iter 6: 85%|████████▍ | 1082/1274 [00:12<00:01, 145.81it/s]
Iter 6: 87%|████████▋ | 1103/1274 [00:13<00:01, 160.84it/s]
Iter 6: 88%|████████▊ | 1120/1274 [00:13<00:01, 149.57it/s]
Iter 6: 89%|████████▉ | 1137/1274 [00:13<00:00, 152.35it/s]
Iter 6: 91%|█████████ | 1154/1274 [00:13<00:00, 152.49it/s]
Iter 6: 92%|█████████▏| 1170/1274 [00:13<00:00, 143.73it/s]
Iter 6: 93%|█████████▎| 1185/1274 [00:13<00:00, 143.41it/s]
Iter 6: 95%|█████████▍| 1204/1274 [00:13<00:00, 155.23it/s]
Iter 6: 96%|█████████▌| 1220/1274 [00:13<00:00, 150.83it/s]
Iter 6: 97%|█████████▋| 1236/1274 [00:13<00:00, 140.51it/s]
Iter 6: 98%|█████████▊| 1251/1274 [00:14<00:00, 128.70it/s]
Iter 6: 100%|█████████▉| 1269/1274 [00:14<00:00, 138.03it/s]
Iter 6: 100%|██████████| 1274/1274 [00:15<00:00, 82.16it/s]
ADLoc locates 1022 events outof 1274 events using 22306 picks outof 24681 picks
Using 17 CPUs for ADLoc inversion.
Iter 7: 0%| | 0/1274 [00:00<?, ?it/s]
Iter 7: 0%| | 1/1274 [00:02<47:51, 2.26s/it]
Iter 7: 0%| | 3/1274 [00:02<14:32, 1.46it/s]
Iter 7: 0%| | 4/1274 [00:02<12:42, 1.67it/s]
Iter 7: 0%| | 5/1274 [00:03<11:07, 1.90it/s]
Iter 7: 0%| | 6/1274 [00:04<14:59, 1.41it/s]
Iter 7: 1%| | 9/1274 [00:04<06:44, 3.13it/s]
Iter 7: 1%|▏ | 16/1274 [00:04<02:29, 8.43it/s]
Iter 7: 2%|▏ | 24/1274 [00:04<01:20, 15.50it/s]
Iter 7: 3%|▎ | 35/1274 [00:04<00:46, 26.74it/s]
Iter 7: 4%|▎ | 47/1274 [00:04<00:30, 40.56it/s]
Iter 7: 5%|▍ | 59/1274 [00:05<00:22, 54.48it/s]
Iter 7: 5%|▌ | 69/1274 [00:05<00:19, 62.20it/s]
Iter 7: 6%|▌ | 79/1274 [00:05<00:17, 70.20it/s]
Iter 7: 7%|▋ | 91/1274 [00:05<00:14, 81.55it/s]
Iter 7: 8%|▊ | 103/1274 [00:05<00:13, 86.38it/s]
Iter 7: 9%|▉ | 113/1274 [00:05<00:13, 84.59it/s]
Iter 7: 10%|▉ | 123/1274 [00:05<00:13, 87.57it/s]
Iter 7: 11%|█ | 136/1274 [00:05<00:11, 96.49it/s]
Iter 7: 12%|█▏ | 148/1274 [00:05<00:11, 100.38it/s]
Iter 7: 13%|█▎ | 162/1274 [00:06<00:10, 107.60it/s]
Iter 7: 14%|█▎ | 175/1274 [00:06<00:10, 104.21it/s]
Iter 7: 15%|█▍ | 188/1274 [00:06<00:10, 108.57it/s]
Iter 7: 16%|█▌ | 200/1274 [00:06<00:09, 111.65it/s]
Iter 7: 17%|█▋ | 212/1274 [00:06<00:09, 107.52it/s]
Iter 7: 18%|█▊ | 223/1274 [00:06<00:10, 103.35it/s]
Iter 7: 18%|█▊ | 234/1274 [00:06<00:10, 103.33it/s]
Iter 7: 19%|█▉ | 246/1274 [00:06<00:11, 92.38it/s]
Iter 7: 20%|██ | 261/1274 [00:07<00:09, 105.26it/s]
Iter 7: 21%|██▏ | 272/1274 [00:07<00:09, 105.31it/s]
Iter 7: 22%|██▏ | 284/1274 [00:07<00:09, 107.53it/s]
Iter 7: 23%|██▎ | 295/1274 [00:07<00:09, 104.34it/s]
Iter 7: 24%|██▍ | 310/1274 [00:07<00:08, 116.58it/s]
Iter 7: 25%|██▌ | 322/1274 [00:07<00:08, 112.78it/s]
Iter 7: 26%|██▌ | 334/1274 [00:07<00:08, 113.53it/s]
Iter 7: 28%|██▊ | 353/1274 [00:07<00:07, 130.75it/s]
Iter 7: 29%|██▉ | 373/1274 [00:07<00:06, 148.42it/s]
Iter 7: 31%|███ | 393/1274 [00:08<00:05, 160.36it/s]
Iter 7: 32%|███▏ | 410/1274 [00:08<00:05, 160.84it/s]
Iter 7: 34%|███▍ | 431/1274 [00:08<00:05, 165.74it/s]
Iter 7: 35%|███▌ | 450/1274 [00:08<00:04, 167.31it/s]
Iter 7: 37%|███▋ | 471/1274 [00:08<00:04, 179.09it/s]
Iter 7: 38%|███▊ | 489/1274 [00:08<00:04, 174.10it/s]
Iter 7: 40%|████ | 513/1274 [00:08<00:03, 191.41it/s]
Iter 7: 42%|████▏ | 536/1274 [00:08<00:03, 200.59it/s]
Iter 7: 44%|████▍ | 558/1274 [00:08<00:03, 203.22it/s]
Iter 7: 46%|████▌ | 581/1274 [00:08<00:03, 205.26it/s]
Iter 7: 47%|████▋ | 602/1274 [00:09<00:03, 200.42it/s]
Iter 7: 49%|████▉ | 623/1274 [00:09<00:03, 203.10it/s]
Iter 7: 51%|█████ | 651/1274 [00:09<00:02, 223.42it/s]
Iter 7: 53%|█████▎ | 678/1274 [00:09<00:02, 236.77it/s]
Iter 7: 55%|█████▌ | 702/1274 [00:09<00:02, 236.74it/s]
Iter 7: 57%|█████▋ | 729/1274 [00:09<00:02, 245.63it/s]
Iter 7: 59%|█████▉ | 754/1274 [00:09<00:02, 209.19it/s]
Iter 7: 61%|██████ | 778/1274 [00:09<00:02, 215.52it/s]
Iter 7: 63%|██████▎ | 806/1274 [00:09<00:02, 226.76it/s]
Iter 7: 65%|██████▌ | 830/1274 [00:10<00:02, 218.39it/s]
Iter 7: 67%|██████▋ | 853/1274 [00:10<00:01, 216.96it/s]
Iter 7: 69%|██████▉ | 877/1274 [00:10<00:01, 216.10it/s]
Iter 7: 71%|███████▏ | 908/1274 [00:10<00:01, 241.72it/s]
Iter 7: 73%|███████▎ | 933/1274 [00:10<00:01, 240.42it/s]
Iter 7: 76%|███████▌ | 967/1274 [00:10<00:01, 264.94it/s]
Iter 7: 78%|███████▊ | 994/1274 [00:10<00:01, 259.80it/s]
Iter 7: 80%|████████ | 1025/1274 [00:10<00:00, 272.09it/s]
Iter 7: 83%|████████▎ | 1053/1274 [00:10<00:00, 265.19it/s]
Iter 7: 85%|████████▍ | 1080/1274 [00:11<00:00, 261.17it/s]
Iter 7: 87%|████████▋ | 1110/1274 [00:11<00:00, 272.05it/s]
Iter 7: 89%|████████▉ | 1140/1274 [00:11<00:00, 278.30it/s]
Iter 7: 92%|█████████▏| 1170/1274 [00:11<00:00, 281.42it/s]
Iter 7: 95%|█████████▍| 1205/1274 [00:11<00:00, 299.11it/s]
Iter 7: 97%|█████████▋| 1235/1274 [00:11<00:00, 270.05it/s]
Iter 7: 99%|█████████▉| 1263/1274 [00:11<00:00, 270.35it/s]
Iter 7: 100%|██████████| 1274/1274 [00:12<00:00, 102.83it/s]
ADLoc locates 1025 events outof 1274 events using 22389 picks outof 24681 picks
In [18]:
Copied!
plot_catalog(adloc_catalog, "ADLoc", region=region, config=config)
plot_catalog(adloc_catalog, "ADLoc", region=region, config=config)
Run HypoDD to relocate relative earthquake locations¶
Credit: Felix Waldhauser
In [19]:
Copied!
def run_hypodd(root_path: str = "local", region: str = "demo"):
data_path = f"{region}/adloc"
result_path = f"{region}/hypodd"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
## Station Format
stations = pd.read_csv(f"{root_path}/{data_path}/adloc_stations.csv")
stations.set_index("station_id", inplace=True)
shift_topo = stations["elevation_m"].max() / 1e3
converted_hypoinverse = []
converted_hypodd = {}
for sta, row in stations.iterrows():
network_code, station_code, comp_code, channel_code = sta.split(".")
station_weight = " "
lat_degree = int(row["latitude"])
lat_minute = (row["latitude"] - lat_degree) * 60
north = "N" if lat_degree >= 0 else "S"
lng_degree = int(row["longitude"])
lng_minute = (row["longitude"] - lng_degree) * 60
west = "W" if lng_degree <= 0 else "E"
elevation = row["elevation_m"]
line_hypoinverse = f"{station_code:<5} {network_code:<2} {comp_code[:-1]:<1}{channel_code:<3} {station_weight}{abs(lat_degree):2.0f} {abs(lat_minute):7.4f}{north}{abs(lng_degree):3.0f} {abs(lng_minute):7.4f}{west}{elevation:4.0f}\n"
converted_hypoinverse.append(line_hypoinverse)
# tmp_code = f"{station_code}{channel_code}"
tmp_code = f"{station_code}"
converted_hypodd[tmp_code] = f"{tmp_code:<8s} {row['latitude']:.3f} {row['longitude']:.3f}\n"
with open(f"{root_path}/{result_path}/stations.dat", "w") as f:
for k, v in converted_hypodd.items():
f.write(v)
## Picks Format
picks_csv = f"{data_path}/adloc_picks.csv"
events_csv = f"{data_path}/adloc_events.csv"
picks = pd.read_csv(f"{root_path}/{picks_csv}")
events = pd.read_csv(f"{root_path}/{events_csv}")
picks["phase_time"] = pd.to_datetime(picks["phase_time"], format="mixed")
events["time"] = pd.to_datetime(events["time"])
# events["magnitude"] = 1.0
events["sigma_time"] = 1.0
# events.sort_values("time", inplace=True)
picks = picks.loc[picks["event_index"].isin(events["event_index"])]
lines = []
picks_by_event = picks.groupby("event_index").groups
for i, event in tqdm(events.iterrows(), desc="Convert catalog", total=len(events)):
# event_time = datetime.strptime(event["time"], "%Y-%m-%dT%H:%M:%S.%f")
event_time = event["time"]
lat = event["latitude"]
lng = event["longitude"]
# dep = event["depth(m)"] / 1e3 + shift_topo
dep = event["depth_km"] + shift_topo
mag = event["magnitude"]
EH = 0
EZ = 0
RMS = event["sigma_time"]
year, month, day, hour, min, sec = (
event_time.year,
event_time.month,
event_time.day,
event_time.hour,
event_time.minute,
float(event_time.strftime("%S.%f")),
)
event_line = f"# {year:4d} {month:2d} {day:2d} {hour:2d} {min:2d} {sec:5.2f} {lat:7.4f} {lng:9.4f} {dep:5.2f} {mag:5.2f} {EH:5.2f} {EZ:5.2f} {RMS:5.2f} {event['event_index']:9d}\n"
lines.append(event_line)
picks_idx = picks_by_event[event["event_index"]]
for j in picks_idx:
# pick = picks.iloc[j]
pick = picks.loc[j]
network_code, station_code, comp_code, channel_code = pick["station_id"].split(".")
phase_type = pick["phase_type"].upper()
phase_score = pick["phase_score"]
# pick_time = (datetime.strptime(pick["phase_time"], "%Y-%m-%dT%H:%M:%S.%f") - event_time).total_seconds()
pick_time = (pick["phase_time"] - event_time).total_seconds()
tmp_code = f"{station_code}"
pick_line = f"{tmp_code:<7s} {pick_time:6.3f} {phase_score:5.4f} {phase_type}\n"
lines.append(pick_line)
with open(f"{root_path}/{result_path}/phase.txt", "w") as fp:
fp.writelines(lines)
## Run Hypodd
print(f"Running Hypodd:")
os.system(f"bash run_hypodd_ct.sh {root_path} {region}")
## Read catalog
columns = ["ID", "LAT", "LON", "DEPTH", "X", "Y", "Z", "EX", "EY", "EZ", "YR", "MO", "DY", "HR", "MI", "SC", "MAG", "NCCP", "NCCS", "NCTP", "NCTS", "RCC", "RCT", "CID"]
catalog_ct_hypodd = pd.read_csv(f"{root_path}/{region}/hypodd/hypodd_ct.reloc", sep="\s+", header=None, names=columns, dtype=float)
catalog_ct_hypodd["time"] = catalog_ct_hypodd.apply(
lambda x: f'{x["YR"]:04.0f}-{x["MO"]:02.0f}-{x["DY"]:02.0f}T{x["HR"]:02.0f}:{x["MI"]:02.0f}:{np.min([float(x["SC"]), 59.999]):05.3f}',
axis=1,
)
catalog_ct_hypodd["time"] = catalog_ct_hypodd["time"].apply(lambda x: datetime.strptime(x, "%Y-%m-%dT%H:%M:%S.%f"))
catalog_ct_hypodd = catalog_ct_hypodd[catalog_ct_hypodd["DEPTH"] != "*********"]
catalog_ct_hypodd["DEPTH"] = catalog_ct_hypodd["DEPTH"].astype(float)
catalog_ct_hypodd.rename({"ID": "event_index", "LAT": "latitude", "LON": "longitude", "DEPTH": "depth_km", "MAG": "magnitude"}, axis=1, inplace=True)
catalog_ct_hypodd.to_csv(f"{root_path}/{region}/hypodd/hypodd_ct.csv", index=False)
return catalog_ct_hypodd
def run_hypodd(root_path: str = "local", region: str = "demo"):
data_path = f"{region}/adloc"
result_path = f"{region}/hypodd"
if not os.path.exists(f"{root_path}/{result_path}"):
os.makedirs(f"{root_path}/{result_path}")
## Station Format
stations = pd.read_csv(f"{root_path}/{data_path}/adloc_stations.csv")
stations.set_index("station_id", inplace=True)
shift_topo = stations["elevation_m"].max() / 1e3
converted_hypoinverse = []
converted_hypodd = {}
for sta, row in stations.iterrows():
network_code, station_code, comp_code, channel_code = sta.split(".")
station_weight = " "
lat_degree = int(row["latitude"])
lat_minute = (row["latitude"] - lat_degree) * 60
north = "N" if lat_degree >= 0 else "S"
lng_degree = int(row["longitude"])
lng_minute = (row["longitude"] - lng_degree) * 60
west = "W" if lng_degree <= 0 else "E"
elevation = row["elevation_m"]
line_hypoinverse = f"{station_code:<5} {network_code:<2} {comp_code[:-1]:<1}{channel_code:<3} {station_weight}{abs(lat_degree):2.0f} {abs(lat_minute):7.4f}{north}{abs(lng_degree):3.0f} {abs(lng_minute):7.4f}{west}{elevation:4.0f}\n"
converted_hypoinverse.append(line_hypoinverse)
# tmp_code = f"{station_code}{channel_code}"
tmp_code = f"{station_code}"
converted_hypodd[tmp_code] = f"{tmp_code:<8s} {row['latitude']:.3f} {row['longitude']:.3f}\n"
with open(f"{root_path}/{result_path}/stations.dat", "w") as f:
for k, v in converted_hypodd.items():
f.write(v)
## Picks Format
picks_csv = f"{data_path}/adloc_picks.csv"
events_csv = f"{data_path}/adloc_events.csv"
picks = pd.read_csv(f"{root_path}/{picks_csv}")
events = pd.read_csv(f"{root_path}/{events_csv}")
picks["phase_time"] = pd.to_datetime(picks["phase_time"], format="mixed")
events["time"] = pd.to_datetime(events["time"])
# events["magnitude"] = 1.0
events["sigma_time"] = 1.0
# events.sort_values("time", inplace=True)
picks = picks.loc[picks["event_index"].isin(events["event_index"])]
lines = []
picks_by_event = picks.groupby("event_index").groups
for i, event in tqdm(events.iterrows(), desc="Convert catalog", total=len(events)):
# event_time = datetime.strptime(event["time"], "%Y-%m-%dT%H:%M:%S.%f")
event_time = event["time"]
lat = event["latitude"]
lng = event["longitude"]
# dep = event["depth(m)"] / 1e3 + shift_topo
dep = event["depth_km"] + shift_topo
mag = event["magnitude"]
EH = 0
EZ = 0
RMS = event["sigma_time"]
year, month, day, hour, min, sec = (
event_time.year,
event_time.month,
event_time.day,
event_time.hour,
event_time.minute,
float(event_time.strftime("%S.%f")),
)
event_line = f"# {year:4d} {month:2d} {day:2d} {hour:2d} {min:2d} {sec:5.2f} {lat:7.4f} {lng:9.4f} {dep:5.2f} {mag:5.2f} {EH:5.2f} {EZ:5.2f} {RMS:5.2f} {event['event_index']:9d}\n"
lines.append(event_line)
picks_idx = picks_by_event[event["event_index"]]
for j in picks_idx:
# pick = picks.iloc[j]
pick = picks.loc[j]
network_code, station_code, comp_code, channel_code = pick["station_id"].split(".")
phase_type = pick["phase_type"].upper()
phase_score = pick["phase_score"]
# pick_time = (datetime.strptime(pick["phase_time"], "%Y-%m-%dT%H:%M:%S.%f") - event_time).total_seconds()
pick_time = (pick["phase_time"] - event_time).total_seconds()
tmp_code = f"{station_code}"
pick_line = f"{tmp_code:<7s} {pick_time:6.3f} {phase_score:5.4f} {phase_type}\n"
lines.append(pick_line)
with open(f"{root_path}/{result_path}/phase.txt", "w") as fp:
fp.writelines(lines)
## Run Hypodd
print(f"Running Hypodd:")
os.system(f"bash run_hypodd_ct.sh {root_path} {region}")
## Read catalog
columns = ["ID", "LAT", "LON", "DEPTH", "X", "Y", "Z", "EX", "EY", "EZ", "YR", "MO", "DY", "HR", "MI", "SC", "MAG", "NCCP", "NCCS", "NCTP", "NCTS", "RCC", "RCT", "CID"]
catalog_ct_hypodd = pd.read_csv(f"{root_path}/{region}/hypodd/hypodd_ct.reloc", sep="\s+", header=None, names=columns, dtype=float)
catalog_ct_hypodd["time"] = catalog_ct_hypodd.apply(
lambda x: f'{x["YR"]:04.0f}-{x["MO"]:02.0f}-{x["DY"]:02.0f}T{x["HR"]:02.0f}:{x["MI"]:02.0f}:{np.min([float(x["SC"]), 59.999]):05.3f}',
axis=1,
)
catalog_ct_hypodd["time"] = catalog_ct_hypodd["time"].apply(lambda x: datetime.strptime(x, "%Y-%m-%dT%H:%M:%S.%f"))
catalog_ct_hypodd = catalog_ct_hypodd[catalog_ct_hypodd["DEPTH"] != "*********"]
catalog_ct_hypodd["DEPTH"] = catalog_ct_hypodd["DEPTH"].astype(float)
catalog_ct_hypodd.rename({"ID": "event_index", "LAT": "latitude", "LON": "longitude", "DEPTH": "depth_km", "MAG": "magnitude"}, axis=1, inplace=True)
catalog_ct_hypodd.to_csv(f"{root_path}/{region}/hypodd/hypodd_ct.csv", index=False)
return catalog_ct_hypodd
<>:101: SyntaxWarning: invalid escape sequence '\s'
<>:101: SyntaxWarning: invalid escape sequence '\s'
/var/folders/0n/_fgybwj14vxfg5mgnsyxy3f40000gn/T/ipykernel_63161/1322883633.py:101: SyntaxWarning: invalid escape sequence '\s'
catalog_ct_hypodd = pd.read_csv(f"{root_path}/{region}/hypodd/hypodd_ct.reloc", sep="\s+", header=None, names=columns, dtype=float)
In [20]:
Copied!
hypodd_catalog = run_hypodd(region=region)
hypodd_catalog = run_hypodd(region=region)
Convert catalog: 0%| | 0/1025 [00:00<?, ?it/s]
Convert catalog: 27%|██▋ | 281/1025 [00:00<00:00, 2805.93it/s]
Convert catalog: 55%|█████▍ | 562/1025 [00:00<00:00, 2784.45it/s]
Convert catalog: 82%|████████▏ | 841/1025 [00:00<00:00, 2748.90it/s]
Convert catalog: 100%|██████████| 1025/1025 [00:00<00:00, 2794.91it/s]
+ '[' 2 -eq 2 ']' + root_path=local + region=demo + data_path=local/demo/hypodd + '[' '!' -d local/demo/hypodd ']' + cd local/demo/hypodd + '[' '!' -d HypoDD ']' + cat + cat + ./HypoDD/src/ph2dt/ph2dt ph2dt.inp
Running Hypodd: starting ph2dt (v2.1b - 08/2012)... Mon Oct 19 22:20:34 2026 reading data ... > stations = 15 > events total = 1025 > events selected = 995 > phases = 22578 forming dtimes...
> stations selected = 14 > P-phase pairs total = 276958 > S-phase pairs total = 307371 > outliers = 21129 ( 3 %) > phases at stations not in station list = 0 > phases at distances larger than MAXDIST = 0 > P-phase pairs selected = 246654 ( 89 %) > S-phase pairs selected = 281990 ( 91 %) > weakly linked events = 39 ( 3 %) > linked event pairs = 30537 > average links per pair = 17 > average offset (km) betw. linked events = 2.59763598 > average offset (km) betw. strongly linked events = 2.59763598 > maximum offset (km) betw. strongly linked events = 9.99354076 Done. Mon Oct 19 22:20:35 2026 Output files: dt.ct; event.dat; event.sel; station.sel; ph2dt.log ph2dt parameters were: (minwght,maxdist,maxsep,maxngh,minlnk,minobs,maxobs) 0.00 200.000 10.000 50 8 8 100 starting hypoDD (v2.1beta - 06/15/2016)... Mon Oct 19 22:20:35 2026 INPUT FILES: cross dtime data: catalog dtime data: dt.ct events: event.sel stations: stations.dat OUTPUT FILES: initial locations: hypodd_ct.loc relocated events: hypodd_ct.reloc event pair residuals: hypodd.res station residuals: hypodd.sta source parameters: hypodd.src Relocate all clusters Relocate all events Use local layered 1D model. Reading data ... Mon Oct 19 22:20:35 2026 # events = 995 # stations < maxdist = 15 # stations w/ neg. elevation (set to 0) = 0
+ ./HypoDD/src/hypoDD/hypoDD ct.inp
# catalog P dtimes = 246654 # catalog S dtimes = 281990 # dtimes total = 528644 # events after dtime match = 985 # stations = 14 no clustering performed. RELOCATION OF CLUSTER: 1 Mon Oct 19 22:20:35 2026 ---------------------- Initial trial sources = 985 1D ray tracing.
IT EV CT RMSCT RMSST DX DY DZ DT OS AQ CND
% % ms % ms m m m ms m
1 100 99 160 -11.7 0 245 201 1993 68 0 26 255
2 1 97 96 157 -1.8 324 242 193 1630 64 270 0 251
3 97 95 137 -12.8 324 159 118 675 34 270 7 253
4 2 97 94 133 -2.3 289 162 117 658 33 134 0 253
5 97 93 118 -11.5 289 119 70 587 20 134 11 247
6 3 96 91 116 -1.5 260 114 67 527 20 518 0 241
7 96 91 111 -5.1 260 96 47 554 15 518 7 234
8 4 95 90 110 -0.8 250 92 46 442 15 842 0 234
9 93 68 97 -11.9 250 67 59 435 15 842 1 240
10 5 93 67 91 -5.7 196 68 59 393 15 723 0 237
11 93 65 86 -5.9 196 41 35 243 10 723 1 233
12 6 93 65 85 -1.3 185 42 36 234 10 846 0 232
13 92 63 82 -3.0 185 31 25 164 6 846 1 232
14 7 92 63 81 -1.0 178 32 26 129 6 960 0 232
15 92 62 80 -1.7 178 23 19 104 4 960 4 224
16 8 92 62 80 -0.5 174 23 19 97 4 977 0 206
17 9 89 27 67 -15.1 125 42 36 117 8 1048 0 133
18 88 25 59 -13.2 125 28 22 104 6 1048 2 130
19 88 25 57 -2.9 125 29 22 93 6 1048 1 129
20 10 88 25 57 -0.6 102 30 23 87 6 1073 0 129
21 87 24 54 -3.9 102 21 16 238 4 1073 1 124
22 11 87 24 54 -0.9 97 21 16 233 4 910 0 124
23 86 24 52 -2.9 97 14 11 71 3 910 1 119
24 12 86 23 52 -0.9 94 15 12 59 3 1073 0 119
25 86 22 42 -17.9 94 19 15 52 3 1073 1 117
26 13 86 21 40 -5.8 68 20 16 56 3 1068 0 118
27 14 86 21 38 -6.3 64 15 12 53 3 1058 0 113
28 15 85 20 36 -5.2 61 13 10 43 2 1083 0 113
29 16 85 19 34 -4.4 59 11 9 36 2 1095 0 107 writing out results ... Program hypoDD finished. Mon Oct 19 22:20:57 2026
In [21]:
Copied!
plot_catalog(hypodd_catalog, "HypoDD", region=region, config=config)
plot_catalog(hypodd_catalog, "HypoDD", region=region, config=config)
