1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import numpy as num
7from pyrocko import model
8from pyrocko.guts import Int, String, List, Float, Bool
10from pyrocko.model.station import load_stations
11from pyrocko.io import stationxml
12from pyrocko.orthodrome import distance_accurate50m_numpy, \
13 geographic_midpoint_locations
15from .base import LocationGenerator
17guts_prefix = 'pf.scenario'
18km = 1e3
19d2r = num.pi / 180.
22class StationGenerator(LocationGenerator):
24 def __init__(self, *args, **kwargs):
25 super().__init__(*args, **kwargs)
26 self._stations = None
28 def get_stations(self):
29 raise NotImplementedError
31 @property
32 def nstations(self):
33 return len(self.get_stations())
35 def has_stations(self):
36 if not self.get_stations():
37 return False
38 return True
40 def clear(self):
41 super().clear()
42 self._stations = None
44 def get_distance_range(self, sources):
45 dists = []
46 for source in sources:
47 for station in self.get_stations():
48 dists.append(
49 source.distance_to(station))
51 return num.min(dists), num.max(dists)
53 def ensure_data(self, engine, sources, path, tmin=None, tmax=None):
54 return []
56 def add_map_artists(self, engine, sources, automap):
57 automap.add_stations(self.get_stations())
60class ImportStationGenerator(StationGenerator):
62 stations_paths = List.T(
63 optional=True,
64 help='List of files with station coordinates in Pyrocko format.')
66 stations_stationxml_paths = List.T(
67 optional=True,
68 help='List of files with station coordinates in StationXML format.')
70 stations = List.T(
71 model.Station.T(),
72 optional=True,
73 help='List of Pyrocko stations.')
75 def get_center_latlon(self):
76 stations = self.get_stations()
77 if not stations:
78 return self._parent.get_center_latlon()
80 return geographic_midpoint_locations(self.get_stations())
82 def get_radius(self):
83 stations = self.get_stations()
84 if not stations:
85 return self._parent.get_radius()
87 clat, clon = self.get_center_latlon()
88 radii = distance_accurate50m_numpy(
89 clat, clon,
90 [st.effective_lat for st in stations],
91 [st.effective_lon for st in stations])
93 return float(radii.max())
95 def get_stations(self):
96 if self._stations is None:
98 stations = []
100 if self.stations_paths:
101 for filename in self.stations_paths:
102 stations.extend(
103 load_stations(filename))
105 if self.stations_stationxml_paths:
106 for filename in self.stations_stationxml_paths:
107 sxml = stationxml.load_xml(filename=filename)
108 stations.extend(
109 sxml.get_pyrocko_stations())
111 if self.stations:
112 stations.extend(self.stations)
114 self._stations = stations
116 return self._stations
118 def nsl(self, istation):
119 stations = self.get_stations()
120 return stations[istation].nsl()
123class CircleStationGenerator(StationGenerator):
125 radius = Float.T(
126 default=100*km,
127 help='Radius of the station circle in km.')
129 azi_start = Float.T(
130 default=0.,
131 help='Start azimuth of circle [deg]. '
132 'Default is a full circle, 0 - 360 deg')
134 azi_end = Float.T(
135 default=360.,
136 help='End azimuth of circle [deg]. '
137 'Default is a full circle, 0 - 360 deg')
139 nstations = Int.T(
140 default=10,
141 help='Number of evenly spaced stations.')
143 network_name = String.T(
144 default='CI',
145 help='Network name.')
147 channels = List.T(
148 default=['BHE', 'BHN', 'BHZ'],
149 help='Seismic channels to generate. Default: BHN, BHE, BHZ.')
151 shift_circle = Bool.T(
152 default=False,
153 help='Rotate circle away by half a station distance.')
155 def get_stations(self):
156 if self._stations is None:
157 if self.channels:
158 channels = [model.station.Channel(c) for c in self.channels]
159 else:
160 channels = None
162 azimuths = num.linspace(
163 self.azi_start*d2r, self.azi_end*d2r,
164 self.nstations, endpoint=False)
166 if self.shift_circle:
167 swath = (self.azi_end - self.azi_start)*d2r
168 azimuths += swath / self.nstations / 2.
170 lat, lon = self.get_center_latlon()
172 stations = []
173 for istation, azi in enumerate(azimuths):
174 net, sta, loc = self.nsl(istation)
176 station = model.Station(
177 net, sta, loc,
178 lat=lat,
179 lon=lon,
180 north_shift=num.cos(azi) * self.radius,
181 east_shift=num.sin(azi) * self.radius,
182 channels=channels)
184 stations.append(station)
186 self._stations = stations
188 return self._stations
190 def nsl(self, istation):
191 return self.network_name, 'S%03i' % (istation + 1), ''
193 def get_radius(self):
194 return self.radius
196 def clear(self):
197 super().clear()
198 self._stations = None
201class RandomStationGenerator(StationGenerator):
203 nstations = Int.T(
204 default=10,
205 help='Number of randomly distributed stations.')
207 network_name = String.T(
208 default='CO',
209 help='Network name')
211 channels = List.T(
212 optional=True,
213 default=['BHE', 'BHN', 'BHZ'],
214 help='Seismic channels to generate. Default: BHN, BHE, BHZ')
216 def nsl(self, istation):
217 return self.network_name, 'S%03i' % (istation + 1), '',
219 def get_stations(self):
220 if self._stations is None:
222 if self.channels:
223 channels = [model.station.Channel(c) for c in self.channels]
224 else:
225 channels = None
227 stations = []
228 for istation in range(self.nstations):
229 lat, lon, north_shift, east_shift, depth = map(
230 float, self.get_coordinates(istation))
232 net, sta, loc = self.nsl(istation)
233 station = model.Station(
234 net, sta, loc,
235 lat=lat,
236 lon=lon,
237 north_shift=north_shift,
238 east_shift=east_shift,
239 depth=depth,
240 channels=channels)
242 stations.append(station)
244 self._stations = stations
246 return self._stations