Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/io/enhanced_sacpz.py: 85%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Reader for enhanced `SAC 

8<http://ds.iris.edu/ds/nodes/dmc/software/downloads/sac/>`_ pole-zero files. 

9 

10Enhanced SAC pole-zero files are SAC pole-zero files with additional channel 

11attributes in the file header comments as they are output e.g. by `rdseed 

12<https://ds.iris.edu/ds/nodes/dmc/software/downloads/rdseed/>`_. 

13''' 

14 

15 

16import time 

17from collections import defaultdict 

18 

19from pyrocko import util, pz, response 

20from pyrocko.io import io_common 

21from pyrocko.io import stationxml as sxml 

22from pyrocko.guts import Object, Tuple, String, Timestamp, Float 

23 

24 

25guts_prefix = 'pf' 

26 

27 

28class EnhancedSacPzError(io_common.FileLoadError): 

29 pass 

30 

31 

32class EnhancedSacPzResponse(Object): 

33 codes = Tuple.T(4, String.T()) 

34 tmin = Timestamp.T(optional=True) 

35 tmax = Timestamp.T(optional=True) 

36 lat = Float.T() 

37 lon = Float.T() 

38 elevation = Float.T() 

39 depth = Float.T() 

40 dip = Float.T() 

41 azimuth = Float.T() 

42 input_unit = String.T() 

43 output_unit = String.T() 

44 response = response.PoleZeroResponse.T() 

45 

46 def spans(self, *args): 

47 if len(args) == 0: 

48 return True 

49 elif len(args) == 1: 

50 return ((self.tmin is None or 

51 self.tmin <= args[0]) and 

52 (self.tmax is None or 

53 args[0] <= self.tmax)) 

54 

55 elif len(args) == 2: 

56 return ((self.tmin is None or 

57 args[1] >= self.tmin) and 

58 (self.tmax is None or 

59 self.tmax >= args[0])) 

60 

61 

62def make_stationxml_response(presponse, input_unit, output_unit): 

63 return sxml.Response.from_pyrocko_pz_response( 

64 presponse, input_unit, output_unit, 1.0) 

65 

66 

67this_year = time.gmtime()[0] 

68 

69 

70def dummy_aware_str_to_time(s, time_format='%Y-%m-%dT%H:%M:%S'): 

71 try: 

72 util.str_to_time(s, format=time_format) 

73 except util.TimeStrError: 

74 year = int(s[:4]) 

75 if year > this_year + 100: 

76 return None # StationXML contained a dummy end date 

77 

78 raise 

79 

80 

81def iload_fh(f, time_format='%Y-%m-%dT%H:%M:%S'): 

82 zeros, poles, constant, comments = pz.read_sac_zpk(file=f, 

83 get_comments=True) 

84 d = {} 

85 for line in comments: 

86 toks = line.split(':', 1) 

87 if len(toks) == 2: 

88 temp = toks[0].strip('* \t') 

89 for k in ('network', 'station', 'location', 'channel', 'start', 

90 'end', 'latitude', 'longitude', 'depth', 'elevation', 

91 'dip', 'azimuth', 'input unit', 'output unit'): 

92 

93 if temp.lower().startswith(k): 

94 d[k] = toks[1].strip() 

95 

96 resp = response.PoleZeroResponse(zeros, poles, constant) 

97 

98 try: 

99 yield EnhancedSacPzResponse( 

100 codes=(d['network'], d['station'], d['location'], d['channel']), 

101 tmin=util.str_to_time(d['start'], format=time_format), 

102 tmax=dummy_aware_str_to_time(d['end']), 

103 lat=float(d['latitude']), 

104 lon=float(d['longitude']), 

105 elevation=float(d['elevation']), 

106 depth=float(d['depth']), 

107 dip=float(d['dip']), 

108 azimuth=float(d['azimuth']), 

109 input_unit=d['input unit'], 

110 output_unit=d['output unit'], 

111 response=resp) 

112 except KeyError as e: 

113 raise EnhancedSacPzError( 

114 'cannot get all required information "%s"' % e.args[0]) 

115 

116 

117iload_filename, iload_dirname, iload_glob, iload = util.make_iload_family( 

118 iload_fh, 'SACPZ', ':py:class:`EnhancedSacPzResponse`') 

119 

120 

121def make_stationxml(responses, inconsistencies='warn'): 

122 ''' 

123 Create stationxml from "enhanced" SACPZ information. 

124 

125 :param responses: iterable yielding 

126 :py:class:`EnhancedSacPzResponse` objects 

127 :returns: :py:class:`pyrocko.io.stationxml.FDSNStationXML` object 

128 ''' 

129 

130 networks = {} 

131 stations = {} 

132 

133 station_coords = defaultdict(list) 

134 station_channels = defaultdict(list) 

135 for cr in responses: 

136 net, sta, loc, cha = cr.codes 

137 station_coords[net, sta].append( 

138 (cr.codes, cr.lat, cr.lon, cr.elevation)) 

139 

140 station_channels[net, sta].append(sxml.Channel( 

141 code=cha, 

142 location_code=loc, 

143 start_date=cr.tmin, 

144 end_date=cr.tmax, 

145 latitude=sxml.Latitude(cr.lat), 

146 longitude=sxml.Longitude(cr.lon), 

147 elevation=sxml.Distance(cr.elevation), 

148 depth=sxml.Distance(cr.depth), 

149 response=make_stationxml_response( 

150 cr.response, cr.input_unit, cr.output_unit))) 

151 

152 for (net, sta), v in sorted(station_coords.items()): 

153 lat, lon, elevation = util.consistency_merge( 

154 v, 

155 message='channel lat/lon/elevation values differ', 

156 error=inconsistencies) 

157 

158 if net not in networks: 

159 networks[net] = sxml.Network(code=net) 

160 

161 stations[net, sta] = sxml.Station( 

162 code=sta, 

163 latitude=sxml.Latitude(lat), 

164 longitude=sxml.Longitude(lon), 

165 elevation=sxml.Distance(elevation), 

166 channel_list=sorted( 

167 station_channels[net, sta], 

168 key=lambda c: (c.location_code, c.code))) 

169 

170 networks[net].station_list.append(stations[net, sta]) 

171 

172 return sxml.FDSNStationXML( 

173 source='Converted from "enhanced" SACPZ information', 

174 created=time.time(), 

175 network_list=[networks[net_] for net_ in sorted(networks.keys())]) 

176 

177 

178if __name__ == '__main__': 

179 

180 import sys 

181 

182 util.setup_logging(__name__) 

183 

184 if len(sys.argv) < 2: 

185 sys.exit('usage: python -m pyrocko.station.enhanced_sacpz <sacpz> ...') 

186 

187 sxml_in = make_stationxml(iload(sys.argv[1:])) 

188 print(sxml_in.dump_xml())