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

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import logging 

5import math 

6import numpy as num 

7 

8from pyrocko import gf, trace, weeding, util 

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

10 Timestamp, List, Dict) 

11from pyrocko.guts_array import Array 

12from pyrocko.gui.snuffler.marker import PhaseMarker 

13 

14from grond.dataset import NotFound 

15from grond.meta import GrondError, store_t, nslcs_to_patterns 

16 

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

18from grond.meta import has_get_plot_classes 

19 

20from pyrocko import crust2x2 

21from string import Template 

22 

23guts_prefix = 'grond' 

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

25 

26g_differentiation_response_1 = trace.DifferentiationResponse(1) 

27g_differentiation_response_2 = trace.DifferentiationResponse(2) 

28 

29 

30class StoreIDSelectorError(GrondError): 

31 pass 

32 

33 

34class StoreIDSelector(Object): 

35 ''' 

36 Base class for GF store selectors. 

37 

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

39 station location, source location or other characteristics. 

40 ''' 

41 

42 pass 

43 

44 

45class Crust2StoreIDSelector(StoreIDSelector): 

46 ''' 

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

48 ''' 

49 

50 template = String.T( 

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

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

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

54 

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

56 s = Template(self.template) 

57 return s.substitute(id=( 

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

59 

60 

61class StationDictStoreIDSelector(StoreIDSelector): 

62 ''' 

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

64 ''' 

65 

66 mapping = Dict.T( 

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

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

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

70 

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

72 try: 

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

74 except KeyError: 

75 try: 

76 store_id = self.mapping['others'] 

77 except KeyError: 

78 raise StoreIDSelectorError( 

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

80 st.network, st.station)) 

81 

82 return store_id 

83 

84 

85class DepthRangeToStoreID(Object): 

86 depth_min = Float.T() 

87 depth_max = Float.T() 

88 store_id = gf.StringID.T() 

89 

90 

91class StationDepthStoreIDSelector(StoreIDSelector): 

92 ''' 

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

94 ''' 

95 

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

97 

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

99 for r in self.depth_ranges: 

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

101 return r.store_id 

102 

103 raise StoreIDSelectorError( 

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

105 st.network, st.station, st.depth)) 

106 

107 

108class DomainChoice(StringChoice): 

109 choices = [ 

110 'time_domain', 

111 'frequency_domain', 

112 'log_frequency_domain', 

113 'envelope', 

114 'absolute', 

115 'cc_max_norm'] 

116 

117 

118class WaveformMisfitConfig(MisfitConfig): 

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

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

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

122 ffactor = Float.T(default=1.5) 

123 tmin = gf.Timing.T( 

124 optional=True, 

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

126 tmax = gf.Timing.T( 

127 optional=True, 

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

129 tfade = Float.T( 

130 optional=True, 

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

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

133 pick_synthetic_traveltime = gf.Timing.T( 

134 optional=True, 

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

136 'and synthetic traces.') 

137 pick_phasename = String.T( 

138 optional=True, 

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

140 'traces.') 

141 domain = DomainChoice.T( 

142 default='time_domain', 

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

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

145 for s in DomainChoice.choices)) 

146 norm_exponent = Int.T( 

147 default=2, 

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

149 tautoshift_max = Float.T( 

150 default=0.0, 

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

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

153 autoshift_penalty_max = Float.T( 

154 default=0.0, 

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

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

157 '``autoshift_penalty_max * normalization_factor * tautoshift**2 ' 

158 '/ tautoshift_max**2``') 

159 

160 ranges = {} 

161 

162 def get_full_frequency_range(self): 

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

164 

165 

166def log_exclude(target, reason): 

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

168 target.string_id(), reason)) 

169 

170 

171class WaveformTargetGroup(TargetGroup): 

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

173 ''' 

174 distance_min = Float.T( 

175 optional=True, 

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

177 distance_max = Float.T( 

178 optional=True, 

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

180 distance_3d_min = Float.T( 

181 optional=True, 

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

183 distance_3d_max = Float.T( 

184 optional=True, 

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

186 depth_min = Float.T( 

187 optional=True, 

188 help='excludes targets with smaller depths') 

189 depth_max = Float.T( 

190 optional=True, 

191 help='excludes targets with larger depths') 

192 include = List.T( 

193 String.T(), 

194 optional=True, 

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

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

197 exclude = List.T( 

198 String.T(), 

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

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

201 only_use_stations_with_picks = Bool.T( 

202 optional=True, 

203 default=False, 

204 help='Choose if targets without arrival picks should be excluded ' 

205 '(True) or should be used without time shifts (False).') 

206 limit = Int.T(optional=True) 

207 channels = List.T( 

208 String.T(), 

209 optional=True, 

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

211 misfit_config = WaveformMisfitConfig.T() 

212 store_id_selector = StoreIDSelector.T( 

213 optional=True, 

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

215 

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

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

218 origin = event 

219 targets = [] 

220 

221 stations = ds.get_stations() 

222 if len(stations) == 0: 

223 logger.warning( 

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

225 

226 for st in ds.get_stations(): 

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

228 for cha in self.channels: 

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

230 

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

232 

233 if self.store_id_selector: 

234 store_id = self.store_id_selector.get_store_id( 

235 event, st, cha) 

236 else: 

237 store_id = self.store_id 

238 

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

240 

241 if self.misfit_config.quantity in [ 

242 'displacement', 

243 'velocity', 

244 'acceleration']: 

245 

246 quantity = 'displacement' 

247 

248 elif self.misfit_config.quantity in [ 

249 'rotation_displacement', 

250 'rotation_velocity', 

251 'rotation_acceleration']: 

252 

253 quantity = 'rotation_displacement' 

254 

255 else: 

256 quantity = self.misfit_config.quantity 

257 

258 target = WaveformMisfitTarget( 

259 quantity=quantity, 

260 codes=nslc, 

261 lat=st.lat, 

262 lon=st.lon, 

263 north_shift=st.north_shift, 

264 east_shift=st.east_shift, 

265 depth=st.depth, 

266 interpolation=self.interpolation, 

267 store_id=store_id, 

268 misfit_config=self.misfit_config, 

269 manual_weight=self.weight, 

270 normalisation_family=self.normalisation_family, 

271 path=self.path or default_path, 

272 only_use_stations_with_picks=self.only_use_stations_with_picks) # noqa 

273 

274 if ds.is_excluded(nslc): 

275 log_exclude(target, 'excluded by dataset') 

276 continue 

277 

278 if util.match_nslc( 

279 nslcs_to_patterns(self.exclude), nslc): 

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

281 continue 

282 

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

284 nslcs_to_patterns(self.include), nslc): 

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

286 continue 

287 

288 if self.distance_min is not None and \ 

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

290 log_exclude(target, 'distance < distance_min') 

291 continue 

292 

293 if self.distance_max is not None and \ 

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

295 log_exclude(target, 'distance > distance_max') 

296 continue 

297 

298 if self.distance_3d_min is not None and \ 

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

300 log_exclude(target, 'distance_3d < distance_3d_min') 

301 continue 

302 

303 if self.distance_3d_max is not None and \ 

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

305 log_exclude(target, 'distance_3d > distance_3d_max') 

306 continue 

307 

308 if self.depth_min is not None and \ 

309 target.depth < self.depth_min: 

310 log_exclude(target, 'depth < depth_min') 

311 continue 

312 

313 if self.depth_max is not None and \ 

314 target.depth > self.depth_max: 

315 log_exclude(target, 'depth > depth_max') 

316 continue 

317 

318 if self.only_use_stations_with_picks: 

319 all_picks_nsl_ids = set( 

320 p.nslc_ids[0][:3] 

321 for p in ds.pick_markers 

322 if isinstance(p, PhaseMarker)) 

323 if st.nsl() not in all_picks_nsl_ids: 

324 log_exclude(target, 'no pick for nsl') 

325 continue 

326 

327 azi, _ = target.azibazi_to(origin) 

328 if cha == 'R': 

329 target.azimuth = azi - 180. 

330 target.dip = 0. 

331 elif cha == 'T': 

332 target.azimuth = azi - 90. 

333 target.dip = 0. 

334 elif cha == 'Z': 

335 target.azimuth = 0. 

336 target.dip = -90. 

337 

338 target.set_dataset(ds) 

339 targets.append(target) 

340 

341 if self.limit: 

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

343 else: 

344 return targets 

345 

346 def is_gf_store_appropriate(self, store, depth_range): 

347 ok = TargetGroup.is_gf_store_appropriate( 

348 self, store, depth_range) 

349 ok &= self._is_gf_store_appropriate_check_extent( 

350 store, depth_range) 

351 ok &= self._is_gf_store_appropriate_check_sample_rate( 

352 store, depth_range) 

353 return ok 

354 

355 

356class TraceSpectrum(Object): 

357 network = String.T() 

358 station = String.T() 

359 location = String.T() 

360 channel = String.T() 

361 deltaf = Float.T(default=1.0) 

362 fmin = Float.T(default=0.0) 

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

364 

365 def get_ydata(self): 

366 return self.ydata 

367 

368 def get_xdata(self): 

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

370 

371 

372class WaveformPiggybackSubtarget(Object): 

373 piggy_id = Int.T() 

374 

375 _next_piggy_id = 0 

376 

377 @classmethod 

378 def new_piggy_id(cls): 

379 piggy_id = WaveformPiggybackSubtarget._next_piggy_id 

380 WaveformPiggybackSubtarget._next_piggy_id += 1 

381 return piggy_id 

382 

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

384 if piggy_id is None: 

385 piggy_id = self.new_piggy_id() 

386 

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

388 

389 def evaluate( 

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

391 

392 raise NotImplementedError() 

393 

394 

395class WaveformPiggybackSubresult(Object): 

396 piggy_id = Int.T() 

397 

398 

399class WaveformMisfitResult(gf.Result, MisfitResult): 

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

401 

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

403 ''' 

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

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

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

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

408 spectrum_obs = TraceSpectrum.T(optional=True) 

409 spectrum_syn = TraceSpectrum.T(optional=True) 

410 

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

412 tobs_shift = Float.T(optional=True) 

413 tsyn_pick = Timestamp.T(optional=True) 

414 tshift = Float.T(optional=True) 

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

416 

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

418 

419 

420@has_get_plot_classes 

421class WaveformMisfitTarget(gf.Target, MisfitTarget): 

422 flip_norm = Bool.T(default=False) 

423 misfit_config = WaveformMisfitConfig.T() 

424 only_use_stations_with_picks = Bool.T(default=False, optional=True) 

425 

426 can_bootstrap_weights = True 

427 

428 def __init__(self, **kwargs): 

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

430 MisfitTarget.__init__(self, **kwargs) 

431 self._piggyback_subtargets = [] 

432 

433 def string_id(self): 

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

435 

436 @classmethod 

437 def get_plot_classes(cls): 

438 from . import plot 

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

440 plots.extend(plot.get_plot_classes()) 

441 return plots 

442 

443 def get_combined_weight(self): 

444 if self._combined_weight is None: 

445 w = self.manual_weight 

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

447 w *= analyser.weight 

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

449 return self._combined_weight 

450 

451 def get_taper_params(self, engine, source): 

452 store = engine.get_store(self.store_id) 

453 config = self.misfit_config 

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

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

456 if config.fmin > 0.0: 

457 tfade = 1.0/config.fmin 

458 else: 

459 tfade = 1.0/config.fmax 

460 

461 if config.tfade is None: 

462 tfade_taper = tfade 

463 else: 

464 tfade_taper = config.tfade 

465 

466 return tmin_fit, tmax_fit, tfade, tfade_taper 

467 

468 def get_backazimuth_for_waveform(self): 

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

470 

471 @property 

472 def backazimuth(self): 

473 return self.azimuth - 180. 

474 

475 def get_freqlimits(self): 

476 config = self.misfit_config 

477 

478 return ( 

479 config.fmin/config.ffactor, 

480 config.fmin, config.fmax, 

481 config.fmax*config.ffactor) 

482 

483 def get_pick_shift(self, engine, source): 

484 config = self.misfit_config 

485 tobs = None 

486 tsyn = None 

487 ds = self.get_dataset() 

488 

489 if config.pick_synthetic_traveltime and config.pick_phasename: 

490 store = engine.get_store(self.store_id) 

491 tsyn = source.time + store.t( 

492 config.pick_synthetic_traveltime, source, self) 

493 

494 marker = ds.get_pick( 

495 source.name, 

496 self.codes[:3], 

497 config.pick_phasename) 

498 

499 if marker: 

500 tobs = marker.tmin 

501 

502 return tobs, tsyn 

503 

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

505 

506 if self.misfit_config.fmin > 0: 

507 tinc_obs = 1.0 / self.misfit_config.fmin 

508 else: 

509 tinc_obs = 10.0 / self.misfit_config.fmax 

510 

511 tmin_obs = (math.floor( 

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

513 tmax_obs = (math.ceil( 

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

515 

516 return tmin_obs, tmax_obs 

517 

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

519 

520 tr_syn = tr_syn.pyrocko_trace() 

521 nslc = self.codes 

522 

523 config = self.misfit_config 

524 

525 tmin_fit, tmax_fit, tfade, tfade_taper = \ 

526 self.get_taper_params(engine, source) 

527 

528 ds = self.get_dataset() 

529 

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

531 

532 if None not in (tobs, tsyn): 

533 tobs_shift = tobs - tsyn 

534 else: 

535 if self.only_use_stations_with_picks: 

536 err = 'No pick available for target: %s' % self.path 

537 logger.debug(err) 

538 raise gf.SeismosizerError(err) 

539 else: 

540 tobs_shift = 0.0 

541 

542 tr_syn.extend( 

543 tmin_fit - tfade * 2.0, 

544 tmax_fit + tfade * 2.0, 

545 fillmethod='repeat') 

546 

547 freqlimits = self.get_freqlimits() 

548 

549 if config.quantity in ('displacement', 'rotation_displacement'): 

550 syn_resp = None 

551 elif config.quantity in ('velocity', 'rotation_velocity'): 

552 syn_resp = g_differentiation_response_1 

553 elif config.quantity in ('acceleration', 'rotation_acceleration'): 

554 syn_resp = g_differentiation_response_2 

555 else: 

556 syn_resp = None 

557 

558 tr_syn = tr_syn.transfer( 

559 freqlimits=freqlimits, 

560 tfade=tfade, 

561 transfer_function=syn_resp) 

562 

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

564 

565 tmin_obs, tmax_obs = self.get_cutout_timespan( 

566 tmin_fit+tobs_shift, tmax_fit+tobs_shift, tfade) 

567 

568 try: 

569 tr_obs = ds.get_waveform( 

570 nslc, 

571 quantity=config.quantity, 

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

573 tmin=tmin_fit+tobs_shift-tfade, 

574 tmax=tmax_fit+tobs_shift+tfade, 

575 tfade=tfade, 

576 freqlimits=freqlimits, 

577 deltat=tr_syn.deltat, 

578 cache=True, 

579 backazimuth=self.get_backazimuth_for_waveform()) 

580 

581 if tobs_shift != 0.0: 

582 tr_obs = tr_obs.copy() 

583 tr_obs.shift(-tobs_shift) 

584 

585 mr = misfit( 

586 tr_obs, tr_syn, 

587 taper=trace.CosTaper( 

588 tmin_fit - tfade_taper, 

589 tmin_fit, 

590 tmax_fit, 

591 tmax_fit + tfade_taper), 

592 domain=config.domain, 

593 exponent=config.norm_exponent, 

594 flip=self.flip_norm, 

595 result_mode=self._result_mode, 

596 tautoshift_max=config.tautoshift_max, 

597 autoshift_penalty_max=config.autoshift_penalty_max, 

598 subtargets=self._piggyback_subtargets) 

599 

600 self._piggyback_subtargets = [] 

601 

602 mr.tobs_shift = float(tobs_shift) 

603 mr.tsyn_pick = float_or_none(tsyn) 

604 

605 return mr 

606 

607 except (NotFound, util.UnavailableDecimation) as e: 

608 logger.debug(str(e)) 

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

610 

611 def get_plain_targets(self, engine, source): 

612 d = dict( 

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

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

615 

616 def add_piggyback_subtarget(self, subtarget): 

617 self._piggyback_subtargets.append(subtarget) 

618 

619 

620def misfit( 

621 tr_obs, tr_syn, taper, domain, exponent, tautoshift_max, 

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

623 

624 ''' 

625 Calculate misfit between observed and synthetic trace. 

626 

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

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

629 :param taper: taper applied in timedomain as 

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

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

632 :param exponent: exponent of Lx type norms 

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

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

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

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

637 ``autoshift_penalty_max * normalization_factor * \ 

638tautoshift**2 / tautoshift_max**2`` 

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

640 computed against *tr_syn* rather than *tr_obs* 

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

642 include only misfit and normalization factor in result 

643 

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

645 ''' 

646 

647 trace.assert_same_sampling_rate(tr_obs, tr_syn) 

648 deltat = tr_obs.deltat 

649 tmin, tmax = taper.time_span() 

650 

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

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

653 

654 piggyback_results = [] 

655 for subtarget in subtargets: 

656 piggyback_results.append( 

657 subtarget.evaluate( 

658 tr_proc_obs, trspec_proc_obs, tr_proc_syn, trspec_proc_syn)) 

659 

660 tshift = None 

661 ctr = None 

662 deltat = tr_proc_obs.deltat 

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

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

665 if flip: 

666 b, a = a, b 

667 

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

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

670 

671 if nshift_max == 0: 

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

673 else: 

674 mns = [] 

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

676 if ishift < 0: 

677 a_cut = a[-ishift:] 

678 b_cut = b[:ishift] 

679 elif ishift == 0: 

680 a_cut = a 

681 b_cut = b 

682 elif ishift > 0: 

683 a_cut = a[:-ishift] 

684 b_cut = b[ishift:] 

685 

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

687 

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

689 

690 iarg = num.argmin(ms) 

691 tshift = (iarg-nshift_max)*deltat 

692 

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

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

695 

696 elif domain == 'cc_max_norm': 

697 

698 ctr = trace.correlate( 

699 tr_proc_syn, 

700 tr_proc_obs, 

701 mode='same', 

702 normalization='normal') 

703 

704 tshift, cc_max = ctr.max() 

705 m = 0.5 - 0.5 * cc_max 

706 n = 0.5 

707 

708 elif domain == 'frequency_domain': 

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

710 if flip: 

711 b, a = a, b 

712 

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

714 

715 elif domain == 'log_frequency_domain': 

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

717 if flip: 

718 b, a = a, b 

719 

720 a = num.abs(a) 

721 b = num.abs(b) 

722 

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

724 if eps == 0.0: 

725 eps = 1e-7 

726 

727 a = num.log(a + eps) 

728 b = num.log(b + eps) 

729 

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

731 

732 if result_mode == 'full': 

733 result = WaveformMisfitResult( 

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

735 processed_obs=tr_proc_obs, 

736 processed_syn=tr_proc_syn, 

737 filtered_obs=tr_obs.copy(), 

738 filtered_syn=tr_syn, 

739 spectrum_obs=trspec_proc_obs, 

740 spectrum_syn=trspec_proc_syn, 

741 taper=taper, 

742 tshift=tshift, 

743 cc=ctr) 

744 

745 elif result_mode == 'sparse': 

746 result = WaveformMisfitResult( 

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

748 else: 

749 assert False 

750 

751 result.piggyback_subresults = piggyback_results 

752 

753 return result 

754 

755 

756def _extend_extract(tr, tmin, tmax): 

757 deltat = tr.deltat 

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

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

760 nframe = itmax_frame - itmin_frame + 1 

761 n = tr.data_len() 

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

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

764 itmax_tr = itmin_tr + n 

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

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

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

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

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

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

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

772 tr = tr.copy(data=False) 

773 tr.tmin = itmin_frame * deltat 

774 tr.set_ydata(a) 

775 return tr 

776 

777 

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

779 tr_proc = _extend_extract(tr, tmin, tmax) 

780 tr_proc.taper(taper) 

781 

782 df = None 

783 trspec_proc = None 

784 

785 if domain == 'envelope': 

786 tr_proc = tr_proc.envelope(inplace=False) 

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

788 

789 elif domain == 'absolute': 

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

791 

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

793 ndata = tr_proc.ydata.size 

794 nfft = trace.nextpow2(ndata) 

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

796 padded[:ndata] = tr_proc.ydata