Coverage for /usr/local/lib/python3.13/dist-packages/pyrocko/model/event.py: 65%

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

4# ---|P------/S----------~Lg----------

6'''

7Simple representation of a seismic event.

8'''

10import re

11import logging

12import numpy as num

13import hashlib

14import base64

16from pyrocko import util, moment_tensor

18from pyrocko.guts import Float, String, Timestamp, Unicode, \

19 StringPattern, List, Dict, Any, Object

20from .location import Location

22logger = logging.getLogger('pyrocko.model.event')

24guts_prefix = 'pf'

26d2r = num.pi / 180.

27km = 1000.

30def cmp(a, b):

31 return (a > b) - (a < b)

34def ehash(s):

35 return str(base64.urlsafe_b64encode(

36 hashlib.sha1(s.encode('utf8')).digest()).decode('ascii'))

39def float_or_none_to_str(x, prec=9):

40 return 'None' if x is None else '{:.{prec}e}'.format(x, prec=prec)

43class FileParseError(Exception):

44 pass

47class EventExtrasDumpError(Exception):

48 pass

51class EOF(Exception):

52 pass

55class EmptyEvent(Exception):

56 pass

59class Tag(StringPattern):

60 pattern = r'^([A-Za-z][A-Za-z0-9._]{0,128})(:([A-Za-z0-9._-]*))?$'

63def opportunistic_cast(v):

64 try:

65 return int(v)

66 except ValueError:

67 pass

69 try:

70 return float(v)

71 except ValueError:

72 pass

74 return v

77def in_range(vmin, vmax, v):

78 return (vmin is None or (v is not None and vmin <= v)) \

79 and (vmax is None or (v is not None and vmax >= v))

82def mult_none(a, b):

83 if None in (a, b):

84 return None

85 else:

86 return a*b

89class EventFilter(Object):

90 '''

91 Filter to select events by given criteria.

92 '''

94 magnitude_min = Float.T(

95 optional=True,

96 help='Minimum magnitude.')

97 magnitude_max = Float.T(

98 optional=True,

99 help='Maximum magnitude.')

100 depth_min = Float.T(

101 optional=True,

102 help='Minimum event depth [m].')

103 depth_max = Float.T(

104 optional=True,

105 help='Maximum event depth [m].')

106

107 @classmethod

108 def setup_argparse(cls, parser):

109

110 parser.add_argument(

111 '--magnitude-min',

112 dest='magnitude_min',

113 metavar='FLOAT',

114 type=float,

115 help='Minimum magnitude for event filter.')

116

117 parser.add_argument(

118 '--magnitude-max',

119 dest='magnitude_max',

120 metavar='FLOAT',

121 type=float,

122 help='Maximum magnitude for event filter.')

123

124 parser.add_argument(

125 '--depth-min',

126 dest='depth_min_km',

127 metavar='FLOAT',

128 type=float,

129 help='Minimum depth for event filter [km].')

130

131 parser.add_argument(

132 '--depth-max',

133 dest='depth_max_km',

134 metavar='FLOAT',

135 type=float,

136 help='Maximum depth for event filter [km].')

137

138 @classmethod

139 def from_argparse(cls, args):

140 return cls(

141 magnitude_min=args.magnitude_min,

142 magnitude_max=args.magnitude_max,

143 depth_min=mult_none(args.depth_min_km, km),

144 depth_max=mult_none(args.depth_max_km, km))

145

146 def get_filter(self):

147 def filter(ev):

148 return (

149 in_range(self.magnitude_min, self.magnitude_max, ev.magnitude)

150 and in_range(self.depth_min, self.depth_max, ev.depth))

151

152 return filter

153

154

155class Event(Location):

156 '''

157 Representation of a seismic event.

158 '''

159

160 time = Timestamp.T(

161 default=Timestamp.D('1970-01-01 00:00:00'),

162 help='Origin time (UTC system timestamp) [s].')

163 depth = Float.T(

164 optional=True,

165 help='Depth below surface [m].')

166 name = String.T(

167 default='',

168 optional=True,

169 yamlstyle="'",

170 help='Event identifier.')

171 magnitude = Float.T(

172 optional=True,

173 help='Magnitude of the event.')

174 magnitude_type = String.T(

175 optional=True,

176 yamlstyle="'",

177 help='Magnitude type :py:gattr:`magnitude` is given in.')

178 region = Unicode.T(

179 optional=True,

180 yamlstyle="'",

181 help='Source region.')

182 catalog = String.T(

183 optional=True,

184 yamlstyle="'",

185 help='Name of catalog that lists this event.')

186 moment_tensor = moment_tensor.MomentTensor.T(

187 optional=True,

188 help='Moment tensor of the event.')

189 duration = Float.T(

190 optional=True,

191 help='Source duration [s].')

192 tags = List.T(

193 Tag.T(),

194 default=[],

195 help='Auxiliary tags.')

196 extras = Dict.T(

197 String.T(),

198 Any.T(),

199 default={},

200 help='Additional user defined event attributes. The given values must '

201 'be YAML-serializable.')

202

203 def __init__(

204 self, lat=0., lon=0., north_shift=0., east_shift=0., time=0.,

205 name='', depth=None, elevation=None,

206 magnitude=None, magnitude_type=None, region=None, load=None,

207 loadf=None, catalog=None, moment_tensor=None, duration=None,

208 tags=None, extras=None):

209

210 if tags is None:

211 tags = []

212

213 if extras is None:

214 extras = {}

215

216 vals = None

217 if load is not None:

218 vals = Event.oldload(load)

219 elif loadf is not None:

220 vals = Event.oldloadf(loadf)

221

222 if vals:

223 lat, lon, north_shift, east_shift, time, name, depth, magnitude, \

224 magnitude_type, region, catalog, moment_tensor, duration, \

225 tags = vals

226

227 Location.__init__(

228 self, lat=lat, lon=lon,

229 north_shift=north_shift, east_shift=east_shift,

230 time=time, name=name, depth=depth,

231 elevation=elevation,

232 magnitude=magnitude, magnitude_type=magnitude_type,

233 region=region, catalog=catalog,

234 moment_tensor=moment_tensor, duration=duration, tags=tags,

235 extras=extras)

236

237 def tags_as_dict(self):

238 d = {}

239 for tag in self.tags:

240 m = re.match(Tag.pattern, tag)

241 if m:

242 k, v = m.group(1), opportunistic_cast(m.group(3))

243 d[k] = None if m.group(2) == '' else v

244 else:

245 logger.warning('Invalid event tag: %s' % tag)

246

247 return d

248

249 def time_as_string(self):

250 return util.time_to_str(self.time)

251

252 def set_name(self, name):

253 self.name = name

254

255 def olddump(self, filename):

256 file = open(filename, 'w')

257 self.olddumpf(file)

258 file.close()

259

260 def olddumpf(self, file):

261 if self.extras:

262 raise EventExtrasDumpError(

263 'Event user-defined extras attributes cannot be dumped in the '

264 '"basic" event file format. Use '

265 'dump_events(..., format="yaml").')

266

267 file.write('name = %s\n' % self.name)

268 file.write('time = %s\n' % util.time_to_str(self.time))

269

270 if self.lat != 0.0:

271 file.write('latitude = %.12g\n' % self.lat)

272 if self.lon != 0.0:

273 file.write('longitude = %.12g\n' % self.lon)

274

275 if self.north_shift != 0.0:

276 file.write('north_shift = %.12g\n' % self.north_shift)

277 if self.east_shift != 0.0:

278 file.write('east_shift = %.12g\n' % self.east_shift)

279

280 if self.magnitude is not None:

281 file.write('magnitude = %g\n' % self.magnitude)

282 file.write('moment = %g\n' %

283 moment_tensor.magnitude_to_moment(self.magnitude))

284 if self.magnitude_type is not None:

285 file.write('magnitude_type = %s\n' % self.magnitude_type)

286 if self.depth is not None:

287 file.write('depth = %.10g\n' % self.depth)

288 if self.region is not None:

289 file.write('region = %s\n' % self.region)

290 if self.catalog is not None:

291 file.write('catalog = %s\n' % self.catalog)

292 if self.moment_tensor is not None:

293 m = self.moment_tensor.m()

294 sdr1, sdr2 = self.moment_tensor.both_strike_dip_rake()

295 file.write((

296 'mnn = %g\nmee = %g\nmdd = %g\nmne = %g\nmnd = %g\nmed = %g\n'

297 'strike1 = %g\ndip1 = %g\nrake1 = %g\n'

298 'strike2 = %g\ndip2 = %g\nrake2 = %g\n') % (

299 (m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]) +

300 sdr1 + sdr2))

301

302 if self.duration is not None:

303 file.write('duration = %g\n' % self.duration)

304

305 if self.tags:

306 file.write('tags = %s\n' % ', '.join(self.tags))

307

308 @staticmethod

309 def unique(events, deltat=10., group_cmp=(lambda a, b:

310 cmp(a.catalog, b.catalog))):

311 groups = Event.grouped(events, deltat)

312

313 events = []

314 for group in groups:

315 if group:

316 group.sort(group_cmp)

317 events.append(group[-1])

318

319 return events

320

321 @staticmethod

322 def grouped(events, deltat=10.):

323 events = list(events)

324 groups = []

325 for ia, a in enumerate(events):

326 groups.append([])

327 haveit = False

328 for ib, b in enumerate(events[:ia]):

329 if abs(b.time - a.time) < deltat:

330 groups[ib].append(a)

331 haveit = True

332 break

333

334 if not haveit:

335 groups[ia].append(a)

336

337 groups = [g for g in groups if g]

338 groups.sort(key=lambda g: sum(e.time for e in g) // len(g))

339 return groups

340

341 @staticmethod

342 def dump_catalog(events, filename=None, stream=None):

343 if filename is not None:

344 file = open(filename, 'w')

345 else:

346 file = stream

347 try:

348 i = 0

349 for ev in events:

350

351 ev.olddumpf(file)

352

353 file.write('--------------------------------------------\n')

354 i += 1

355

356 finally:

357 if filename is not None:

358 file.close()

359

360 @staticmethod

361 def oldload(filename):

362 with open(filename, 'r') as file:

363 return Event.oldloadf(file)

364

365 @staticmethod

366 def oldloadf(file):

367 d = {}

368 try:

369 for line in file:

370 if line.lstrip().startswith('#'):

371 continue

372

373 toks = line.split(' = ', 1)

374 if len(toks) == 2:

375 k, v = toks[0].strip(), toks[1].strip()

376 if k in ('name', 'region', 'catalog', 'magnitude_type'):

377 d[k] = v

378 if k in (('latitude longitude magnitude depth duration '

379 'north_shift east_shift '

380 'mnn mee mdd mne mnd med strike1 dip1 rake1 '

381 'strike2 dip2 rake2 duration').split()):

382 d[k] = float(v)

383 if k == 'time':

384 d[k] = util.str_to_time(v)

385 if k == 'tags':

386 d[k] = [x.strip() for x in v.split(',')]

387

388 if line.startswith('---'):

389 d['have_separator'] = True

390 break

391

392 except Exception as e:

393 raise FileParseError(e)

394

395 if not d:

396 raise EOF()

397

398 if 'have_separator' in d and len(d) == 1:

399 raise EmptyEvent()

400

401 mt = None

402 m6 = [d[x] for x in 'mnn mee mdd mne mnd med'.split() if x in d]

403 if len(m6) == 6:

404 mt = moment_tensor.MomentTensor(m=moment_tensor.symmat6(*m6))

405 else:

406 sdr = [d[x] for x in 'strike1 dip1 rake1'.split() if x in d]

407 if len(sdr) == 3:

408 moment = 1.0

409 if 'moment' in d:

410 moment = d['moment']

411 elif 'magnitude' in d:

412 moment = moment_tensor.magnitude_to_moment(d['magnitude'])

413

414 mt = moment_tensor.MomentTensor(

415 strike=sdr[0], dip=sdr[1], rake=sdr[2],

416 scalar_moment=moment)

417

418 return (

419 d.get('latitude', 0.0),

420 d.get('longitude', 0.0),

421 d.get('north_shift', 0.0),

422 d.get('east_shift', 0.0),

423 d.get('time', 0.0),

424 d.get('name', ''),

425 d.get('depth', None),

426 d.get('magnitude', None),

427 d.get('magnitude_type', None),

428 d.get('region', None),

429 d.get('catalog', None),

430 mt,

431 d.get('duration', None),

432 d.get('tags', []))

433

434 @staticmethod

435 def load_catalog(filename):

436

437 with open(filename, 'r') as file:

438 try:

439 while True:

440 try:

441 ev = Event(loadf=file)

442 yield ev

443 except EmptyEvent:

444 pass

445

446 except EOF:

447 pass

448

449 def get_hash(self):

450 '''

451 Get a pseudo-unique hash over the main attributes of the event.

452

453 The following attributes are hashed: :py:gattr:`time`,

454 :py:gattr:`~pyrocko.model.location.Location.lat`,

455 :py:gattr:`~pyrocko.model.location.Location.lon`, :py:gattr:`depth`,

456 :py:gattr:`magnitude`, :py:gattr:`catalog`, :py:gattr:`name`,

457 :py:gattr:`region`.

458

459 :returns:

460 URL-safe base64 encoded SHA1 hash.

461 '''

462 e = self

463 if isinstance(e.time, float):

464 stime = util.time_to_str(e.time, format='%Y-%m-%d %H:%M:%S.3FRAC')

465 else:

466 stime = util.time_to_str(e.time, format='%Y-%m-%d %H:%M:%S.6FRAC')

467

468 s = float_or_none_to_str

469

470 to_hash = ', '.join((

471 stime,

472 s(e.lat), s(e.lon), s(e.depth),

473 float_or_none_to_str(e.magnitude, 5),

474 str(e.catalog), str(e.name or ''),

475 str(e.region)))

476

477 return ehash(to_hash)

478

479 def human_str(self):

480 s = [

481 'Latitude [deg]: %g' % self.lat,

482 'Longitude [deg]: %g' % self.lon,

483 'Time [UTC]: %s' % util.time_to_str(self.time)]

484

485 if self.name:

486 s.append('Name: %s' % self.name)

487

488 if self.depth is not None:

489 s.append('Depth [km]: %g' % (self.depth / 1000.))

490

491 if self.magnitude is not None:

492 s.append('Magnitude [%s]: %3.1f' % (

493 self.magnitude_type or 'M?', self.magnitude))

494

495 if self.region:

496 s.append('Region: %s' % self.region)

497

498 if self.catalog:

499 s.append('Catalog: %s' % self.catalog)

500

501 if self.moment_tensor:

502 s.append(str(self.moment_tensor))

503

504 return '\n'.join(s)

505

506 @property

507 def summary(self):

508 return '%s: %s, %s, %s, %s' % (

509 self.__class__.__name__,

510 self.name,

511 util.time_to_str(self.time),

512 '%-3s %3.1f' % (self.magnitude_type or ' ', self.magnitude)

513 if self.magnitude is not None else 'M ---',

514 self.region)

515

516

517def detect_format(filename):

518 with open(filename, 'r') as f:

519 for line in f:

520 line = line.strip()

521 if not line or line.startswith('#') or line.startswith('%'):

522 continue

523 if line.startswith('--- !pf.Event'):

524 return 'yaml'

525 else:

526 return 'basic'

527

528 return 'basic'

529

530

531def load_events(filename, format='detect'):

532 '''

533 Read events file.

534

535 :param filename: name of file as str

536 :param format: file format: ``'detect'``, ``'basic'``, or ``'yaml'``

537 :returns: list of :py:class:`Event` objects

538 '''

539

540 if filename.startswith('http://') or filename.startswith('https://'):

541 import tempfile

542 with tempfile.NamedTemporaryFile() as fp:

543 util.download_file(filename, fp.name)

544 return load_events(fp.name, format=format)

545

546 if format == 'detect':

547 format = detect_format(filename)

548

549 if format == 'yaml':

550 from pyrocko import guts

551 events = [

552 ev for ev in guts.load_all(filename=filename)

553 if isinstance(ev, Event)]

554

555 return events

556 elif format == 'basic':

557 return list(Event.load_catalog(filename))

558 else:

559 from pyrocko.io.io_common import FileLoadError

560 raise FileLoadError('unknown event file format: %s' % format)

561

562

563class OneEventRequired(Exception):

564 pass

565

566

567def load_one_event(filename, format='detect'):

568 events = load_events(filename)

569 if len(events) != 1:

570 raise OneEventRequired(

571 'exactly one event is required in "%s"' % filename)

572

573 return events[0]

574

575

576def dump_events(events, filename=None, stream=None, format='basic'):

577 '''

578 Write events file.

579

580 :param events: list of :py:class:`Event` objects

581 :param filename: name of file as str

582 :param format: file format: ``'basic'``, or ``'yaml'``

583 '''

584

585 if format == 'basic':

586 Event.dump_catalog(events, filename=filename, stream=stream)

587

588 elif format == 'yaml':

589 from pyrocko import guts

590 events = [ev for ev in events if isinstance(ev, Event)]

591 guts.dump_all(events, filename=filename, stream=None)

592

593 else:

594 from pyrocko.io.io_common import FileSaveError

595 raise FileSaveError('unknown event file format: %s' % format)

596

597

598def load_kps_event_list(filename):

599 elist = []

600 f = open(filename, 'r')

601 for line in f:

602 toks = line.split()

603 if len(toks) < 7:

604 continue

605

606 tim = util.to_time_float(util.ctimegm(toks[0]+' '+toks[1]))

607 lat, lon, depth, magnitude = [float(x) for x in toks[2:6]]

608 duration = float(toks[10])

609 region = toks[-1]

610 name = util.gmctime_fn(tim)

611 e = Event(

612 lat, lon, tim,

613 name=name,

614 depth=depth,

615 magnitude=magnitude,

616 duration=duration,

617 region=region)

618

619 elist.append(e)

620

621 f.close()

622 return elist

623

624

625def load_gfz_event_list(filename):

626 from pyrocko import catalog

627 cat = catalog.Geofon()

628

629 elist = []

630 f = open(filename, 'r')

631 for line in f:

632 e = cat.get_event(line.strip())

633 elist.append(e)

634

635 f.close()

636 return elist