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.

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

Conversation

casarotti
Copy link

I think there is a problem with the SAC header construction.

here

pysep/pysep/pysep.py

Line 1862 in 92c4c62

self.st = append_sac_headers(self.st, self.event, self.inv)

pysep constructs the SAC header (before the trim), which means that it calculates the B value from the stream start time
here

pysep/pysep/pysep.py

Line 1868 in 92c4c62

self.st = self.preprocess()

pysep calls the preprocess routine and trims the waves, setting the new stream start time.
Then it writes the SAC files (with the new start time and the previously defined header).

When the SAC file is read back, OBSPY takes into account the start time (defined by the trim) and the B value, and adjusts nzsec, nzmsec.
which will be used later to create the origin time.

So the origin time is modified in the sac files

I modified pysep.py moving the preprocess (where the trim occurs) before the sac header.

the proprocess routine assumed that a sac header was already written but only to write the scale_factor value, so I move this writing out of the routine after the sac header construction

move sac header writer after preprocess (trim rewrite b, e, and origin time)
reinstate the nzmsec conversion
@@ -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

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

@@ -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
    }

@@ -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

@bch0w
Copy link
Member

bch0w commented Feb 10, 2025

Thank you @casarotti! Great catch on the incorrect SAC header construction, indeed b and e were not constructed correctly, both because of the trimming but also because e did not account for b. I just added two comments on making the changes slightly more explicit if you wouldn't mind implementing. Happy to merge after that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants