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

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

20from .location import Location

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

24guts_prefix = 'pf'

26d2r = num.pi / 180.

29def cmp(a, b):

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

33def ehash(s):

34 return str(base64.urlsafe_b64encode(

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

38def float_or_none_to_str(x, prec=9):

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

42class FileParseError(Exception):

43 pass

46class EventExtrasDumpError(Exception):

47 pass

50class EOF(Exception):

51 pass

54class EmptyEvent(Exception):

55 pass

58class Tag(StringPattern):

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

62def opportunistic_cast(v):

63 try:

64 return int(v)

65 except ValueError:

66 pass

68 try:

69 return float(v)

70 except ValueError:

71 pass

73 return v

76class Event(Location):

77 '''

78 Representation of a seismic event.

79 '''

81 time = Timestamp.T(

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

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

84 depth = Float.T(

85 optional=True,

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

87 name = String.T(

88 default='',

89 optional=True,

90 yamlstyle="'",

91 help='Event identifier.')

92 magnitude = Float.T(

93 optional=True,

94 help='Magnitude of the event.')

95 magnitude_type = String.T(

96 optional=True,

97 yamlstyle="'",

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

99 region = Unicode.T(

100 optional=True,

101 yamlstyle="'",

102 help='Source region.')

103 catalog = String.T(

104 optional=True,

105 yamlstyle="'",

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

107 moment_tensor = moment_tensor.MomentTensor.T(

108 optional=True,

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

110 duration = Float.T(

111 optional=True,

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

113 tags = List.T(

114 Tag.T(),

115 default=[],

116 help='Auxiliary tags.')

117 extras = Dict.T(

118 String.T(),

119 Any.T(),

120 default={},

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

122 'be YAML-serializable.')

123

124 def __init__(

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

126 name='', depth=None, elevation=None,

127 magnitude=None, magnitude_type=None, region=None, load=None,

128 loadf=None, catalog=None, moment_tensor=None, duration=None,

129 tags=None, extras=None):

130

131 if tags is None:

132 tags = []

133

134 if extras is None:

135 extras = {}

136

137 vals = None

138 if load is not None:

139 vals = Event.oldload(load)

140 elif loadf is not None:

141 vals = Event.oldloadf(loadf)

142

143 if vals:

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

145 magnitude_type, region, catalog, moment_tensor, duration, \

146 tags = vals

147

148 Location.__init__(

149 self, lat=lat, lon=lon,

150 north_shift=north_shift, east_shift=east_shift,

151 time=time, name=name, depth=depth,

152 elevation=elevation,

153 magnitude=magnitude, magnitude_type=magnitude_type,

154 region=region, catalog=catalog,

155 moment_tensor=moment_tensor, duration=duration, tags=tags,

156 extras=extras)

157

158 def tags_as_dict(self):

159 d = {}

160 for tag in self.tags:

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

162 if m:

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

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

165 else:

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

167

168 return d

169

170 def time_as_string(self):

171 return util.time_to_str(self.time)

172

173 def set_name(self, name):

174 self.name = name

175

176 def olddump(self, filename):

177 file = open(filename, 'w')

178 self.olddumpf(file)

179 file.close()

180

181 def olddumpf(self, file):

182 if self.extras:

183 raise EventExtrasDumpError(

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

185 '"basic" event file format. Use '

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

187

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

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

190

191 if self.lat != 0.0:

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

193 if self.lon != 0.0:

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

195

196 if self.north_shift != 0.0:

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

198 if self.east_shift != 0.0:

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

200

201 if self.magnitude is not None:

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

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

204 moment_tensor.magnitude_to_moment(self.magnitude))

205 if self.magnitude_type is not None:

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

207 if self.depth is not None:

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

209 if self.region is not None:

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

211 if self.catalog is not None:

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

213 if self.moment_tensor is not None:

214 m = self.moment_tensor.m()

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

216 file.write((

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

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

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

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

221 sdr1 + sdr2))

222

223 if self.duration is not None:

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

225

226 if self.tags:

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

228

229 @staticmethod

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

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

232 groups = Event.grouped(events, deltat)

233

234 events = []

235 for group in groups:

236 if group:

237 group.sort(group_cmp)

238 events.append(group[-1])

239

240 return events

241

242 @staticmethod

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

244 events = list(events)

245 groups = []

246 for ia, a in enumerate(events):

247 groups.append([])

248 haveit = False

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

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

251 groups[ib].append(a)

252 haveit = True

253 break

254

255 if not haveit:

256 groups[ia].append(a)

257

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

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

260 return groups

261

262 @staticmethod

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

264 if filename is not None:

265 file = open(filename, 'w')

266 else:

267 file = stream

268 try:

269 i = 0

270 for ev in events:

271

272 ev.olddumpf(file)

273

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

275 i += 1

276

277 finally:

278 if filename is not None:

279 file.close()

280

281 @staticmethod

282 def oldload(filename):

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

284 return Event.oldloadf(file)

285

286 @staticmethod

287 def oldloadf(file):

288 d = {}

289 try:

290 for line in file:

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

292 continue

293

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

295 if len(toks) == 2:

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

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

298 d[k] = v

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

300 'north_shift east_shift '

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

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

303 d[k] = float(v)

304 if k == 'time':

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

306 if k == 'tags':

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

308

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

310 d['have_separator'] = True

311 break

312

313 except Exception as e:

314 raise FileParseError(e)

315

316 if not d:

317 raise EOF()

318

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

320 raise EmptyEvent()

321

322 mt = None

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

324 if len(m6) == 6:

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

326 else:

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

328 if len(sdr) == 3:

329 moment = 1.0

330 if 'moment' in d:

331 moment = d['moment']

332 elif 'magnitude' in d:

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

334

335 mt = moment_tensor.MomentTensor(

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

337 scalar_moment=moment)

338

339 return (

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

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

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

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

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

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

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

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

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

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

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

351 mt,

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

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

354

355 @staticmethod

356 def load_catalog(filename):

357

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

359 try:

360 while True:

361 try:

362 ev = Event(loadf=file)

363 yield ev

364 except EmptyEvent:

365 pass

366

367 except EOF:

368 pass

369

370 def get_hash(self):

371 '''

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

373

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

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

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

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

378 :py:gattr:`region`.

379

380 :returns:

381 URL-safe base64 encoded SHA1 hash.

382 '''

383 e = self

384 if isinstance(e.time, float):

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

386 else:

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

388

389 s = float_or_none_to_str

390

391 to_hash = ', '.join((

392 stime,

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

394 float_or_none_to_str(e.magnitude, 5),

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

396 str(e.region)))

397

398 return ehash(to_hash)

399

400 def human_str(self):

401 s = [

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

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

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

405

406 if self.name:

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

408

409 if self.depth is not None:

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

411

412 if self.magnitude is not None:

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

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

415

416 if self.region:

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

418

419 if self.catalog:

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

421

422 if self.moment_tensor:

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

424

425 return '\n'.join(s)

426

427 @property

428 def summary(self):

429 return '%s: %s, %s, %s, %s' % (

430 self.__class__.__name__,

431 self.name,

432 util.time_to_str(self.time),

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

434 if self.magnitude else 'M ---',

435 self.region)

436

437

438def detect_format(filename):

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

440 for line in f:

441 line = line.strip()

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

443 continue

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

445 return 'yaml'

446 else:

447 return 'basic'

448

449 return 'basic'

450

451

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

453 '''

454 Read events file.

455

456 :param filename: name of file as str

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

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

459 '''

460

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

462 import tempfile

463 with tempfile.NamedTemporaryFile() as fp:

464 util.download_file(filename, fp.name)

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

466

467 if format == 'detect':

468 format = detect_format(filename)

469

470 if format == 'yaml':

471 from pyrocko import guts

472 events = [

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

474 if isinstance(ev, Event)]

475

476 return events

477 elif format == 'basic':

478 return list(Event.load_catalog(filename))

479 else:

480 from pyrocko.io.io_common import FileLoadError

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

482

483

484class OneEventRequired(Exception):

485 pass

486

487

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

489 events = load_events(filename)

490 if len(events) != 1:

491 raise OneEventRequired(

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

493

494 return events[0]

495

496

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

498 '''

499 Write events file.

500

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

502 :param filename: name of file as str

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

504 '''

505

506 if format == 'basic':

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

508

509 elif format == 'yaml':

510 from pyrocko import guts

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

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

513

514 else:

515 from pyrocko.io.io_common import FileSaveError

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

517

518

519def load_kps_event_list(filename):

520 elist = []

521 f = open(filename, 'r')

522 for line in f:

523 toks = line.split()

524 if len(toks) < 7:

525 continue

526

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

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

529 duration = float(toks[10])

530 region = toks[-1]

531 name = util.gmctime_fn(tim)

532 e = Event(

533 lat, lon, tim,

534 name=name,

535 depth=depth,

536 magnitude=magnitude,

537 duration=duration,

538 region=region)

539

540 elist.append(e)

541

542 f.close()

543 return elist

544

545

546def load_gfz_event_list(filename):

547 from pyrocko import catalog

548 cat = catalog.Geofon()

549

550 elist = []

551 f = open(filename, 'r')

552 for line in f:

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

554 elist.append(e)

555

556 f.close()

557 return elist