Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/phase_pick/target.py: 80%
153 statements
« prev ^ index » next coverage.py v6.5.0, created at 2024-11-27 15:15 +0000
1# https://pyrocko.org/grond - GPLv3
2#
3# The Grond Developers, 21st Century
4import logging
6import numpy as num
8from pyrocko.guts import Float, Tuple, String, Int, List
9from pyrocko import gf, util
10from pyrocko.plot import beachball
12from ..base import MisfitTarget, TargetGroup, MisfitResult
13from grond.meta import has_get_plot_classes, GrondError, nslcs_to_patterns
15guts_prefix = 'grond'
16logger = logging.getLogger('grond.targets.phase_pick.target')
19def log_exclude(target, reason):
20 logger.debug('Excluding potential target %s: %s' % (
21 target.string_id(), reason))
24class PhasePickTargetGroup(TargetGroup):
26 '''
27 Generate targets to fit phase arrival times.
28 '''
30 distance_min = Float.T(optional=True)
31 distance_max = Float.T(optional=True)
32 distance_3d_min = Float.T(optional=True)
33 distance_3d_max = Float.T(optional=True)
34 depth_min = Float.T(optional=True)
35 depth_max = Float.T(optional=True)
36 store_id = gf.StringID.T(optional=True)
37 include = List.T(
38 String.T(),
39 optional=True,
40 help='If not None, list of stations/components to include according '
41 'to their STA, NET.STA, NET.STA.LOC, or NET.STA.LOC.CHA codes.')
42 exclude = List.T(
43 String.T(),
44 help='Stations/components to be excluded according to their STA, '
45 'NET.STA, NET.STA.LOC, or NET.STA.LOC.CHA codes.')
46 limit = Int.T(optional=True)
48 pick_synthetic_traveltime = gf.Timing.T(
49 help='Synthetic phase arrival definition.')
50 pick_phasename = String.T(
51 help='Name of phase in pick file.')
53 weight_traveltime = Float.T(default=1.0)
54 weight_polarity = Float.T(default=1.0)
56 def get_targets(self, ds, event, default_path='none'):
57 logger.debug('Selecting phase pick targets...')
58 origin = event
59 targets = []
61 stations = ds.get_stations()
62 if len(stations) == 0:
63 logger.warning(
64 'No stations found to create waveform target group.')
66 for st in stations:
67 nslc = st.nsl() + ('',)
68 target = PhasePickTarget(
69 codes=nslc[:3],
70 lat=st.lat,
71 lon=st.lon,
72 north_shift=st.north_shift,
73 east_shift=st.east_shift,
74 depth=st.depth,
75 store_id=self.store_id,
76 manual_weight=self.weight,
77 normalisation_family=self.normalisation_family,
78 path=self.path or default_path,
79 pick_synthetic_traveltime=self.pick_synthetic_traveltime,
80 pick_phasename=self.pick_phasename)
82 if util.match_nslc(
83 nslcs_to_patterns(self.exclude), nslc):
84 log_exclude(target, 'excluded by target group')
85 continue
87 if self.include is not None and not util.match_nslc(
88 nslcs_to_patterns(self.include), nslc):
89 log_exclude(target, 'excluded by target group')
90 continue
92 if self.distance_min is not None and \
93 target.distance_to(origin) < self.distance_min:
94 log_exclude(target, 'distance < distance_min')
95 continue
97 if self.distance_max is not None and \
98 target.distance_to(origin) > self.distance_max:
99 log_exclude(target, 'distance > distance_max')
100 continue
102 if self.distance_3d_min is not None and \
103 target.distance_3d_to(origin) < self.distance_3d_min:
104 log_exclude(target, 'distance_3d < distance_3d_min')
105 continue
107 if self.distance_3d_max is not None and \
108 target.distance_3d_to(origin) > self.distance_3d_max:
109 log_exclude(target, 'distance_3d > distance_3d_max')
110 continue
112 if self.depth_min is not None and \
113 target.depth < self.depth_min:
114 log_exclude(target, 'depth < depth_min')
115 continue
117 if self.depth_max is not None and \
118 target.depth > self.depth_max:
119 log_exclude(target, 'depth > depth_max')
120 continue
122 marker = ds.get_pick(
123 event.name,
124 target.codes[:3],
125 target.pick_phasename)
127 if not marker:
128 log_exclude(target, 'no pick available')
129 continue
131 target.set_dataset(ds)
132 targets.append(target)
134 return targets
136 def is_gf_store_appropriate(self, store, depth_range):
137 ok = TargetGroup.is_gf_store_appropriate(
138 self, store, depth_range)
139 ok &= self._is_gf_store_appropriate_check_extent(
140 store, depth_range)
141 return ok
144class PhasePickResult(MisfitResult):
145 tobs = Float.T(optional=True)
146 tsyn = Float.T(optional=True)
147 polarity_obs = Float.T(optional=True)
148 polarity_syn = Float.T(optional=True)
149 azimuth = Float.T(optional=True)
150 takeoff_angle = Float.T(optional=True)
153@has_get_plot_classes
154class PhasePickTarget(gf.Location, MisfitTarget):
156 '''
157 Target to fit phase arrival times.
158 '''
160 codes = Tuple.T(
161 3, String.T(),
162 help='network, station, location codes.')
164 store_id = gf.StringID.T(
165 help='ID of Green\'s function store (only used for earth model).')
167 pick_synthetic_traveltime = gf.Timing.T(
168 help='Synthetic phase arrival definition.')
170 pick_phasename = String.T(
171 help='Name of phase in pick file.')
173 can_bootstrap_weights = True
175 weight_traveltime = Float.T(default=1.0)
176 weight_polarity = Float.T(default=1.0)
178 def __init__(self, **kwargs):
179 gf.Location.__init__(self, **kwargs)
180 MisfitTarget.__init__(self, **kwargs)
181 self._tobs_cache = {}
183 @property
184 def nmisfits(self):
185 return 2
187 @classmethod
188 def get_plot_classes(cls):
189 from . import plot
190 plots = super(PhasePickTarget, cls).get_plot_classes()
191 plots.extend(plot.get_plot_classes())
192 return plots
194 def string_id(self):
195 return '.'.join(x for x in (self.path,) + self.codes)
197 def set_dataset(self, ds):
198 MisfitTarget.set_dataset(self, ds)
199 self._obs_cache = {}
201 def get_plain_targets(self, engine, source):
202 return self.prepare_modelling(engine, source, None)
204 def prepare_modelling(self, engine, source, targets):
205 return []
207 def get_times_polarities(self, engine, source):
208 tsyn = None
210 k = source.name
211 if k not in self._obs_cache:
212 ds = self.get_dataset()
213 tobs = None
214 marker = ds.get_pick(
215 source.name,
216 self.codes[:3],
217 self.pick_phasename)
219 if marker:
220 polarity_obs = marker.get_polarity()
221 if polarity_obs not in (1, -1, None):
222 raise GrondError(
223 'Polarity of pick %s.%s.%s.%s must be 1 or -1 or None'
224 % marker.one_nslc())
225 self._obs_cache[k] = marker.tmin, polarity_obs
226 else:
227 self._obs_cache[k] = None
229 tobs, polarity_obs = self._obs_cache[k]
231 store = engine.get_store(self.store_id)
233 use_polarities = self.weight_polarity > 0. and polarity_obs is not None
234 use_traveltimes = self.weight_traveltime > 0. and tobs is not None
236 if use_polarities:
237 traveltime, takeoff_angle = store.t(
238 self.pick_synthetic_traveltime, source, self,
239 attributes=['takeoff_angle'])
240 else:
241 traveltime = store.t(
242 self.pick_synthetic_traveltime, source, self)
243 takeoff_angle = None
245 if use_traveltimes:
246 tsyn = source.time + traveltime
247 else:
248 tsyn = None
250 if use_polarities:
251 mt = source.pyrocko_moment_tensor()
252 azimuth = source.azibazi_to(self)[0]
253 amp = beachball.amplitudes(mt, [azimuth], [takeoff_angle])[0]
254 polarity_syn = -1.0 if amp < 0. else 1.0
255 else:
256 polarity_syn = None
257 azimuth = None
259 return tobs, tsyn, polarity_obs, polarity_syn, azimuth, takeoff_angle
261 def finalize_modelling(
262 self, engine, source, modelling_targets, modelling_results):
264 ds = self.get_dataset() # noqa
266 tobs, tsyn, polarity_obs, polarity_syn, azimuth, takeoff_angle = \
267 self.get_times_polarities(engine, source)
269 misfits = num.full((2, 2), num.nan)
270 misfits[:, 1] = 1.0
272 if self.weight_traveltime > 0. and None not in (tobs, tsyn):
273 misfits[0, 0] = self.weight_traveltime * abs(tobs - tsyn)
275 if self.weight_polarity > 0. \
276 and None not in (polarity_obs, polarity_syn):
277 misfits[1, 0] = self.weight_polarity \
278 * 0.5 * abs(polarity_obs - polarity_syn)
280 result = PhasePickResult(
281 tobs=tobs,
282 tsyn=tsyn,
283 polarity_obs=polarity_obs,
284 polarity_syn=polarity_syn,
285 azimuth=azimuth,
286 takeoff_angle=takeoff_angle,
287 misfits=misfits)
289 return result
292__all__ = '''
293 PhasePickTargetGroup
294 PhasePickTarget
295 PhasePickResult
296'''.split()