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( 

364 shape=(None,), 

365 dtype=num.complex128, 

366 serialize_dtype='<c16', 

367 serialize_as='base64+meta') 

368 

369 def get_ydata(self): 

370 return self.ydata 

371 

372 def get_xdata(self): 

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

374 

375 

376class WaveformPiggybackSubtarget(Object): 

377 piggy_id = Int.T() 

378 

379 _next_piggy_id = 0 

380 

381 @classmethod 

382 def new_piggy_id(cls): 

383 piggy_id = WaveformPiggybackSubtarget._next_piggy_id 

384 WaveformPiggybackSubtarget._next_piggy_id += 1 

385 return piggy_id 

386 

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

388 if piggy_id is None: 

389 piggy_id = self.new_piggy_id() 

390 

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

392 

393 def evaluate( 

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

395 

396 raise NotImplementedError() 

397 

398 

399class WaveformPiggybackSubresult(Object): 

400 piggy_id = Int.T() 

401 

402 

403class WaveformMisfitResult(gf.Result, MisfitResult): 

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

405 

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

407 ''' 

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

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

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

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

412 spectrum_obs = TraceSpectrum.T(optional=True) 

413 spectrum_syn = TraceSpectrum.T(optional=True) 

414 

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

416 tobs_shift = Float.T(optional=True) 

417 tsyn_pick = Timestamp.T(optional=True) 

418 tshift = Float.T(optional=True) 

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

420 

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

422 

423 

424@has_get_plot_classes 

425class WaveformMisfitTarget(gf.Target, MisfitTarget): 

426 flip_norm = Bool.T(default=False) 

427 misfit_config = WaveformMisfitConfig.T() 

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

429 

430 can_bootstrap_weights = True 

431 

432 def __init__(self, **kwargs): 

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

434 MisfitTarget.__init__(self, **kwargs) 

435 self._piggyback_subtargets = [] 

436 

437 def string_id(self): 

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

439 

440 @classmethod 

441 def get_plot_classes(cls): 

442 from . import plot 

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

444 plots.extend(plot.get_plot_classes()) 

445 return plots 

446 

447 def get_combined_weight(self): 

448 if self._combined_weight is None: 

449 w = self.manual_weight 

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

451 w *= analyser.weight 

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

453 return self._combined_weight 

454 

455 def get_taper_params(self, engine, source): 

456 store = engine.get_store(self.store_id) 

457 config = self.misfit_config 

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

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

460 if config.fmin > 0.0: 

461 tfade = 1.0/config.fmin 

462 else: 

463 tfade = 1.0/config.fmax 

464 

465 if config.tfade is None: 

466 tfade_taper = tfade 

467 else: 

468 tfade_taper = config.tfade 

469 

470 return tmin_fit, tmax_fit, tfade, tfade_taper 

471 

472 def get_backazimuth_for_waveform(self): 

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

474 

475 @property 

476 def backazimuth(self): 

477 return self.azimuth - 180. 

478 

479 def get_freqlimits(self): 

480 config = self.misfit_config 

481 

482 return ( 

483 config.fmin/config.ffactor, 

484 config.fmin, config.fmax, 

485 config.fmax*config.ffactor) 

486 

487 def get_pick_shift(self, engine, source): 

488 config = self.misfit_config 

489 tobs = None 

490 tsyn = None 

491 ds = self.get_dataset() 

492 

493 if config.pick_synthetic_traveltime and config.pick_phasename: 

494 store = engine.get_store(self.store_id) 

495 tsyn = source.time + store.t( 

496 config.pick_synthetic_traveltime, source, self) 

497 

498 marker = ds.get_pick( 

499 source.name, 

500 self.codes[:3], 

501 config.pick_phasename) 

502 

503 if marker: 

504 tobs = marker.tmin 

505 

506 return tobs, tsyn 

507 

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

509 

510 if self.misfit_config.fmin > 0: 

511 tinc_obs = 1.0 / self.misfit_config.fmin 

512 else: 

513 tinc_obs = 10.0 / self.misfit_config.fmax 

514 

515 tmin_obs = (math.floor( 

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

517 tmax_obs = (math.ceil( 

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

519 

520 return tmin_obs, tmax_obs 

521 

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

523 

524 tr_syn = tr_syn.pyrocko_trace() 

525 nslc = self.codes 

526 

527 config = self.misfit_config 

528 

529 tmin_fit, tmax_fit, tfade, tfade_taper = \ 

530 self.get_taper_params(engine, source) 

531 

532 ds = self.get_dataset() 

533 

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

535 

536 if None not in (tobs, tsyn): 

537 tobs_shift = tobs - tsyn 

538 else: 

539 if self.only_use_stations_with_picks: 

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

541 logger.debug(err) 

542 raise gf.SeismosizerError(err) 

543 else: 

544 tobs_shift = 0.0 

545 

546 tr_syn.extend( 

547 tmin_fit - tfade * 2.0, 

548 tmax_fit + tfade * 2.0, 

549 fillmethod='repeat') 

550 

551 freqlimits = self.get_freqlimits() 

552 

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

554 syn_resp = None 

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

556 syn_resp = g_differentiation_response_1 

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

558 syn_resp = g_differentiation_response_2 

559 else: 

560 syn_resp = None 

561 

562 tr_syn = tr_syn.transfer( 

563 freqlimits=freqlimits, 

564 tfade=tfade, 

565 transfer_function=syn_resp) 

566 

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

568 

569 tmin_obs, tmax_obs = self.get_cutout_timespan( 

570 tmin_fit+tobs_shift, tmax_fit+tobs_shift, tfade) 

571 

572 try: 

573 tr_obs = ds.get_waveform( 

574 nslc, 

575 quantity=config.quantity, 

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

577 tmin=tmin_fit+tobs_shift-tfade, 

578 tmax=tmax_fit+tobs_shift+tfade, 

579 tfade=tfade, 

580 freqlimits=freqlimits, 

581 deltat=tr_syn.deltat, 

582 cache=True, 

583 backazimuth=self.get_backazimuth_for_waveform()) 

584 

585 if tobs_shift != 0.0: 

586 tr_obs = tr_obs.copy() 

587 tr_obs.shift(-tobs_shift) 

588 

589 mr = misfit( 

590 tr_obs, tr_syn, 

591 taper=trace.CosTaper( 

592 tmin_fit - tfade_taper, 

593 tmin_fit, 

594 tmax_fit, 

595 tmax_fit + tfade_taper), 

596 domain=config.domain, 

597 exponent=config.norm_exponent, 

598 flip=self.flip_norm, 

599 result_mode=self._result_mode, 

600 tautoshift_max=config.tautoshift_max, 

601 autoshift_penalty_max=config.autoshift_penalty_max, 

602 subtargets=self._piggyback_subtargets) 

603 

604 self._piggyback_subtargets = [] 

605 

606 mr.tobs_shift = float(tobs_shift) 

607 mr.tsyn_pick = float_or_none(tsyn) 

608 

609 return mr 

610 

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

612 logger.debug(str(e)) 

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

614 

615 def get_plain_targets(self, engine, source): 

616 d = dict( 

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

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

619 

620 def add_piggyback_subtarget(self, subtarget): 

621 self._piggyback_subtargets.append(subtarget) 

622 

623 

624def misfit( 

625 tr_obs, tr_syn, taper, domain, exponent, tautoshift_max, 

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

627 

628 ''' 

629 Calculate misfit between observed and synthetic trace. 

630 

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

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

633 :param taper: taper applied in timedomain as 

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

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

636 :param exponent: exponent of Lx type norms 

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

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

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

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

641 ``autoshift_penalty_max * normalization_factor * \ 

642tautoshift**2 / tautoshift_max**2`` 

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

644 computed against *tr_syn* rather than *tr_obs* 

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

646 include only misfit and normalization factor in result 

647 

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

649 ''' 

650 

651 trace.assert_same_sampling_rate(tr_obs, tr_syn) 

652 deltat = tr_obs.deltat 

653 tmin, tmax = taper.time_span() 

654 

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

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

657 

658 piggyback_results = [] 

659 for subtarget in subtargets: 

660 piggyback_results.append( 

661 subtarget.evaluate( 

662 tr_proc_obs, trspec_proc_obs, tr_proc_syn, trspec_proc_syn)) 

663 

664 tshift = None 

665 ctr = None 

666 deltat = tr_proc_obs.deltat 

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

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

669 if flip: 

670 b, a = a, b 

671 

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

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

674 

675 if nshift_max == 0: 

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

677 else: 

678 mns = [] 

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

680 if ishift < 0: 

681 a_cut = a[-ishift:] 

682 b_cut = b[:ishift] 

683 elif ishift == 0: 

684 a_cut = a 

685 b_cut = b 

686 elif ishift > 0: 

687 a_cut = a[:-ishift] 

688 b_cut = b[ishift:] 

689 

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

691 

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

693 

694 iarg = num.argmin(ms) 

695 tshift = (iarg-nshift_max)*deltat 

696 

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

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

699 

700 elif domain == 'cc_max_norm': 

701 

702 ctr = trace.correlate( 

703 tr_proc_syn, 

704 tr_proc_obs, 

705 mode='same', 

706 normalization='normal') 

707 

708 tshift, cc_max = ctr.max() 

709 m = 0.5 - 0.5 * cc_max 

710 n = 0.5 

711 

712 elif domain == 'frequency_domain': 

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

714 if flip: 

715 b, a = a, b 

716 

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

718 

719 elif domain == 'log_frequency_domain': 

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

721 if flip: 

722 b, a = a, b 

723 

724 a = num.abs(a) 

725 b = num.abs(b) 

726 

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

728 if eps == 0.0: 

729 eps = 1e-7 

730 

731 a = num.log(a + eps) 

732 b = num.log(b + eps) 

733 

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

735 

736 if result_mode == 'full': 

737 result = WaveformMisfitResult( 

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

739 processed_obs=tr_proc_obs, 

740 processed_syn=tr_proc_syn, 

741 filtered_obs=tr_obs.copy(), 

742 filtered_syn=tr_syn, 

743 spectrum_obs=trspec_proc_obs, 

744 spectrum_syn=trspec_proc_syn, 

745 taper=taper, 

746 tshift=tshift, 

747 cc=ctr) 

748 

749 elif result_mode == 'sparse': 

750 result = WaveformMisfitResult( 

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

752 else: 

753 assert False 

754 

755 result.piggyback_subresults = piggyback_results 

756 

757 return result 

758 

759 

760def _extend_extract(tr, tmin, tmax): 

761 deltat = tr.deltat 

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

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

764 nframe = itmax_frame - itmin_frame + 1 

765 n = tr.data_len() 

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

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

768 itmax_tr = itmin_tr + n 

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

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

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

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

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

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

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

776 tr = tr.copy(data=False) 

777 tr.tmin = itmin_frame * deltat 

778 tr.set_ydata(a) 

779 return tr 

780 

781 

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

783 tr_proc = _extend_extract(tr, tmin, tmax) 

784 tr_proc.taper(taper) 

785 

786 df = None 

787 trspec_proc = None 

788 

789 if domain == 'envelope': 

790 tr_proc = tr_proc.envelope(inplace=False) 

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

792 

793 elif domain == 'absolute': 

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

795 

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

797 ndata = tr_proc.ydata.size 

798 nfft = trace.nextpow2(ndata)