1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import math
8import copy
9import logging
10import numpy as num
12from pyrocko import orthodrome
13from pyrocko.orthodrome import wrap
14from pyrocko.guts import Object, Float, String, List, dump_all
16from .location import Location
18logger = logging.getLogger('pyrocko.model.station')
20guts_prefix = 'pf'
22d2r = num.pi / 180.
25class ChannelsNotOrthogonal(Exception):
26 pass
29def guess_azimuth_from_name(channel_name):
30 if channel_name.endswith('N'):
31 return 0.
32 elif channel_name.endswith('E'):
33 return 90.
34 elif channel_name.endswith('Z'):
35 return 0.
37 return None
40def guess_dip_from_name(channel_name):
41 if channel_name.endswith('N'):
42 return 0.
43 elif channel_name.endswith('E'):
44 return 0.
45 elif channel_name.endswith('Z'):
46 return -90.
48 return None
51def guess_azimuth_dip_from_name(channel_name):
52 return guess_azimuth_from_name(channel_name), \
53 guess_dip_from_name(channel_name)
56def mkvec(x, y, z):
57 return num.array([x, y, z], dtype=float)
60def are_orthogonal(enus, eps=0.05):
61 return all(abs(x) < eps for x in [
62 num.dot(enus[0], enus[1]),
63 num.dot(enus[1], enus[2]),
64 num.dot(enus[2], enus[0])])
67def fill_orthogonal(enus):
69 nmiss = sum(x is None for x in enus)
71 if nmiss == 1:
72 for ic in range(len(enus)):
73 if enus[ic] is None:
74 enus[ic] = num.cross(enus[(ic-2) % 3], enus[(ic-1) % 3])
76 if nmiss == 2:
77 for ic in range(len(enus)):
78 if enus[ic] is not None:
79 xenu = enus[ic] + mkvec(1, 1, 1)
80 enus[(ic+1) % 3] = num.cross(enus[ic], xenu)
81 enus[(ic+2) % 3] = num.cross(enus[ic], enus[(ic+1) % 3])
83 if nmiss == 3:
84 # holy camoly..
85 enus[0] = mkvec(1, 0, 0)
86 enus[1] = mkvec(0, 1, 0)
87 enus[2] = mkvec(0, 0, 1)
90class Channel(Object):
91 name = String.T()
92 azimuth = Float.T(optional=True)
93 dip = Float.T(optional=True)
94 gain = Float.T(default=1.0)
96 def __init__(self, name, azimuth=None, dip=None, gain=1.0):
97 if azimuth is None:
98 azimuth = guess_azimuth_from_name(name)
99 if dip is None:
100 dip = guess_dip_from_name(name)
102 Object.__init__(
103 self,
104 name=name,
105 azimuth=float_or_none(azimuth),
106 dip=float_or_none(dip),
107 gain=float(gain))
109 @property
110 def ned(self):
111 if None in (self.azimuth, self.dip):
112 return None
114 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
115 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
116 d = math.sin(self.dip*d2r)
117 return mkvec(n, e, d)
119 @property
120 def enu(self):
121 if None in (self.azimuth, self.dip):
122 return None
124 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
125 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
126 d = math.sin(self.dip*d2r)
127 return mkvec(e, n, -d)
129 def __str__(self):
130 return '%s %f %f %g' % (self.name, self.azimuth, self.dip, self.gain)
133class Station(Location):
134 network = String.T()
135 station = String.T()
136 location = String.T()
137 name = String.T(default='')
138 channels = List.T(Channel.T())
140 def __init__(self, network='', station='', location='',
141 lat=0.0, lon=0.0,
142 elevation=0.0, depth=0.0,
143 north_shift=0.0, east_shift=0.0,
144 name='', channels=None):
146 Location.__init__(
147 self,
148 network=network, station=station, location=location,
149 lat=float(lat), lon=float(lon),
150 elevation=float(elevation),
151 depth=float(depth),
152 north_shift=float(north_shift),
153 east_shift=float(east_shift),
154 name=name or '',
155 channels=channels or [])
157 self.dist_deg = None
158 self.dist_m = None
159 self.azimuth = None
160 self.backazimuth = None
162 def copy(self):
163 return copy.deepcopy(self)
165 def set_event_relative_data(self, event, distance_3d=False):
166 surface_dist = self.distance_to(event)
167 if distance_3d:
168 if event.depth is None:
169 logger.warning('No event depth given: using 0 m.')
170 dd = 0.0 - self.depth
171 else:
172 dd = event.depth - self.depth
174 self.dist_m = math.sqrt(dd**2 + surface_dist**2)
175 else:
176 self.dist_m = surface_dist
178 self.dist_deg = surface_dist / orthodrome.earthradius_equator * \
179 orthodrome.r2d
181 self.azimuth, self.backazimuth = event.azibazi_to(self)
183 def set_channels_by_name(self, *args):
184 self.set_channels([])
185 for name in args:
186 self.add_channel(Channel(name))
188 def set_channels(self, channels):
189 self.channels = []
190 for ch in channels:
191 self.add_channel(ch)
193 def get_channels(self):
194 return list(self.channels)
196 def get_channel_names(self):
197 return set(ch.name for ch in self.channels)
199 def remove_channel_by_name(self, name):
200 todel = [ch for ch in self.channels if ch.name == name]
201 for ch in todel:
202 self.channels.remove(ch)
204 def add_channel(self, channel):
205 self.remove_channel_by_name(channel.name)
206 self.channels.append(channel)
207 self.channels.sort(key=lambda ch: ch.name)
209 def get_channel(self, name):
210 for ch in self.channels:
211 if ch.name == name:
212 return ch
214 return None
216 def rotation_ne_to_rt(self, in_channel_names, out_channel_names):
218 angle = wrap(self.backazimuth + 180., -180., 180.)
219 in_channels = [self.get_channel(name) for name in in_channel_names]
220 out_channels = [
221 Channel(out_channel_names[0],
222 wrap(self.backazimuth+180., -180., 180.), 0., 1.),
223 Channel(out_channel_names[1],
224 wrap(self.backazimuth+270., -180., 180.), 0., 1.)]
225 return angle, in_channels, out_channels
227 def _projection_to(
228 self, to, in_channel_names, out_channel_names, use_gains=False):
230 in_channels = [self.get_channel(name) for name in in_channel_names]
232 # create orthogonal vectors for missing components, such that this
233 # won't break projections when components are missing.
235 vecs = []
236 for ch in in_channels:
237 if ch is None:
238 vecs.append(None)
239 else:
240 vec = getattr(ch, to)
241 if use_gains:
242 vec /= ch.gain
243 vecs.append(vec)
245 fill_orthogonal(vecs)
246 if not are_orthogonal(vecs):
247 raise ChannelsNotOrthogonal(
248 'components are not orthogonal: station %s.%s.%s, '
249 'channels %s, %s, %s'
250 % (self.nsl() + tuple(in_channel_names)))
252 m = num.hstack([vec2[:, num.newaxis] for vec2 in vecs])
254 m = num.where(num.abs(m) < num.max(num.abs(m))*1e-16, 0., m)
256 if to == 'ned':
257 out_channels = [
258 Channel(out_channel_names[0], 0., 0., 1.),
259 Channel(out_channel_names[1], 90., 0., 1.),
260 Channel(out_channel_names[2], 0., 90., 1.)]
262 elif to == 'enu':
263 out_channels = [
264 Channel(out_channel_names[0], 90., 0., 1.),
265 Channel(out_channel_names[1], 0., 0., 1.),
266 Channel(out_channel_names[2], 0., -90., 1.)]
268 return m, in_channels, out_channels
270 def guess_channel_groups(self):
271 cg = {}
272 for channel in self.get_channels():
273 if len(channel.name) >= 1:
274 kind = channel.name[:-1]
275 if kind not in cg:
276 cg[kind] = []
277 cg[kind].append(channel.name[-1])
279 def allin(a, b):
280 return all(x in b for x in a)
282 out_groups = []
283 for kind, components in cg.items():
284 for sys in ('ENZ', '12Z', 'XYZ', 'RTZ', '123'):
285 if allin(sys, components):
286 out_groups.append(tuple([kind+c for c in sys]))
288 return out_groups
290 def guess_projections_to_enu(self, out_channels=('E', 'N', 'U'), **kwargs):
291 proj = []
292 for cg in self.guess_channel_groups():
293 try:
294 proj.append(self.projection_to_enu(
295 cg, out_channels=out_channels, **kwargs))
297 except ChannelsNotOrthogonal as e:
298 logger.warning(str(e))
300 return proj
302 def guess_projections_to_rtu(
303 self, out_channels=('R', 'T', 'U'), backazimuth=None, **kwargs):
305 if backazimuth is None:
306 backazimuth = self.backazimuth
307 out_channels_ = [
308 Channel(
309 out_channels[0], wrap(backazimuth+180., -180., 180.), 0., 1.),
310 Channel(
311 out_channels[1], wrap(backazimuth+270., -180., 180.), 0., 1.),
312 Channel(
313 out_channels[2], 0., -90., 1.)]
315 proj = []
316 for (m, in_channels, _) in self.guess_projections_to_enu(**kwargs):
317 phi = (backazimuth + 180.)*d2r
318 r = num.array([[math.sin(phi), math.cos(phi), 0.0],
319 [math.cos(phi), -math.sin(phi), 0.0],
320 [0.0, 0.0, 1.0]])
321 proj.append((num.dot(r, m), in_channels, out_channels_))
323 return proj
325 def projection_to_enu(
326 self,
327 in_channels,
328 out_channels=('E', 'N', 'U'),
329 **kwargs):
331 return self._projection_to('enu', in_channels, out_channels, **kwargs)
333 def projection_to_ned(
334 self,
335 in_channels,
336 out_channels=('N', 'E', 'D'),
337 **kwargs):
339 return self._projection_to('ned', in_channels, out_channels, **kwargs)
341 def projection_from_enu(
342 self,
343 in_channels=('E', 'N', 'U'),
344 out_channels=('X', 'Y', 'Z'),
345 **kwargs):
347 m, out_channels, in_channels = self._projection_to(
348 'enu', out_channels, in_channels, **kwargs)
350 return num.linalg.inv(m), in_channels, out_channels
352 def projection_from_ned(
353 self,
354 in_channels=('N', 'E', 'D'),
355 out_channels=('X', 'Y', 'Z'),
356 **kwargs):
358 m, out_channels, in_channels = self._projection_to(
359 'ned', out_channels, in_channels, **kwargs)
361 return num.linalg.inv(m), in_channels, out_channels
363 def nsl_string(self):
364 return '.'.join((self.network, self.station, self.location))
366 def nsl(self):
367 return self.network, self.station, self.location
369 def cannot_handle_offsets(self):
370 if self.north_shift != 0.0 or self.east_shift != 0.0:
371 logger.warning(
372 'Station %s.%s.%s has a non-zero Cartesian offset. Such '
373 'offsets cannot be saved in the basic station file format. '
374 'Effective lat/lons are saved only. Please save the stations '
375 'in YAML format to preserve the reference-and-offset '
376 'coordinates.' % self.nsl())
378 def oldstr(self):
379 self.cannot_handle_offsets()
380 nsl = '%s.%s.%s' % (self.network, self.station, self.location)
381 s = '%-15s %14.5f %14.5f %14.1f %14.1f %s' % (
382 nsl, self.effective_lat, self.effective_lon, self.elevation,
383 self.depth, self.name)
384 return s
387def dump_stations(stations, filename):
388 '''
389 Write stations file.
391 :param stations: list of :py:class:`Station` objects
392 :param filename: filename as str
393 '''
394 f = open(filename, 'w')
395 for sta in stations:
396 f.write(sta.oldstr()+'\n')
397 for cha in sta.get_channels():
398 azimuth = 'NaN'
399 if cha.azimuth is not None:
400 azimuth = '%7g' % cha.azimuth
402 dip = 'NaN'
403 if cha.dip is not None:
404 dip = '%7g' % cha.dip
406 f.write('%5s %14s %14s %14g\n' % (
407 cha.name, azimuth, dip, cha.gain))
409 f.close()
412def dump_stations_yaml(stations, filename):
413 '''
414 Write stations file in YAML format.
416 :param stations: list of :py:class:`Station` objects
417 :param filename: filename as str
418 '''
420 dump_all(stations, filename=filename)
423def float_or_none(s):
424 if s is None:
425 return None
426 elif isinstance(s, str) and s.lower() == 'nan':
427 return None
428 else:
429 return float(s)
432def detect_format(filename):
433 with open(filename, 'r') as f:
434 for line in f:
435 line = line.strip()
436 if not line or line.startswith('#') or line.startswith('%'):
437 continue
438 if line.startswith('--- !pf.Station'):
439 return 'yaml'
440 else:
441 return 'basic'
443 return 'basic'
446def load_stations(filename, format='detect'):
447 '''
448 Read stations file.
450 :param filename: filename
451 :returns: list of :py:class:`Station` objects
452 '''
454 if format == 'detect':
455 format = detect_format(filename)
457 if format == 'yaml':
458 from pyrocko import guts
459 stations = [
460 st for st in guts.load_all(filename=filename)
461 if isinstance(st, Station)]
463 return stations
465 elif format == 'basic':
466 stations = []
467 f = open(filename, 'r')
468 station = None
469 channel_names = []
470 for (iline, line) in enumerate(f):
471 toks = line.split(None, 5)
472 if line.strip().startswith('#') or line.strip() == '':
473 continue
475 if len(toks) == 5 or len(toks) == 6:
476 net, sta, loc = toks[0].split('.')
477 lat, lon, elevation, depth = [float(x) for x in toks[1:5]]
478 if len(toks) == 5:
479 name = ''
480 else:
481 name = toks[5].rstrip()
483 station = Station(
484 net, sta, loc, lat, lon,
485 elevation=elevation, depth=depth, name=name)
487 stations.append(station)
488 channel_names = []
490 elif len(toks) == 4 and station is not None:
491 name, azi, dip, gain = (
492 toks[0],
493 float_or_none(toks[1]),
494 float_or_none(toks[2]),
495 float(toks[3]))
496 if name in channel_names:
497 logger.warning(
498 'redefined channel! (line: %i, file: %s)' %
499 (iline + 1, filename))
500 else:
501 channel_names.append(name)
503 channel = Channel(name, azimuth=azi, dip=dip, gain=gain)
504 station.add_channel(channel)
506 else:
507 logger.warning('skipping invalid station/channel definition '
508 '(line: %i, file: %s' % (iline + 1, filename))
510 f.close()
511 return stations
513 else:
514 from pyrocko.io.io_common import FileLoadError
515 raise FileLoadError('unknown event file format: %s' % format)
518def dump_kml(objects, filename):
519 station_template = '''
520 <Placemark>
521 <name>%(network)s.%(station)s.%(location)s</name>
522 <description></description>
523 <styleUrl>#msn_S</styleUrl>
524 <Point>
525 <coordinates>%(elon)f,%(elat)f,%(elevation)f</coordinates>
526 </Point>
527 </Placemark>
528'''
530 f = open(filename, 'w')
531 f.write('<?xml version="1.0" encoding="UTF-8"?>\n')
532 f.write('<kml xmlns="http://www.opengis.net/kml/2.2">\n')
533 f.write('<Document>\n')
534 f.write(''' <Style id="sh_S">
535 <IconStyle>
536 <scale>1.3</scale>
537 <Icon>
538 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
539 </Icon>
540 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
541 </IconStyle>
542 <ListStyle>
543 <ItemIcon>
544 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
545 </ItemIcon>
546 </ListStyle>
547 </Style>
548 <Style id="sn_S">
549 <IconStyle>
550 <scale>1.1</scale>
551 <Icon>
552 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
553 </Icon>
554 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
555 </IconStyle>
556 <ListStyle>
557 <ItemIcon>
558 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
559 </ItemIcon>
560 </ListStyle>
561 </Style>
562 <StyleMap id="msn_S">
563 <Pair>
564 <key>normal</key>
565 <styleUrl>#sn_S</styleUrl>
566 </Pair>
567 <Pair>
568 <key>highlight</key>
569 <styleUrl>#sh_S</styleUrl>
570 </Pair>
571 </StyleMap>
572''')
573 for obj in objects:
575 if isinstance(obj, Station):
576 d = obj.__dict__.copy()
577 d['elat'], d['elon'] = obj.effective_latlon
578 f.write(station_template % d)
580 f.write('</Document>')
581 f.write('</kml>\n')
582 f.close()