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'''
10import logging
11import re
12from xml.parsers.expat import ParserCreate
14from pyrocko import util, pz, model
15from pyrocko.util import Request, HTTPError, urlopen, urlencode
18logger = logging.getLogger('pyrocko.client.iris')
20base_url = 'http://service.iris.edu/irisws'
23def tdate(s):
24 return util.str_to_time(s, '%Y-%m-%d')
27def sdate(t):
28 return util.time_to_str(t, '%Y-%m-%d')
31def tdatetime(s):
32 return util.str_to_time(s, format='%Y-%m-%dT%H:%M:%S')
35def sdatetime(t):
36 return util.time_to_str(t, format='%Y-%m-%dT%H:%M:%S')
39class Element(object):
41 def __init__(self, name, depth, attrs):
42 self.name = name
43 self.depth = depth
44 self.attrs = attrs
46 def __getattr__(self, k):
47 return self.attrs[k]
49 def __str__(self):
50 return '%s:\n ' % self.name + '\n '.join(
51 '%s : %s' % (k, v) for (k, v) in self.attrs.items())
54class XMLZipHandler(object):
56 def __init__(self, watch):
57 self.watch = watch
58 self.stack = []
59 self.wstack = []
60 self.outstack = []
62 def startElement(self, name, attrs):
63 self.stack.append((name, []))
64 if len(self.wstack) < len(self.watch) and \
65 name == self.watch[len(self.wstack)]:
67 el = Element(name, len(self.stack), dict(attrs.items()))
68 self.wstack.append(el)
70 def characters(self, content):
71 if self.wstack and len(self.stack) == self.wstack[-1].depth + 1:
72 if content.strip():
73 self.stack[-1][1].append(content)
75 def endElement(self, name):
76 if self.wstack:
77 if len(self.stack) == self.wstack[-1].depth + 1 and \
78 self.stack[-1][1]:
80 self.wstack[-1].attrs[name] = ''.join(self.stack[-1][1])
82 if name == self.watch[len(self.wstack)-1]:
83 if len(self.wstack) == len(self.watch):
84 self.outstack.append(list(self.wstack))
86 self.wstack.pop()
88 self.stack.pop()
90 def getQueuedElements(self):
91 outstack = self.outstack
92 self.outstack = []
93 return outstack
96def xmlzip(source, watch, bufsize=10000):
97 parser = ParserCreate()
98 handler = XMLZipHandler(watch)
100 parser.StartElementHandler = handler.startElement
101 parser.EndElementHandler = handler.endElement
102 parser.CharacterDataHandler = handler.characters
104 while True:
105 data = source.read(bufsize)
106 if not data:
107 break
109 parser.Parse(data, False)
110 for elements in handler.getQueuedElements():
111 yield elements
113 parser.Parse('', True)
114 for elements in handler.getQueuedElements():
115 yield elements
118class NotFound(Exception):
119 '''
120 Raised when the server sends an 404 response.
121 '''
122 def __init__(self, url):
123 Exception.__init__(self)
124 self._url = url
126 def __str__(self):
127 return 'No results for request %s' % self._url
130def ws_request(url, post=False, **kwargs):
131 url_values = urlencode(kwargs)
132 url = url + '?' + url_values
133 logger.debug('Accessing URL %s' % url)
135 req = Request(url)
136 if post:
137 req.add_data(post)
139 req.add_header('Accept', '*/*')
141 try:
142 return urlopen(req)
144 except HTTPError as e:
145 if e.code == 404:
146 raise NotFound(url)
147 else:
148 raise e
151def ws_virtualnetwork(**kwargs):
152 '''
153 Query IRIS virtualnetwork web service.
154 '''
156 for k in 'starttime', 'endtime':
157 if k in kwargs:
158 kwargs[k] = sdate(kwargs[k])
160 if 'timewindow' in kwargs:
161 tmin, tmax = kwargs.pop('timewindow')
162 kwargs['starttime'] = sdate(tmin)
163 kwargs['endtime'] = sdate(tmax)
165 return ws_request(base_url + '/virtualnetwork/1/query', **kwargs)
168def ws_sacpz(
169 network=None,
170 station=None,
171 location=None,
172 channel=None,
173 time=None,
174 tmin=None,
175 tmax=None):
176 '''
177 Query IRIS sacpz web service.
178 '''
180 d = {}
181 if network:
182 d['network'] = network
183 if station:
184 d['station'] = station
185 if location:
186 d['location'] = location
187 else:
188 d['location'] = '--'
189 if channel:
190 d['channel'] = channel
192 if tmin is not None and tmax is not None:
193 d['starttime'] = sdatetime(tmin)
194 d['endtime'] = sdatetime(tmax)
195 elif time is not None:
196 d['time'] = sdatetime(time)
198 return ws_request(base_url + '/sacpz/1/query', **d)
201def ws_resp(
202 network=None,
203 station=None,
204 location=None,
205 channel=None,
206 time=None,
207 tmin=None,
208 tmax=None):
209 '''
210 Query IRIS resp web service.
211 '''
213 d = {}
214 if network:
215 d['network'] = network
216 if station:
217 d['station'] = station
218 if location:
219 d['location'] = location
220 else:
221 d['location'] = '--'
222 if channel:
223 d['channel'] = channel
225 if tmin is not None and tmax is not None:
226 d['starttime'] = sdatetime(tmin)
227 d['endtime'] = sdatetime(tmax)
228 elif time is not None:
229 d['time'] = sdatetime(time)
231 return ws_request(base_url + '/resp/1/query', **d)
234class ChannelInfo(object):
235 def __init__(
236 self, network, station, location, channel, start, end, azimuth,
237 dip, elevation, depth, latitude, longitude, sample, input, output,
238 zpk):
240 self.network = network
241 self.station = station
242 self.location = location
243 self.channel = channel
244 self.start = start
245 self.end = end
246 self.azimuth = azimuth
247 self.dip = dip
248 self.elevation = elevation
249 self.depth = depth
250 self.latitude = latitude
251 self.longitude = longitude
252 self.sample = sample
253 self.input = input
254 self.output = output
255 self.zpk = zpk
257 def __str__(self):
258 return '%s.%s.%s.%s' % (
259 self.network, self.station, self.location, self.channel)
262def nslc(x):
263 return x.network, x.station, x.location, x.channel
266def grok_sacpz(data):
267 pzlines = []
268 d = {}
269 responses = []
270 float_keys = ('latitude', 'longitude', 'elevation', 'depth', 'dip',
271 'azimuth', 'sample')
272 string_keys = ('input', 'output', 'network', 'station', 'location',
273 'channel')
274 time_keys = ('start', 'end')
275 for line in data.splitlines():
276 line = line.strip()
277 if line.startswith('*'):
278 if pzlines:
279 if any(pzlines):
280 d['zpk'] = pz.read_sac_zpk(string='\n'.join(pzlines))
281 responses.append(d)
282 d = {}
283 pzlines = []
285 m = re.match(r'^\* ([A-Z]+)[^:]*:(.*)$', line)
286 if m:
287 k, v = m.group(1).lower(), m.group(2).strip()
288 if k in d:
289 assert False, 'duplicate entry? %s' % k
291 if k in float_keys:
292 d[k] = float(v)
293 elif k in string_keys:
294 d[k] = v
295 elif k in time_keys:
296 d[k] = tdatetime(v)
298 else:
299 pzlines.append(line)
301 if pzlines and any(pzlines):
302 d['zpk'] = pz.read_sac_zpk(string='\n'.join(pzlines))
303 responses.append(d)
305 cis = {}
306 for kwargs in responses:
307 try:
308 for k in float_keys + string_keys + time_keys:
309 if k not in kwargs:
310 logger.error('Missing entry: %s' % k)
311 raise Exception()
313 ci = ChannelInfo(**kwargs)
315 cis[nslc(ci)] = ci
317 except Exception:
318 logger.error('Error while parsing SACPZ data')
320 return cis
323def grok_station_xml(data, tmin, tmax):
325 stations = {}
327 for (sta, sta_epo, cha, cha_epo) in xmlzip(data, (
328 'Station', 'StationEpoch', 'Channel', 'Epoch')):
330 sta_beg, sta_end, cha_beg, cha_end = [tdatetime(x) for x in (
331 sta_epo.StartDate, sta_epo.EndDate, cha_epo.StartDate,
332 cha_epo.EndDate)]
334 if not (sta_beg <= tmin and tmax <= sta_end and
335 cha_beg <= tmin and tmax <= cha_end):
337 continue
339 nslc = tuple([str(x.strip()) for x in (
340 sta.net_code, sta.sta_code, cha.loc_code, cha.chan_code)])
342 lat, lon, ele, dep, azi, dip = [
343 float(cha_epo.attrs[x])
344 for x in 'Lat Lon Elevation Depth Azimuth Dip'.split()]
346 nsl = nslc[:3]
347 if nsl not in stations:
348 stations[nsl] = model.Station(
349 nsl[0], nsl[1], nsl[2], lat, lon, ele, dep)
351 stations[nsl].add_channel(model.station.Channel(nslc[-1], azi, dip))
353 return list(stations.values())
356def grok_virtualnet_xml(data):
357 net_sta = set()
359 for network, station in xmlzip(data, ('network', 'station')):
360 net_sta.add((network.code, station.code))
362 return net_sta
365def data_selection(stations, tmin, tmax, channel_prio=[
366 ['BHZ', 'HHZ'],
367 ['BH1', 'BHN', 'HH1', 'HHN'],
368 ['BH2', 'BHE', 'HH2', 'HHE']]):
370 selection = []
371 for station in stations:
372 wanted = []
373 for group in channel_prio:
374 gchannels = []
375 for channel in station.get_channels():
376 if channel.name in group:
377 gchannels.append(channel)
378 if gchannels:
379 gchannels.sort(key=lambda a: group.index(a.name))
380 wanted.append(gchannels[0])
382 if wanted:
383 for channel in wanted:
384 selection.append((
385 station.network, station.station, station.location,
386 channel.name, tmin, tmax))
388 return selection