Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/phase_pick/target.py: 82%
149 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-06-12 12:25 +0000
1import logging
3import numpy as num
5from pyrocko.guts import Float, Tuple, String, Int, List
6from pyrocko import gf, util
7from pyrocko.plot import beachball
9from ..base import MisfitTarget, TargetGroup, MisfitResult
10from grond.meta import has_get_plot_classes, GrondError, nslcs_to_patterns
12guts_prefix = 'grond'
13logger = logging.getLogger('grond.targets.phase_pick.target')
16def log_exclude(target, reason):
17 logger.debug('Excluding potential target %s: %s' % (
18 target.string_id(), reason))
21class PhasePickTargetGroup(TargetGroup):
23 '''
24 Generate targets to fit phase arrival times.
25 '''
27 distance_min = Float.T(optional=True)
28 distance_max = Float.T(optional=True)
29 distance_3d_min = Float.T(optional=True)
30 distance_3d_max = Float.T(optional=True)
31 depth_min = Float.T(optional=True)
32 depth_max = Float.T(optional=True)
33 store_id = gf.StringID.T(optional=True)
34 include = List.T(
35 String.T(),
36 optional=True,
37 help='If not None, list of stations/components to include according '
38 'to their STA, NET.STA, NET.STA.LOC, or NET.STA.LOC.CHA codes.')
39 exclude = List.T(
40 String.T(),
41 help='Stations/components to be excluded according to their STA, '
42 'NET.STA, NET.STA.LOC, or NET.STA.LOC.CHA codes.')
43 limit = Int.T(optional=True)
45 pick_synthetic_traveltime = gf.Timing.T(
46 help='Synthetic phase arrival definition.')
47 pick_phasename = String.T(
48 help='Name of phase in pick file.')
50 weight_traveltime = Float.T(default=1.0)
51 weight_polarity = Float.T(default=1.0)
53 def get_targets(self, ds, event, default_path='none'):
54 logger.debug('Selecting phase pick targets...')
55 origin = event
56 targets = []
58 stations = ds.get_stations()
59 if len(stations) == 0:
60 logger.warning(
61 'No stations found to create waveform target group.')
63 for st in stations:
64 nslc = st.nsl() + ('',)
65 target = PhasePickTarget(
66 codes=nslc[:3],
67 lat=st.lat,
68 lon=st.lon,
69 north_shift=st.north_shift,
70 east_shift=st.east_shift,
71 depth=st.depth,
72 store_id=self.store_id,
73 manual_weight=self.weight,
74 normalisation_family=self.normalisation_family,
75 path=self.path or default_path,
76 pick_synthetic_traveltime=self.pick_synthetic_traveltime,
77 pick_phasename=self.pick_phasename)
79 if util.match_nslc(
80 nslcs_to_patterns(self.exclude), nslc):
81 log_exclude(target, 'excluded by target group')
82 continue
84 if self.include is not None and not util.match_nslc(
85 nslcs_to_patterns(self.include), nslc):
86 log_exclude(target, 'excluded by target group')
87 continue
89 if self.distance_min is not None and \
90 target.distance_to(origin) < self.distance_min:
91 log_exclude(target, 'distance < distance_min')
92 continue
94 if self.distance_max is not None and \
95 target.distance_to(origin) > self.distance_max:
96 log_exclude(target, 'distance > distance_max')
97 continue
99 if self.distance_3d_min is not None and \
100 target.distance_3d_to(origin) < self.distance_3d_min:
101 log_exclude(target, 'distance_3d < distance_3d_min')
102 continue
104 if self.distance_3d_max is not None and \
105 target.distance_3d_to(origin) > self.distance_3d_max:
106 log_exclude(target, 'distance_3d > distance_3d_max')
107 continue
109 if self.depth_min is not None and \
110 target.depth < self.depth_min:
111 log_exclude(target, 'depth < depth_min')
112 continue
114 if self.depth_max is not None and \
115 target.depth > self.depth_max:
116 log_exclude(target, 'depth > depth_max')
117 continue
119 marker = ds.get_pick(
120 event.name,
121 target.codes[:3],
122 target.pick_phasename)
124 if not marker:
125 log_exclude(target, 'no pick available')
126 continue
128 target.set_dataset(ds)
129 targets.append(target)
131 return targets
134class PhasePickResult(MisfitResult):
135 tobs = Float.T(optional=True)
136 tsyn = Float.T(optional=True)
137 polarity_obs = Float.T(optional=True)
138 polarity_syn = Float.T(optional=True)
139 azimuth = Float.T(optional=True)
140 takeoff_angle = Float.T(optional=True)
143@has_get_plot_classes
144class PhasePickTarget(gf.Location, MisfitTarget):
146 '''
147 Target to fit phase arrival times.
148 '''
150 codes = Tuple.T(
151 3, String.T(),
152 help='network, station, location codes.')
154 store_id = gf.StringID.T(
155 help='ID of Green\'s function store (only used for earth model).')
157 pick_synthetic_traveltime = gf.Timing.T(
158 help='Synthetic phase arrival definition.')
160 pick_phasename = String.T(
161 help='Name of phase in pick file.')
163 can_bootstrap_weights = True
165 weight_traveltime = Float.T(default=1.0)
166 weight_polarity = Float.T(default=1.0)
168 def __init__(self, **kwargs):
169 gf.Location.__init__(self, **kwargs)
170 MisfitTarget.__init__(self, **kwargs)
171 self._tobs_cache = {}
173 @property
174 def nmisfits(self):
175 return 2
177 @classmethod
178 def get_plot_classes(cls):
179 from . import plot
180 plots = super(PhasePickTarget, cls).get_plot_classes()
181 plots.extend(plot.get_plot_classes())
182 return plots
184 def string_id(self):
185 return '.'.join(x for x in (self.path,) + self.codes)
187 def set_dataset(self, ds):
188 MisfitTarget.set_dataset(self, ds)
189 self._obs_cache = {}
191 def get_plain_targets(self, engine, source):
192 return self.prepare_modelling(engine, source, None)
194 def prepare_modelling(self, engine, source, targets):
195 return []
197 def get_times_polarities(self, engine, source):
198 tsyn = None
200 k = source.name
201 if k not in self._obs_cache:
202 ds = self.get_dataset()
203 tobs = None
204 marker = ds.get_pick(
205 source.name,
206 self.codes[:3],
207 self.pick_phasename)
209 if marker:
210 polarity_obs = marker.get_polarity()
211 if polarity_obs not in (1, -1, None):
212 raise GrondError(
213 'Polarity of pick %s.%s.%s.%s must be 1 or -1 or None'
214 % marker.one_nslc())
215 self._obs_cache[k] = marker.tmin, polarity_obs
216 else:
217 self._obs_cache[k] = None
219 tobs, polarity_obs = self._obs_cache[k]
221 store = engine.get_store(self.store_id)
223 use_polarities = self.weight_polarity > 0. and polarity_obs is not None
224 use_traveltimes = self.weight_traveltime > 0. and tobs is not None
226 if use_polarities:
227 traveltime, takeoff_angle = store.t(
228 self.pick_synthetic_traveltime, source, self,
229 attributes=['takeoff_angle'])
230 else:
231 traveltime = store.t(
232 self.pick_synthetic_traveltime, source, self)
233 takeoff_angle = None
235 if use_traveltimes:
236 tsyn = source.time + traveltime
237 else:
238 tsyn = None
240 if use_polarities:
241 mt = source.pyrocko_moment_tensor()
242 azimuth = source.azibazi_to(self)[0]
243 amp = beachball.amplitudes(mt, [azimuth], [takeoff_angle])[0]
244 polarity_syn = -1.0 if amp < 0. else 1.0
245 else:
246 polarity_syn = None
247 azimuth = None
249 return tobs, tsyn, polarity_obs, polarity_syn, azimuth, takeoff_angle
251 def finalize_modelling(
252 self, engine, source, modelling_targets, modelling_results):
254 ds = self.get_dataset() # noqa
256 tobs, tsyn, polarity_obs, polarity_syn, azimuth, takeoff_angle = \
257 self.get_times_polarities(engine, source)
259 misfits = num.full((2, 2), num.nan)
260 misfits[:, 1] = 1.0
262 if self.weight_traveltime > 0. and None not in (tobs, tsyn):
263 misfits[0, 0] = self.weight_traveltime * abs(tobs - tsyn)
265 if self.weight_polarity > 0. \
266 and None not in (polarity_obs, polarity_syn):
267 misfits[1, 0] = self.weight_polarity \
268 * 0.5 * abs(polarity_obs - polarity_syn)
270 result = PhasePickResult(
271 tobs=tobs,
272 tsyn=tsyn,
273 polarity_obs=polarity_obs,
274 polarity_syn=polarity_syn,
275 azimuth=azimuth,
276 takeoff_angle=takeoff_angle,
277 misfits=misfits)
279 return result
282__all__ = '''
283 PhasePickTargetGroup
284 PhasePickTarget
285 PhasePickResult
286'''.split()