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.
36 elif channel_name.endswith('U'):
37 return 0.
39 return None
42def guess_dip_from_name(channel_name):
43 if channel_name.endswith('N'):
44 return 0.
45 elif channel_name.endswith('E'):
46 return 0.
47 elif channel_name.endswith('Z'):
48 return -90.
49 elif channel_name.endswith('U'):
50 return -90.
52 return None
55def guess_azimuth_dip_from_name(channel_name):
56 return guess_azimuth_from_name(channel_name), \
57 guess_dip_from_name(channel_name)
60def mkvec(x, y, z):
61 return num.array([x, y, z], dtype=float)
64def are_orthogonal(enus, eps=0.05):
65 return all(abs(x) < eps for x in [
66 num.dot(enus[0], enus[1]),
67 num.dot(enus[1], enus[2]),
68 num.dot(enus[2], enus[0])])
71def fill_orthogonal(enus):
73 nmiss = sum(x is None for x in enus)
75 if nmiss == 1:
76 for ic in range(len(enus)):
77 if enus[ic] is None:
78 enus[ic] = num.cross(enus[(ic-2) % 3], enus[(ic-1) % 3])
80 if nmiss == 2:
81 for ic in range(len(enus)):
82 if enus[ic] is not None:
83 xenu = enus[ic] + mkvec(1, 1, 1)
84 enus[(ic+1) % 3] = num.cross(enus[ic], xenu)
85 enus[(ic+2) % 3] = num.cross(enus[ic], enus[(ic+1) % 3])
87 if nmiss == 3:
88 # holy camoly..
89 enus[0] = mkvec(1, 0, 0)
90 enus[1] = mkvec(0, 1, 0)
91 enus[2] = mkvec(0, 0, 1)
94class Channel(Object):
95 name = String.T()
96 azimuth = Float.T(optional=True)
97 dip = Float.T(optional=True)
98 gain = Float.T(default=1.0)
100 def __init__(self, name, azimuth=None, dip=None, gain=1.0):
101 if azimuth is None:
102 azimuth = guess_azimuth_from_name(name)
103 if dip is None:
104 dip = guess_dip_from_name(name)
106 Object.__init__(
107 self,
108 name=name,
109 azimuth=float_or_none(azimuth),
110 dip=float_or_none(dip),
111 gain=float(gain))
113 @property
114 def ned(self):
115 if None in (self.azimuth, self.dip):
116 return None
118 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
119 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
120 d = math.sin(self.dip*d2r)
121 return mkvec(n, e, d)
123 @property
124 def enu(self):
125 if None in (self.azimuth, self.dip):
126 return None
128 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
129 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
130 d = math.sin(self.dip*d2r)
131 return mkvec(e, n, -d)
133 def __str__(self):
134 return '%s %f %f %g' % (self.name, self.azimuth, self.dip, self.gain)
137class Station(Location):
138 network = String.T()
139 station = String.T()
140 location = String.T()
141 name = String.T(default='')
142 channels = List.T(Channel.T())
144 def __init__(self, network='', station='', location='',
145 lat=0.0, lon=0.0,
146 elevation=0.0, depth=0.0,
147 north_shift=0.0, east_shift=0.0,
148 name='', channels=None):
150 Location.__init__(
151 self,
152 network=network, station=station, location=location,
153 lat=float(lat), lon=float(lon),
154 elevation=float(elevation),
155 depth=float(depth),
156 north_shift=float(north_shift),
157 east_shift=float(east_shift),
158 name=name or '',
159 channels=channels or [])
161 self.dist_deg = None
162 self.dist_m = None
163 self.azimuth = None
164 self.backazimuth = None
166 def copy(self):
167 return copy.deepcopy(self)
169 def set_event_relative_data(self, event, distance_3d=False):
170 surface_dist = self.distance_to(event)
171 if distance_3d:
172 if event.depth is None:
173 logger.warning('No event depth given: using 0 m.')
174 dd = 0.0 - self.depth
175 else:
176 dd = event.depth - self.depth
178 self.dist_m = math.sqrt(dd**2 + surface_dist**2)
179 else:
180 self.dist_m = surface_dist
182 self.dist_deg = surface_dist / orthodrome.earthradius_equator * \
183 orthodrome.r2d
185 self.azimuth, self.backazimuth = event.azibazi_to(self)
187 def set_channels_by_name(self, *args):
188 self.set_channels([])
189 for name in args:
190 self.add_channel(Channel(name))
192 def set_channels(self, channels):
193 self.channels = []
194 for ch in channels:
195 self.add_channel(ch)
197 def get_channels(self):
198 return list(self.channels)
200 def get_channel_names(self):
201 return set(ch.name for ch in self.channels)
203 def remove_channel_by_name(self, name):
204 todel = [ch for ch in self.channels if ch.name == name]
205 for ch in todel:
206 self.channels.remove(ch)
208 def add_channel(self, channel):
209 self.remove_channel_by_name(channel.name)
210 self.channels.append(channel)
211 self.channels.sort(key=lambda ch: ch.name)
213 def get_channel(self, name):
214 for ch in self.channels:
215 if ch.name == name:
216 return ch
218 return None
220 def rotation_ne_to_rt(self, in_channel_names, out_channel_names):
222 angle = wrap(self.backazimuth + 180., -180., 180.)
223 in_channels = [self.get_channel(name) for name in in_channel_names]
224 out_channels = [
225 Channel(out_channel_names[0],
226 wrap(self.backazimuth+180., -180., 180.), 0., 1.),
227 Channel(out_channel_names[1],
228 wrap(self.backazimuth+270., -180., 180.), 0., 1.)]
229 return angle, in_channels, out_channels
231 def _projection_to(
232 self, to, in_channel_names, out_channel_names, use_gains=False):
234 in_channels = [self.get_channel(name) for name in in_channel_names]
236 # create orthogonal vectors for missing components, such that this
237 # won't break projections when components are missing.
239 vecs = []
240 for ch in in_channels:
241 if ch is None:
242 vecs.append(None)
243 else:
244 vec = getattr(ch, to)
245 if use_gains:
246 vec /= ch.gain
247 vecs.append(vec)
249 fill_orthogonal(vecs)
250 if not are_orthogonal(vecs):
251 raise ChannelsNotOrthogonal(
252 'components are not orthogonal: station %s.%s.%s, '
253 'channels %s, %s, %s'
254 % (self.nsl() + tuple(in_channel_names)))
256 m = num.hstack([vec2[:, num.newaxis] for vec2 in vecs])
258 m = num.where(num.abs(m) < num.max(num.abs(m))*1e-16, 0., m)
260 if to == 'ned':
261 out_channels = [
262 Channel(out_channel_names[0], 0., 0., 1.),
263 Channel(out_channel_names[1], 90., 0., 1.),
264 Channel(out_channel_names[2], 0., 90., 1.)]
266 elif to == 'enu':
267 out_channels = [
268 Channel(out_channel_names[0], 90., 0., 1.),
269 Channel(out_channel_names[1], 0., 0., 1.),
270 Channel(out_channel_names[2], 0., -90., 1.)]
272 return m, in_channels, out_channels
274 def guess_channel_groups(self):
275 cg = {}
276 for channel in self.get_channels():
277 if len(channel.name) >= 1:
278 kind = channel.name[:-1]
279 if kind not in cg:
280 cg[kind] = []
281 cg[kind].append(channel.name[-1])
283 def allin(a, b):
284 return all(x in b for x in a)
286 out_groups = []
287 for kind, components in cg.items():
288 for sys in ('ENU', 'ENZ', '12Z', 'XYZ', 'RTZ', '123'):
289 if allin(sys, components):
290 out_groups.append(tuple([kind+c for c in sys]))
292 return out_groups
294 def guess_projections_to_enu(self, out_channels=('E', 'N', 'U'), **kwargs):
295 proj = []
296 for cg in self.guess_channel_groups():
297 try:
298 proj.append(self.projection_to_enu(
299 cg, out_channels=out_channels, **kwargs))
301 except ChannelsNotOrthogonal as e:
302 logger.warning(str(e))
304 return proj
306 def guess_projections_to_rtu(
307 self, out_channels=('R', 'T', 'U'), backazimuth=None, **kwargs):
309 if backazimuth is None:
310 backazimuth = self.backazimuth
311 out_channels_ = [
312 Channel(
313 out_channels[0], wrap(backazimuth+180., -180., 180.), 0., 1.),
314 Channel(
315 out_channels[1], wrap(backazimuth+270., -180., 180.), 0., 1.),
316 Channel(
317 out_channels[2], 0., -90., 1.)]
319 proj = []
320 for (m, in_channels, _) in self.guess_projections_to_enu(**kwargs):
321 phi = (backazimuth + 180.)*d2r
322 r = num.array([[math.sin(phi), math.cos(phi), 0.0],
323 [math.cos(phi), -math.sin(phi), 0.0],
324 [0.0, 0.0, 1.0]])
325 proj.append((num.dot(r, m), in_channels, out_channels_))
327 return proj
329 def projection_to_enu(
330 self,
331 in_channels,
332 out_channels=('E', 'N', 'U'),
333 **kwargs):
335 return self._projection_to('enu', in_channels, out_channels, **kwargs)
337 def projection_to_ned(
338 self,
339 in_channels,
340 out_channels=('N', 'E', 'D'),
341 **kwargs):
343 return self._projection_to('ned', in_channels, out_channels, **kwargs)
345 def projection_from_enu(
346 self,
347 in_channels=('E', 'N', 'U'),
348 out_channels=('X', 'Y', 'Z'),
349 **kwargs):
351 m, out_channels, in_channels = self._projection_to(
352 'enu', out_channels, in_channels, **kwargs)
354 return num.linalg.inv(m), in_channels, out_channels
356 def projection_from_ned(
357 self,
358 in_channels=('N', 'E', 'D'),
359 out_channels=('X', 'Y', 'Z'),
360 **kwargs):
362 m, out_channels, in_channels = self._projection_to(
363 'ned', out_channels, in_channels, **kwargs)
365 return num.linalg.inv(m), in_channels, out_channels
367 def nsl_string(self):
368 return '.'.join((self.network, self.station, self.location))
370 def nsl(self):
371 return self.network, self.station, self.location
373 def cannot_handle_offsets(self):
374 if self.north_shift != 0.0 or self.east_shift != 0.0:
375 logger.warning(
376 'Station %s.%s.%s has a non-zero Cartesian offset. Such '
377 'offsets cannot be saved in the basic station file format. '
378 'Effective lat/lons are saved only. Please save the stations '
379 'in YAML format to preserve the reference-and-offset '
380 'coordinates.' % self.nsl())
382 def oldstr(self):
383 self.cannot_handle_offsets()
384 nsl = '%s.%s.%s' % (self.network, self.station, self.location)
385 s = '%-15s %14.5f %14.5f %14.1f %14.1f %s' % (
386 nsl, self.effective_lat, self.effective_lon, self.elevation,
387 self.depth, self.name)
388 return s
391def dump_stations(stations, filename):
392 '''
393 Write stations file.
395 :param stations: list of :py:class:`Station` objects
396 :param filename: filename as str
397 '''
398 f = open(filename, 'w')
399 for sta in stations:
400 f.write(sta.oldstr()+'\n')
401 for cha in sta.get_channels():
402 azimuth = 'NaN'
403 if cha.azimuth is not None:
404 azimuth = '%7g' % cha.azimuth
406 dip = 'NaN'
407 if cha.dip is not None:
408 dip = '%7g' % cha.dip
410 f.write('%5s %14s %14s %14g\n' % (
411 cha.name, azimuth, dip, cha.gain))
413 f.close()
416def dump_stations_yaml(stations, filename):
417 '''
418 Write stations file in YAML format.
420 :param stations: list of :py:class:`Station` objects
421 :param filename: filename as str
422 '''
424 dump_all(stations, filename=filename)
427def float_or_none(s):
428 if s is None:
429 return None
430 elif isinstance(s, str) and s.lower() == 'nan':
431 return None
432 else:
433 return float(s)
436def detect_format(filename):
437 with open(filename, 'r') as f:
438 for line in f:
439 line = line.strip()
440 if not line or line.startswith('#') or line.startswith('%'):
441 continue
442 if line.startswith('--- !pf.Station'):
443 return 'yaml'
444 else:
445 return 'basic'
447 return 'basic'
450def load_stations(filename, format='detect'):
451 '''
452 Read stations file.
454 :param filename: filename
455 :returns: list of :py:class:`Station` objects
456 '''
458 if format == 'detect':
459 format = detect_format(filename)
461 if format == 'yaml':
462 from pyrocko import guts
463 stations = [
464 st for st in guts.load_all(filename=filename)
465 if isinstance(st, Station)]
467 return stations
469 elif format == 'basic':
470 stations = []
471 f = open(filename, 'r')
472 station = None
473 channel_names = []
474 for (iline, line) in enumerate(f):
475 toks = line.split(None, 5)
476 if line.strip().startswith('#') or line.strip() == '':
477 continue
479 if len(toks) == 5 or len(toks) == 6:
480 net, sta, loc = toks[0].split('.')
481 lat, lon, elevation, depth = [float(x) for x in toks[1:5]]
482 if len(toks) == 5:
483 name = ''
484 else:
485 name = toks[5].rstrip()
487 station = Station(
488 net, sta, loc, lat, lon,
489 elevation=elevation, depth=depth, name=name)
491 stations.append(station)
492 channel_names = []
494 elif len(toks) == 4 and station is not None:
495 name, azi, dip, gain = (
496 toks[0],
497 float_or_none(toks[1]),
498 float_or_none(toks[2]),
499 float(toks[3]))
500 if name in channel_names:
501 logger.warning(
502 'redefined channel! (line: %i, file: %s)' %
503 (iline + 1, filename))
504 else:
505 channel_names.append(name)
507 channel = Channel(name, azimuth=azi, dip=dip, gain=gain)
508 station.add_channel(channel)
510 else:
511 logger.warning('skipping invalid station/channel definition '
512 '(line: %i, file: %s' % (iline + 1, filename))
514 f.close()
515 return stations
517 else:
518 from pyrocko.io.io_common import FileLoadError
519 raise FileLoadError('unknown event file format: %s' % format)
522def dump_kml(objects, filename):
523 station_template = '''
524 <Placemark>
525 <name>%(network)s.%(station)s.%(location)s</name>
526 <description></description>
527 <styleUrl>#msn_S</styleUrl>
528 <Point>
529 <coordinates>%(elon)f,%(elat)f,%(elevation)f</coordinates>
530 </Point>
531 </Placemark>
532'''
534 f = open(filename, 'w')
535 f.write('<?xml version="1.0" encoding="UTF-8"?>\n')
536 f.write('<kml xmlns="http://www.opengis.net/kml/2.2">\n')
537 f.write('<Document>\n')
538 f.write(''' <Style id="sh_S">
539 <IconStyle>
540 <scale>1.3</scale>
541 <Icon>
542 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
543 </Icon>
544 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
545 </IconStyle>
546 <ListStyle>
547 <ItemIcon>
548 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
549 </ItemIcon>
550 </ListStyle>
551 </Style>
552 <Style id="sn_S">
553 <IconStyle>
554 <scale>1.1</scale>
555 <Icon>
556 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
557 </Icon>
558 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
559 </IconStyle>
560 <ListStyle>
561 <ItemIcon>
562 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
563 </ItemIcon>
564 </ListStyle>
565 </Style>
566 <StyleMap id="msn_S">
567 <Pair>
568 <key>normal</key>
569 <styleUrl>#sn_S</styleUrl>
570 </Pair>
571 <Pair>
572 <key>highlight</key>
573 <styleUrl>#sh_S</styleUrl>
574 </Pair>
575 </StyleMap>
576''')
577 for obj in objects:
579 if isinstance(obj, Station):
580 d = obj.__dict__.copy()
581 d['elat'], d['elon'] = obj.effective_latlon
582 f.write(station_template % d)
584 f.write('</Document>')
585 f.write('</kml>\n')
586 f.close()