Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/waveform/target.py: 77%

1from __future__ import print_function 

2 

3import logging 

4import math 

5import numpy as num 

6 

7from pyrocko import gf, trace, weeding, util 

8from pyrocko.guts import (Object, String, Float, Bool, Int, StringChoice, 

9 Timestamp, List, Dict) 

10from pyrocko.guts_array import Array 

11 

12from grond.dataset import NotFound 

13from grond.meta import GrondError, store_t, nslcs_to_patterns 

14 

15from ..base import (MisfitConfig, MisfitTarget, MisfitResult, TargetGroup) 

16from grond.meta import has_get_plot_classes 

17 

18from pyrocko import crust2x2 

19from string import Template 

20 

21guts_prefix = 'grond' 

22logger = logging.getLogger('grond.targets.waveform.target') 

23 

24 

25class StoreIDSelectorError(GrondError): 

26 pass 

27 

28 

29class StoreIDSelector(Object): 

30 ''' 

31 Base class for GF store selectors. 

32 

33 GF store selectors can be implemented to select different stores, based on 

34 station location, source location or other characteristics. 

35 ''' 

36 

37 pass 

38 

39 

40class Crust2StoreIDSelector(StoreIDSelector): 

41 ''' 

42 Store ID selector picking CRUST 2.0 model based on event location. 

43 ''' 

44 

45 template = String.T( 

46 help="Template for the GF store ID. For example ``'crust2_${id}'`` " 

47 "where ``'${id}'`` will be replaced with the corresponding CRUST " 

48 "2.0 profile identifier for the source location.") 

49 

50 def get_store_id(self, event, st, cha): 

51 s = Template(self.template) 

52 return s.substitute(id=( 

53 crust2x2.get_profile(event.lat, event.lon)._ident).lower()) 

54 

55 

56class StationDictStoreIDSelector(StoreIDSelector): 

57 ''' 

58 Store ID selector using a manual station to store ID mapping. 

59 ''' 

60 

61 mapping = Dict.T( 

62 String.T(), gf.StringID.T(), 

63 help='Dictionary with station to store ID pairs, keys are NET.STA. ' 

64 "Add a fallback store ID under the key ``'others'``.") 

65 

66 def get_store_id(self, event, st, cha): 

67 try: 

68 store_id = self.mapping['%s.%s' % (st.network, st.station)] 

69 except KeyError: 

70 try: 

71 store_id = self.mapping['others'] 

72 except KeyError: 

73 raise StoreIDSelectorError( 

74 'No store ID found for station "%s.%s".' % ( 

75 st.network, st.station)) 

76 

77 return store_id 

78 

79 

80class DepthRangeToStoreID(Object): 

81 depth_min = Float.T() 

82 depth_max = Float.T() 

83 store_id = gf.StringID.T() 

84 

85 

86class StationDepthStoreIDSelector(StoreIDSelector): 

87 ''' 

88 Store ID selector using a mapping from station depth range to store ID. 

89 ''' 

90 

91 depth_ranges = List.T(DepthRangeToStoreID.T()) 

92 

93 def get_store_id(self, event, st, cha): 

94 for r in self.depth_ranges: 

95 if r.depth_min <= st.depth < r.depth_max: 

96 return r.store_id 

97 

98 raise StoreIDSelectorError( 

99 'No store ID found for station "%s.%s" at %g m depth.' % ( 

100 st.network, st.station, st.depth)) 

101 

102 

103class DomainChoice(StringChoice): 

104 choices = [ 

105 'time_domain', 

106 'frequency_domain', 

107 'log_frequency_domain', 

108 'envelope', 

109 'absolute', 

110 'cc_max_norm'] 

111 

112 

113class WaveformMisfitConfig(MisfitConfig): 

114 quantity = gf.QuantityType.T(default='displacement') 

115 fmin = Float.T(default=0.0, help='minimum frequency of bandpass filter') 

116 fmax = Float.T(help='maximum frequency of bandpass filter') 

117 ffactor = Float.T(default=1.5) 

118 tmin = gf.Timing.T( 

119 optional=True, 

120 help='Start of main time window used for waveform fitting.') 

121 tmax = gf.Timing.T( 

122 optional=True, 

123 help='End of main time window used for waveform fitting.') 

124 tfade = Float.T( 

125 optional=True, 

126 help='Decay time of taper prepended and appended to main time window ' 

127 'used for waveform fitting [s].') 

128 pick_synthetic_traveltime = gf.Timing.T( 

129 optional=True, 

130 help='Synthetic phase arrival definition for alignment of observed ' 

131 'and synthetic traces.') 

132 pick_phasename = String.T( 

133 optional=True, 

134 help='Name of picked phase for alignment of observed and synthetic ' 

135 'traces.') 

136 domain = DomainChoice.T( 

137 default='time_domain', 

138 help='Type of data characteristic to be fitted.\n\nAvailable choices ' 

139 'are: %s' % ', '.join("``'%s'``" % s 

140 for s in DomainChoice.choices)) 

141 norm_exponent = Int.T( 

142 default=2, 

143 help='Exponent to use in norm (1: L1-norm, 2: L2-norm)') 

144 tautoshift_max = Float.T( 

145 default=0.0, 

146 help='If non-zero, allow synthetic and observed traces to be shifted ' 

147 'against each other by up to +/- the given value [s].') 

148 autoshift_penalty_max = Float.T( 

149 default=0.0, 

150 help='If non-zero, a penalty misfit is added for non-zero shift ' 

151 'values.\n\nThe penalty value is computed as ' 

152 '``autoshift_penalty_max * normalization_factor * tautoshift**2 ' 

153 '/ tautoshift_max**2``') 

154 

155 ranges = {} 

156 

157 def get_full_frequency_range(self): 

158 return self.fmin / self.ffactor, self.fmax * self.ffactor 

159 

160 

161def log_exclude(target, reason): 

162 logger.debug('Excluding potential target %s: %s' % ( 

163 target.string_id(), reason)) 

164 

165 

166class WaveformTargetGroup(TargetGroup): 

167 '''Handles seismogram targets or other targets of dynamic ground motion. 

168 ''' 

169 distance_min = Float.T( 

170 optional=True, 

171 help='excludes targets nearer to source, along a great circle') 

172 distance_max = Float.T( 

173 optional=True, 

174 help='excludes targets farther from source, along a great circle') 

175 distance_3d_min = Float.T( 

176 optional=True, 

177 help='excludes targets nearer from source (direct distance)') 

178 distance_3d_max = Float.T( 

179 optional=True, 

180 help='excludes targets farther from source (direct distance)') 

181 depth_min = Float.T( 

182 optional=True, 

183 help='excludes targets with smaller depths') 

184 depth_max = Float.T( 

185 optional=True, 

186 help='excludes targets with larger depths') 

187 include = List.T( 

188 String.T(), 

189 optional=True, 

190 help='If not None, list of stations/components to include according ' 

191 'to their STA, NET.STA, NET.STA.LOC, or NET.STA.LOC.CHA codes.') 

192 exclude = List.T( 

193 String.T(), 

194 help='Stations/components to be excluded according to their STA, ' 

195 'NET.STA, NET.STA.LOC, or NET.STA.LOC.CHA codes.') 

196 limit = Int.T(optional=True) 

197 channels = List.T( 

198 String.T(), 

199 optional=True, 

200 help="set channels to include, e.g. ['Z', 'T']") 

201 misfit_config = WaveformMisfitConfig.T() 

202 store_id_selector = StoreIDSelector.T( 

203 optional=True, 

204 help='select GF store based on event-station geometry.') 

205 

206 def get_targets(self, ds, event, default_path='none'): 

207 logger.debug('Selecting waveform targets...') 

208 origin = event 

209 targets = [] 

210 

211 stations = ds.get_stations() 

212 if len(stations) == 0: 

213 logger.warning( 

214 'No stations found to create waveform target group.') 

215 

216 for st in ds.get_stations(): 

217 logger.debug('Selecting waveforms for station %s.%s.%s' % st.nsl()) 

218 for cha in self.channels: 

219 nslc = st.nsl() + (cha,) 

220 

221 logger.debug('Selecting waveforms for %s.%s.%s.%s' % nslc) 

222 

223 if self.store_id_selector: 

224 store_id = self.store_id_selector.get_store_id( 

225 event, st, cha) 

226 else: 

227 store_id = self.store_id 

228 

229 logger.debug('Selecting waveforms for %s.%s.%s.%s' % nslc) 

230 

231 target = WaveformMisfitTarget( 

232 quantity='displacement', 

233 codes=nslc, 

234 lat=st.lat, 

235 lon=st.lon, 

236 north_shift=st.north_shift, 

237 east_shift=st.east_shift, 

238 depth=st.depth, 

239 interpolation=self.interpolation, 

240 store_id=store_id, 

241 misfit_config=self.misfit_config, 

242 manual_weight=self.weight, 

243 normalisation_family=self.normalisation_family, 

244 path=self.path or default_path) 

245 

246 if ds.is_blacklisted(nslc): 

247 log_exclude(target, 'excluded by dataset') 

248 continue 

249 

250 if util.match_nslc( 

251 nslcs_to_patterns(self.exclude), nslc): 

252 log_exclude(target, 'excluded by target group') 

253 continue 

254 

255 if self.include is not None and not util.match_nslc( 

256 nslcs_to_patterns(self.include), nslc): 

257 log_exclude(target, 'excluded by target group') 

258 continue 

259 

260 if self.distance_min is not None and \ 

261 target.distance_to(origin) < self.distance_min: 

262 log_exclude(target, 'distance < distance_min') 

263 continue 

264 

265 if self.distance_max is not None and \ 

266 target.distance_to(origin) > self.distance_max: 

267 log_exclude(target, 'distance > distance_max') 

268 continue 

269 

270 if self.distance_3d_min is not None and \ 

271 target.distance_3d_to(origin) < self.distance_3d_min: 

272 log_exclude(target, 'distance_3d < distance_3d_min') 

273 continue 

274 

275 if self.distance_3d_max is not None and \ 

276 target.distance_3d_to(origin) > self.distance_3d_max: 

277 log_exclude(target, 'distance_3d > distance_3d_max') 

278 continue 

279 

280 if self.depth_min is not None and \ 

281 target.depth < self.depth_min: 

282 log_exclude(target, 'depth < depth_min') 

283 continue 

284 

285 if self.depth_max is not None and \ 

286 target.depth > self.depth_max: 

287 log_exclude(target, 'depth > depth_max') 

288 continue 

289 

290 azi, _ = target.azibazi_to(origin) 

291 if cha == 'R': 

292 target.azimuth = azi - 180. 

293 target.dip = 0. 

294 elif cha == 'T': 

295 target.azimuth = azi - 90. 

296 target.dip = 0. 

297 elif cha == 'Z': 

298 target.azimuth = 0. 

299 target.dip = -90. 

300 

301 target.set_dataset(ds) 

302 targets.append(target) 

303 

304 if self.limit: 

305 return weed(origin, targets, self.limit)[0] 

306 else: 

307 return targets 

308 

309 

310class TraceSpectrum(Object): 

311 network = String.T() 

312 station = String.T() 

313 location = String.T() 

314 channel = String.T() 

315 deltaf = Float.T(default=1.0) 

316 fmin = Float.T(default=0.0) 

317 ydata = Array.T(shape=(None,), dtype=num.complex128, serialize_as='list') 

318 

319 def get_ydata(self): 

320 return self.ydata 

321 

322 def get_xdata(self): 

323 return self.fmin + num.arange(self.ydata.size) * self.deltaf 

324 

325 

326class WaveformPiggybackSubtarget(Object): 

327 piggy_id = Int.T() 

328 

329 _next_piggy_id = 0 

330 

331 @classmethod 

332 def new_piggy_id(cls): 

333 piggy_id = WaveformPiggybackSubtarget._next_piggy_id 

334 WaveformPiggybackSubtarget._next_piggy_id += 1 

335 return piggy_id 

336 

337 def __init__(self, piggy_id=None, **kwargs): 

338 if piggy_id is None: 

339 piggy_id = self.new_piggy_id() 

340 

341 Object.__init__(self, piggy_id=piggy_id, **kwargs) 

342 

343 def evaluate( 

344 self, tr_proc_obs, trspec_proc_obs, tr_proc_syn, trspec_proc_syn): 

345 

346 raise NotImplementedError() 

347 

348 

349class WaveformPiggybackSubresult(Object): 

350 piggy_id = Int.T() 

351 

352 

353class WaveformMisfitResult(gf.Result, MisfitResult): 

354 '''Carries the observations for a target and corresponding synthetics. 

355 

356 A number of different waveform or phase representations are possible. 

357 ''' 

358 processed_obs = trace.Trace.T(optional=True) 

359 processed_syn = trace.Trace.T(optional=True) 

360 filtered_obs = trace.Trace.T(optional=True) 

361 filtered_syn = trace.Trace.T(optional=True) 

362 spectrum_obs = TraceSpectrum.T(optional=True) 

363 spectrum_syn = TraceSpectrum.T(optional=True) 

364 

365 taper = trace.Taper.T(optional=True) 

366 tobs_shift = Float.T(optional=True) 

367 tsyn_pick = Timestamp.T(optional=True) 

368 tshift = Float.T(optional=True) 

369 cc = trace.Trace.T(optional=True) 

370 

371 piggyback_subresults = List.T(WaveformPiggybackSubresult.T()) 

372 

373 

374@has_get_plot_classes 

375class WaveformMisfitTarget(gf.Target, MisfitTarget): 

376 flip_norm = Bool.T(default=False) 

377 misfit_config = WaveformMisfitConfig.T() 

378 

379 can_bootstrap_weights = True 

380 

381 def __init__(self, **kwargs): 

382 gf.Target.__init__(self, **kwargs) 

383 MisfitTarget.__init__(self, **kwargs) 

384 self._piggyback_subtargets = [] 

385 

386 def string_id(self): 

387 return '.'.join(x for x in (self.path,) + self.codes) 

388 

389 @classmethod 

390 def get_plot_classes(cls): 

391 from . import plot 

392 plots = super(WaveformMisfitTarget, cls).get_plot_classes() 

393 plots.extend(plot.get_plot_classes()) 

394 return plots 

395 

396 def get_combined_weight(self): 

397 if self._combined_weight is None: 

398 w = self.manual_weight 

399 for analyser in self.analyser_results.values(): 

400 w *= analyser.weight 

401 self._combined_weight = num.array([w], dtype=float) 

402 return self._combined_weight 

403 

404 def get_taper_params(self, engine, source): 

405 store = engine.get_store(self.store_id) 

406 config = self.misfit_config 

407 tmin_fit = source.time + store_t(store, config.tmin, source, self) 

408 tmax_fit = source.time + store_t(store, config.tmax, source, self) 

409 if config.fmin > 0.0: 

410 tfade = 1.0/config.fmin 

411 else: 

412 tfade = 1.0/config.fmax 

413 

414 if config.tfade is None: 

415 tfade_taper = tfade 

416 else: 

417 tfade_taper = config.tfade 

418 

419 return tmin_fit, tmax_fit, tfade, tfade_taper 

420 

421 def get_backazimuth_for_waveform(self): 

422 return backazimuth_for_waveform(self.azimuth, self.codes) 

423 

424 @property 

425 def backazimuth(self): 

426 return self.azimuth - 180. 

427 

428 def get_freqlimits(self): 

429 config = self.misfit_config 

430 

431 return ( 

432 config.fmin/config.ffactor, 

433 config.fmin, config.fmax, 

434 config.fmax*config.ffactor) 

435 

436 def get_pick_shift(self, engine, source): 

437 config = self.misfit_config 

438 tobs = None 

439 tsyn = None 

440 ds = self.get_dataset() 

441 

442 if config.pick_synthetic_traveltime and config.pick_phasename: 

443 store = engine.get_store(self.store_id) 

444 tsyn = source.time + store.t( 

445 config.pick_synthetic_traveltime, source, self) 

446 

447 marker = ds.get_pick( 

448 source.name, 

449 self.codes[:3], 

450 config.pick_phasename) 

451 

452 if marker: 

453 tobs = marker.tmin 

454 

455 return tobs, tsyn 

456 

457 def get_cutout_timespan(self, tmin, tmax, tfade): 

458 

459 if self.misfit_config.fmin > 0: 

460 tinc_obs = 1.0 / self.misfit_config.fmin 

461 else: 

462 tinc_obs = 10.0 / self.misfit_config.fmax 

463 

464 tmin_obs = (math.floor( 

465 (tmin - tfade) / tinc_obs) - 1.0) * tinc_obs 

466 tmax_obs = (math.ceil( 

467 (tmax + tfade) / tinc_obs) + 1.0) * tinc_obs 

468 

469 return tmin_obs, tmax_obs 

470 

471 def post_process(self, engine, source, tr_syn): 

472 

473 tr_syn = tr_syn.pyrocko_trace() 

474 nslc = self.codes 

475 

476 config = self.misfit_config 

477 

478 tmin_fit, tmax_fit, tfade, tfade_taper = \ 

479 self.get_taper_params(engine, source) 

480 

481 ds = self.get_dataset() 

482 

483 tobs, tsyn = self.get_pick_shift(engine, source) 

484 if None not in (tobs, tsyn): 

485 tobs_shift = tobs - tsyn 

486 else: 

487 tobs_shift = 0.0 

488 

489 tr_syn.extend( 

490 tmin_fit - tfade * 2.0, 

491 tmax_fit + tfade * 2.0, 

492 fillmethod='repeat') 

493 

494 freqlimits = self.get_freqlimits() 

495 

496 if config.quantity == 'displacement': 

497 syn_resp = None 

498 elif config.quantity == 'velocity': 

499 syn_resp = trace.DifferentiationResponse(1) 

500 elif config.quantity == 'acceleration': 

501 syn_resp = trace.DifferentiationResponse(2) 

502 else: 

503 GrondError('Unsupported quantity: %s' % config.quantity) 

504 

505 tr_syn = tr_syn.transfer( 

506 freqlimits=freqlimits, 

507 tfade=tfade, 

508 transfer_function=syn_resp) 

509 

510 tr_syn.chop(tmin_fit - 2*tfade, tmax_fit + 2*tfade) 

511 

512 tmin_obs, tmax_obs = self.get_cutout_timespan( 

513 tmin_fit+tobs_shift, tmax_fit+tobs_shift, tfade) 

514 

515 try: 

516 tr_obs = ds.get_waveform( 

517 nslc, 

518 quantity=config.quantity, 

519 tinc_cache=1.0/(config.fmin or 0.1*config.fmax), 

520 tmin=tmin_fit+tobs_shift-tfade, 

521 tmax=tmax_fit+tobs_shift+tfade, 

522 tfade=tfade, 

523 freqlimits=freqlimits, 

524 deltat=tr_syn.deltat, 

525 cache=True, 

526 backazimuth=self.get_backazimuth_for_waveform()) 

527 

528 if tobs_shift != 0.0: 

529 tr_obs = tr_obs.copy() 

530 tr_obs.shift(-tobs_shift) 

531 

532 mr = misfit( 

533 tr_obs, tr_syn, 

534 taper=trace.CosTaper( 

535 tmin_fit - tfade_taper, 

536 tmin_fit, 

537 tmax_fit, 

538 tmax_fit + tfade_taper), 

539 domain=config.domain, 

540 exponent=config.norm_exponent, 

541 flip=self.flip_norm, 

542 result_mode=self._result_mode, 

543 tautoshift_max=config.tautoshift_max, 

544 autoshift_penalty_max=config.autoshift_penalty_max, 

545 subtargets=self._piggyback_subtargets) 

546 

547 self._piggyback_subtargets = [] 

548 

549 mr.tobs_shift = float(tobs_shift) 

550 mr.tsyn_pick = float_or_none(tsyn) 

551 

552 return mr 

553 

554 except NotFound as e: 

555 logger.debug(str(e)) 

556 raise gf.SeismosizerError('No waveform data: %s' % str(e)) 

557 

558 def get_plain_targets(self, engine, source): 

559 d = dict( 

560 (k, getattr(self, k)) for k in gf.Target.T.propnames) 

561 return [gf.Target(**d)] 

562 

563 def add_piggyback_subtarget(self, subtarget): 

564 self._piggyback_subtargets.append(subtarget) 

565 

566 

567def misfit( 

568 tr_obs, tr_syn, taper, domain, exponent, tautoshift_max, 

569 autoshift_penalty_max, flip, result_mode='sparse', subtargets=[]): 

570 

571 ''' 

572 Calculate misfit between observed and synthetic trace. 

573 

574 :param tr_obs: observed trace as :py:class:`pyrocko.trace.Trace` 

575 :param tr_syn: synthetic trace as :py:class:`pyrocko.trace.Trace` 

576 :param taper: taper applied in timedomain as 

577 :py:class:`pyrocko.trace.Taper` 

578 :param domain: how to calculate difference, see :py:class:`DomainChoice` 

579 :param exponent: exponent of Lx type norms 

580 :param tautoshift_max: if non-zero, return lowest misfit when traces are 

581 allowed to shift against each other by up to +/- ``tautoshift_max`` 

582 :param autoshift_penalty_max: if non-zero, a penalty misfit is added for 

583 for non-zero shift values. The penalty value is 

584 ``autoshift_penalty_max * normalization_factor * \ 

585tautoshift**2 / tautoshift_max**2`` 

586 :param flip: ``bool``, if set to ``True``, normalization factor is 

587 computed against *tr_syn* rather than *tr_obs* 

588 :param result_mode: ``'full'``, include traces and spectra or ``'sparse'``, 

589 include only misfit and normalization factor in result 

590 

591 :returns: object of type :py:class:`WaveformMisfitResult` 

592 ''' 

593 

594 trace.assert_same_sampling_rate(tr_obs, tr_syn) 

595 deltat = tr_obs.deltat 

596 tmin, tmax = taper.time_span() 

597 

598 tr_proc_obs, trspec_proc_obs = _process(tr_obs, tmin, tmax, taper, domain) 

599 tr_proc_syn, trspec_proc_syn = _process(tr_syn, tmin, tmax, taper, domain) 

600 

601 piggyback_results = [] 

602 for subtarget in subtargets: 

603 piggyback_results.append( 

604 subtarget.evaluate( 

605 tr_proc_obs, trspec_proc_obs, tr_proc_syn, trspec_proc_syn)) 

606 

607 tshift = None 

608 ctr = None 

609 deltat = tr_proc_obs.deltat 

610 if domain in ('time_domain', 'envelope', 'absolute'): 

611 a, b = tr_proc_syn.ydata, tr_proc_obs.ydata 

612 if flip: 

613 b, a = a, b 

614 

615 nshift_max = max(0, min(a.size-1, 

616 int(math.floor(tautoshift_max / deltat)))) 

617 

618 if nshift_max == 0: 

619 m, n = trace.Lx_norm(a, b, norm=exponent) 

620 else: 

621 mns = [] 

622 for ishift in range(-nshift_max, nshift_max+1): 

623 if ishift < 0: 

624 a_cut = a[-ishift:] 

625 b_cut = b[:ishift] 

626 elif ishift == 0: 

627 a_cut = a 

628 b_cut = b 

629 elif ishift > 0: 

630 a_cut = a[:-ishift] 

631 b_cut = b[ishift:] 

632 

633 mns.append(trace.Lx_norm(a_cut, b_cut, norm=exponent)) 

634 

635 ms, ns = num.array(mns).T 

636 

637 iarg = num.argmin(ms) 

638 tshift = (iarg-nshift_max)*deltat 

639 

640 m, n = ms[iarg], ns[iarg] 

641 m += autoshift_penalty_max * n * tshift**2 / tautoshift_max**2 

642 

643 elif domain == 'cc_max_norm': 

644 

645 ctr = trace.correlate( 

646 tr_proc_syn, 

647 tr_proc_obs, 

648 mode='same', 

649 normalization='normal') 

650 

651 tshift, cc_max = ctr.max() 

652 m = 0.5 - 0.5 * cc_max 

653 n = 0.5 

654 

655 elif domain == 'frequency_domain': 

656 a, b = trspec_proc_syn.ydata, trspec_proc_obs.ydata 

657 if flip: 

658 b, a = a, b 

659 

660 m, n = trace.Lx_norm(num.abs(a), num.abs(b), norm=exponent) 

661 

662 elif domain == 'log_frequency_domain': 

663 a, b = trspec_proc_syn.ydata, trspec_proc_obs.ydata 

664 if flip: 

665 b, a = a, b 

666 

667 a = num.abs(a) 

668 b = num.abs(b) 

669 

670 eps = (num.mean(a) + num.mean(b)) * 1e-7 

671 if eps == 0.0: 

672 eps = 1e-7 

673 

674 a = num.log(a + eps) 

675 b = num.log(b + eps) 

676 

677 m, n = trace.Lx_norm(a, b, norm=exponent) 

678 

679 if result_mode == 'full': 

680 result = WaveformMisfitResult( 

681 misfits=num.array([[m, n]], dtype=float), 

682 processed_obs=tr_proc_obs, 

683 processed_syn=tr_proc_syn, 

684 filtered_obs=tr_obs.copy(), 

685 filtered_syn=tr_syn, 

686 spectrum_obs=trspec_proc_obs, 

687 spectrum_syn=trspec_proc_syn, 

688 taper=taper, 

689 tshift=tshift, 

690 cc=ctr) 

691 

692 elif result_mode == 'sparse': 

693 result = WaveformMisfitResult( 

694 misfits=num.array([[m, n]], dtype=float)) 

695 else: 

696 assert False 

697 

698 result.piggyback_subresults = piggyback_results 

699 

700 return result 

701 

702 

703def _extend_extract(tr, tmin, tmax): 

704 deltat = tr.deltat 

705 itmin_frame = int(math.floor(tmin/deltat)) 

706 itmax_frame = int(math.ceil(tmax/deltat)) 

707 nframe = itmax_frame - itmin_frame + 1 

708 n = tr.data_len() 

709 a = num.empty(nframe, dtype=float) 

710 itmin_tr = int(round(tr.tmin / deltat)) 

711 itmax_tr = itmin_tr + n 

712 icut1 = min(max(0, itmin_tr - itmin_frame), nframe) 

713 icut2 = min(max(0, itmax_tr - itmin_frame), nframe) 

714 icut1_tr = min(max(0, icut1 + itmin_frame - itmin_tr), n) 

715 icut2_tr = min(max(0, icut2 + itmin_frame - itmin_tr), n) 

716 a[:icut1] = tr.ydata[0] 

717 a[icut1:icut2] = tr.ydata[icut1_tr:icut2_tr] 

718 a[icut2:] = tr.ydata[-1] 

719 tr = tr.copy(data=False) 

720 tr.tmin = itmin_frame * deltat 

721 tr.set_ydata(a) 

722 return tr 

723 

724 

725def _process(tr, tmin, tmax, taper, domain): 

726 tr_proc = _extend_extract(tr, tmin, tmax) 

727 tr_proc.taper(taper) 

728 

729 df = None 

730 trspec_proc = None 

731 

732 if domain == 'envelope': 

733 tr_proc = tr_proc.envelope(inplace=False) 

734 tr_proc.set_ydata(num.abs(tr_proc.get_ydata())) 

735 

736 elif domain == 'absolute': 

737 tr_proc.set_ydata(num.abs(tr_proc.get_ydata())) 

738 

739 elif domain in ('frequency_domain', 'log_frequency_domain'): 

740 ndata = tr_proc.ydata.size 

741 nfft = trace.nextpow2(ndata) 

742 padded = num.zeros(nfft, dtype=float) 

743 padded[:ndata] = tr_proc.ydata 

744 spectrum = num.fft.rfft(padded) 

745 df = 1.0 / (tr_proc.deltat * nfft) 

746 

747 trspec_proc = TraceSpectrum( 

748 network=tr_proc.network, 

749 station=tr_proc.station, 

750 location=tr_proc.location, 

751 channel=tr_proc.channel, 

752 deltaf=df, 

753 fmin=0.0, 

754 ydata=spectrum) 

755 

756 return tr_proc, trspec_proc 

757 

758 

759def backazimuth_for_waveform(azimuth, nslc): 

760 if nslc[-1] == 'R': 

761 backazimuth = azimuth + 180. 

762 elif nslc[-1] == 'T': 

763 backazimuth = azimuth + 90. 

764 else: 

765 backazimuth = None 

766 

767 return backazimuth 

768 

769 

770def float_or_none(x): 

771 if x is None: 

772 return x 

773 else: 

774 return float(x) 

775 

776 

777def weed(origin, targets, limit, neighborhood=3): 

778 azimuths = num.zeros(len(targets)) 

779 dists = num.zeros(len(targets)) 

780 for i, target in enumerate(targets): 

781 _, azimuths[i] = target.azibazi_to(origin) 

782 dists[i] = target.distance_to(origin) 

783 

784 badnesses = num.ones(len(targets), dtype=float) 

785 deleted, meandists_kept = weeding.weed( 

786 azimuths, dists, badnesses, 

787 nwanted=limit, 

788 neighborhood=neighborhood) 

789 

790 targets_weeded = [ 

791 target for (delete, target) in zip(deleted, targets) if not delete]