Skip to content

Commit

Permalink
More verbose logs in bulk_station_orientations.py
Browse files Browse the repository at this point in the history
  • Loading branch information
geojunky committed Dec 13, 2022
1 parent 1aebe58 commit 74d5560
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 16 deletions.
13 changes: 6 additions & 7 deletions seismic/rf_station_orientations.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,11 +98,9 @@ def analyze_station_orientations(ned, curation_opts=DEFAULT_CURATION_OPTS,
if len(ned) == 0:
return {}

logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

# Determine limiting date range per station
results = defaultdict(dict)
full_code = None
for sta, db_evid in ned.by_station():
start_time = end_time = None
full_code = '.'.join([ned.network, sta])
Expand All @@ -121,6 +119,9 @@ def analyze_station_orientations(ned, curation_opts=DEFAULT_CURATION_OPTS,
results[full_code]['date_range'] = [str(start_time), str(end_time)]
# end for

logger = logging.getLogger(__name__ + ':' + full_code)
logger.setLevel(logging.INFO)

evids_orig = set([evid for _, evid, _ in ned])

# Trim streams to time window
Expand Down Expand Up @@ -159,7 +160,7 @@ def analyze_station_orientations(ned, curation_opts=DEFAULT_CURATION_OPTS,
mask = np.isfinite(y)
n_valid = np.sum(mask)
if n_valid < 5:
logger.info('{}: Insufficient data ({} valid points)'.format(sta, n_valid))
logger.info('Insufficient data ({} valid points)'.format(n_valid))
continue
# end if
x_valid = x[mask]
Expand All @@ -182,9 +183,7 @@ def analyze_station_orientations(ned, curation_opts=DEFAULT_CURATION_OPTS,
angle_max = angles_fine[np.argmax(y_fitted)]
code = '.'.join([ned.network, sta])
results[code]['azimuth_correction'] = angle_max
logger.info('{}: {:2.3f}°, stddev {:2.3f}° (N = {:3d})'.format(
sta, angle_max, ph_uncertainty, N))

logger.info('corr {:2.3f}°, stddev {:2.3f}° (N = {:3d})'.format(angle_max, ph_uncertainty, N))

if save_plots_path is not None:
_f = plt.figure(figsize=(16, 9))
Expand Down
25 changes: 16 additions & 9 deletions seismic/swp_station_orientations.py
Original file line number Diff line number Diff line change
Expand Up @@ -290,14 +290,17 @@ def org(T,nameconv):
return s
# end func

def compute_phis(ned, grv_dict):
def compute_phis(ned, grv_dict, logger=None):
nameconv = 1
hrs = 60 * 60 * 4

R1phi = None
R1cc = None
R2phi = None
R2cc = None

discarded = 0
sta = None
# Loop over stations
for i, (sta, evts) in enumerate(ned.by_station()): # ned contains 1 station
nevts = len(evts)
Expand All @@ -323,6 +326,7 @@ def compute_phis(ned, grv_dict):
st = org(st.copy(), nameconv)

if checklen(st, hrs):
discarded += 1
continue
# end if

Expand Down Expand Up @@ -445,6 +449,8 @@ def compute_phis(ned, grv_dict):
#break
# end for
# end for
if(logger):
logger.info("Discarded {}/{} events".format(discarded, len(ned.db_sta[sta])))

return R1cc, R1phi, R2cc, R2phi
# end func
Expand Down Expand Up @@ -543,12 +549,14 @@ def fcalc1(phi,cc,lim,R1cc,R2cc):
# end func


def summary_calculations(R1cc, R1phi, R2cc, R2phi):
def summary_calculations(R1cc, R1phi, R2cc, R2phi, logger=None):
"""
Final Orientation Calculation File
A. Doran and G. Laske
"""

if(logger): logger.info('Calculating summary')

#####################
## ANGLE CALCULATION PARAMETERS
#####################
Expand Down Expand Up @@ -649,9 +657,6 @@ def analyze_station_orientations(ned, grv_dict, save_plots_path=None, ax=None):
if len(ned) == 0:
return {}

logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

# Determine limiting date range per station
results = defaultdict(dict)
full_code = None
Expand All @@ -673,10 +678,13 @@ def analyze_station_orientations(ned, grv_dict, save_plots_path=None, ax=None):
results[full_code]['date_range'] = [str(start_time), str(end_time)]
# end for

logger = logging.getLogger(__name__ + ':' + full_code)
logger.setLevel(logging.INFO)

logger.info('Analysing arrivals')

r1cc, r1phi, r2cc, r2phi = compute_phis(ned, grv_dict)
corr, err, ndata, nevents = summary_calculations(r1cc, r1phi, r2cc, r2phi)
r1cc, r1phi, r2cc, r2phi = compute_phis(ned, grv_dict, logger)
corr, err, ndata, nevents = summary_calculations(r1cc, r1phi, r2cc, r2phi, logger)

corr *= -1 # converting to azimuth correction
while (corr > 180): corr -= 360
Expand All @@ -685,8 +693,7 @@ def analyze_station_orientations(ned, grv_dict, save_plots_path=None, ax=None):
results[full_code]['azimuth_correction'] = corr
results[full_code]['uncertainty'] = err

logger.info('{}: {:2.3f}°, stddev {:2.3f}° (N = {:3d})'.format(
full_code, -corr, err, int(nevents)))
logger.info('corr {:2.3f}°, stddev {:2.3f}° (N = {:3d})'.format(-corr, err, int(nevents)))

CEN = corr
c = np.matlib.repmat([40, 35, 30, 25, 20, 15, 10], r1cc.shape[0], 1)
Expand Down

0 comments on commit 74d5560

Please sign in to comment.