1 | |
---|
2 | # file EventDsc.py |
---|
3 | # =========== |
---|
4 | # |
---|
5 | # version 1, 14-Nov-2007 |
---|
6 | # |
---|
7 | # Python interface to SHM evt files. |
---|
8 | # K. Stammler, 14-Nov-2007 |
---|
9 | |
---|
10 | |
---|
11 | import sys |
---|
12 | import os |
---|
13 | import datetime |
---|
14 | |
---|
15 | # ------------------------------------------------------------------------------ |
---|
16 | |
---|
17 | |
---|
18 | class Event: |
---|
19 | |
---|
20 | "Event type class." |
---|
21 | |
---|
22 | def __init__( self, evtname=None ): |
---|
23 | self.phase = [] |
---|
24 | self.epi = {} |
---|
25 | if evtname: |
---|
26 | self.LoadEvt( evtname ) |
---|
27 | |
---|
28 | def LoadEvt( self, evtname ): |
---|
29 | "Reads in evt file." |
---|
30 | currphase = {} |
---|
31 | for line in file(evtname): |
---|
32 | for k in typedict.keys(): |
---|
33 | if line.startswith(k): |
---|
34 | val = line.lstrip(k).strip().lstrip(':').lstrip() |
---|
35 | if typedict[k][1] == 'i': |
---|
36 | val = int(val) |
---|
37 | elif typedict[k][1] == 'f': |
---|
38 | val = float(val) |
---|
39 | if typedict[k][0] == 'e': |
---|
40 | self.epi[k] = val |
---|
41 | elif typedict[k][0] == 'p': |
---|
42 | currphase[k] = val |
---|
43 | break |
---|
44 | else: |
---|
45 | if line.startswith(EvITEM_END_OF_EVENT): |
---|
46 | self.phase.append( currphase ) |
---|
47 | currphase = {} |
---|
48 | |
---|
49 | if len(currphase.keys()) > 0: |
---|
50 | print 'lost phase info:', currphase |
---|
51 | |
---|
52 | |
---|
53 | # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
---|
54 | |
---|
55 | |
---|
56 | def WriteLocsatFiles( self, outprefix, tabprefix='tab', autostartloc=True,\ |
---|
57 | fixed_depth=True, conflevel=0.9, damp=-1, stderr=1.0, degfree=8,\ |
---|
58 | maxiter=20, initdepth=10.0 ): |
---|
59 | |
---|
60 | "Writes LocSAT input files." |
---|
61 | |
---|
62 | if len(self.phase) < 4: |
---|
63 | print 'too small number of phases' |
---|
64 | sys.exit() |
---|
65 | |
---|
66 | # write control file |
---|
67 | fp = open( outprefix+"_c.dat", "w" ) |
---|
68 | tpath = os.environ['SH_LOCSAT']+'/tables' |
---|
69 | if len(tpath) > 30: |
---|
70 | print 'SH_LOCSAT path too long.' |
---|
71 | sys.exit() |
---|
72 | fp.write( "%s %s\n" % (tpath,tabprefix) ) |
---|
73 | fp.write( "%s/scorr TT\n" % os.environ['SH_LOCSAT'] ) |
---|
74 | ac = bc = cc = 'y' |
---|
75 | if not autostartloc: ac = 'n' |
---|
76 | if not fixed_depth: bc = 'n' |
---|
77 | fp.write( " %s %s %s %5.2f %6.2f %5.2f %5d %3d\n" \ |
---|
78 | % (ac,bc,cc,conflevel,damp,stderr, degfree, maxiter) ) |
---|
79 | fp.close() |
---|
80 | |
---|
81 | # first arrival time |
---|
82 | olist = [Sh2Datetime(p[EvITEM_ONSET_TIME]) for p in self.phase] |
---|
83 | olist.sort() |
---|
84 | minonset = olist[0] |
---|
85 | |
---|
86 | # station with first onset |
---|
87 | firststat = None |
---|
88 | for p in self.phase: |
---|
89 | if Sh2Datetime(p[EvITEM_ONSET_TIME]) == minonset: |
---|
90 | firststat = p[EvITEM_STATION] |
---|
91 | break |
---|
92 | |
---|
93 | # exclude phases with quality smaller than 2 |
---|
94 | numphases = 0 |
---|
95 | for p in self.phase: |
---|
96 | if EvITEM_QUALITY in p.keys() and p[EvITEM_QUALITY] > 1: |
---|
97 | numphases += 1 |
---|
98 | |
---|
99 | # inital location |
---|
100 | initlat = initlon = initorig = None |
---|
101 | if EvITEM_LATITUDE in self.epi.keys() and self.epi[EvITEM_LATITUDE] != None: |
---|
102 | initlat = self.epi[EvITEM_LATITUDE] |
---|
103 | if EvITEM_LONGITUDE in self.epi.keys() and self.epi[EvITEM_LONGITUDE] != None: |
---|
104 | initlon = self.epi[EvITEM_LONGITUDE] |
---|
105 | if EvITEM_ORIGIN_TIME in self.epi.keys() and self.epi[EvITEM_ORIGIN_TIME] != None: |
---|
106 | initorig = Sh2Datetime( self.epi[EvITEM_ORIGIN_TIME] ) |
---|
107 | if initlat == None or initlon == None: |
---|
108 | initlat, initlon = StationCoords( firststat ) |
---|
109 | if initlat == None or initlon == None: |
---|
110 | print "station info not found for", firststat, " Abort." |
---|
111 | return |
---|
112 | if initorig == None: |
---|
113 | initorig = minonset |
---|
114 | |
---|
115 | # write station info file |
---|
116 | fp = open( outprefix+"_s.dat", "w" ) |
---|
117 | for p in self.phase: |
---|
118 | code = p[EvITEM_STATION] |
---|
119 | (slat, slon, selev) = StationCoords( code ) |
---|
120 | if slat == None: continue |
---|
121 | if selev == None: |
---|
122 | selev = 0.0 |
---|
123 | else: |
---|
124 | selev /= 1000.0 |
---|
125 | fp.write( "%6s%9.4lf%9.4lf%9.4f\n" % (code,slat,slon,selev) ) |
---|
126 | fp.close() |
---|
127 | |
---|
128 | # write phase data file |
---|
129 | year = initorig.year |
---|
130 | if year > 1900: year -= 1900 |
---|
131 | if year >= 100: year -= 100 |
---|
132 | arrival_id = 0 |
---|
133 | fp = open( outprefix+"_d.dat", "w" ) |
---|
134 | fp.write( "%02d%02d%02d %02d%02d%3d.%02d" % (year, \ |
---|
135 | initorig.month, initorig.day, initorig.hour, initorig.minute, \ |
---|
136 | initorig.second, initorig.microsecond/10000) ) |
---|
137 | fp.write( " %f %f %f 0.0 %d\n" % (initlat, initlon, initdepth, \ |
---|
138 | numphases) ) |
---|
139 | for p in self.phase: |
---|
140 | if EvITEM_QUALITY in p.keys() and p[EvITEM_QUALITY] > 1: |
---|
141 | arrival_id += 1 |
---|
142 | td = Sh2Datetime(p[EvITEM_ONSET_TIME]) - minonset |
---|
143 | tdiff = float(td.seconds) + float(td.microseconds)/1000000.0 |
---|
144 | fp.write( "%8d %6s %8s %4s%2s %f %f\n" % (arrival_id, \ |
---|
145 | p[EvITEM_STATION], p[EvITEM_PHASE], "t", "d", tdiff, 1.0 ) ) |
---|
146 | fp.close() |
---|
147 | |
---|
148 | self.minonset = minonset |
---|
149 | |
---|
150 | |
---|
151 | |
---|
152 | # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
---|
153 | |
---|
154 | |
---|
155 | def LocsatResult( self, prefix ): |
---|
156 | |
---|
157 | "Calls LocSAT and reads the result." |
---|
158 | |
---|
159 | os.system( \ |
---|
160 | "$SH_UTIL/LocSAT -s %s_s.dat -d %s_d.dat -c %s_c.dat -o %s_r.dat"\ |
---|
161 | % (prefix,prefix,prefix,prefix) ) |
---|
162 | |
---|
163 | fp = open( prefix+"_r.dat", "r" ) |
---|
164 | valid = False |
---|
165 | for line in fp.readlines(): |
---|
166 | if line.endswith("Converged!\n"): |
---|
167 | valid = True |
---|
168 | if not valid: continue |
---|
169 | if line.startswith(" Latitude:"): |
---|
170 | x = line.split() |
---|
171 | self.epi[EvITEM_LATITUDE] = float(x[1]) |
---|
172 | sign = x[3] |
---|
173 | if sign == "S": self.epi[EvITEM_LATITUDE] *= -1.0 |
---|
174 | self.epi[EvITEM_ERR_LAT] = float(x[5]) |
---|
175 | elif line.startswith(" Longitude:"): |
---|
176 | x = line.split() |
---|
177 | self.epi[EvITEM_LONGITUDE] = float(x[1]) |
---|
178 | sign = x[3] |
---|
179 | if sign == "W": self.epi[EvITEM_LONGITUDE] *= -1.0 |
---|
180 | self.epi[EvITEM_ERR_LON] = float(x[5]) |
---|
181 | elif line.startswith(" Depth:"): |
---|
182 | x = line.strip().lstrip("Depth:").strip().split() |
---|
183 | self.epi[EvITEM_DEPTH] = float(x[0]) |
---|
184 | self.epi[EvITEM_ERR_DEPTH] = float(x[3]) |
---|
185 | elif line.startswith(" Relative O.T.:"): |
---|
186 | x = line.strip().lstrip("Relative O.T.:").strip().split() |
---|
187 | relorig = float(x[0]) |
---|
188 | tmp = self.minonset + datetime.timedelta(seconds=relorig) |
---|
189 | millisec = tmp.microsecond / 1000 |
---|
190 | msstr = "%03d" % millisec |
---|
191 | self.epi[EvITEM_ORIGIN_TIME] = tmp.strftime("%d-%h-%Y_%T")+'.'+msstr |
---|
192 | self.epi[EvITEM_ERR_ORIGIN] = float(x[3]) |
---|
193 | fp.close() |
---|
194 | |
---|
195 | |
---|
196 | |
---|
197 | #------------------------------------------------------------------------------- |
---|
198 | |
---|
199 | |
---|
200 | |
---|
201 | def StationCoords( code ): |
---|
202 | "Returns station coordinates or (None,None)." |
---|
203 | res = os.popen( "$SH_UTIL/statinf "+code ).readline().split() |
---|
204 | if len(res) < 3: |
---|
205 | return (None,None) |
---|
206 | for i in range(len(res)): |
---|
207 | if res[i].startswith("elev:"): |
---|
208 | elev = res[i].lstrip("elev:") |
---|
209 | if elev == "": |
---|
210 | if len(res)>i+1: |
---|
211 | elev = float(res[i+1]) |
---|
212 | else: |
---|
213 | elev = None |
---|
214 | else: |
---|
215 | elev = float( elev ) |
---|
216 | break |
---|
217 | else: |
---|
218 | elev = None |
---|
219 | return (float(res[1]),float(res[2]),elev) |
---|
220 | |
---|
221 | #------------------------------------------------------------------------------- |
---|
222 | |
---|
223 | |
---|
224 | |
---|
225 | def Sh2Datetime( shtime ): |
---|
226 | """Convert SH time string to datetime. |
---|
227 | |
---|
228 | shtime: SH time string |
---|
229 | """ |
---|
230 | |
---|
231 | x = os.popen( "$SH_UTIL/timename time_to_int "+shtime ).readlines() |
---|
232 | if len(x) < 1: |
---|
233 | print 'Sh2Datetime failed on', '>', shtime, '<' |
---|
234 | return None |
---|
235 | y = x[0].rstrip() |
---|
236 | res = y.split( " " ) |
---|
237 | if (len(res) < 7): |
---|
238 | print "Syntax error in shtime:", shtime |
---|
239 | sys.exit() |
---|
240 | year, month, day, hour, minute, sec, ms = res |
---|
241 | s_time = datetime.datetime( int(year), int(month), int(day), int(hour), \ |
---|
242 | int(minute), int(sec), int(ms)*1000 ) |
---|
243 | return s_time |
---|
244 | |
---|
245 | |
---|
246 | |
---|
247 | # ------------------------------------------------------------------------------ |
---|
248 | |
---|
249 | # string definitions of item descriptors |
---|
250 | EvITEM_ARRAY = "Array name" |
---|
251 | EvITEM_STATION = "Station code" |
---|
252 | EvITEM_ONSET_TIME = "Onset time" |
---|
253 | EvITEM_ONSET_TYPE = "Onset type" |
---|
254 | EvITEM_ONSET_COUNT = "Following Onsets" |
---|
255 | EvITEM_PHASE = "Phase name" |
---|
256 | EvITEM_SIGN = "Sign" |
---|
257 | EvITEM_COMPONENT = "Component" |
---|
258 | EvITEM_PERIOD = "Period (sec)" |
---|
259 | EvITEM_AMPLITUDE = "Amplitude (nm)" |
---|
260 | EvITEM_AMPLITUDE_TIME = "Amplitude Time (sec)" |
---|
261 | EvITEM_AMPLITUDE_VEL = "Vel. Amplitude (nm/sec)" |
---|
262 | EvITEM_LP_COMPONENT = "LP component" |
---|
263 | EvITEM_LP_PERIOD = "LP period (sec)" |
---|
264 | EvITEM_LP_AMPLITUDE = "LP amplitude" |
---|
265 | EvITEM_BB_AMPLITUDE = "BB amplitude (nm/sec)" |
---|
266 | EvITEM_BB_PERIOD = "BB period (sec)" |
---|
267 | EvITEM_B_SLOWNESS = "Beam-Slowness (sec/deg)" |
---|
268 | EvITEM_B_AZIMUTH = "Beam-Azimuth (deg)" |
---|
269 | EvITEM_L_SLOWNESS = "Epi-Slowness (sec/deg)" |
---|
270 | EvITEM_L_AZIMUTH = "Epi-Azimuth (deg)" |
---|
271 | EvITEM_THEO_AZIM = "Theo. Azimuth (deg)" |
---|
272 | EvITEM_THEO_BACK_AZIM = "Theo. Backazimuth (deg)" |
---|
273 | EvITEM_DISTANCE_DEG = "Distance (deg)" |
---|
274 | EvITEM_DISTANCE_KM = "Distance (km)" |
---|
275 | EvITEM_QUALITY = "Quality number" |
---|
276 | EvITEM_MS = "Magnitude ms" |
---|
277 | EvITEM_MB = "Magnitude mb" |
---|
278 | EvITEM_ML = "Magnitude ml" |
---|
279 | EvITEM_MW = "Magnitude mw" |
---|
280 | EvITEM_MBB = "Broadband Magnitude" |
---|
281 | EvITEM_MU = "User Magnitude" |
---|
282 | EvITEM_MEAN_MS = "Mean Magnitude ms" |
---|
283 | EvITEM_MEAN_MB = "Mean Magnitude mb" |
---|
284 | EvITEM_MEAN_ML = "Mean Magnitude ml" |
---|
285 | EvITEM_MEAN_MW = "Mean Magnitude mw" |
---|
286 | EvITEM_MEAN_MBB = "Mean BB Magnitude" |
---|
287 | EvITEM_MEAN_MU = "Mean User Magnitude" |
---|
288 | EvITEM_MU_DESCR = "User Magn. Description" |
---|
289 | EvITEM_REGION = "Source region" |
---|
290 | EvITEM_COMMENT = "Comment" |
---|
291 | EvITEM_LATITUDE = "Latitude" |
---|
292 | EvITEM_LONGITUDE = "Longitude" |
---|
293 | EvITEM_LOC_METHOD = "Location method" |
---|
294 | EvITEM_LOC_QUALITY = "Location quality" |
---|
295 | EvITEM_LOC_VELMOD = "Velocity Model" |
---|
296 | EvITEM_LOC_ADDPAR = "Location Input Params" |
---|
297 | EvITEM_FILTER = "Applied filter" |
---|
298 | EvITEM_DEPTH = "Depth (km)" |
---|
299 | EvITEM_DEPTH_TYPE = "Depth type" |
---|
300 | EvITEM_ORIGIN_TIME = "Origin time" |
---|
301 | EvITEM_EVID = "Event ID" |
---|
302 | EvITEM_WEIGHT = "Weight" |
---|
303 | EvITEM_ONSET_ACC = "Onset Accuracy" |
---|
304 | EvITEM_REF_LATITUDE = "Reference Latitude" |
---|
305 | EvITEM_REF_LONGITUDE = "Reference Longitude" |
---|
306 | EvITEM_REF_NAME = "Reference Location Name" |
---|
307 | EvITEM_ANALYST = "Analyst" |
---|
308 | EvITEM_STATIONS_USED = "No. of Stations used" |
---|
309 | EvITEM_REGION_TABLE = "Region Table" |
---|
310 | EvITEM_REGION_ID = "Region ID" |
---|
311 | EvITEM_EVENT_TYPE = "Event Type" |
---|
312 | EvITEM_SOURCE = "Source of Information" |
---|
313 | EvITEM_AP_SOURCE = "Ampl&Period Source" |
---|
314 | EvITEM_ONSET_WDW_L = "Onset Window Left" |
---|
315 | EvITEM_ONSET_WDW_R = "Onset Window Right" |
---|
316 | EvITEM_ERR_LAT = "Error in Latitude (km)" |
---|
317 | EvITEM_ERR_LON = "Error in Longitude (km)" |
---|
318 | EvITEM_ERR_DEPTH = "Error in Depth (km)" |
---|
319 | EvITEM_ERR_ORIGIN = "Error in Origin Time" |
---|
320 | EvITEM_ERR_SMAJOR = "Error Ellipse Major" |
---|
321 | EvITEM_ERR_SMINOR = "Error Ellipse Minor" |
---|
322 | EvITEM_ERR_MAJSTRIKE = "Error Ellipse Strike" |
---|
323 | EvITEM_AZIM_MAX_GAP = "Max Azimuthal Gap (deg)" |
---|
324 | EvITEM_RESID_RMS = "RMS of Residuals (sec)" |
---|
325 | EvITEM_PHASE_FLAGS = "Phase Flags" |
---|
326 | EvITEM_SIGNOISE = "Signal/Noise" |
---|
327 | EvITEM_RESIDUAL = "Residual Time" |
---|
328 | EvITEM_RESID_CORR = "Residual Correction" |
---|
329 | EvITEM_PICK_TYPE = "Pick Type" |
---|
330 | EvITEM_CORNERFREQ = "Corner Frequency" |
---|
331 | EvITEM_LOWFREQLEVEL = "Low Frequency Level" |
---|
332 | EvITEM_MOMTEN = "Moment Tensor Elements" |
---|
333 | EvITEM_MOMTEN_DESCR = "Moment Tensor Descr." |
---|
334 | EvITEM_M0 = "Scalar Moment" |
---|
335 | EvITEM_FPS_ANGLES = "Fault Plane Solution" |
---|
336 | EvITEM_FPS_DESCR = "FPS Description" |
---|
337 | EvITEM_END_OF_EVENT = "--- End of Phase ---" |
---|
338 | |
---|
339 | typedict = { |
---|
340 | EvITEM_ARRAY : "ps", |
---|
341 | EvITEM_STATION : "ps", |
---|
342 | EvITEM_ONSET_TIME : "ps", |
---|
343 | EvITEM_ONSET_TYPE : "ps", |
---|
344 | EvITEM_ONSET_COUNT : "pi", |
---|
345 | EvITEM_PHASE : "ps", |
---|
346 | EvITEM_SIGN : "ps", |
---|
347 | EvITEM_COMPONENT : "ps", |
---|
348 | EvITEM_PERIOD : "pf", |
---|
349 | EvITEM_AMPLITUDE : "pf", |
---|
350 | EvITEM_AMPLITUDE_TIME : "pf", |
---|
351 | EvITEM_AMPLITUDE_VEL : "pf", |
---|
352 | EvITEM_LP_COMPONENT : "ps", |
---|
353 | EvITEM_LP_PERIOD : "pf", |
---|
354 | EvITEM_LP_AMPLITUDE : "pf", |
---|
355 | EvITEM_BB_AMPLITUDE : "pf", |
---|
356 | EvITEM_BB_PERIOD : "pf", |
---|
357 | EvITEM_B_SLOWNESS : "pf", |
---|
358 | EvITEM_B_AZIMUTH : "pf", |
---|
359 | EvITEM_L_SLOWNESS : "pf", |
---|
360 | EvITEM_L_AZIMUTH : "pf", |
---|
361 | EvITEM_THEO_AZIM : "pf", |
---|
362 | EvITEM_THEO_BACK_AZIM : "pf", |
---|
363 | EvITEM_DISTANCE_DEG : "pf", |
---|
364 | EvITEM_DISTANCE_KM : "pf", |
---|
365 | EvITEM_QUALITY : "pi", |
---|
366 | EvITEM_MS : "pf", |
---|
367 | EvITEM_MB : "pf", |
---|
368 | EvITEM_ML : "pf", |
---|
369 | EvITEM_MW : "pf", |
---|
370 | EvITEM_MBB : "pf", |
---|
371 | EvITEM_MU : "pf", |
---|
372 | EvITEM_MEAN_MS : "ef", |
---|
373 | EvITEM_MEAN_MB : "ef", |
---|
374 | EvITEM_MEAN_ML : "ef", |
---|
375 | EvITEM_MEAN_MW : "ef", |
---|
376 | EvITEM_MEAN_MBB : "ef", |
---|
377 | EvITEM_MEAN_MU : "ef", |
---|
378 | EvITEM_MU_DESCR : "es", |
---|
379 | EvITEM_REGION : "es", |
---|
380 | EvITEM_COMMENT : "es", |
---|
381 | EvITEM_LATITUDE : "ef", |
---|
382 | EvITEM_LONGITUDE : "ef", |
---|
383 | EvITEM_LOC_METHOD : "es", |
---|
384 | EvITEM_LOC_QUALITY : "es", |
---|
385 | EvITEM_LOC_VELMOD : "es", |
---|
386 | EvITEM_LOC_ADDPAR : "es", |
---|
387 | EvITEM_FILTER : "ps", |
---|
388 | EvITEM_DEPTH : "ef", |
---|
389 | EvITEM_DEPTH_TYPE : "es", |
---|
390 | EvITEM_ORIGIN_TIME : "es", |
---|
391 | EvITEM_EVID : "ei", |
---|
392 | EvITEM_WEIGHT : "pi", |
---|
393 | EvITEM_ONSET_ACC : "pf", |
---|
394 | EvITEM_REF_LATITUDE : "ef", |
---|
395 | EvITEM_REF_LONGITUDE : "ef", |
---|
396 | EvITEM_REF_NAME : "es", |
---|
397 | EvITEM_ANALYST : "es", |
---|
398 | EvITEM_STATIONS_USED : "ei", |
---|
399 | EvITEM_REGION_TABLE : "es", |
---|
400 | EvITEM_REGION_ID : "ei", |
---|
401 | EvITEM_EVENT_TYPE : "es", |
---|
402 | EvITEM_SOURCE : "es", |
---|
403 | EvITEM_AP_SOURCE : "ps", |
---|
404 | EvITEM_ONSET_WDW_L : "pf", |
---|
405 | EvITEM_ONSET_WDW_R : "pf", |
---|
406 | EvITEM_ERR_LAT : "ef", |
---|
407 | EvITEM_ERR_LON : "ef", |
---|
408 | EvITEM_ERR_DEPTH : "ef", |
---|
409 | EvITEM_ERR_ORIGIN : "ef", |
---|
410 | EvITEM_ERR_SMAJOR : "ef", |
---|
411 | EvITEM_ERR_SMINOR : "ef", |
---|
412 | EvITEM_ERR_MAJSTRIKE : "ef", |
---|
413 | EvITEM_AZIM_MAX_GAP : "ef", |
---|
414 | EvITEM_RESID_RMS : "ef", |
---|
415 | EvITEM_PHASE_FLAGS : "es", |
---|
416 | EvITEM_SIGNOISE : "pf", |
---|
417 | EvITEM_RESIDUAL : "pf", |
---|
418 | EvITEM_RESID_CORR : "pf", |
---|
419 | EvITEM_PICK_TYPE : "ps", |
---|
420 | EvITEM_CORNERFREQ : "ef", |
---|
421 | EvITEM_LOWFREQLEVEL : "ef", |
---|
422 | EvITEM_MOMTEN : "es", |
---|
423 | EvITEM_MOMTEN_DESCR : "es", |
---|
424 | EvITEM_M0 : "ef", |
---|
425 | EvITEM_FPS_ANGLES : "es", |
---|
426 | EvITEM_FPS_DESCR : "es", |
---|
427 | } |
---|