A seismology toolkit for Python
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

483 lines
14 KiB

  1. # http://pyrocko.org - GPLv3
  2. #
  3. # The Pyrocko Developers, 21st Century
  4. # ---|P------/S----------~Lg----------
  5. from __future__ import absolute_import, division
  6. import logging
  7. import numpy as num
  8. import hashlib
  9. import base64
  10. from pyrocko import util, moment_tensor
  11. from pyrocko.guts import Float, String, Timestamp, Unicode, \
  12. StringPattern, List, Dict, Any
  13. from .location import Location
  14. logger = logging.getLogger('pyrocko.model.event')
  15. guts_prefix = 'pf'
  16. d2r = num.pi / 180.
  17. def cmp(a, b):
  18. return (a > b) - (a < b)
  19. def ehash(s):
  20. return str(base64.urlsafe_b64encode(
  21. hashlib.sha1(s.encode('utf8')).digest()).decode('ascii'))
  22. def float_or_none_to_str(x, prec=9):
  23. return 'None' if x is None else '{:.{prec}e}'.format(x, prec=prec)
  24. class FileParseError(Exception):
  25. pass
  26. class EventExtrasDumpError(Exception):
  27. pass
  28. class EOF(Exception):
  29. pass
  30. class EmptyEvent(Exception):
  31. pass
  32. class Tag(StringPattern):
  33. pattern = r'^[A-Za-z][A-Za-z0-9._]{0,128}(:[A-Za-z0-9._-]*)?$'
  34. class Event(Location):
  35. '''Seismic event representation
  36. :param lat: latitude of hypocenter (default 0.0)
  37. :param lon: longitude of hypocenter (default 0.0)
  38. :param time: origin time as float in seconds after '1970-01-01 00:00:00
  39. :param name: event identifier as string (optional)
  40. :param depth: source depth (optional)
  41. :param magnitude: magnitude of event (optional)
  42. :param region: source region (optional)
  43. :param catalog: name of catalog that lists this event (optional)
  44. :param moment_tensor: moment tensor as
  45. :py:class:`moment_tensor.MomentTensor` instance (optional)
  46. :param duration: source duration as float (optional)
  47. :param tags: list of tags describing event (optional)
  48. :param extras: dictionary for user defined event attributes (optional).
  49. Keys must be strings, values must be YAML serializable.
  50. '''
  51. time = Timestamp.T(default=util.str_to_time('1970-01-01 00:00:00'))
  52. depth = Float.T(optional=True)
  53. name = String.T(default='', optional=True, yamlstyle="'")
  54. magnitude = Float.T(optional=True)
  55. magnitude_type = String.T(optional=True, yamlstyle="'")
  56. region = Unicode.T(optional=True, yamlstyle="'")
  57. catalog = String.T(optional=True, yamlstyle="'")
  58. moment_tensor = moment_tensor.MomentTensor.T(optional=True)
  59. duration = Float.T(optional=True)
  60. tags = List.T(Tag.T(), default=[])
  61. extras = Dict.T(String.T(), Any.T(), default={})
  62. def __init__(
  63. self, lat=0., lon=0., north_shift=0., east_shift=0., time=0.,
  64. name='', depth=None, elevation=None,
  65. magnitude=None, magnitude_type=None, region=None, load=None,
  66. loadf=None, catalog=None, moment_tensor=None, duration=None,
  67. tags=None, extras=None):
  68. if tags is None:
  69. tags = []
  70. if extras is None:
  71. extras = {}
  72. vals = None
  73. if load is not None:
  74. vals = Event.oldload(load)
  75. elif loadf is not None:
  76. vals = Event.oldloadf(loadf)
  77. if vals:
  78. lat, lon, north_shift, east_shift, time, name, depth, magnitude, \
  79. magnitude_type, region, catalog, moment_tensor, duration, \
  80. tags = vals
  81. Location.__init__(
  82. self, lat=lat, lon=lon,
  83. north_shift=north_shift, east_shift=east_shift,
  84. time=time, name=name, depth=depth,
  85. elevation=elevation,
  86. magnitude=magnitude, magnitude_type=magnitude_type,
  87. region=region, catalog=catalog,
  88. moment_tensor=moment_tensor, duration=duration, tags=tags,
  89. extras=extras)
  90. def time_as_string(self):
  91. return util.time_to_str(self.time)
  92. def set_name(self, name):
  93. self.name = name
  94. def olddump(self, filename):
  95. file = open(filename, 'w')
  96. self.olddumpf(file)
  97. file.close()
  98. def olddumpf(self, file):
  99. if self.extras:
  100. raise EventExtrasDumpError(
  101. 'Event user-defined extras attributes cannot be dumped in the '
  102. '"basic" event file format. Use '
  103. 'dump_events(..., format="yaml").')
  104. file.write('name = %s\n' % self.name)
  105. file.write('time = %s\n' % util.time_to_str(self.time))
  106. if self.lat != 0.0:
  107. file.write('latitude = %.12g\n' % self.lat)
  108. if self.lon != 0.0:
  109. file.write('longitude = %.12g\n' % self.lon)
  110. if self.north_shift != 0.0:
  111. file.write('north_shift = %.12g\n' % self.north_shift)
  112. if self.east_shift != 0.0:
  113. file.write('east_shift = %.12g\n' % self.east_shift)
  114. if self.magnitude is not None:
  115. file.write('magnitude = %g\n' % self.magnitude)
  116. file.write('moment = %g\n' %
  117. moment_tensor.magnitude_to_moment(self.magnitude))
  118. if self.magnitude_type is not None:
  119. file.write('magnitude_type = %s\n' % self.magnitude_type)
  120. if self.depth is not None:
  121. file.write('depth = %.10g\n' % self.depth)
  122. if self.region is not None:
  123. file.write('region = %s\n' % self.region)
  124. if self.catalog is not None:
  125. file.write('catalog = %s\n' % self.catalog)
  126. if self.moment_tensor is not None:
  127. m = self.moment_tensor.m()
  128. sdr1, sdr2 = self.moment_tensor.both_strike_dip_rake()
  129. file.write((
  130. 'mnn = %g\nmee = %g\nmdd = %g\nmne = %g\nmnd = %g\nmed = %g\n'
  131. 'strike1 = %g\ndip1 = %g\nrake1 = %g\n'
  132. 'strike2 = %g\ndip2 = %g\nrake2 = %g\n') % (
  133. (m[0, 0], m[1, 1], m[2, 2], m[0, 1], m[0, 2], m[1, 2]) +
  134. sdr1 + sdr2))
  135. if self.duration is not None:
  136. file.write('duration = %g\n' % self.duration)
  137. if self.tags:
  138. file.write('tags = %s\n' % ', '.join(self.tags))
  139. @staticmethod
  140. def unique(events, deltat=10., group_cmp=(lambda a, b:
  141. cmp(a.catalog, b.catalog))):
  142. groups = Event.grouped(events, deltat)
  143. events = []
  144. for group in groups:
  145. if group:
  146. group.sort(group_cmp)
  147. events.append(group[-1])
  148. return events
  149. @staticmethod
  150. def grouped(events, deltat=10.):
  151. events = list(events)
  152. groups = []
  153. for ia, a in enumerate(events):
  154. groups.append([])
  155. haveit = False
  156. for ib, b in enumerate(events[:ia]):
  157. if abs(b.time - a.time) < deltat:
  158. groups[ib].append(a)
  159. haveit = True
  160. break
  161. if not haveit:
  162. groups[ia].append(a)
  163. groups = [g for g in groups if g]
  164. groups.sort(key=lambda g: sum(e.time for e in g) // len(g))
  165. return groups
  166. @staticmethod
  167. def dump_catalog(events, filename=None, stream=None):
  168. if filename is not None:
  169. file = open(filename, 'w')
  170. else:
  171. file = stream
  172. try:
  173. i = 0
  174. for ev in events:
  175. ev.olddumpf(file)
  176. file.write('--------------------------------------------\n')
  177. i += 1
  178. finally:
  179. if filename is not None:
  180. file.close()
  181. @staticmethod
  182. def oldload(filename):
  183. with open(filename, 'r') as file:
  184. return Event.oldloadf(file)
  185. @staticmethod
  186. def oldloadf(file):
  187. d = {}
  188. try:
  189. for line in file:
  190. if line.lstrip().startswith('#'):
  191. continue
  192. toks = line.split(' = ', 1)
  193. if len(toks) == 2:
  194. k, v = toks[0].strip(), toks[1].strip()
  195. if k in ('name', 'region', 'catalog', 'magnitude_type'):
  196. d[k] = v
  197. if k in (('latitude longitude magnitude depth duration '
  198. 'north_shift east_shift '
  199. 'mnn mee mdd mne mnd med strike1 dip1 rake1 '
  200. 'strike2 dip2 rake2 duration').split()):
  201. d[k] = float(v)
  202. if k == 'time':
  203. d[k] = util.str_to_time(v)
  204. if k == 'tags':
  205. d[k] = [x.strip() for x in v.split(',')]
  206. if line.startswith('---'):
  207. d['have_separator'] = True
  208. break
  209. except Exception as e:
  210. raise FileParseError(e)
  211. if not d:
  212. raise EOF()
  213. if 'have_separator' in d and len(d) == 1:
  214. raise EmptyEvent()
  215. mt = None
  216. m6 = [d[x] for x in 'mnn mee mdd mne mnd med'.split() if x in d]
  217. if len(m6) == 6:
  218. mt = moment_tensor.MomentTensor(m=moment_tensor.symmat6(*m6))
  219. else:
  220. sdr = [d[x] for x in 'strike1 dip1 rake1'.split() if x in d]
  221. if len(sdr) == 3:
  222. moment = 1.0
  223. if 'moment' in d:
  224. moment = d['moment']
  225. elif 'magnitude' in d:
  226. moment = moment_tensor.magnitude_to_moment(d['magnitude'])
  227. mt = moment_tensor.MomentTensor(
  228. strike=sdr[0], dip=sdr[1], rake=sdr[2],
  229. scalar_moment=moment)
  230. return (
  231. d.get('latitude', 0.0),
  232. d.get('longitude', 0.0),
  233. d.get('north_shift', 0.0),
  234. d.get('east_shift', 0.0),
  235. d.get('time', 0.0),
  236. d.get('name', ''),
  237. d.get('depth', None),
  238. d.get('magnitude', None),
  239. d.get('magnitude_type', None),
  240. d.get('region', None),
  241. d.get('catalog', None),
  242. mt,
  243. d.get('duration', None),
  244. d.get('tags', []))
  245. @staticmethod
  246. def load_catalog(filename):
  247. file = open(filename, 'r')
  248. try:
  249. while True:
  250. try:
  251. ev = Event(loadf=file)
  252. yield ev
  253. except EmptyEvent:
  254. pass
  255. except EOF:
  256. pass
  257. file.close()
  258. def get_hash(self):
  259. e = self
  260. if isinstance(e.time, util.hpfloat):
  261. stime = util.time_to_str(e.time, format='%Y-%m-%d %H:%M:%S.6FRAC')
  262. else:
  263. stime = util.time_to_str(e.time, format='%Y-%m-%d %H:%M:%S.3FRAC')
  264. s = float_or_none_to_str
  265. to_hash = ', '.join((
  266. stime,
  267. s(e.lat), s(e.lon), s(e.depth),
  268. float_or_none_to_str(e.magnitude, 5),
  269. str(e.catalog), str(e.name or ''),
  270. str(e.region)))
  271. return ehash(to_hash)
  272. def human_str(self):
  273. s = [
  274. 'Latitude [deg]: %g' % self.lat,
  275. 'Longitude [deg]: %g' % self.lon,
  276. 'Time [UTC]: %s' % util.time_to_str(self.time)]
  277. if self.name:
  278. s.append('Name: %s' % self.name)
  279. if self.depth is not None:
  280. s.append('Depth [km]: %g' % (self.depth / 1000.))
  281. if self.magnitude is not None:
  282. s.append('Magnitude [%s]: %3.1f' % (
  283. self.magnitude_type or 'M?', self.magnitude))
  284. if self.region:
  285. s.append('Region: %s' % self.region)
  286. if self.catalog:
  287. s.append('Catalog: %s' % self.catalog)
  288. if self.moment_tensor:
  289. s.append(str(self.moment_tensor))
  290. return '\n'.join(s)
  291. def detect_format(filename):
  292. with open(filename, 'r') as f:
  293. for line in f:
  294. line = line.strip()
  295. if not line or line.startswith('#') or line.startswith('%'):
  296. continue
  297. if line.startswith('--- !pf.Event'):
  298. return 'yaml'
  299. else:
  300. return 'basic'
  301. return 'basic'
  302. def load_events(filename, format='detect'):
  303. '''Read events file.
  304. :param filename: name of file as str
  305. :param format: file format: ``'detect'``, ``'basic'``, or ``'yaml'``
  306. :returns: list of :py:class:`Event` objects
  307. '''
  308. if format == 'detect':
  309. format = detect_format(filename)
  310. if format == 'yaml':
  311. from pyrocko import guts
  312. events = [
  313. ev for ev in guts.load_all(filename=filename)
  314. if isinstance(ev, Event)]
  315. return events
  316. elif format == 'basic':
  317. return list(Event.load_catalog(filename))
  318. else:
  319. from pyrocko.io.io_common import FileLoadError
  320. raise FileLoadError('unknown event file format: %s' % format)
  321. class OneEventRequired(Exception):
  322. pass
  323. def load_one_event(filename, format='detect'):
  324. events = load_events(filename)
  325. if len(events) != 1:
  326. raise OneEventRequired(
  327. 'exactly one event is required in "%s"' % filename)
  328. return events[0]
  329. def dump_events(events, filename=None, stream=None, format='basic'):
  330. '''Write events file.
  331. :param events: list of :py:class:`Event` objects
  332. :param filename: name of file as str
  333. :param format: file format: ``'basic'``, or ``'yaml'``
  334. '''
  335. if format == 'basic':
  336. Event.dump_catalog(events, filename=filename, stream=stream)
  337. elif format == 'yaml':
  338. from pyrocko import guts
  339. events = [ev for ev in events if isinstance(ev, Event)]
  340. guts.dump_all(object=events, filename=filename, stream=None)
  341. else:
  342. from pyrocko.io.io_common import FileSaveError
  343. raise FileSaveError('unknown event file format: %s' % format)
  344. def load_kps_event_list(filename):
  345. elist = []
  346. f = open(filename, 'r')
  347. for line in f:
  348. toks = line.split()
  349. if len(toks) < 7:
  350. continue
  351. tim = util.ctimegm(toks[0]+' '+toks[1])
  352. lat, lon, depth, magnitude = [float(x) for x in toks[2:6]]
  353. duration = float(toks[10])
  354. region = toks[-1]
  355. name = util.gmctime_fn(tim)
  356. e = Event(
  357. lat, lon, tim,
  358. name=name,
  359. depth=depth,
  360. magnitude=magnitude,
  361. duration=duration,
  362. region=region)
  363. elist.append(e)
  364. f.close()
  365. return elist
  366. def load_gfz_event_list(filename):
  367. from pyrocko import catalog
  368. cat = catalog.Geofon()
  369. elist = []
  370. f = open(filename, 'r')
  371. for line in f:
  372. e = cat.get_event(line.strip())
  373. elist.append(e)
  374. f.close()
  375. return elist