Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rearrange sac header writer #154

Open
wants to merge 2 commits into
base: devel
Choose a base branch
from
Open

rearrange sac header writer #154

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0d1328f
rearrange sac header writer
casarotti Feb 3, 2025
610bf18
update cap_sac.py with the nzmsec conversion
casarotti Feb 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 18 additions & 11 deletions pysep/pysep.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
Python Seismogram Extraction and Processing (PySEP)

Download, pre-process, and organize seismic waveforms, event and station
metadata. Save waveforms as SAC files for use in moment tensor inversion and
metadata. Save waveforms as SAC files for use in moment tensor inversion and
adjoint tomography codes.
"""
import argparse
Expand Down Expand Up @@ -513,7 +513,7 @@ def __init__(self, config_file=None, event_selection="default",
self.demean = bool(demean)
self.detrend = bool(detrend)
self.taper_percentage = taper_percentage
self.rotate = rotate
self.rotate = rotate
self.remove_response = bool(remove_response)
self.output_unit = output_unit
self.water_level = water_level
Expand Down Expand Up @@ -1065,7 +1065,7 @@ def _bulk_query_waveforms_from_client(self):
bulk.append((net.code, sta.code, self.locations,
self.channels, t1, t2))

# Catch edge case where len(bulk)==0 which will cause ObsPy to fail
# Catch edge case where len(bulk)==0 which will cause ObsPy to fail
assert bulk, (
f"station curtailing has removed any stations to query data for. "
f"please check your `distance` and `azimuth` curtailing criteria "
Expand Down Expand Up @@ -1255,7 +1255,7 @@ def preprocess(self):
logger.info(f"applying amplitude scale factor: {self.scale_factor}")
for tr in st_out:
tr.data = tr.data * self.scale_factor
tr.stats.sac["scale"] = self.scale_factor
#tr.stats.sac["scale"] = self.scale_factor
Copy link
Member

@bch0w bch0w Feb 10, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can probably remove this

if self.client == "LLNL":
# This won't do anything if we don't have any 'LL' network codes
st_out = scale_llnl_waveform_amplitudes(st_out)
Expand Down Expand Up @@ -1554,7 +1554,7 @@ def write(self, write_files=None, _return_filenames=False, _subset=None,
logger.info("writing stations file")
logger.debug(fid)
write_pysep_station_file(
self.inv, self.event, fid,
self.inv, self.event, fid,
order_station_list_by=order_station_list_by
)

Expand Down Expand Up @@ -1660,8 +1660,8 @@ def write_config(self, fid=None, overwrite=False):
:type fid: str
:param fid: name of the file to write. defaults to config.yaml
:type overwrite: bool
:param overwrite: if True and `fid` already exists, save a new config
file with the same name, overwriting the old file. if False
:param overwrite: if True and `fid` already exists, save a new config
file with the same name, overwriting the old file. if False
(default), throws a warning if encountering existing `fid` and does
not write config file
"""
Expand Down Expand Up @@ -1774,7 +1774,7 @@ def _set_log_file(self, mode):

if not save_log:
return

# Make the output directory that the log file will be saved in
if not os.path.exists(self._output_dir):
os.makedirs(self._output_dir)
Expand Down Expand Up @@ -1859,13 +1859,20 @@ def run(self, event=None, inv=None, st=None, **kwargs):
# user wants it, it should have been written out
del self.st_raw

## Waveform preprocessing and standardization
self.st = self.preprocess()

self.st = append_sac_headers(self.st, self.event, self.inv)
if self.taup_model is not None:
self.st = format_sac_header_w_taup_traveltimes(self.st,
self.st = format_sac_header_w_taup_traveltimes(self.st,
self.taup_model,
self.phase_list)
# Waveform preprocessing and standardization
self.st = self.preprocess()

## Waveform preprocessing and standardization
#self.st = self.preprocess()
if self.scale_factor:
for tr in self.st:
tr.stats.sac["scale"] = self.scale_factor

# Rotation to various orientations. The output stream will have ALL
# components, both rotated and non-rotated
Expand Down
10 changes: 6 additions & 4 deletions pysep/utils/cap_sac.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,8 @@ def _append_sac_headers_trace(tr, event, inv):
"lcalda": 1, # 1 if DIST, AZ, BAZ, GCARC to be calc'd from metadata
}

sac_header["e"]=sac_header["e"]+sac_header["b"]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it might be better to construct these explicitly before the dict construction, can we do something like:

    otime = event.preferred_origin().time
    begin_time = tr.stats.starttime - otime
    end_time = begin_time + (tr.stats.npts * tr.stats.delta)

    sac_header = {
        "iztype": 9,  # Ref time equivalence, IB (9): Begin time
        "b": begin_time,
        "e": end_time,
        "evla": event.preferred_origin().latitude,
        "evlo": event.preferred_origin().longitude,
        "stla": sta.latitude,
        "stlo": sta.longitude,
        "stel": sta.elevation / 1E3,  # elevation in km
        "kevnm": format_event_tag_legacy(event),  # only take date code
        "nzyear": otime.year,
        "nzjday": otime.julday,
        "nzhour": otime.hour,
        "nzmin": otime.minute,
        "nzsec": otime.second,
        "nzmsec": int(f"{otime.microsecond:0>6}"[:3]),  # micros->ms see #152
        "dist": dist_km,
        "az": az,  # degrees
        "baz": baz,  # degrees
        "gcarc": dist_deg,  # degrees
        "lpspol": 0,  # 1 if left-hand polarity (usually no in passive seis)
        "lcalda": 1,  # 1 if DIST, AZ, BAZ, GCARC to be calc'd from metadata
    }


# Some Inventory objects will not go all the way to channel, only to station
try:
cha = sta[0]
Expand Down Expand Up @@ -235,7 +237,7 @@ def _append_sac_headers_trace(tr, event, inv):
_warn_about.append(key)
if _warn_about:
logger.warning(f"no SAC header values found for: {_warn_about}")

# Append SAC header and include back azimuth for rotation
tr.stats.sac = sac_header
tr.stats.back_azimuth = baz
Expand Down Expand Up @@ -317,7 +319,7 @@ def _append_sac_headers_cartesian_trace(tr, event, rcv_x, rcv_y):
dist_deg = kilometer2degrees(dist_km) # spherical earth approximation
azimuth = np.degrees(np.arctan2((rcv_x - src_x), (rcv_y - src_y))) % 360
backazimuth = (azimuth - 180) % 360

sac_header = {
"iztype": 9, # Ref time equivalence, IB (9): Begin time
"b": tr.stats.starttime - otime, # begin time
Expand All @@ -341,9 +343,9 @@ def _append_sac_headers_cartesian_trace(tr, event, rcv_x, rcv_y):
"lpspol": 0, # 1 if left-hand polarity (usually no in passive seis)
"lcalda": 1, # 1 if DIST, AZ, BAZ, GCARC to be calc'd from metadata
}
sac_header["e"]=sac_header["e"]+sac_header["b"]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above


tr.stats.sac = sac_header

return tr


Expand Down Expand Up @@ -382,7 +384,7 @@ def format_sac_header_w_taup_traveltimes(st, model="ak135",
for tr in st_out[:]:
net_sta = ".".join(tr.get_id().split(".")[:2]) # NN.SSS

# Missing SAC header values may cause certain or all stations to not
# Missing SAC header values may cause certain or all stations to not
# be present in the `phase_dict`
if net_sta not in phase_dict:
logger.warning(f"{tr.get_id()} not present in TauP arrivals, cant "
Expand Down