 Timestamp:
 02/08/16 22:52:16 (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

SHX/trunk/SeismicHandler/modules/seismics.py
r1167 r1168 72 72 else: 73 73 origin = xphase.picktime  ttrav  residtime 74 75 print "dbg: depthphases", determineDepthFromDepthPhases( elat, elon ) 74 76 75 77 return (origin,elat,elon,center_lat,center_lon) … … 133 135 134 136 137 def findDepthPhases(): 138 """ 139 Returns list of depth phases together with the time distance to the main 140 phase and station name, e.g. [('pP',123.4,'GR.GRA1'),('pP',120.4,'GR.GRA2')] 141 """ 142 phaselist = PhaseList() 143 manual = phaselist.translatePicktype( 'manual' ) 144 auto = phaselist.translatePicktype( 'auto' ) 145 retlist = [] 146 for station in phaselist.getStations(): 147 depthphases = [] 148 for phase in phaselist.getPhaseList(station): 149 if phase.picktype not in (manual,auto): 150 continue 151 if phase.name.startswith('p') or phase.name.startswith('s'): 152 depthphases.append( phase ) 153 for phase in phaselist.getPhaseList(station): 154 if phase.picktype not in (manual,auto): 155 continue 156 for dphase in depthphases: 157 mainphasename = dphase.name[1:] 158 if phase.name == mainphasename: 159 retlist.append( (dphase.name, 160 (dphase.picktimephase.picktime),station) ) 161 return retlist 162 163 164 # 165 166 167 def determineDepthFromDepthPhases( epilat, epilon ): 168 dphases = findDepthPhases() 169 stations = Stations() 170 if len(dphases) == 0: 171 return None 172 deplist = [] 173 for phase,ttdist,station in dphases: 174 try: 175 r = stations[(station,None)] 176 lat = float( r.latitude ) 177 lon = float( r.longitude ) 178 except: 179 continue 180 dist, az, baz = gps2DistAzimuth( epilat, epilon, lat, lon ) 181 dist /= (1000.*DEG_TO_KM) 182 dep, acc = fit_depth( phase, ttdist, dist ) 183 deplist.append( dep ) 184 return (sum(deplist)/float(len(deplist))) 185 186 187 # 188 189 135 190 def absToRelPositions( lons, lats ): 136 191 "Compute relative distances in km from lon,lat pairs." … … 376 431 hislo = midslo 377 432 acc = midsloslo 433 if abs(acc) < minacc: 434 return (mid,acc) 435 return (mid,acc) 436 437 438 # 439 440 441 def fit_depth( phase, ttdist, distance, minacc=0.01 ): 442 "Fit depth to pair of main phase and depth phase." 443 depth_range = (1.0,750.0) 444 mainphase = phase[1:] 445 model = TauPyModel(model="iasp91") 446 lo, hi = depth_range 447 arr = model.get_travel_times( source_depth_in_km=lo, 448 distance_in_degree=distance, phase_list=[mainphase,phase] ) 449 lotime = arr[1].time  arr[0].time 450 arr = model.get_travel_times( source_depth_in_km=hi, 451 distance_in_degree=distance, phase_list=[mainphase,phase] ) 452 hitime = arr[1].time  arr[0].time 453 if ttdist < lotime: 454 return None 455 if ttdist > hitime: 456 return None 457 for i in range(50): 458 mid = (lo+hi)/2. 459 arr = model.get_travel_times( source_depth_in_km=mid, 460 distance_in_degree=distance, phase_list=[mainphase,phase] ) 461 midtime = arr[1].time  arr[0].time 462 #print loslo, midslo, hislo, '', lo, hi 463 if ttdist < midtime: 464 hi = mid 465 hitime = midtime 466 else: 467 lo = mid 468 lotime = midtime 469 acc = midtimettdist 378 470 if abs(acc) < minacc: 379 471 return (mid,acc)
Note: See TracChangeset
for help on using the changeset viewer.