Changeset 1168 for SHX


Ignore:
Timestamp:
08.02.2016 22:52:16 (4 years ago)
Author:
klaus
Message:

fit depth using depth phases

File:
1 edited

Legend:

Unmodified
Added
Removed
  • SHX/trunk/SeismicHandler/modules/seismics.py

    r1167 r1168  
    7272    else: 
    7373        origin = xphase.picktime - ttrav - residtime 
     74     
     75    print "dbg: depthphases", determineDepthFromDepthPhases( elat, elon ) 
    7476 
    7577    return (origin,elat,elon,center_lat,center_lon) 
     
    133135 
    134136 
     137def 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.picktime-phase.picktime),station) ) 
     161    return retlist 
     162 
     163 
     164#------------------------------------------------------------------------------- 
     165 
     166 
     167def 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 
    135190def absToRelPositions( lons, lats ): 
    136191    "Compute relative distances in km from lon,lat pairs." 
     
    376431                hislo = midslo 
    377432        acc = midslo-slo 
     433        if abs(acc) < minacc: 
     434            return (mid,acc) 
     435    return (mid,acc) 
     436 
     437 
     438#------------------------------------------------------------------------------- 
     439 
     440 
     441def 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 = midtime-ttdist 
    378470        if abs(acc) < minacc: 
    379471            return (mid,acc) 
Note: See TracChangeset for help on using the changeset viewer.