 Timestamp:
 02/05/16 21:35:37 (4 years ago)
 Location:
 SHX/trunk/SeismicHandler/commands
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

SHX/trunk/SeismicHandler/commands/corrpick.py
r1160 r1161 6 6 7 7 """ 8 Correlation picker. 8 Smart correlation picker. Computes correlation of the selected time window with 9 all other traces. Takes the N (N=3) largest correlation values and takes the 10 best fitting a plane wave. If residuals stay over a tolerance value, no pick 11 is set on this trace. 9 12 """ 10 13 … … 31 34 """ 32 35 # Define number and type of parameters, possible qualifiers. 33 numberOfParameters = [4 ]36 numberOfParameters = [4,5] 34 37 parameterQueries = [ 35 38 { … … 55 58 ] 56 59 known_qualifiers = [ 60 "CORRMAXIMA", 61 "MAXTOL", 62 "MAXLOOP", 57 63 ] 58 64 redraw = True … … 67 73 #@timeit 68 74 def run(self): 69 " Correlation picker."75 "Smart correlation picker." 70 76 71 77 nummax = 3 78 if self.qualifiers["CORRMAXIMA"]: 79 nummax = int( self.qualifers["CORRMAXIMA"] ) 80 maxtol = 2. 81 if self.qualifiers["MAXTOL"]: 82 maxtol = float( self.qualifers["MAXTOL"] ) 83 maxloop = 100 84 if self.qualifiers["MAXLOOP"]: 85 maxloop = int( self.qualifers["MAXLOOP"] ) 72 86 self.modifdict = {} 73 87 self.maxmodif = 2*nummax … … 78 92 self.phaselist = PhaseList() 79 93 trace = traces[0] 94 # compute time windows for correlation 80 95 _start = float( self.parameters[1] ) 81 96 _end = float( self.parameters[2] ) 82 97 _shift = float( self.parameters[3] ) 98 if len(self.parameters) > 4: 99 outvar = self.parameters[4][1:] 100 else: 101 outvar = None 83 102 wdwlen = _end  _start 84 103 trace_sname = "%s.%s" % (trace.stats.network,trace.stats.station) 85 # find phase 104 # find phase in selected window 86 105 absstart = trace.stats.starttime + _start  trace.get_info('torigin') 87 106 absend = absstart + wdwlen … … 96 115 shiftend = trace.stats.endtime 97 116 _shift = ((shiftend  shiftstart)  wdwlen)/2. 117 # Prepare traces ad correlate. 118 # Results are in corrinfo, which is a dictionary on station names. 119 # Each station has a list of pairs (reltime,relcorrvalue), reltime is 120 # the time of the correlation maxima, relcorrvalue is the relative 121 # correlation coefficient (ratio to the largest coefficient). 98 122 corrtrc = trace.slice( shiftstart, shiftend ) 99 123 shiftlength = int( _shift / trace.stats.delta )/2  1 … … 127 151 except: 128 152 pass 129 corrinfo = self.findOptShifts( corrinfo ) 153 # Modify the maxima to fit a plane wave. 154 corrinfo, res = self.findOptShifts( corrinfo, maxtol, maxloop ) 155 # Define phases as given in the first list entry in each station. 130 156 for sname in corrinfo.keys(): 131 reltime = corrinfo[sname][0][0] 157 if len(corrinfo[sname]) > 0: 158 reltime = corrinfo[sname][0][0] 159 else: 160 reltime = None 132 161 self.phaselist.clearPhasenameStation( phase, sname, auto ) 133 162 if reltime == None: … … 135 164 self.phaselist.addPhasePick( sname, phasetime+reltime, name=phase, 136 165 picktype=auto ) 166 if outvar: 167 self.symbols.set( outvar, "%g" % res ) 137 168 138 169 def findPhase( self, trc, start, end ): 170 "Determin name and time of phase in the selected time window." 139 171 sname = "%s.%s" % (trc.stats.network,trc.stats.station) 140 172 manual = self.phaselist.translatePicktype( 'manual' ) … … 150 182 return (None,None) 151 183 152 def findOptShifts( self, corrinfo ): 184 def findOptShifts( self, corrinfo, maxtol=2., maxloop=100 ): 185 """ 186 Modify the correlation maxima (cyclic exchange or delete) to match 187 a plane wave. 188 """ 153 189 # find station locations 154 190 optcorrinfo = None … … 168 204 lats.append( lat ) 169 205 lons.append( lon ) 206 # recompute to relative coos in km 170 207 relx, rely = absToRelPositions( lons, lats ) 171 for i in range(50): 208 # Several tries to modify corrinfo. 209 nonecnt = 0 210 for i in range(maxloop): 211 # Collect rel. times from correlation maxima in the current order. 212 # 'None' values at stations with too large persistent residuals. 172 213 reltimes = [] 173 214 for s in names: … … 176 217 else: 177 218 reltimes.append( None ) 219 # Get residuals for each entry. 178 220 resid = self.planeWaveResiduals( reltimes, relx, rely ) 179 print "dbg: names", names 180 print "dbg: resid", resid 181 print "dbg: modif", self.modifdict 221 #print "dbg: names", names 222 #print "dbg: resid", resid 223 #print "dbg: modif", self.modifdict 224 # Compute quality (norm of resid vector, excluding Nones) 182 225 try: 183 226 vecnorm = np.sqrt( (resid*resid).sum() ) … … 189 232 continue 190 233 sqsum += r*r 191 vecnomr = np.sqrt( sqsum ) 234 vecnorm = np.sqrt( sqsum ) 235 # Store config with minimum residuals. 192 236 if minvecnorm == None or vecnorm < minvecnorm: 193 237 minvecnorm = vecnorm 194 238 optcorrinfo = deepcopy( corrinfo ) 195 print "dbg: copied corrinfo at", vecnorm 239 #print "dbg: copied corrinfo at", vecnorm 240 # Evaluate improvement to last step. If result is much worse than 241 # before, delete this maximum. 196 242 if lastvecnorm == None: 197 243 improvement = 0. … … 199 245 improvement = lastvecnorm  vecnorm 200 246 if improvement < 5. and modsta != None: 201 # much worse, throw away this corr247 # Much worse, throw away this corr. 202 248 if len(corrinfo[modsta]) > 1: 203 249 corrinfo[modsta] = corrinfo[modsta][1:] 204 print "dbg: deleted entry at", modsta250 #print "dbg: deleted entry at", modsta 205 251 lastvecnorm = vecnorm 206 print "dbg: resid norm", vecnorm, improvement 207 # put first correlation max to last position 208 modsta = self.modifyCorrinfo( corrinfo, resid, names ) 209 return optcorrinfo 252 #print "dbg: resid norm", vecnorm, improvement 253 # Put first correlation max to last position and reiterate. 254 modsta = self.modifyCorrinfo( corrinfo, resid, names, maxtol ) 255 #print "dbg: findOptShift counter", i 256 if modsta == None: 257 nonecnt += 1 258 else: 259 nonecnt = 0 260 if nonecnt > 3: 261 break 262 #print "dbg: minimum norm", minvecnorm 263 return (optcorrinfo,minvecnorm) 210 264 211 265 def planeWaveResiduals( self, reltimes, relx, rely ): 266 "Compute residuals of given shift times to a plane wave." 267 # Bit complicated due to possible None values. 212 268 if None in reltimes: 213 269 rreltimes = [] … … 226 282 az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes), 227 283 np.array(rrelx), np.array(rrely) ) 284 # Now we have slowness and azimuth, let's compute residuals. 228 285 reltimes_mean = sum(rreltimes)/float(len(rreltimes)) 229 # compute beam shifts (as in command 'beam' 286 # compute beam shifts (as in command 'beam') 230 287 radaz = (90.0  az) / (180./np.pi) 231 288 kmslo = slo / DEG_TO_KM … … 237 294 bshifts_mean = sum(bshifts)/float(len(bshifts)) 238 295 resids = [] 296 # Residuals consists of sum (difference) of correlation maxima and beam. 239 297 for i,relt in enumerate(reltimes): 240 298 if relt == None: … … 245 303 return np.array(resids) 246 304 247 def modifyCorrinfo( self, corrinfo, resid, names, maxtol=2. ): 248 for i in range(50): 305 def modifyCorrinfo( self, corrinfo, resid, names, maxtol=2., loopcnt=50 ): 306 """ 307 Modify corrinfo for the next iteration. Returns name of changed station. 308 resid and names contain the residuals and station names (in the same 309 order). 310 """ 311 for i in range(loopcnt): 312 # Determine maximum residual, corresponding station 249 313 amax = resid.argmax() 250 314 maxresid = resid[amax] 251 315 maxsta = names[amax] 316 #print "dbg: argmax", amax, maxsta, len(corrinfo[maxsta]), repr(maxresid) 317 # Globally count modifications at this station. 252 318 if not maxsta in self.modifdict: 253 319 self.modifdict[maxsta] = 0 254 # delete entry if max counter reach ed and still over tol 255 if self.modifdict[maxsta] >= self.maxmodif \ 256 and len(corrinfo[maxsta]) > 0 and maxresid > maxtol: 320 # Delete current (=first) entry if max counter reached or only one 321 # entry left and still over tolerance. 322 if (self.modifdict[maxsta] >= self.maxmodif \ 323 and len(corrinfo[maxsta]) > 0) or len(corrinfo[maxsta]) == 1 \ 324 and maxresid > maxtol: 257 325 if len(corrinfo[maxsta]) > 1: 326 # More than one, reduce by one. 258 327 corrinfo[maxsta] = corrinfo[maxsta][1:] 328 #print "dbg: delete entry at", maxsta 259 329 else: 330 # Delete the last one. 260 331 corrinfo[maxsta] = [] 261 print "dbg: *** empty", maxsta332 #print "dbg: *** empty", maxsta 262 333 self.modifdict[maxsta] += 1 263 334 return maxsta 335 # Change the station with a large residual below maximum count or 336 # a station with a smaller residual and comparable correlation 337 # maximum. 264 338 if self.modifdict[maxsta] < self.maxmodif \ 265 339 and len(corrinfo[maxsta]) > 1 \ 266 and ( corrinfo[maxsta][0][0]> maxtol \340 and (maxresid > maxtol \ 267 341 or corrinfo[maxsta][0][1] > 0.95): 268 print "dbg: rotating", maxsta, corrinfo[maxsta][0]342 #print "dbg: rotating", maxsta, corrinfo[maxsta][0] 269 343 corrinfo[maxsta] = corrinfo[maxsta][1:] + [corrinfo[maxsta][0]] 270 344 self.modifdict[maxsta] += 1 271 345 return maxsta 272 346 resid[amax] = resid[resid.argmin()] 273 print "dbg: failed to modify corrinfo" 347 #print "dbg: modifyCorrinfo counter", i 348 #print "dbg: failed to modify corrinfo" 274 349 return None 275 350 
SHX/trunk/SeismicHandler/commands/phase.py
r1158 r1161 183 183 elif subcmd == 'slowness': 184 184 pdict = getPhaseSlownesses() 185 print "dbg:", pdict186 185 pname = self.parameters[1] 187 186 if pname not in pdict.keys():
Note: See TracChangeset
for help on using the changeset viewer.