1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import copy
8import logging
9import numpy as num
11from pyrocko import orthodrome
12from pyrocko.orthodrome import wrap
13from pyrocko.guts import Object, Float, String, List, dump_all
15from .location import Location
17logger = logging.getLogger('pyrocko.model.station')
19guts_prefix = 'pf'
21d2r = num.pi / 180.
24class ChannelsNotOrthogonal(Exception):
25 pass
28def guess_azimuth_from_name(channel_name):
29 if channel_name.endswith('N'):
30 return 0.
31 elif channel_name.endswith('E'):
32 return 90.
33 elif channel_name.endswith('Z'):
34 return 0.
36 return None
39def guess_dip_from_name(channel_name):
40 if channel_name.endswith('N'):
41 return 0.
42 elif channel_name.endswith('E'):
43 return 0.
44 elif channel_name.endswith('Z'):
45 return -90.
47 return None
50def guess_azimuth_dip_from_name(channel_name):
51 return guess_azimuth_from_name(channel_name), \
52 guess_dip_from_name(channel_name)
55def mkvec(x, y, z):
56 return num.array([x, y, z], dtype=float)
59def are_orthogonal(enus, eps=0.05):
60 return all(abs(x) < eps for x in [
61 num.dot(enus[0], enus[1]),
62 num.dot(enus[1], enus[2]),
63 num.dot(enus[2], enus[0])])
66def fill_orthogonal(enus):
68 nmiss = sum(x is None for x in enus)
70 if nmiss == 1:
71 for ic in range(len(enus)):
72 if enus[ic] is None:
73 enus[ic] = num.cross(enus[(ic-2) % 3], enus[(ic-1) % 3])
75 if nmiss == 2:
76 for ic in range(len(enus)):
77 if enus[ic] is not None:
78 xenu = enus[ic] + mkvec(1, 1, 1)
79 enus[(ic+1) % 3] = num.cross(enus[ic], xenu)
80 enus[(ic+2) % 3] = num.cross(enus[ic], enus[(ic+1) % 3])
82 if nmiss == 3:
83 # holy camoly..
84 enus[0] = mkvec(1, 0, 0)
85 enus[1] = mkvec(0, 1, 0)
86 enus[2] = mkvec(0, 0, 1)
89class Channel(Object):
90 name = String.T()
91 azimuth = Float.T(optional=True)
92 dip = Float.T(optional=True)
93 gain = Float.T(default=1.0)
95 def __init__(self, name, azimuth=None, dip=None, gain=1.0):
96 if azimuth is None:
97 azimuth = guess_azimuth_from_name(name)
98 if dip is None:
99 dip = guess_dip_from_name(name)
101 Object.__init__(
102 self,
103 name=name,
104 azimuth=float_or_none(azimuth),
105 dip=float_or_none(dip),
106 gain=float(gain))
108 @property
109 def ned(self):
110 if None in (self.azimuth, self.dip):
111 return None
113 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
114 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
115 d = math.sin(self.dip*d2r)
116 return mkvec(n, e, d)
118 @property
119 def enu(self):
120 if None in (self.azimuth, self.dip):
121 return None
123 n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
124 e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
125 d = math.sin(self.dip*d2r)
126 return mkvec(e, n, -d)
128 def __str__(self):
129 return '%s %f %f %g' % (self.name, self.azimuth, self.dip, self.gain)
132class Station(Location):
133 network = String.T()
134 station = String.T()
135 location = String.T()
136 name = String.T(default='')
137 channels = List.T(Channel.T())
139 def __init__(self, network='', station='', location='',
140 lat=0.0, lon=0.0,
141 elevation=0.0, depth=0.0,
142 north_shift=0.0, east_shift=0.0,
143 name='', channels=None):
145 Location.__init__(
146 self,
147 network=network, station=station, location=location,
148 lat=float(lat), lon=float(lon),
149 elevation=float(elevation),
150 depth=float(depth),
151 north_shift=float(north_shift),
152 east_shift=float(east_shift),
153 name=name or '',
154 channels=channels or [])
156 self.dist_deg = None
157 self.dist_m = None
158 self.azimuth = None
159 self.backazimuth = None
161 def copy(self):
162 return copy.deepcopy(self)
164 def set_event_relative_data(self, event, distance_3d=False):
165 surface_dist = self.distance_to(event)
166 if distance_3d:
167 if event.depth is None:
168 logger.warning('No event depth given: using 0 m.')
169 dd = 0.0 - self.depth
170 else:
171 dd = event.depth - self.depth
173 self.dist_m = math.sqrt(dd**2 + surface_dist**2)
174 else:
175 self.dist_m = surface_dist
177 self.dist_deg = surface_dist / orthodrome.earthradius_equator * \
178 orthodrome.r2d
180 self.azimuth, self.backazimuth = event.azibazi_to(self)
182 def set_channels_by_name(self, *args):
183 self.set_channels([])
184 for name in args:
185 self.add_channel(Channel(name))
187 def set_channels(self, channels):
188 self.channels = []
189 for ch in channels:
190 self.add_channel(ch)
192 def get_channels(self):
193 return list(self.channels)
195 def get_channel_names(self):
196 return set(ch.name for ch in self.channels)
198 def remove_channel_by_name(self, name):
199 todel = [ch for ch in self.channels if ch.name == name]
200 for ch in todel:
201 self.channels.remove(ch)
203 def add_channel(self, channel):
204 self.remove_channel_by_name(channel.name)
205 self.channels.append(channel)
206 self.channels.sort(key=lambda ch: ch.name)
208 def get_channel(self, name):
209 for ch in self.channels:
210 if ch.name == name:
211 return ch
213 return None
215 def rotation_ne_to_rt(self, in_channel_names, out_channel_names):
217 angle = wrap(self.backazimuth + 180., -180., 180.)
218 in_channels = [self.get_channel(name) for name in in_channel_names]
219 out_channels = [
220 Channel(out_channel_names[0],
221 wrap(self.backazimuth+180., -180., 180.), 0., 1.),
222 Channel(out_channel_names[1],
223 wrap(self.backazimuth+270., -180., 180.), 0., 1.)]
224 return angle, in_channels, out_channels
226 def _projection_to(
227 self, to, in_channel_names, out_channel_names, use_gains=False):
229 in_channels = [self.get_channel(name) for name in in_channel_names]
231 # create orthogonal vectors for missing components, such that this
232 # won't break projections when components are missing.
234 vecs = []
235 for ch in in_channels:
236 if ch is None:
237 vecs.append(None)
238 else:
239 vec = getattr(ch, to)
240 if use_gains:
241 vec /= ch.gain
242 vecs.append(vec)
244 fill_orthogonal(vecs)
245 if not are_orthogonal(vecs):
246 raise ChannelsNotOrthogonal(
247 'components are not orthogonal: station %s.%s.%s, '
248 'channels %s, %s, %s'
249 % (self.nsl() + tuple(in_channel_names)))
251 m = num.hstack([vec2[:, num.newaxis] for vec2 in vecs])
253 m = num.where(num.abs(m) < num.max(num.abs(m))*1e-16, 0., m)
255 if to == 'ned':
256 out_channels = [
257 Channel(out_channel_names[0], 0., 0., 1.),
258 Channel(out_channel_names[1], 90., 0., 1.),
259 Channel(out_channel_names[2], 0., 90., 1.)]
261 elif to == 'enu':
262 out_channels = [
263 Channel(out_channel_names[0], 90., 0., 1.),
264 Channel(out_channel_names[1], 0., 0., 1.),
265 Channel(out_channel_names[2], 0., -90., 1.)]
267 return m, in_channels, out_channels
269 def guess_channel_groups(self):
270 cg = {}
271 for channel in self.get_channels():
272 if len(channel.name) >= 1:
273 kind = channel.name[:-1]
274 if kind not in cg:
275 cg[kind] = []
276 cg[kind].append(channel.name[-1])
278 def allin(a, b):
279 return all(x in b for x in a)
281 out_groups = []
282 for kind, components in cg.items():
283 for sys in ('ENZ', '12Z', 'XYZ', 'RTZ', '123'):
284 if allin(sys, components):
285 out_groups.append(tuple([kind+c for c in sys]))
287 return out_groups
289 def guess_projections_to_enu(self, out_channels=('E', 'N', 'U'), **kwargs):
290 proj = []
291 for cg in self.guess_channel_groups():
292 try:
293 proj.append(self.projection_to_enu(
294 cg, out_channels=out_channels, **kwargs))
296 except ChannelsNotOrthogonal as e:
297 logger.warning(str(e))
299 return proj
301 def guess_projections_to_rtu(
302 self, out_channels=('R', 'T', 'U'), backazimuth=None, **kwargs):
304 if backazimuth is None:
305 backazimuth = self.backazimuth
306 out_channels_ = [
307 Channel(
308 out_channels[0], wrap(backazimuth+180., -180., 180.), 0., 1.),
309 Channel(
310 out_channels[1], wrap(backazimuth+270., -180., 180.), 0., 1.),
311 Channel(
312 out_channels[2], 0., -90., 1.)]
314 proj = []
315 for (m, in_channels, _) in self.guess_projections_to_enu(**kwargs):
316 phi = (backazimuth + 180.)*d2r
317 r = num.array([[math.sin(phi), math.cos(phi), 0.0],
318 [math.cos(phi), -math.sin(phi), 0.0],
319 [0.0, 0.0, 1.0]])
320 proj.append((num.dot(r, m), in_channels, out_channels_))
322 return proj
324 def projection_to_enu(
325 self,
326 in_channels,
327 out_channels=('E', 'N', 'U'),
328 **kwargs):
330 return self._projection_to('enu', in_channels, out_channels, **kwargs)
332 def projection_to_ned(
333 self,
334 in_channels,
335 out_channels=('N', 'E', 'D'),
336 **kwargs):
338 return self._projection_to('ned', in_channels, out_channels, **kwargs)
340 def projection_from_enu(
341 self,
342 in_channels=('E', 'N', 'U'),
343 out_channels=('X', 'Y', 'Z'),
344 **kwargs):
346 m, out_channels, in_channels = self._projection_to(
347 'enu', out_channels, in_channels, **kwargs)
349 return num.linalg.inv(m), in_channels, out_channels
351 def projection_from_ned(
352 self,
353 in_channels=('N', 'E', 'D'),
354 out_channels=('X', 'Y', 'Z'),
355 **kwargs):
357 m, out_channels, in_channels = self._projection_to(
358 'ned', out_channels, in_channels, **kwargs)
360 return num.linalg.inv(m), in_channels, out_channels
362 def nsl_string(self):
363 return '.'.join((self.network, self.station, self.location))
365 def nsl(self):
366 return self.network, self.station, self.location
368 def cannot_handle_offsets(self):
369 if self.north_shift != 0.0 or self.east_shift != 0.0:
370 logger.warning(
371 'Station %s.%s.%s has a non-zero Cartesian offset. Such '
372 'offsets cannot be saved in the basic station file format. '
373 'Effective lat/lons are saved only. Please save the stations '
374 'in YAML format to preserve the reference-and-offset '
375 'coordinates.' % self.nsl())
377 def oldstr(self):
378 self.cannot_handle_offsets()
379 nsl = '%s.%s.%s' % (self.network, self.station, self.location)
380 s = '%-15s %14.5f %14.5f %14.1f %14.1f %s' % (
381 nsl, self.effective_lat, self.effective_lon, self.elevation,
382 self.depth, self.name)
383 return s
386def dump_stations(stations, filename):
387 '''
388 Write stations file.
390 :param stations: list of :py:class:`Station` objects
391 :param filename: filename as str
392 '''
393 f = open(filename, 'w')
394 for sta in stations:
395 f.write(sta.oldstr()+'\n')
396 for cha in sta.get_channels():
397 azimuth = 'NaN'
398 if cha.azimuth is not None:
399 azimuth = '%7g' % cha.azimuth
401 dip = 'NaN'
402 if cha.dip is not None:
403 dip = '%7g' % cha.dip
405 f.write('%5s %14s %14s %14g\n' % (
406 cha.name, azimuth, dip, cha.gain))
408 f.close()
411def dump_stations_yaml(stations, filename):
412 '''
413 Write stations file in YAML format.
415 :param stations: list of :py:class:`Station` objects
416 :param filename: filename as str
417 '''
419 dump_all(stations, filename=filename)
422def float_or_none(s):
423 if s is None:
424 return None
425 elif isinstance(s, str) and s.lower() == 'nan':
426 return None
427 else:
428 return float(s)
431def detect_format(filename):
432 with open(filename, 'r') as f:
433 for line in f:
434 line = line.strip()
435 if not line or line.startswith('#') or line.startswith('%'):
436 continue
437 if line.startswith('--- !pf.Station'):
438 return 'yaml'
439 else:
440 return 'basic'
442 return 'basic'
445def load_stations(filename, format='detect'):
446 '''
447 Read stations file.
449 :param filename: filename
450 :returns: list of :py:class:`Station` objects
451 '''
453 if format == 'detect':
454 format = detect_format(filename)
456 if format == 'yaml':
457 from pyrocko import guts
458 stations = [
459 st for st in guts.load_all(filename=filename)
460 if isinstance(st, Station)]
462 return stations
464 elif format == 'basic':
465 stations = []
466 f = open(filename, 'r')
467 station = None
468 channel_names = []
469 for (iline, line) in enumerate(f):
470 toks = line.split(None, 5)
471 if line.strip().startswith('#') or line.strip() == '':
472 continue
474 if len(toks) == 5 or len(toks) == 6:
475 net, sta, loc = toks[0].split('.')
476 lat, lon, elevation, depth = [float(x) for x in toks[1:5]]
477 if len(toks) == 5:
478 name = ''
479 else:
480 name = toks[5].rstrip()
482 station = Station(
483 net, sta, loc, lat, lon,
484 elevation=elevation, depth=depth, name=name)
486 stations.append(station)
487 channel_names = []
489 elif len(toks) == 4 and station is not None:
490 name, azi, dip, gain = (
491 toks[0],
492 float_or_none(toks[1]),
493 float_or_none(toks[2]),
494 float(toks[3]))
495 if name in channel_names:
496 logger.warning(
497 'redefined channel! (line: %i, file: %s)' %
498 (iline + 1, filename))
499 else:
500 channel_names.append(name)
502 channel = Channel(name, azimuth=azi, dip=dip, gain=gain)
503 station.add_channel(channel)
505 else:
506 logger.warning('skipping invalid station/channel definition '
507 '(line: %i, file: %s' % (iline + 1, filename))
509 f.close()
510 return stations
512 else:
513 from pyrocko.io.io_common import FileLoadError
514 raise FileLoadError('unknown event file format: %s' % format)
517def dump_kml(objects, filename):
518 station_template = '''
519 <Placemark>
520 <name>%(network)s.%(station)s.%(location)s</name>
521 <description></description>
522 <styleUrl>#msn_S</styleUrl>
523 <Point>
524 <coordinates>%(elon)f,%(elat)f,%(elevation)f</coordinates>
525 </Point>
526 </Placemark>
527'''
529 f = open(filename, 'w')
530 f.write('<?xml version="1.0" encoding="UTF-8"?>\n')
531 f.write('<kml xmlns="http://www.opengis.net/kml/2.2">\n')
532 f.write('<Document>\n')
533 f.write(''' <Style id="sh_S">
534 <IconStyle>
535 <scale>1.3</scale>
536 <Icon>
537 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
538 </Icon>
539 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
540 </IconStyle>
541 <ListStyle>
542 <ItemIcon>
543 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
544 </ItemIcon>
545 </ListStyle>
546 </Style>
547 <Style id="sn_S">
548 <IconStyle>
549 <scale>1.1</scale>
550 <Icon>
551 <href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
552 </Icon>
553 <hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
554 </IconStyle>
555 <ListStyle>
556 <ItemIcon>
557 <href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
558 </ItemIcon>
559 </ListStyle>
560 </Style>
561 <StyleMap id="msn_S">
562 <Pair>
563 <key>normal</key>
564 <styleUrl>#sn_S</styleUrl>
565 </Pair>
566 <Pair>
567 <key>highlight</key>
568 <styleUrl>#sh_S</styleUrl>
569 </Pair>
570 </StyleMap>
571''')
572 for obj in objects:
574 if isinstance(obj, Station):
575 d = obj.__dict__.copy()
576 d['elat'], d['elon'] = obj.effective_latlon
577 f.write(station_template % d)
579 f.write('</Document>')
580 f.write('</kml>\n')
581 f.close()