1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import
7import os
8import sys
9import subprocess
10import tempfile
11import calendar
12import time
13import re
14import logging
15import shutil
17from pyrocko import ExternalProgramMissing
18from pyrocko import trace, pile, model, util
19from pyrocko.io import eventdata
21pjoin = os.path.join
22logger = logging.getLogger('pyrocko.io.rdseed')
25def cmp(a, b):
26 return (a > b) - (a < b)
29def read_station_header_file(fn):
31 def m(i, *args):
32 if len(ltoks) >= i + len(args) \
33 and (tuple(ltoks[i:i+len(args)]) == args):
34 return True
35 return False
37 def f(i):
38 return float(toks[i])
40 def s(i):
41 if len(toks) > i:
42 return toks[i]
43 else:
44 return ''
46 with open(fn, 'r') as fi:
48 stations = []
49 atsec, station, channel = None, None, None
51 for line in fi:
52 toks = line.split()
53 ltoks = line.lower().split()
54 if m(2, 'station', 'header'):
55 atsec = 'station'
56 station = {'channels': []}
57 stations.append(station)
58 continue
60 if m(2, 'station') and m(5, 'channel'):
61 atsec = 'channel'
62 channel = {}
63 station['channels'].append(channel)
64 continue
66 if atsec == 'station':
67 if m(1, 'station', 'code:'):
68 station['station'] = s(3)
70 elif m(1, 'network', 'code:'):
71 station['network'] = s(3)
73 elif m(1, 'name:'):
74 station['name'] = ' '.join(toks[2:])
76 if atsec == 'channel':
77 if m(1, 'channel:'):
78 channel['channel'] = s(2)
80 elif m(1, 'location:'):
81 channel['location'] = s(2)
83 elif m(1, 'latitude:'):
84 station['lat'] = f(2)
86 elif m(1, 'longitude:'):
87 station['lon'] = f(2)
89 elif m(1, 'elevation:'):
90 station['elevation'] = f(2)
92 elif m(1, 'local', 'depth:'):
93 channel['depth'] = f(3)
95 elif m(1, 'azimuth:'):
96 channel['azimuth'] = f(2)
98 elif m(1, 'dip:'):
99 channel['dip'] = f(2)
101 nsl_stations = {}
102 for station in stations:
103 for channel in station['channels']:
104 def cs(k, default=None):
105 return channel.get(k, station.get(k, default))
107 nsl = station['network'], station['station'], channel['location']
108 if nsl not in nsl_stations:
109 nsl_stations[nsl] = model.Station(
110 network=station['network'],
111 station=station['station'],
112 location=channel['location'],
113 lat=cs('lat'),
114 lon=cs('lon'),
115 elevation=cs('elevation'),
116 depth=cs('depth', None),
117 name=station['name'])
119 nsl_stations[nsl].add_channel(model.Channel(
120 channel['channel'],
121 azimuth=channel['azimuth'],
122 dip=channel['dip']))
124 return list(nsl_stations.values())
127def cmp_version(a, b):
128 ai = [int(x) for x in a.split('.')]
129 bi = [int(x) for x in b.split('.')]
130 return cmp(ai, bi)
133def dumb_parser(data):
135 (in_ws, in_kw, in_str) = (1, 2, 3)
137 state = in_ws
139 rows = []
140 cols = []
141 accu = ''
142 for c in data:
143 if state == in_ws:
144 if c == '"':
145 new_state = in_str
147 elif c not in (' ', '\t', '\n', '\r'):
148 new_state = in_kw
150 if state == in_kw:
151 if c in (' ', '\t', '\n', '\r'):
152 cols.append(accu)
153 accu = ''
154 if c in ('\n', '\r'):
155 rows.append(cols)
156 cols = []
157 new_state = in_ws
159 if state == in_str:
160 if c == '"':
161 accu += c
162 cols.append(accu[1:-1])
163 accu = ''
164 if c in ('\n', '\r'):
165 rows.append(cols)
166 cols = []
167 new_state = in_ws
169 state = new_state
171 if state in (in_kw, in_str):
172 accu += c
174 if len(cols) != 0:
175 rows.append(cols)
177 return rows
180class Programs(object):
181 rdseed = 'rdseed'
182 avail = None
184 @staticmethod
185 def check():
186 if Programs.avail is not None:
187 return Programs.avail
189 else:
190 try:
191 rdseed_proc = subprocess.Popen(
192 [Programs.rdseed],
193 stdin=subprocess.PIPE,
194 stdout=subprocess.PIPE,
195 stderr=subprocess.PIPE)
197 (out, err) = rdseed_proc.communicate()
199 except OSError as e:
200 if e.errno == 2:
201 reason = "Could not find executable: '%s'." \
202 % Programs.rdseed
203 else:
204 reason = str(e)
206 logging.debug('Failed to run rdseed program. %s' % reason)
207 Programs.avail = False
208 return Programs.avail
210 ms = [re.search(
211 r'Release (\d+(\.\d+(\.\d+)?)?)', s.decode())
212 for s in (err, out)]
213 ms = list(filter(bool, ms))
214 if not ms:
215 logger.error('Cannot determine rdseed version number.')
216 else:
217 version = ms[0].group(1)
218 if cmp_version('4.7.5', version) == 1 \
219 or cmp_version(version, '5.1') == 1:
221 logger.warning(
222 'Module pyrocko.rdseed has not been tested with '
223 'version %s of rdseed.' % version)
225 Programs.avail = True
226 return Programs.avail
229class SeedVolumeNotFound(Exception):
230 pass
233class SeedVolumeAccess(eventdata.EventDataAccess):
235 def __init__(self, seedvolume, datapile=None):
237 '''
238 Create new SEED Volume access object.
240 :param seedvolume: filename of seed volume
241 :param datapile: if not ``None``, this should be a
242 :py:class:`pile.Pile` object with data traces which are then used
243 instead of the data provided by the SEED volume.
244 (This is useful for dataless SEED volumes.)
245 '''
247 eventdata.EventDataAccess.__init__(self, datapile=datapile)
248 self.tempdir = None
249 Programs.check()
251 self.tempdir = None
252 self.seedvolume = seedvolume
253 if not os.path.isfile(self.seedvolume):
254 raise SeedVolumeNotFound()
256 self.tempdir = tempfile.mkdtemp("", "SeedVolumeAccess-")
257 self.station_headers_file = os.path.join(
258 self.tempdir, 'station_header_infos')
259 self._unpack()
260 self.shutil = shutil
262 def __del__(self):
263 if self.tempdir:
264 self.shutil.rmtree(self.tempdir)
266 def get_pile(self):
267 if self._pile is None:
268 fns = util.select_files([self.tempdir], regex=r'\.SAC$')
269 self._pile = pile.Pile()
270 self._pile.load_files(fns, fileformat='sac')
272 return self._pile
274 def get_resp_file(self, tr):
275 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % tr.nslc_id)
276 if not os.path.exists(respfile):
277 raise eventdata.NoRestitution(
278 'no response information available for trace %s.%s.%s.%s'
279 % tr.nslc_id)
281 return respfile
283 def get_stationxml(self):
284 stations = self.get_pyrocko_stations()
285 respfiles = []
286 for station in stations:
287 for channel in station.get_channels():
288 nslc = station.nsl() + (channel.name,)
289 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % nslc)
290 respfiles.append(respfile)
292 from pyrocko.fdsn import resp
293 sxml = resp.make_stationxml(stations, resp.iload(respfiles))
294 return sxml
296 def get_restitution(self, tr, allowed_methods):
297 if 'evalresp' in allowed_methods:
298 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % tr.nslc_id)
299 if not os.path.exists(respfile):
300 raise eventdata.NoRestitution(
301 'no response information available for trace %s.%s.%s.%s'
302 % tr.nslc_id)
304 trans = trace.InverseEvalresp(respfile, tr)
305 return trans
306 else:
307 raise eventdata.NoRestitution(
308 'no response information available for trace %s.%s.%s.%s'
309 % tr.nslc_id)
311 def get_pyrocko_response(self, tr, target):
312 '''
313 Extract the frequency response as :py:class:`trace.EvalResp`
314 instance for *tr*.
316 :param tr: :py:class:`trace.Trace` instance
317 '''
318 respfile = self.get_resp_file(tr)
319 return trace.Evalresp(respfile, tr, target)
321 def _unpack(self):
322 input_fn = self.seedvolume
323 output_dir = self.tempdir
325 def strerr(s):
326 return '\n'.join(['rdseed: '+line.decode()
327 for line in s.splitlines()])
328 try:
330 # seismograms:
331 if self._pile is None:
332 rdseed_proc = subprocess.Popen(
333 [Programs.rdseed, '-f', input_fn, '-d', '-z', '3', '-o',
334 '1', '-p', '-R', '-q', output_dir],
335 stdin=subprocess.PIPE,
336 stdout=subprocess.PIPE,
337 stderr=subprocess.PIPE)
339 (out, err) = rdseed_proc.communicate()
340 logging.info(strerr(err))
341 else:
342 rdseed_proc = subprocess.Popen(
343 [Programs.rdseed, '-f', input_fn, '-z', '3', '-p', '-R',
344 '-q', output_dir],
345 stdin=subprocess.PIPE,
346 stdout=subprocess.PIPE,
347 stderr=subprocess.PIPE)
349 (out, err) = rdseed_proc.communicate()
350 logging.info(strerr(err))
352 # event data:
353 rdseed_proc = subprocess.Popen(
354 [Programs.rdseed, '-f', input_fn, '-e', '-q', output_dir],
355 stdin=subprocess.PIPE,
356 stdout=subprocess.PIPE,
357 stderr=subprocess.PIPE)
359 (out, err) = rdseed_proc.communicate()
360 logging.info(strerr(err))
362 # station summary information:
363 rdseed_proc = subprocess.Popen(
364 [Programs.rdseed, '-f', input_fn, '-S', '-q', output_dir],
365 stdin=subprocess.PIPE,
366 stdout=subprocess.PIPE,
367 stderr=subprocess.PIPE)
369 (out, err) = rdseed_proc.communicate()
370 logging.info(strerr(err))
372 # station headers:
373 rdseed_proc = subprocess.Popen(
374 [Programs.rdseed, '-f', input_fn, '-s', '-q', output_dir],
375 stdin=subprocess.PIPE,
376 stdout=subprocess.PIPE,
377 stderr=subprocess.PIPE)
379 (out, err) = rdseed_proc.communicate()
380 with open(self.station_headers_file, 'w') as fout:
381 fout.write(out.decode())
382 logging.info(strerr(err))
384 except OSError as e:
385 if e.errno == 2:
386 reason = "Could not find executable: '%s'." % Programs.rdseed
387 else:
388 reason = str(e)
390 logging.fatal('Failed to unpack SEED volume. %s' % reason)
391 sys.exit(1)
393 def _get_events_from_file(self):
394 rdseed_event_file = os.path.join(self.tempdir, 'rdseed.events')
395 if not os.path.isfile(rdseed_event_file):
396 return []
398 with open(rdseed_event_file, 'r') as f:
399 events = []
400 for line in f:
401 toks = line.split(', ')
402 if len(toks) == 9:
403 datetime = toks[1].split('.')[0]
404 format = '%Y/%m/%d %H:%M:%S'
405 secs = util.to_time_float(calendar.timegm(
406 time.strptime(datetime, format)))
407 e = model.Event(
408 lat=float(toks[2]),
409 lon=float(toks[3]),
410 depth=float(toks[4])*1000.,
411 magnitude=float(toks[8]),
412 time=secs)
414 events.append(e)
415 else:
416 raise Exception('Event description in unrecognized format')
418 return events
420 def _get_stations_from_file(self):
421 stations = read_station_header_file(self.station_headers_file)
422 return stations
424 def _insert_channel_descriptions(self, stations):
425 # this is done beforehand in this class
426 pass
429def check_have_rdseed():
430 if not Programs.check():
431 raise ExternalProgramMissing(
432 'rdseed is not installed or cannot be found')