1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import os
7import sys
8import subprocess
9import tempfile
10import calendar
11import time
12import re
13import logging
14import shutil
16from pyrocko import ExternalProgramMissing
17from pyrocko import pile, model, util, response
18from pyrocko.io import eventdata
20pjoin = os.path.join
21logger = logging.getLogger('pyrocko.io.rdseed')
24def cmp(a, b):
25 return (a > b) - (a < b)
28def read_station_header_file(fn):
30 def m(i, *args):
31 if len(ltoks) >= i + len(args) \
32 and (tuple(ltoks[i:i+len(args)]) == args):
33 return True
34 return False
36 def f(i):
37 return float(toks[i])
39 def s(i):
40 if len(toks) > i:
41 return toks[i]
42 else:
43 return ''
45 with open(fn, 'r') as fi:
47 stations = []
48 atsec, station, channel = None, None, None
50 for line in fi:
51 toks = line.split()
52 ltoks = line.lower().split()
53 if m(2, 'station', 'header'):
54 atsec = 'station'
55 station = {'channels': []}
56 stations.append(station)
57 continue
59 if m(2, 'station') and m(5, 'channel'):
60 atsec = 'channel'
61 channel = {}
62 station['channels'].append(channel)
63 continue
65 if atsec == 'station':
66 if m(1, 'station', 'code:'):
67 station['station'] = s(3)
69 elif m(1, 'network', 'code:'):
70 station['network'] = s(3)
72 elif m(1, 'name:'):
73 station['name'] = ' '.join(toks[2:])
75 if atsec == 'channel':
76 if m(1, 'channel:'):
77 channel['channel'] = s(2)
79 elif m(1, 'location:'):
80 channel['location'] = s(2)
82 elif m(1, 'latitude:'):
83 station['lat'] = f(2)
85 elif m(1, 'longitude:'):
86 station['lon'] = f(2)
88 elif m(1, 'elevation:'):
89 station['elevation'] = f(2)
91 elif m(1, 'local', 'depth:'):
92 channel['depth'] = f(3)
94 elif m(1, 'azimuth:'):
95 channel['azimuth'] = f(2)
97 elif m(1, 'dip:'):
98 channel['dip'] = f(2)
100 nsl_stations = {}
101 for station in stations:
102 for channel in station['channels']:
103 def cs(k, default=None):
104 return channel.get(k, station.get(k, default))
106 nsl = station['network'], station['station'], channel['location']
107 if nsl not in nsl_stations:
108 nsl_stations[nsl] = model.Station(
109 network=station['network'],
110 station=station['station'],
111 location=channel['location'],
112 lat=cs('lat'),
113 lon=cs('lon'),
114 elevation=cs('elevation'),
115 depth=cs('depth', None),
116 name=station['name'])
118 nsl_stations[nsl].add_channel(model.Channel(
119 channel['channel'],
120 azimuth=channel['azimuth'],
121 dip=channel['dip']))
123 return list(nsl_stations.values())
126def cmp_version(a, b):
127 ai = [int(x) for x in a.split('.')]
128 bi = [int(x) for x in b.split('.')]
129 return cmp(ai, bi)
132def dumb_parser(data):
134 (in_ws, in_kw, in_str) = (1, 2, 3)
136 state = in_ws
138 rows = []
139 cols = []
140 accu = ''
141 for c in data:
142 if state == in_ws:
143 if c == '"':
144 new_state = in_str
146 elif c not in (' ', '\t', '\n', '\r'):
147 new_state = in_kw
149 if state == in_kw:
150 if c in (' ', '\t', '\n', '\r'):
151 cols.append(accu)
152 accu = ''
153 if c in ('\n', '\r'):
154 rows.append(cols)
155 cols = []
156 new_state = in_ws
158 if state == in_str:
159 if c == '"':
160 accu += c
161 cols.append(accu[1:-1])
162 accu = ''
163 if c in ('\n', '\r'):
164 rows.append(cols)
165 cols = []
166 new_state = in_ws
168 state = new_state
170 if state in (in_kw, in_str):
171 accu += c
173 if len(cols) != 0:
174 rows.append(cols)
176 return rows
179class Programs(object):
180 rdseed = 'rdseed'
181 avail = None
183 @staticmethod
184 def check():
185 if Programs.avail is not None:
186 return Programs.avail
188 else:
189 try:
190 rdseed_proc = subprocess.Popen(
191 [Programs.rdseed],
192 stdin=subprocess.PIPE,
193 stdout=subprocess.PIPE,
194 stderr=subprocess.PIPE)
196 (out, err) = rdseed_proc.communicate()
198 except OSError as e:
199 if e.errno == 2:
200 reason = "Could not find executable: '%s'." \
201 % Programs.rdseed
202 else:
203 reason = str(e)
205 logging.debug('Failed to run rdseed program. %s' % reason)
206 Programs.avail = False
207 return Programs.avail
209 ms = [re.search(
210 r'Release (\d+(\.\d+(\.\d+)?)?)', s.decode())
211 for s in (err, out)]
212 ms = list(filter(bool, ms))
213 if not ms:
214 logger.error('Cannot determine rdseed version number.')
215 else:
216 version = ms[0].group(1)
217 if cmp_version('4.7.5', version) == 1 \
218 or cmp_version(version, '5.1') == 1:
220 logger.warning(
221 'Module pyrocko.rdseed has not been tested with '
222 'version %s of rdseed.' % version)
224 Programs.avail = True
225 return Programs.avail
228class SeedVolumeNotFound(Exception):
229 pass
232class SeedVolumeAccess(eventdata.EventDataAccess):
234 def __init__(self, seedvolume, datapile=None):
236 '''
237 Create new SEED Volume access object.
239 :param seedvolume: filename of seed volume
240 :param datapile: if not ``None``, this should be a
241 :py:class:`pile.Pile` object with data traces which are then used
242 instead of the data provided by the SEED volume.
243 (This is useful for dataless SEED volumes.)
244 '''
246 eventdata.EventDataAccess.__init__(self, datapile=datapile)
247 self.tempdir = None
248 Programs.check()
250 self.tempdir = None
251 self.seedvolume = seedvolume
252 if not os.path.isfile(self.seedvolume):
253 raise SeedVolumeNotFound()
255 self.tempdir = tempfile.mkdtemp('', 'SeedVolumeAccess-')
256 self.station_headers_file = os.path.join(
257 self.tempdir, 'station_header_infos')
258 self._unpack()
259 self.shutil = shutil
261 def __del__(self):
262 if self.tempdir:
263 self.shutil.rmtree(self.tempdir)
265 def get_pile(self):
266 if self._pile is None:
267 fns = util.select_files([self.tempdir], include=r'\.SAC$')
268 self._pile = pile.Pile()
269 self._pile.load_files(fns, fileformat='sac')
271 return self._pile
273 def get_resp_file(self, tr):
274 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % tr.nslc_id)
275 if not os.path.exists(respfile):
276 raise eventdata.NoRestitution(
277 'no response information available for trace %s.%s.%s.%s'
278 % tr.nslc_id)
280 return respfile
282 def get_stationxml(self):
283 stations = self.get_pyrocko_stations()
284 respfiles = []
285 for station in stations:
286 for channel in station.get_channels():
287 nslc = station.nsl() + (channel.name,)
288 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % nslc)
289 respfiles.append(respfile)
291 from pyrocko.fdsn import resp
292 sxml = resp.make_stationxml(stations, resp.iload(respfiles))
293 return sxml
295 def get_restitution(self, tr, allowed_methods):
296 if 'evalresp' in allowed_methods:
297 respfile = pjoin(self.tempdir, 'RESP.%s.%s.%s.%s' % tr.nslc_id)
298 if not os.path.exists(respfile):
299 raise eventdata.NoRestitution(
300 'no response information available for trace %s.%s.%s.%s'
301 % tr.nslc_id)
303 trans = response.InverseEvalresp(respfile, tr)
304 return trans
305 else:
306 raise eventdata.NoRestitution(
307 'no response information available for trace %s.%s.%s.%s'
308 % tr.nslc_id)
310 def get_pyrocko_response(self, tr, target):
311 '''
312 Extract the frequency response as :py:class:`response.EvalResp`
313 instance for *tr*.
315 :param tr: :py:class:`trace.Trace` instance
316 '''
317 respfile = self.get_resp_file(tr)
318 return response.Evalresp(respfile, tr, target)
320 def _unpack(self):
321 input_fn = self.seedvolume
322 output_dir = self.tempdir
324 def strerr(s):
325 return '\n'.join(['rdseed: '+line.decode()
326 for line in s.splitlines()])
327 try:
329 # seismograms:
330 if self._pile is None:
331 rdseed_proc = subprocess.Popen(
332 [Programs.rdseed, '-f', input_fn, '-d', '-z', '3', '-o',
333 '1', '-p', '-R', '-q', output_dir],
334 stdin=subprocess.PIPE,
335 stdout=subprocess.PIPE,
336 stderr=subprocess.PIPE)
338 (out, err) = rdseed_proc.communicate()
339 logging.info(strerr(err))
340 else:
341 rdseed_proc = subprocess.Popen(
342 [Programs.rdseed, '-f', input_fn, '-z', '3', '-p', '-R',
343 '-q', output_dir],
344 stdin=subprocess.PIPE,
345 stdout=subprocess.PIPE,
346 stderr=subprocess.PIPE)
348 (out, err) = rdseed_proc.communicate()
349 logging.info(strerr(err))
351 # event data:
352 rdseed_proc = subprocess.Popen(
353 [Programs.rdseed, '-f', input_fn, '-e', '-q', output_dir],
354 stdin=subprocess.PIPE,
355 stdout=subprocess.PIPE,
356 stderr=subprocess.PIPE)
358 (out, err) = rdseed_proc.communicate()
359 logging.info(strerr(err))
361 # station summary information:
362 rdseed_proc = subprocess.Popen(
363 [Programs.rdseed, '-f', input_fn, '-S', '-q', output_dir],
364 stdin=subprocess.PIPE,
365 stdout=subprocess.PIPE,
366 stderr=subprocess.PIPE)
368 (out, err) = rdseed_proc.communicate()
369 logging.info(strerr(err))
371 # station headers:
372 rdseed_proc = subprocess.Popen(
373 [Programs.rdseed, '-f', input_fn, '-s', '-q', output_dir],
374 stdin=subprocess.PIPE,
375 stdout=subprocess.PIPE,
376 stderr=subprocess.PIPE)
378 (out, err) = rdseed_proc.communicate()
379 with open(self.station_headers_file, 'w') as fout:
380 fout.write(out.decode())
381 logging.info(strerr(err))
383 except OSError as e:
384 if e.errno == 2:
385 reason = "Could not find executable: '%s'." % Programs.rdseed
386 else:
387 reason = str(e)
389 logging.fatal('Failed to unpack SEED volume. %s' % reason)
390 sys.exit(1)
392 def _get_events_from_file(self):
393 rdseed_event_file = os.path.join(self.tempdir, 'rdseed.events')
394 if not os.path.isfile(rdseed_event_file):
395 return []
397 with open(rdseed_event_file, 'r') as f:
398 events = []
399 for line in f:
400 toks = line.split(', ')
401 if len(toks) == 9:
402 datetime = toks[1].split('.')[0]
403 format = '%Y/%m/%d %H:%M:%S'
404 secs = util.to_time_float(calendar.timegm(
405 time.strptime(datetime, format)))
406 e = model.Event(
407 lat=float(toks[2]),
408 lon=float(toks[3]),
409 depth=float(toks[4])*1000.,
410 magnitude=float(toks[8]),
411 time=secs)
413 events.append(e)
414 else:
415 raise Exception('Event description in unrecognized format')
417 return events
419 def _get_stations_from_file(self):
420 stations = read_station_header_file(self.station_headers_file)
421 return stations
423 def _insert_channel_descriptions(self, stations):
424 # this is done beforehand in this class
425 pass
428def check_have_rdseed():
429 if not Programs.check():
430 raise ExternalProgramMissing(
431 'rdseed is not installed or cannot be found')
434if __name__ == '__main__':
435 print(SeedVolumeAccess(sys.argv[1]).get_stationxml().dump_xml())