1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6'''
7Access to some IRIS (non-FDSN) web services.
8'''
10from __future__ import absolute_import
13import logging
14import re
15from xml.parsers.expat import ParserCreate
17from pyrocko import util, pz, model
18from pyrocko.util import Request, HTTPError, urlopen, urlencode
21logger = logging.getLogger('pyrocko.client.iris')
23base_url = 'http://service.iris.edu/irisws'
26def tdate(s):
27 return util.str_to_time(s, '%Y-%m-%d')
30def sdate(t):
31 return util.time_to_str(t, '%Y-%m-%d')
34def tdatetime(s):
35 return util.str_to_time(s, format='%Y-%m-%dT%H:%M:%S')
38def sdatetime(t):
39 return util.time_to_str(t, format='%Y-%m-%dT%H:%M:%S')
42class Element(object):
44 def __init__(self, name, depth, attrs):
45 self.name = name
46 self.depth = depth
47 self.attrs = attrs
49 def __getattr__(self, k):
50 return self.attrs[k]
52 def __str__(self):
53 return '%s:\n ' % self.name + '\n '.join(
54 '%s : %s' % (k, v) for (k, v) in self.attrs.items())
57class XMLZipHandler(object):
59 def __init__(self, watch):
60 self.watch = watch
61 self.stack = []
62 self.wstack = []
63 self.outstack = []
65 def startElement(self, name, attrs):
66 self.stack.append((name, []))
67 if len(self.wstack) < len(self.watch) and \
68 name == self.watch[len(self.wstack)]:
70 el = Element(name, len(self.stack), dict(attrs.items()))
71 self.wstack.append(el)
73 def characters(self, content):
74 if self.wstack and len(self.stack) == self.wstack[-1].depth + 1:
75 if content.strip():
76 self.stack[-1][1].append(content)
78 def endElement(self, name):
79 if self.wstack:
80 if len(self.stack) == self.wstack[-1].depth + 1 and \
81 self.stack[-1][1]:
83 self.wstack[-1].attrs[name] = ''.join(self.stack[-1][1])
85 if name == self.watch[len(self.wstack)-1]:
86 if len(self.wstack) == len(self.watch):
87 self.outstack.append(list(self.wstack))
89 self.wstack.pop()
91 self.stack.pop()
93 def getQueuedElements(self):
94 outstack = self.outstack
95 self.outstack = []
96 return outstack
99def xmlzip(source, watch, bufsize=10000):
100 parser = ParserCreate()
101 handler = XMLZipHandler(watch)
103 parser.StartElementHandler = handler.startElement
104 parser.EndElementHandler = handler.endElement
105 parser.CharacterDataHandler = handler.characters
107 while True:
108 data = source.read(bufsize)
109 if not data:
110 break
112 parser.Parse(data, False)
113 for elements in handler.getQueuedElements():
114 yield elements
116 parser.Parse('', True)
117 for elements in handler.getQueuedElements():
118 yield elements
121class NotFound(Exception):
122 '''
123 Raised when the server sends an 404 response.
124 '''
125 def __init__(self, url):
126 Exception.__init__(self)
127 self._url = url
129 def __str__(self):
130 return 'No results for request %s' % self._url
133def ws_request(url, post=False, **kwargs):
134 url_values = urlencode(kwargs)
135 url = url + '?' + url_values
136 logger.debug('Accessing URL %s' % url)
138 req = Request(url)
139 if post:
140 req.add_data(post)
142 req.add_header('Accept', '*/*')
144 try:
145 return urlopen(req)
147 except HTTPError as e:
148 if e.code == 404:
149 raise NotFound(url)
150 else:
151 raise e
154def ws_virtualnetwork(**kwargs):
155 '''
156 Query IRIS virtualnetwork web service.
157 '''
159 for k in 'starttime', 'endtime':
160 if k in kwargs:
161 kwargs[k] = sdate(kwargs[k])
163 if 'timewindow' in kwargs:
164 tmin, tmax = kwargs.pop('timewindow')
165 kwargs['starttime'] = sdate(tmin)
166 kwargs['endtime'] = sdate(tmax)
168 return ws_request(base_url + '/virtualnetwork/1/query', **kwargs)
171def ws_sacpz(
172 network=None,
173 station=None,
174 location=None,
175 channel=None,
176 time=None,
177 tmin=None,
178 tmax=None):
179 '''
180 Query IRIS sacpz web service.
181 '''
183 d = {}
184 if network:
185 d['network'] = network
186 if station:
187 d['station'] = station
188 if location:
189 d['location'] = location
190 else:
191 d['location'] = '--'
192 if channel:
193 d['channel'] = channel
195 if tmin is not None and tmax is not None:
196 d['starttime'] = sdatetime(tmin)
197 d['endtime'] = sdatetime(tmax)
198 elif time is not None:
199 d['time'] = sdatetime(time)
201 return ws_request(base_url + '/sacpz/1/query', **d)
204def ws_resp(
205 network=None,
206 station=None,
207 location=None,
208 channel=None,
209 time=None,
210 tmin=None,
211 tmax=None):
212 '''
213 Query IRIS resp web service.
214 '''
216 d = {}
217 if network:
218 d['network'] = network
219 if station:
220 d['station'] = station
221 if location:
222 d['location'] = location
223 else:
224 d['location'] = '--'
225 if channel:
226 d['channel'] = channel
228 if tmin is not None and tmax is not None:
229 d['starttime'] = sdatetime(tmin)
230 d['endtime'] = sdatetime(tmax)
231 elif time is not None:
232 d['time'] = sdatetime(time)
234 return ws_request(base_url + '/resp/1/query', **d)
237class ChannelInfo(object):
238 def __init__(
239 self, network, station, location, channel, start, end, azimuth,
240 dip, elevation, depth, latitude, longitude, sample, input, output,
241 zpk):
243 self.network = network
244 self.station = station
245 self.location = location
246 self.channel = channel
247 self.start = start
248 self.end = end
249 self.azimuth = azimuth
250 self.dip = dip
251 self.elevation = elevation
252 self.depth = depth
253 self.latitude = latitude
254 self.longitude = longitude
255 self.sample = sample
256 self.input = input
257 self.output = output
258 self.zpk = zpk
260 def __str__(self):
261 return '%s.%s.%s.%s' % (
262 self.network, self.station, self.location, self.channel)
265def nslc(x):
266 return x.network, x.station, x.location, x.channel
269def grok_sacpz(data):
270 pzlines = []
271 d = {}
272 responses = []
273 float_keys = ('latitude', 'longitude', 'elevation', 'depth', 'dip',
274 'azimuth', 'sample')
275 string_keys = ('input', 'output', 'network', 'station', 'location',
276 'channel')
277 time_keys = ('start', 'end')
278 for line in data.splitlines():
279 line = line.strip()
280 if line.startswith('*'):
281 if pzlines:
282 if any(pzlines):
283 d['zpk'] = pz.read_sac_zpk(string='\n'.join(pzlines))
284 responses.append(d)
285 d = {}
286 pzlines = []
288 m = re.match(r'^\* ([A-Z]+)[^:]*:(.*)$', line)
289 if m:
290 k, v = m.group(1).lower(), m.group(2).strip()
291 if k in d:
292 assert False, 'duplicate entry? %s' % k
294 if k in float_keys:
295 d[k] = float(v)
296 elif k in string_keys:
297 d[k] = v
298 elif k in time_keys:
299 d[k] = tdatetime(v)
301 else:
302 pzlines.append(line)
304 if pzlines and any(pzlines):
305 d['zpk'] = pz.read_sac_zpk(string='\n'.join(pzlines))
306 responses.append(d)
308 cis = {}
309 for kwargs in responses:
310 try:
311 for k in float_keys + string_keys + time_keys:
312 if k not in kwargs:
313 logger.error('Missing entry: %s' % k)
314 raise Exception()
316 ci = ChannelInfo(**kwargs)
318 cis[nslc(ci)] = ci
320 except Exception:
321 logger.error('Error while parsing SACPZ data')
323 return cis
326def grok_station_xml(data, tmin, tmax):
328 stations = {}
330 for (sta, sta_epo, cha, cha_epo) in xmlzip(data, (
331 'Station', 'StationEpoch', 'Channel', 'Epoch')):
333 sta_beg, sta_end, cha_beg, cha_end = [tdatetime(x) for x in (
334 sta_epo.StartDate, sta_epo.EndDate, cha_epo.StartDate,
335 cha_epo.EndDate)]
337 if not (sta_beg <= tmin and tmax <= sta_end and
338 cha_beg <= tmin and tmax <= cha_end):
340 continue
342 nslc = tuple([str(x.strip()) for x in (
343 sta.net_code, sta.sta_code, cha.loc_code, cha.chan_code)])
345 lat, lon, ele, dep, azi, dip = [
346 float(cha_epo.attrs[x])
347 for x in 'Lat Lon Elevation Depth Azimuth Dip'.split()]
349 nsl = nslc[:3]
350 if nsl not in stations:
351 stations[nsl] = model.Station(
352 nsl[0], nsl[1], nsl[2], lat, lon, ele, dep)
354 stations[nsl].add_channel(model.station.Channel(nslc[-1], azi, dip))
356 return list(stations.values())
359def grok_virtualnet_xml(data):
360 net_sta = set()
362 for network, station in xmlzip(data, ('network', 'station')):
363 net_sta.add((network.code, station.code))
365 return net_sta
368def data_selection(stations, tmin, tmax, channel_prio=[
369 ['BHZ', 'HHZ'],
370 ['BH1', 'BHN', 'HH1', 'HHN'],
371 ['BH2', 'BHE', 'HH2', 'HHE']]):
373 selection = []
374 for station in stations:
375 wanted = []
376 for group in channel_prio:
377 gchannels = []
378 for channel in station.get_channels():
379 if channel.name in group:
380 gchannels.append(channel)
381 if gchannels:
382 gchannels.sort(key=lambda a: group.index(a.name))
383 wanted.append(gchannels[0])
385 if wanted:
386 for channel in wanted:
387 selection.append((
388 station.network, station.station, station.location,
389 channel.name, tmin, tmax))
391 return selection