Source code for pyrocko.scenario.station

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
from __future__ import absolute_import, division, print_function

import numpy as num
from pyrocko import model
from pyrocko.guts import Int, String, List, Float, Bool

from pyrocko.model.station import load_stations
from pyrocko.io import stationxml
from pyrocko.orthodrome import distance_accurate50m_numpy, \
    geographic_midpoint_locations

from .base import LocationGenerator

guts_prefix = 'pf.scenario'
km = 1e3
d2r = num.pi / 180.


[docs]class StationGenerator(LocationGenerator): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) self._stations = None def get_stations(self): raise NotImplementedError @property def nstations(self): return len(self.get_stations()) def has_stations(self): if not self.get_stations(): return False return True def clear(self): super().clear() self._stations = None def get_distance_range(self, sources): dists = [] for source in sources: for station in self.get_stations(): dists.append( source.distance_to(station)) return num.min(dists), num.max(dists) def ensure_data(self, engine, sources, path, tmin=None, tmax=None): return [] def add_map_artists(self, engine, sources, automap): automap.add_stations(self.get_stations())
[docs]class ImportStationGenerator(StationGenerator): stations_paths = List.T( optional=True, help='List of files with station coordinates in Pyrocko format.') stations_stationxml_paths = List.T( optional=True, help='List of files with station coordinates in StationXML format.') stations = List.T( model.Station.T(), optional=True, help='List of Pyrocko stations.') def get_center_latlon(self): stations = self.get_stations() if not stations: return self._parent.get_center_latlon() return geographic_midpoint_locations(self.get_stations()) def get_radius(self): stations = self.get_stations() if not stations: return self._parent.get_radius() clat, clon = self.get_center_latlon() radii = distance_accurate50m_numpy( clat, clon, [st.effective_lat for st in stations], [st.effective_lon for st in stations]) return float(radii.max()) def get_stations(self): if self._stations is None: stations = [] if self.stations_paths: for filename in self.stations_paths: stations.extend( load_stations(filename)) if self.stations_stationxml_paths: for filename in self.stations_stationxml_paths: sxml = stationxml.load_xml(filename=filename) stations.extend( sxml.cko_stations()) if self.stations: stations.extend(self.stations) self._stations = stations return self._stations def nsl(self, istation): stations = self.get_stations() return stations[istation].nsl()
[docs]class CircleStationGenerator(StationGenerator): radius = Float.T( default=100*km, help='Radius of the station circle in km.') azi_start = Float.T( default=0., help='Start azimuth of circle [deg]. ' 'Default is a full circle, 0 - 360 deg') azi_end = Float.T( default=360., help='End azimuth of circle [deg]. ' 'Default is a full circle, 0 - 360 deg') nstations = Int.T( default=10, help='Number of evenly spaced stations.') network_name = String.T( default='CI', help='Network name.') channels = List.T( default=['BHE', 'BHN', 'BHZ'], help='Seismic channels to generate. Default: BHN, BHE, BHZ.') shift_circle = Bool.T( default=False, help='Rotate circle away by half a station distance.') def get_stations(self): if self._stations is None: if self.channels: channels = [model.station.Channel(c) for c in self.channels] else: channels = None azimuths = num.linspace( self.azi_start*d2r, self.azi_end*d2r, self.nstations, endpoint=False) if self.shift_circle: swath = (self.azi_end - self.azi_start)*d2r azimuths += swath / self.nstations / 2. lat, lon = self.get_center_latlon() stations = [] for istation, azi in enumerate(azimuths): net, sta, loc = self.nsl(istation) station = model.Station( net, sta, loc, lat=lat, lon=lon, north_shift=num.cos(azi) * self.radius, east_shift=num.sin(azi) * self.radius, channels=channels) stations.append(station) self._stations = stations return self._stations def nsl(self, istation): return self.network_name, 'S%03i' % (istation + 1), '' def get_radius(self): return self.radius def clear(self): super().clear() self._stations = None
[docs]class RandomStationGenerator(StationGenerator): nstations = Int.T( default=10, help='Number of randomly distributed stations.') network_name = String.T( default='CO', help='Network name') channels = List.T( optional=True, default=['BHE', 'BHN', 'BHZ'], help='Seismic channels to generate. Default: BHN, BHE, BHZ') def nsl(self, istation): return self.network_name, 'S%03i' % (istation + 1), '', def get_stations(self): if self._stations is None: if self.channels: channels = [model.station.Channel(c) for c in self.channels] else: channels = None stations = [] for istation in range(self.nstations): lat, lon, north_shift, east_shift, depth = map( float, self.get_coordinates(istation)) net, sta, loc = self.nsl(istation) station = model.Station( net, sta, loc, lat=lat, lon=lon, north_shift=north_shift, east_shift=east_shift, depth=depth, channels=channels) stations.append(station) self._stations = stations return self._stations