Coverage for /usr/local/lib/python3.7/dist-packages/pyrocko/model/station.py : 76%

# 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()