Data-Data Misfit

The following code snippet quantifies misfit for an M6.0 in northern Alaska. Waveforms are downloaded from two stations in interior Alaska (TA.POKR and TA.J25K).

We show data-data misfit due to the difficulty of generating synthetic data for a small example problem, so one set of data stands in for the synthetics. After misfit quantification, SPECFEM-ready files are output for an adjoint simulation.

Note

The source receiver map generated by the plot function will not show one of the stations as it assumes st_obs and st_syn are the same station, which is not the case in this example.

from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from pyatoa import Config, Manager
from pyatoa.utils.srcrcv import merge_inventories

# First we download event metadata from IRIS
c = Client("IRIS")
event = c.get_events(eventid="10934221")[0]
origintime = event.preferred_origin().time

# Then we grab station metadata, merging the two inventories
inv = merge_inventories(
        inv_a=c.get_stations(network="TA", station="POKR", location="",
                             channel="BH?", starttime=origintime,
                             endtime=origintime + 300, level="response"),
        inv_b=c.get_stations(network="TA", station="J25K", location="",
                             channel="BH?", starttime=origintime,
                             endtime=origintime + 300, level="response")
)

# Then we grab waveforms for both stations
st_1 = c.get_waveforms(network="TA", station="POKR", location="",
                       channel="BH?", starttime=origintime,
                       endtime=origintime + 300)
st_2 = c.get_waveforms(network="TA", station="J25K", location="",
                       channel="BH?", starttime=origintime,
                       endtime=origintime + 300)

# Now we initiate Pyatoa Config object which controls processing
pyflex_cfg = {"stalta_waterlevel": 0.05,
              "tshift_acceptance_level": 25,
              "cc_acceptance_level": 0.6
              }
cfg = Config(min_period=10, max_period=30,
             unit_output="DISP", rotate_to_rtz=True,
             pyflex_preset="default", adj_src_type="cc_traveltime_misfit",
             st_obs_type="obs", st_syn_type="obs", **pyflex_cfg
             )

# Provide Manager with Config and data and let it do the rest
mgmt = Manager(config=cfg, st_obs=st_1, st_syn=st_2,
               inv=inv, event=event
               )
mgmt.flow()
mgmt.plot()

# Output adjoint sources
mgmt.write_adjsrcs(path="./", write_blanks=True)