source: SH_SHM/trunk/util/EventDsc.py @ 92

Revision 18, 13.1 KB checked in by marcus, 15 years ago (diff)

r3 | svn | 2008-01-04 19:36:01 +0100 (Fr, 04 Jan 2008) | 1 line

fixed bug in milliseconds of origin time

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Rev Id Date
Line 
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
11import sys
12import os
13import datetime
14
15# ------------------------------------------------------------------------------
16
17
18class 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
201def 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
225def 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
250EvITEM_ARRAY          = "Array name"
251EvITEM_STATION        = "Station code"
252EvITEM_ONSET_TIME     = "Onset time"
253EvITEM_ONSET_TYPE     = "Onset type"
254EvITEM_ONSET_COUNT    = "Following Onsets"
255EvITEM_PHASE          = "Phase name"
256EvITEM_SIGN           = "Sign"
257EvITEM_COMPONENT      = "Component"
258EvITEM_PERIOD         = "Period (sec)"
259EvITEM_AMPLITUDE      = "Amplitude (nm)"
260EvITEM_AMPLITUDE_TIME = "Amplitude Time (sec)"
261EvITEM_AMPLITUDE_VEL  = "Vel. Amplitude (nm/sec)"
262EvITEM_LP_COMPONENT   = "LP component"
263EvITEM_LP_PERIOD      = "LP period (sec)"
264EvITEM_LP_AMPLITUDE   = "LP amplitude"
265EvITEM_BB_AMPLITUDE   = "BB amplitude (nm/sec)"
266EvITEM_BB_PERIOD      = "BB period (sec)"
267EvITEM_B_SLOWNESS     = "Beam-Slowness (sec/deg)"
268EvITEM_B_AZIMUTH      = "Beam-Azimuth (deg)"
269EvITEM_L_SLOWNESS     = "Epi-Slowness (sec/deg)"
270EvITEM_L_AZIMUTH      = "Epi-Azimuth (deg)"
271EvITEM_THEO_AZIM      = "Theo. Azimuth (deg)"
272EvITEM_THEO_BACK_AZIM = "Theo. Backazimuth (deg)"
273EvITEM_DISTANCE_DEG   = "Distance (deg)"
274EvITEM_DISTANCE_KM    = "Distance (km)"
275EvITEM_QUALITY        = "Quality number"
276EvITEM_MS             = "Magnitude ms"
277EvITEM_MB             = "Magnitude mb"
278EvITEM_ML             = "Magnitude ml"
279EvITEM_MW             = "Magnitude mw"
280EvITEM_MBB            = "Broadband Magnitude"
281EvITEM_MU             = "User Magnitude"
282EvITEM_MEAN_MS        = "Mean Magnitude ms"
283EvITEM_MEAN_MB        = "Mean Magnitude mb"
284EvITEM_MEAN_ML        = "Mean Magnitude ml"
285EvITEM_MEAN_MW        = "Mean Magnitude mw"
286EvITEM_MEAN_MBB       = "Mean BB Magnitude"
287EvITEM_MEAN_MU        = "Mean User Magnitude"
288EvITEM_MU_DESCR       = "User Magn. Description"
289EvITEM_REGION         = "Source region"
290EvITEM_COMMENT        = "Comment"
291EvITEM_LATITUDE       = "Latitude"
292EvITEM_LONGITUDE      = "Longitude"
293EvITEM_LOC_METHOD     = "Location method"
294EvITEM_LOC_QUALITY    = "Location quality"
295EvITEM_LOC_VELMOD     = "Velocity Model"
296EvITEM_LOC_ADDPAR     = "Location Input Params"
297EvITEM_FILTER         = "Applied filter"
298EvITEM_DEPTH          = "Depth (km)"
299EvITEM_DEPTH_TYPE     = "Depth type"
300EvITEM_ORIGIN_TIME    = "Origin time"
301EvITEM_EVID           = "Event ID"
302EvITEM_WEIGHT         = "Weight"
303EvITEM_ONSET_ACC      = "Onset Accuracy"
304EvITEM_REF_LATITUDE   = "Reference Latitude"
305EvITEM_REF_LONGITUDE  = "Reference Longitude"
306EvITEM_REF_NAME       = "Reference Location Name"
307EvITEM_ANALYST        = "Analyst"
308EvITEM_STATIONS_USED  = "No. of Stations used"
309EvITEM_REGION_TABLE   = "Region Table"
310EvITEM_REGION_ID      = "Region ID"
311EvITEM_EVENT_TYPE     = "Event Type"
312EvITEM_SOURCE         = "Source of Information"
313EvITEM_AP_SOURCE      = "Ampl&Period Source"
314EvITEM_ONSET_WDW_L    = "Onset Window Left"
315EvITEM_ONSET_WDW_R    = "Onset Window Right"
316EvITEM_ERR_LAT        = "Error in Latitude (km)"
317EvITEM_ERR_LON        = "Error in Longitude (km)"
318EvITEM_ERR_DEPTH      = "Error in Depth (km)"
319EvITEM_ERR_ORIGIN     = "Error in Origin Time"
320EvITEM_ERR_SMAJOR     = "Error Ellipse Major"
321EvITEM_ERR_SMINOR     = "Error Ellipse Minor"
322EvITEM_ERR_MAJSTRIKE  = "Error Ellipse Strike"
323EvITEM_AZIM_MAX_GAP   = "Max Azimuthal Gap (deg)"
324EvITEM_RESID_RMS      = "RMS of Residuals (sec)"
325EvITEM_PHASE_FLAGS    = "Phase Flags"
326EvITEM_SIGNOISE       = "Signal/Noise"
327EvITEM_RESIDUAL       = "Residual Time"
328EvITEM_RESID_CORR     = "Residual Correction"
329EvITEM_PICK_TYPE      = "Pick Type"
330EvITEM_CORNERFREQ     = "Corner Frequency"
331EvITEM_LOWFREQLEVEL   = "Low Frequency Level"
332EvITEM_MOMTEN         = "Moment Tensor Elements"
333EvITEM_MOMTEN_DESCR   = "Moment Tensor Descr."
334EvITEM_M0             = "Scalar Moment"
335EvITEM_FPS_ANGLES     = "Fault Plane Solution"
336EvITEM_FPS_DESCR      = "FPS Description"
337EvITEM_END_OF_EVENT   = "--- End of Phase ---"
338
339typedict = {
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}
Note: See TracBrowser for help on using the repository browser.