1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division, print_function
7import numpy as num
8from pyrocko import model
9from pyrocko.guts import Int, String, List, Float, Bool
11from pyrocko.model.station import load_stations
12from pyrocko.io import stationxml
13from pyrocko.orthodrome import distance_accurate50m_numpy, \
14 geographic_midpoint_locations
16from .base import LocationGenerator
18guts_prefix = 'pf.scenario'
19km = 1e3
20d2r = num.pi / 180.
23class StationGenerator(LocationGenerator):
25 def __init__(self, *args, **kwargs):
26 super().__init__(*args, **kwargs)
27 self._stations = None
29 def get_stations(self):
30 raise NotImplementedError
32 @property
33 def nstations(self):
34 return len(self.get_stations())
36 def has_stations(self):
37 if not self.get_stations():
38 return False
39 return True
41 def clear(self):
42 super().clear()
43 self._stations = None
45 def get_distance_range(self, sources):
46 dists = []
47 for source in sources:
48 for station in self.get_stations():
49 dists.append(
50 source.distance_to(station))
52 return num.min(dists), num.max(dists)
54 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
55 return []
57 def add_map_artists(self, engine, sources, automap):
58 automap.add_stations(self.get_stations())
61class ImportStationGenerator(StationGenerator):
63 stations_paths = List.T(
64 optional=True,
65 help='List of files with station coordinates in Pyrocko format.')
67 stations_stationxml_paths = List.T(
68 optional=True,
69 help='List of files with station coordinates in StationXML format.')
71 stations = List.T(
72 model.Station.T(),
73 optional=True,
74 help='List of Pyrocko stations.')
76 def get_center_latlon(self):
77 stations = self.get_stations()
78 if not stations:
79 return self._parent.get_center_latlon()
81 return geographic_midpoint_locations(self.get_stations())
83 def get_radius(self):
84 stations = self.get_stations()
85 if not stations:
86 return self._parent.get_radius()
88 clat, clon = self.get_center_latlon()
89 radii = distance_accurate50m_numpy(
90 clat, clon,
91 [st.effective_lat for st in stations],
92 [st.effective_lon for st in stations])
94 return float(radii.max())
96 def get_stations(self):
97 if self._stations is None:
99 stations = []
101 if self.stations_paths:
102 for filename in self.stations_paths:
103 stations.extend(
104 load_stations(filename))
106 if self.stations_stationxml_paths:
107 for filename in self.stations_stationxml_paths:
108 sxml = stationxml.load_xml(filename=filename)
109 stations.extend(
110 sxml.get_pyrocko_stations())
112 if self.stations:
113 stations.extend(self.stations)
115 self._stations = stations
117 return self._stations
119 def nsl(self, istation):
120 stations = self.get_stations()
121 return stations[istation].nsl()
124class CircleStationGenerator(StationGenerator):
126 radius = Float.T(
127 default=100*km,
128 help='Radius of the station circle in km.')
130 azi_start = Float.T(
131 default=0.,
132 help='Start azimuth of circle [deg]. '
133 'Default is a full circle, 0 - 360 deg')
135 azi_end = Float.T(
136 default=360.,
137 help='End azimuth of circle [deg]. '
138 'Default is a full circle, 0 - 360 deg')
140 nstations = Int.T(
141 default=10,
142 help='Number of evenly spaced stations.')
144 network_name = String.T(
145 default='CI',
146 help='Network name.')
148 channels = List.T(
149 default=['BHE', 'BHN', 'BHZ'],
150 help='Seismic channels to generate. Default: BHN, BHE, BHZ.')
152 shift_circle = Bool.T(
153 default=False,
154 help='Rotate circle away by half a station distance.')
156 def get_stations(self):
157 if self._stations is None:
158 if self.channels:
159 channels = [model.station.Channel(c) for c in self.channels]
160 else:
161 channels = None
163 azimuths = num.linspace(
164 self.azi_start*d2r, self.azi_end*d2r,
165 self.nstations, endpoint=False)
167 if self.shift_circle:
168 swath = (self.azi_end - self.azi_start)*d2r
169 azimuths += swath / self.nstations / 2.
171 lat, lon = self.get_center_latlon()
173 stations = []
174 for istation, azi in enumerate(azimuths):
175 net, sta, loc = self.nsl(istation)
177 station = model.Station(
178 net, sta, loc,
179 lat=lat,
180 lon=lon,
181 north_shift=num.cos(azi) * self.radius,
182 east_shift=num.sin(azi) * self.radius,
183 channels=channels)
185 stations.append(station)
187 self._stations = stations
189 return self._stations
191 def nsl(self, istation):
192 return self.network_name, 'S%03i' % (istation + 1), ''
194 def get_radius(self):
195 return self.radius
197 def clear(self):
198 super().clear()
199 self._stations = None
202class RandomStationGenerator(StationGenerator):
204 nstations = Int.T(
205 default=10,
206 help='Number of randomly distributed stations.')
208 network_name = String.T(
209 default='CO',
210 help='Network name')
212 channels = List.T(
213 optional=True,
214 default=['BHE', 'BHN', 'BHZ'],
215 help='Seismic channels to generate. Default: BHN, BHE, BHZ')
217 def nsl(self, istation):
218 return self.network_name, 'S%03i' % (istation + 1), '',
220 def get_stations(self):
221 if self._stations is None:
223 if self.channels:
224 channels = [model.station.Channel(c) for c in self.channels]
225 else:
226 channels = None
228 stations = []
229 for istation in range(self.nstations):
230 lat, lon, north_shift, east_shift, depth = map(
231 float, self.get_coordinates(istation))
233 net, sta, loc = self.nsl(istation)
234 station = model.Station(
235 net, sta, loc,
236 lat=lat,
237 lon=lon,
238 north_shift=north_shift,
239 east_shift=east_shift,
240 depth=depth,
241 channels=channels)
243 stations.append(station)
245 self._stations = stations
247 return self._stations