Skip to content

Commit

Permalink
Miscellaneous Changes
Browse files Browse the repository at this point in the history
* Added support for filtering out unwanted channels and band-codes from
  data exporter for Salvus
* Added station-names as an extra column in the output from write_gmt_data.py
* Replaced pandas.to_csv call with custom formatted csv writer in bulk_rf_report.py
  • Loading branch information
geojunky committed Feb 24, 2023
1 parent a53103b commit aff76f7
Show file tree
Hide file tree
Showing 3 changed files with 101 additions and 3 deletions.
80 changes: 80 additions & 0 deletions seismic/ASDFdatabase/asdf2salvus.py
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,84 @@ def extract_data_for_event(fds:FederatedASDFDataSet,
@return:
"""

def filter_channels(rows, include_sband=False):
"""
@param rows: output of FederatedASDFDataset.get_stations
@param include_sband:
@return: filtered rows
"""

def is_preferred_component(cha: str):
pcomps = ['Z', 'N', 'E'] # preferred components

for pc in pcomps:
if (pc == cha[-1]): return True
# end for

return False

# end func

itype_dict = defaultdict(lambda: defaultdict( \
lambda: defaultdict(lambda: defaultdict(list))))

for irow, row in enumerate(rows):
net, sta, loc, cha, lon, lat, elev = row

if (cha[0] == 'S' and not include_sband): continue
if (not is_preferred_component(cha)): continue

comp = cha[-1]
itype_dict[net][sta][loc][comp].append([cha, irow])
# end for

bands = ['H', 'B', 'S'] if include_sband else ['H', 'B']
orows = []
for nk, nv in itype_dict.items():
for sk, sv in itype_dict[nk].items():
for lk, lv in itype_dict[nk][sk].items():
for ck, cv in itype_dict[nk][sk][lk].items():

# More than one band-codes found for component
if (len(cv) > 1):
# select H, B or S in order of preference
for band in bands:
found = False
for item in cv:
if (band == item[0][0]):
found = True
orows.append(rows[item[1]])
break
# end if
# end for
if (found): break
# end for
# print(nk,sk,lk, ck, cv)
else:
orows.append(rows[cv[0][1]])
# end if
# end for
# end for
# end for

return orows
# end func

def remove_comments(iinv: Inventory) -> Inventory:
oinv = iinv.copy()
for net in oinv.networks:
net.comments = []
for sta in net.stations:
sta.comments = []
for cha in sta.channels:
cha.comments = []
# end for
# end for
# end for

return oinv
# end func

assert len(event) == 1, 'Invalid event dict {}. Must contain a single event. Aborting..'.format(event)

tag = 'raw_recording'
Expand All @@ -140,6 +218,7 @@ def extract_data_for_event(fds:FederatedASDFDataSet,

st, et = otime - seconds_before, otime + seconds_after
rows = fds.get_stations(st, et)
rows = filter_channels(rows) # filter out unwanted channels and band-codes

receivers_dict = defaultdict(dict)
pbar = tqdm(total=len(rows), desc=ename)
Expand Down Expand Up @@ -170,6 +249,7 @@ def extract_data_for_event(fds:FederatedASDFDataSet,
try:
oinv = inventory.select(network=net, station=sta,
location=loc, channel=cha)
oinv = remove_comments(oinv)
ods.add_stationxml(oinv)

ods.add_waveforms(stream, tag)
Expand Down
17 changes: 16 additions & 1 deletion seismic/receiver_fn/bulk_rf_report.py
Original file line number Diff line number Diff line change
Expand Up @@ -685,7 +685,22 @@ def flatten_dict_list(dict_list):
colnames = ['Longitude', 'Latitude'] + list(itertools.chain.from_iterable(colnames))
df.columns = colnames
df.index.name = 'Station'
df.to_csv(fname)

# write formatted csv file (not currently supported by pandas)
hcolnames = ['Station'] + colnames

hformatter = '{:>15s}, ' + ', '.join(['{:>10s}']*len(hcolnames[1:])) + '\n'
lineformatter = '{:>15s}, ' + ', '.join(['{:>10.3f}']*len(hcolnames[1:])) + '\n'

header = hformatter.format(*hcolnames)
with open(fname, 'w+') as fh:
fh.write(header)
for irow in np.arange(len(df)):
line = lineformatter.format(df.index[irow],
*tuple(df.values[irow, :])).replace('nan', ' ')
fh.write(line)
# end for
# end with
# end for

# gather pdf-names, flatten list and merge pdfs
Expand Down
7 changes: 5 additions & 2 deletions seismic/receiver_fn/write_gmt_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,11 @@ def from_params(params):
# Output network/method locations and sample weight * dataset weight
# Total relative weighting can be used as symbol size on maps
for data in params.method_datasets:
formatted_data = np.array((data.lon, data.lat, data.total_weight, data.val)).T
outfile = params.gmt_loc_file.format(data.name)
with open(outfile, 'w') as fw:
np.savetxt(fw, formatted_data, fmt=['%.6f', '%.6f', '%.2f', '%.6f'], delimiter= ' ', header='lon lat weight depth')
fw.write('# lon lat weight depth station\n')
for i in np.arange(len(data.lon)):
fw.write('%.6f %.6f %.2f %.6f %s\n'%(data.lon[i], data.lat[i],
data.total_weight[i], data.val[i],
data.sta[i]))
print(f"Complete! Data files saved to '{params.gmt_dir}'")

0 comments on commit aff76f7

Please sign in to comment.