Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/scenario/station.py: 56%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6''' 

7Implementation of the station generators. 

8''' 

9 

10import numpy as num 

11from pyrocko import model 

12from pyrocko.guts import Int, String, List, Float, Bool 

13 

14from pyrocko.model.station import load_stations 

15from pyrocko.io import stationxml 

16from pyrocko.orthodrome import distance_accurate50m_numpy, \ 

17 geographic_midpoint_locations 

18 

19from .base import LocationGenerator 

20 

21guts_prefix = 'pf.scenario' 

22km = 1e3 

23d2r = num.pi / 180. 

24 

25 

26class StationGenerator(LocationGenerator): 

27 

28 def __init__(self, *args, **kwargs): 

29 super().__init__(*args, **kwargs) 

30 self._stations = None 

31 

32 def get_stations(self): 

33 raise NotImplementedError 

34 

35 @property 

36 def nstations(self): 

37 return len(self.get_stations()) 

38 

39 def has_stations(self): 

40 if not self.get_stations(): 

41 return False 

42 return True 

43 

44 def clear(self): 

45 super().clear() 

46 self._stations = None 

47 

48 def get_distance_range(self, sources): 

49 dists = [] 

50 for source in sources: 

51 for station in self.get_stations(): 

52 dists.append( 

53 source.distance_to(station)) 

54 

55 return num.min(dists), num.max(dists) 

56 

57 def ensure_data(self, engine, sources, path, tmin=None, tmax=None): 

58 return [] 

59 

60 def add_map_artists(self, engine, sources, automap): 

61 automap.add_stations(self.get_stations()) 

62 

63 

64class ImportStationGenerator(StationGenerator): 

65 

66 stations_paths = List.T( 

67 optional=True, 

68 help='List of files with station coordinates in Pyrocko format.') 

69 

70 stations_stationxml_paths = List.T( 

71 optional=True, 

72 help='List of files with station coordinates in StationXML format.') 

73 

74 stations = List.T( 

75 model.Station.T(), 

76 optional=True, 

77 help='List of Pyrocko stations.') 

78 

79 def get_center_latlon(self): 

80 stations = self.get_stations() 

81 if not stations: 

82 return self._parent.get_center_latlon() 

83 

84 return geographic_midpoint_locations(self.get_stations()) 

85 

86 def get_radius(self): 

87 stations = self.get_stations() 

88 if not stations: 

89 return self._parent.get_radius() 

90 

91 clat, clon = self.get_center_latlon() 

92 radii = distance_accurate50m_numpy( 

93 clat, clon, 

94 [st.effective_lat for st in stations], 

95 [st.effective_lon for st in stations]) 

96 

97 return float(radii.max()) 

98 

99 def get_stations(self): 

100 if self._stations is None: 

101 

102 stations = [] 

103 

104 if self.stations_paths: 

105 for filename in self.stations_paths: 

106 stations.extend( 

107 load_stations(filename)) 

108 

109 if self.stations_stationxml_paths: 

110 for filename in self.stations_stationxml_paths: 

111 sxml = stationxml.load_xml(filename=filename) 

112 stations.extend( 

113 sxml.get_pyrocko_stations()) 

114 

115 if self.stations: 

116 stations.extend(self.stations) 

117 

118 self._stations = stations 

119 

120 return self._stations 

121 

122 def nsl(self, istation): 

123 stations = self.get_stations() 

124 return stations[istation].nsl() 

125 

126 

127class CircleStationGenerator(StationGenerator): 

128 

129 radius = Float.T( 

130 default=100*km, 

131 help='Radius of the station circle in km.') 

132 

133 azi_start = Float.T( 

134 default=0., 

135 help='Start azimuth of circle [deg]. ' 

136 'Default is a full circle, 0 - 360 deg') 

137 

138 azi_end = Float.T( 

139 default=360., 

140 help='End azimuth of circle [deg]. ' 

141 'Default is a full circle, 0 - 360 deg') 

142 

143 nstations = Int.T( 

144 default=10, 

145 help='Number of evenly spaced stations.') 

146 

147 network_name = String.T( 

148 default='CI', 

149 help='Network name.') 

150 

151 channels = List.T( 

152 default=['BHE', 'BHN', 'BHZ'], 

153 help='Seismic channels to generate. Default: BHN, BHE, BHZ.') 

154 

155 shift_circle = Bool.T( 

156 default=False, 

157 help='Rotate circle away by half a station distance.') 

158 

159 def get_stations(self): 

160 if self._stations is None: 

161 if self.channels: 

162 channels = [model.station.Channel(c) for c in self.channels] 

163 else: 

164 channels = None 

165 

166 azimuths = num.linspace( 

167 self.azi_start*d2r, self.azi_end*d2r, 

168 self.nstations, endpoint=False) 

169 

170 if self.shift_circle: 

171 swath = (self.azi_end - self.azi_start)*d2r 

172 azimuths += swath / self.nstations / 2. 

173 

174 lat, lon = self.get_center_latlon() 

175 

176 stations = [] 

177 for istation, azi in enumerate(azimuths): 

178 net, sta, loc = self.nsl(istation) 

179 

180 station = model.Station( 

181 net, sta, loc, 

182 lat=lat, 

183 lon=lon, 

184 north_shift=num.cos(azi) * self.radius, 

185 east_shift=num.sin(azi) * self.radius, 

186 channels=channels) 

187 

188 stations.append(station) 

189 

190 self._stations = stations 

191 

192 return self._stations 

193 

194 def nsl(self, istation): 

195 return self.network_name, 'S%03i' % (istation + 1), '' 

196 

197 def get_radius(self): 

198 return self.radius 

199 

200 def clear(self): 

201 super().clear() 

202 self._stations = None 

203 

204 

205class RandomStationGenerator(StationGenerator): 

206 

207 nstations = Int.T( 

208 default=10, 

209 help='Number of randomly distributed stations.') 

210 

211 network_name = String.T( 

212 default='CO', 

213 help='Network name') 

214 

215 channels = List.T( 

216 optional=True, 

217 default=['BHE', 'BHN', 'BHZ'], 

218 help='Seismic channels to generate. Default: BHN, BHE, BHZ') 

219 

220 def nsl(self, istation): 

221 return self.network_name, 'S%03i' % (istation + 1), '', 

222 

223 def get_stations(self): 

224 if self._stations is None: 

225 

226 if self.channels: 

227 channels = [model.station.Channel(c) for c in self.channels] 

228 else: 

229 channels = None 

230 

231 stations = [] 

232 for istation in range(self.nstations): 

233 lat, lon, north_shift, east_shift, depth = map( 

234 float, self.get_coordinates(istation)) 

235 

236 net, sta, loc = self.nsl(istation) 

237 station = model.Station( 

238 net, sta, loc, 

239 lat=lat, 

240 lon=lon, 

241 north_shift=north_shift, 

242 east_shift=east_shift, 

243 depth=depth, 

244 channels=channels) 

245 

246 stations.append(station) 

247 

248 self._stations = stations 

249 

250 return self._stations