Skip to content

Commit de73872

Browse files
committed
fix bugs in back azimuth calculation
1 parent 569332e commit de73872

File tree

3 files changed

+67
-5
lines changed

3 files changed

+67
-5
lines changed

docs/source/tutorial.rst

Lines changed: 64 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,67 @@
11
Tutorial
2-
*******************
2+
========
33

4-
This is tutorial...To be continued...
4+
The data used here should be read in by `obspy`_ as a ``Stream`` or ``Trace``.
55

6+
.. _obspy: https://github.com/obspy/obspy/wiki
7+
8+
1. Signal processing
9+
--------------------
10+
11+
Give observed seismograms, you want to apply signal processing operations to remove the instrument response(stationxml file is required), filter to a certain band, re-sampling and rotate from ``NE`` to ``RT``. You can write your script this way::
12+
13+
from pytomo3d.signal.process import process
14+
from obspy import read, read_inventory
15+
16+
# read in waveform data
17+
obs = read("II.AAK.obs.mseed")
18+
# read in stationxml
19+
inv = read_inventory("II.AAK.xml")
20+
# set up your filter frequency band
21+
pre_filt = [1/150., 1/100., 1/50., 1/40.]
22+
# setup cutting starttime and endtime
23+
starttime = stream[0].stats.starttime + 10 # second
24+
endtime = stream[0].stats.starttime + 3610 # second
25+
new_obs = process(obs, remove_response_flag=True, inventory=inv,
26+
filter_flag=True, pre_filt=pre_filt,
27+
starttime=starttime, endtime=endtime,
28+
resample_flag=True, sampling_rate=1.0,
29+
rotate_flag=True, event_latitude=12.2,
30+
event_longitude=-95.6)
31+
# write out processed stream
32+
new_obs.write("II.AAK.obs.proc.mseed", format="MSEED")
33+
34+
Given an synthetic stream, you want to filter, re-sample and rotate::
35+
36+
new_syn = process(syn, remove_response_flag=False, inventory=inv,
37+
filter_flag=True, pre_filt=pre_filt,
38+
starttime=starttime, endtime=endtime,
39+
resample_flag=True, sampling_rate=1.0,
40+
rotate_flag=True, event_latitude=12.2,
41+
event_longitude=-95.6)
42+
# write out processed stream
43+
new_obs.write("II.AAK.syn.proc.mseed", format="MSEED")
44+
45+
2. Window Selection
46+
-------------------
47+
To make window selections, you need first prepare window selection config dictionary in python::
48+
49+
{
50+
"Z": pyflex.Config,
51+
"R": pyflex.Config,
52+
"T": pyflex.Config
53+
}
54+
55+
And the selection script::
56+
57+
window = window_on_stream(
58+
obsd, synt, config_dict, station=inv, event=event, figure_mode=True,
59+
figure_dir="./figure/", verbose=False)
60+
61+
3. Adjoint Sources
62+
------------------
63+
For adjoint source calculate, the script::
64+
65+
adjsrcs = calcualte_adjsrc_on_stream(
66+
obsd, synt, windows, config, 'multitaper_misfit',
67+
figure_mode=True, figure_dir="./figure")

pytomo3d/adjoint/adjsrc.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -200,7 +200,7 @@ def calculate_baz(elat, elon, slat, slon):
200200
:return: back azimuth
201201
"""
202202

203-
_, baz, _ = gps2DistAzimuth(elat, elon, slat, slon)
203+
_, _, baz = gps2DistAzimuth(elat, elon, slat, slon)
204204

205205
return baz
206206

pytomo3d/signal/rotate.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -315,8 +315,8 @@ def rotate_stream(st, event_latitude, event_longitude,
315315
raise ValueError("Mode %s required inventory(stationxml) "
316316
"information provided" % mode)
317317

318-
_, baz, _ = gps2DistAzimuth(station_latitude, station_longitude,
319-
event_latitude, event_longitude)
318+
_, _, baz = gps2DistAzimuth(event_latitude, event_longitude,
319+
station_latitude, station_longitude)
320320

321321
components = [tr.stats.channel[-1] for tr in st]
322322

0 commit comments

Comments
 (0)