|
# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
from __future__ import absolute_import, division
import math
import copy
import logging
import numpy as num
from pyrocko import orthodrome
from pyrocko.orthodrome import wrap
from pyrocko.guts import Object, Float, String, List, dump_all
from .location import Location
logger = logging.getLogger('pyrocko.model.station')
guts_prefix = 'pf'
d2r = num.pi / 180.
class ChannelsNotOrthogonal(Exception):
pass
def guess_azimuth_from_name(channel_name):
if channel_name.endswith('N'):
return 0.
elif channel_name.endswith('E'):
return 90.
elif channel_name.endswith('Z'):
return 0.
return None
def guess_dip_from_name(channel_name):
if channel_name.endswith('N'):
return 0.
elif channel_name.endswith('E'):
return 0.
elif channel_name.endswith('Z'):
return -90.
return None
def guess_azimuth_dip_from_name(channel_name):
return guess_azimuth_from_name(channel_name), \
guess_dip_from_name(channel_name)
def mkvec(x, y, z):
return num.array([x, y, z], dtype=num.float)
def are_orthogonal(enus, eps=0.05):
return all(abs(x) < eps for x in [
num.dot(enus[0], enus[1]),
num.dot(enus[1], enus[2]),
num.dot(enus[2], enus[0])])
def fill_orthogonal(enus):
nmiss = sum(x is None for x in enus)
if nmiss == 1:
for ic in range(len(enus)):
if enus[ic] is None:
enus[ic] = num.cross(enus[(ic-2) % 3], enus[(ic-1) % 3])
if nmiss == 2:
for ic in range(len(enus)):
if enus[ic] is not None:
xenu = enus[ic] + mkvec(1, 1, 1)
enus[(ic+1) % 3] = num.cross(enus[ic], xenu)
enus[(ic+2) % 3] = num.cross(enus[ic], enus[(ic+1) % 3])
if nmiss == 3:
# holy camoly..
enus[0] = mkvec(1, 0, 0)
enus[1] = mkvec(0, 1, 0)
enus[2] = mkvec(0, 0, 1)
class Channel(Object):
name = String.T()
azimuth = Float.T(optional=True)
dip = Float.T(optional=True)
gain = Float.T(default=1.0)
def __init__(self, name, azimuth=None, dip=None, gain=1.0):
if azimuth is None:
azimuth = guess_azimuth_from_name(name)
if dip is None:
dip = guess_dip_from_name(name)
Object.__init__(
self,
name=name,
azimuth=float_or_none(azimuth),
dip=float_or_none(dip),
gain=float(gain))
@property
def ned(self):
if None in (self.azimuth, self.dip):
return None
n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
d = math.sin(self.dip*d2r)
return mkvec(n, e, d)
@property
def enu(self):
if None in (self.azimuth, self.dip):
return None
n = math.cos(self.azimuth*d2r)*math.cos(self.dip*d2r)
e = math.sin(self.azimuth*d2r)*math.cos(self.dip*d2r)
d = math.sin(self.dip*d2r)
return mkvec(e, n, -d)
def __str__(self):
return '%s %f %f %g' % (self.name, self.azimuth, self.dip, self.gain)
class Station(Location):
network = String.T()
station = String.T()
location = String.T()
name = String.T(default='')
channels = List.T(Channel.T())
def __init__(self, network='', station='', location='',
lat=0.0, lon=0.0,
elevation=0.0, depth=0.0,
north_shift=0.0, east_shift=0.0,
name='', channels=None):
Location.__init__(
self,
network=network, station=station, location=location,
lat=float(lat), lon=float(lon),
elevation=float(elevation),
depth=float(depth),
north_shift=float(north_shift),
east_shift=float(east_shift),
name=name or '',
channels=channels or [])
self.dist_deg = None
self.dist_m = None
self.azimuth = None
self.backazimuth = None
def copy(self):
return copy.deepcopy(self)
def set_event_relative_data(self, event, distance_3d=False):
surface_dist = self.distance_to(event)
if distance_3d:
if event.depth is None:
logger.warn('No event depth given: using 0 m.')
dd = 0.0 - self.depth
else:
dd = event.depth - self.depth
self.dist_m = math.sqrt(dd**2 + surface_dist**2)
else:
self.dist_m = surface_dist
self.dist_deg = surface_dist / orthodrome.earthradius_equator * \
orthodrome.r2d
self.azimuth, self.backazimuth = event.azibazi_to(self)
def set_channels_by_name(self, *args):
self.set_channels([])
for name in args:
self.add_channel(Channel(name))
def set_channels(self, channels):
self.channels = []
for ch in channels:
self.add_channel(ch)
def get_channels(self):
return list(self.channels)
def get_channel_names(self):
return set(ch.name for ch in self.channels)
def remove_channel_by_name(self, name):
todel = [ch for ch in self.channels if ch.name == name]
for ch in todel:
self.channels.remove(ch)
def add_channel(self, channel):
self.remove_channel_by_name(channel.name)
self.channels.append(channel)
self.channels.sort(key=lambda ch: ch.name)
def get_channel(self, name):
for ch in self.channels:
if ch.name == name:
return ch
return None
def rotation_ne_to_rt(self, in_channel_names, out_channel_names):
angle = wrap(self.backazimuth + 180., -180., 180.)
in_channels = [self.get_channel(name) for name in in_channel_names]
out_channels = [
Channel(out_channel_names[0],
wrap(self.backazimuth+180., -180., 180.), 0., 1.),
Channel(out_channel_names[1],
wrap(self.backazimuth+270., -180., 180.), 0., 1.)]
return angle, in_channels, out_channels
def _projection_to(
self, to, in_channel_names, out_channel_names, use_gains=False):
in_channels = [self.get_channel(name) for name in in_channel_names]
# create orthogonal vectors for missing components, such that this
# won't break projections when components are missing.
vecs = []
for ch in in_channels:
if ch is None:
vecs.append(None)
else:
vec = getattr(ch, to)
if use_gains:
vec /= ch.gain
vecs.append(vec)
fill_orthogonal(vecs)
if not are_orthogonal(vecs):
raise ChannelsNotOrthogonal(
'components are not orthogonal: station %s.%s.%s, '
'channels %s, %s, %s'
% (self.nsl() + tuple(in_channel_names)))
m = num.hstack([vec2[:, num.newaxis] for vec2 in vecs])
m = num.where(num.abs(m) < num.max(num.abs(m))*1e-16, 0., m)
if to == 'ned':
out_channels = [
Channel(out_channel_names[0], 0., 0., 1.),
Channel(out_channel_names[1], 90., 0., 1.),
Channel(out_channel_names[2], 0., 90., 1.)]
elif to == 'enu':
out_channels = [
Channel(out_channel_names[0], 90., 0., 1.),
Channel(out_channel_names[1], 0., 0., 1.),
Channel(out_channel_names[2], 0., -90., 1.)]
return m, in_channels, out_channels
def guess_channel_groups(self):
cg = {}
for channel in self.get_channels():
if len(channel.name) >= 1:
kind = channel.name[:-1]
if kind not in cg:
cg[kind] = []
cg[kind].append(channel.name[-1])
def allin(a, b):
return all(x in b for x in a)
out_groups = []
for kind, components in cg.items():
for sys in ('ENZ', '12Z', 'XYZ', 'RTZ', '123'):
if allin(sys, components):
out_groups.append(tuple([kind+c for c in sys]))
return out_groups
def guess_projections_to_enu(self, out_channels=('E', 'N', 'U'), **kwargs):
proj = []
for cg in self.guess_channel_groups():
try:
proj.append(self.projection_to_enu(
cg, out_channels=out_channels, **kwargs))
except ChannelsNotOrthogonal as e:
logger.warning(str(e))
return proj
def guess_projections_to_rtu(
self, out_channels=('R', 'T', 'U'), backazimuth=None, **kwargs):
if backazimuth is None:
backazimuth = self.backazimuth
out_channels_ = [
Channel(
out_channels[0], wrap(backazimuth+180., -180., 180.), 0., 1.),
Channel(
out_channels[1], wrap(backazimuth+270., -180., 180.), 0., 1.),
Channel(
out_channels[2], 0., -90., 1.)]
proj = []
for (m, in_channels, _) in self.guess_projections_to_enu(**kwargs):
phi = (backazimuth + 180.)*d2r
r = num.array([[math.sin(phi), math.cos(phi), 0.0],
[math.cos(phi), -math.sin(phi), 0.0],
[0.0, 0.0, 1.0]])
proj.append((num.dot(r, m), in_channels, out_channels_))
return proj
def projection_to_enu(
self,
in_channels,
out_channels=('E', 'N', 'U'),
**kwargs):
return self._projection_to('enu', in_channels, out_channels, **kwargs)
def projection_to_ned(
self,
in_channels,
out_channels=('N', 'E', 'D'),
**kwargs):
return self._projection_to('ned', in_channels, out_channels, **kwargs)
def projection_from_enu(
self,
in_channels=('E', 'N', 'U'),
out_channels=('X', 'Y', 'Z'),
**kwargs):
m, out_channels, in_channels = self._projection_to(
'enu', out_channels, in_channels, **kwargs)
return num.linalg.inv(m), in_channels, out_channels
def projection_from_ned(
self,
in_channels=('N', 'E', 'D'),
out_channels=('X', 'Y', 'Z'),
**kwargs):
m, out_channels, in_channels = self._projection_to(
'ned', out_channels, in_channels, **kwargs)
return num.linalg.inv(m), in_channels, out_channels
def nsl_string(self):
return '.'.join((self.network, self.station, self.location))
def nsl(self):
return self.network, self.station, self.location
def cannot_handle_offsets(self):
if self.north_shift != 0.0 or self.east_shift != 0.0:
logger.warn(
'Station %s.%s.%s has a non-zero Cartesian offset. Such '
'offsets cannot be saved in the basic station file format. '
'Effective lat/lons are saved only. Please save the stations '
'in YAML format to preserve the reference-and-offset '
'coordinates.' % self.nsl())
def oldstr(self):
self.cannot_handle_offsets()
nsl = '%s.%s.%s' % (self.network, self.station, self.location)
s = '%-15s %14.5f %14.5f %14.1f %14.1f %s' % (
nsl, self.effective_lat, self.effective_lon, self.elevation,
self.depth, self.name)
return s
def dump_stations(stations, filename):
'''Write stations file.
:param stations: list of :py:class:`Station` objects
:param filename: filename as str
'''
f = open(filename, 'w')
for sta in stations:
f.write(sta.oldstr()+'\n')
for cha in sta.get_channels():
azimuth = 'NaN'
if cha.azimuth is not None:
azimuth = '%7g' % cha.azimuth
dip = 'NaN'
if cha.dip is not None:
dip = '%7g' % cha.dip
f.write('%5s %14s %14s %14g\n' % (
cha.name, azimuth, dip, cha.gain))
f.close()
def dump_stations_yaml(stations, filename):
'''Write stations file in YAML format.
:param stations: list of :py:class:`Station` objects
:param filename: filename as str
'''
dump_all(stations, filename=filename)
def float_or_none(s):
if s is None:
return None
elif isinstance(s, str) and s.lower() == 'nan':
return None
else:
return float(s)
def detect_format(filename):
with open(filename, 'r') as f:
for line in f:
line = line.strip()
if not line or line.startswith('#') or line.startswith('%'):
continue
if line.startswith('--- !pf.Station'):
return 'yaml'
else:
return 'basic'
return 'basic'
def load_stations(filename, format='detect'):
'''Read stations file.
:param filename: filename
:returns: list of :py:class:`Station` objects
'''
if format == 'detect':
format = detect_format(filename)
if format == 'yaml':
from pyrocko import guts
stations = [
st for st in guts.load_all(filename=filename)
if isinstance(st, Station)]
return stations
elif format == 'basic':
stations = []
f = open(filename, 'r')
station = None
channel_names = []
for (iline, line) in enumerate(f):
toks = line.split(None, 5)
if line.strip().startswith('#') or line.strip() == '':
continue
if len(toks) == 5 or len(toks) == 6:
net, sta, loc = toks[0].split('.')
lat, lon, elevation, depth = [float(x) for x in toks[1:5]]
if len(toks) == 5:
name = ''
else:
name = toks[5].rstrip()
station = Station(
net, sta, loc, lat, lon,
elevation=elevation, depth=depth, name=name)
stations.append(station)
channel_names = []
elif len(toks) == 4 and station is not None:
name, azi, dip, gain = (
toks[0],
float_or_none(toks[1]),
float_or_none(toks[2]),
float(toks[3]))
if name in channel_names:
logger.warning(
'redefined channel! (line: %i, file: %s)' %
(iline + 1, filename))
else:
channel_names.append(name)
channel = Channel(name, azimuth=azi, dip=dip, gain=gain)
station.add_channel(channel)
else:
logger.warning('skipping invalid station/channel definition '
'(line: %i, file: %s' % (iline + 1, filename))
f.close()
return stations
else:
from pyrocko.io.io_common import FileLoadError
raise FileLoadError('unknown event file format: %s' % format)
def dump_kml(objects, filename):
station_template = '''
<Placemark>
<name>%(network)s.%(station)s.%(location)s</name>
<description></description>
<styleUrl>#msn_S</styleUrl>
<Point>
<coordinates>%(elon)f,%(elat)f,%(elevation)f</coordinates>
</Point>
</Placemark>
'''
f = open(filename, 'w')
f.write('<?xml version="1.0" encoding="UTF-8"?>\n')
f.write('<kml xmlns="http://www.opengis.net/kml/2.2">\n')
f.write('<Document>\n')
f.write(''' <Style id="sh_S">
<IconStyle>
<scale>1.3</scale>
<Icon>
<href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
</Icon>
<hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
</IconStyle>
<ListStyle>
<ItemIcon>
<href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
</ItemIcon>
</ListStyle>
</Style>
<Style id="sn_S">
<IconStyle>
<scale>1.1</scale>
<Icon>
<href>http://maps.google.com/mapfiles/kml/paddle/S.png</href>
</Icon>
<hotSpot x="32" y="1" xunits="pixels" yunits="pixels"/>
</IconStyle>
<ListStyle>
<ItemIcon>
<href>http://maps.google.com/mapfiles/kml/paddle/S-lv.png</href>
</ItemIcon>
</ListStyle>
</Style>
<StyleMap id="msn_S">
<Pair>
<key>normal</key>
<styleUrl>#sn_S</styleUrl>
</Pair>
<Pair>
<key>highlight</key>
<styleUrl>#sh_S</styleUrl>
</Pair>
</StyleMap>
''')
for obj in objects:
if isinstance(obj, Station):
d = obj.__dict__.copy()
d['elat'], d['elon'] = obj.effective_latlon
f.write(station_template % d)
f.write('</Document>')
f.write('</kml>\n')
f.close()
|