1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6from pyrocko import model, response 

7 

8import logging 

9import copy 

10import pickle 

11 

12logger = logging.getLogger('pyrocko.io.eventdata') 

13 

14 

15class NoRestitution(Exception): 

16 pass 

17 

18 

19class FileNotFound(Exception): 

20 

21 def __init__(self, s): 

22 self.s = s 

23 

24 def __str__(self): 

25 return 'File not found: %s' % self.s 

26 

27 

28class Problems(object): 

29 def __init__(self): 

30 self._problems = {} 

31 

32 def add(self, kind, nslct): 

33 if kind not in self._problems: 

34 self._problems[kind] = set() 

35 problem = self._problems[kind] 

36 problem.add(nslct) 

37 

38 def dump(self, fn): 

39 f = open(fn, 'wb') 

40 pickle.dump(self._problems, f) 

41 f.close() 

42 

43 def load(self, fn): 

44 f = open(fn, 'rb') 

45 self._problems = pickle.load(f) 

46 f.close() 

47 

48 def mapped(self, mapping=lambda nslct: nslct[:3]): 

49 p = {} 

50 for kind, problem in self._problems.items(): 

51 nsl = set() 

52 for nslct in problem: 

53 nsl.add(mapping(nslct)) 

54 p[kind] = nsl 

55 

56 return p 

57 

58 

59class EventDataAccess(object): 

60 ''' 

61 Abstract base class for event data access (see rdseed.py). 

62 ''' 

63 

64 def __init__(self, events=None, stations=None, datapile=None): 

65 

66 self._pile = datapile 

67 self._events = events 

68 

69 if stations is None: 

70 self._stations = None 

71 else: 

72 self._stations = {} 

73 for station in stations: 

74 self._stations[station.nsl()] = station 

75 

76 self._problems = Problems() 

77 

78 def get_pile(self): 

79 return self._pile 

80 

81 def get_pyrocko_events(self): 

82 ''' 

83 Extract :py:class:`model.Event` instances from the volume. 

84 ''' 

85 

86 if not self._events: 

87 self._events = self._get_events_from_file() 

88 return self._events 

89 

90 def get_pyrocko_station(self, tr, relative_event=None): 

91 ''' 

92 Get station information for a given trace. 

93 

94 :param tr: :py:class:`trace.Trace` instance 

95 

96 :returns: :py:class:`model.station.Station` objects. 

97 ''' 

98 

99 self._update_stations() 

100 s = copy.deepcopy(self._stations[tr.nslc_id[:3]]) 

101 if relative_event is not None: 

102 s.set_event_relative_data(relative_event) 

103 return s 

104 

105 def get_pyrocko_channel(self, tr): 

106 ''' 

107 Get channel information for a given trace. 

108 

109 :param tr: :py:class:`trace.Trace` instance 

110 

111 :returns: :py:class:`model.station.Channel` objects. 

112 ''' 

113 sta = self.get_station(tr) 

114 return sta.get_channel(tr.channel) 

115 

116 def get_pyrocko_stations(self): 

117 ''' 

118 Exctract a list of :py:class:`model.Station` instances. 

119 ''' 

120 return list(self._get_stations().values()) 

121 

122 def _get_stations(self, relative_event=None): 

123 self._update_stations() 

124 stations = copy.deepcopy(self._stations) 

125 if relative_event is not None: 

126 for s in stations.values(): 

127 s.set_event_relative_data(relative_event) 

128 

129 return stations 

130 

131 def _update_stations(self): 

132 if not self._stations: 

133 self._stations = {} 

134 for station in self._get_stations_from_file(): 

135 self._stations[station.nsl()] = station 

136 self._insert_channel_descriptions(self._stations) 

137 

138 def _insert_channel_descriptions(self, stations): 

139 pile = self.get_pile() 

140 nslc_ids = pile.gather_keys( 

141 lambda tr: (tr.network, tr.station, tr.location, tr.channel)) 

142 

143 for nslc in nslc_ids: 

144 if nslc[:3] not in stations: 

145 logger.warning( 

146 'No station description for trace %s.%s.%s.%s' % nslc) 

147 continue 

148 

149 sta = stations[nslc[:3]] 

150 try: 

151 sta.add_channel(self._get_channel_description_from_file(nslc)) 

152 except FileNotFound: 

153 logger.warning( 

154 'No channel description for trace %s.%s.%s.%s' % nslc) 

155 

156 def _get_channel_description_from_file(self, nslc): 

157 return model.Channel(nslc[3], None, None, 1.) 

158 

159 def iter_traces(self, group_selector=None, trace_selector=None): 

160 

161 for traces in self.get_pile().chopper_grouped( 

162 gather=lambda tr: (tr.network, tr.station, tr.location), 

163 group_selector=group_selector, 

164 trace_selector=trace_selector): 

165 

166 yield traces 

167 

168 def problems(self): 

169 return self._problems 

170 

171 def _redundant_channel_weeder(self, redundant_channel_priorities, nslcs): 

172 

173 if redundant_channel_priorities is None: 

174 return [] 

175 

176 # n=network,s=station,l=location,c=channel 

177 # channels by station 

178 by_nsl = {} 

179 for nslc in nslcs: 

180 nsl = nslc[:3] 

181 if nsl not in by_nsl: 

182 by_nsl[nsl] = [] 

183 

184 by_nsl[nsl].append(nslc) 

185 

186 # figure out what items to remove 

187 to_delete = [] 

188 for ((h1, h2), (l1, l2)) in redundant_channel_priorities: 

189 for nsl, nslcs in by_nsl.items(): 

190 channels = [nslc[3] for nslc in nslcs] 

191 if h1 in channels and \ 

192 h2 in channels and \ 

193 l1 in channels and \ 

194 l2 in channels: 

195 

196 to_delete.append(nslc[:3] + (l1,)) 

197 to_delete.append(nslc[:3] + (l2,)) 

198 

199 return to_delete 

200 

201 def get_restitution(self, tr, allowed_methods): 

202 if 'integration' in allowed_methods: 

203 response.IntegrationResponse() 

204 else: 

205 raise Exception('only "integration" restitution method is allowed')