Changeset 1163 for SHX


Ignore:
Timestamp:
07.02.2016 01:12:13 (4 years ago)
Author:
klaus
Message:

geodesic lib for location; add origin time to location

Location:
SHX/trunk/SeismicHandler
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • SHX/trunk/SeismicHandler/basics/analysispar.py

    r1158 r1163  
    194194            ) ) 
    195195        return ret 
     196     
     197    def findPhase( self, station, phasename, picktypes ): 
     198        for picktype in picktypes: 
     199            for phase in self.getPhaseList(station,picktype): 
     200                if phase.name != phasename: 
     201                    continue 
     202                return phase 
     203        return None 
    196204     
    197205    def setDefaultPhase( self, phase ): 
  • SHX/trunk/SeismicHandler/cmdscripts/SHX_MENU_LOCATE_TELE.SHC

    r1162 r1163  
    11 
     2sdef origin 
    23sdef lat 
    34sdef lon 
     5sdef depth 
     6sdef depth_type 
    47 
    5 locate_tele p 33 &lat &lon 
     8param get source_depth &depth 
     9if "depth nes ;; goto depthok: 
     10    calc r &depth = 33. 
     11    param set source_depth 33. 
     12depthok: 
     13 
     14locate_tele p "depth &origin &lat &lon 
    615param set epi_latitude "lat 
    716param set epi_longitude "lon 
     17param set origin_time "origin 
     18 
     19param get depth_type &depth_type 
     20@IF "DEPTH_TYPE NES unknown GOTO TYPEOK: 
     21    @PARAM SET DEPTH_TYPE fixed 
     22typeok: 
    823 
    924return 
  • SHX/trunk/SeismicHandler/commands/corrpick.py

    r1161 r1163  
    1818from SeismicHandler.modules.stations import Stations 
    1919from SeismicHandler.modules.seismics import absToRelPositions, \ 
    20     computeSlownessFromPicks 
     20    computeSlownessFromPicks, planeWaveResiduals 
    2121from SeismicHandler.basics.messages import log_message 
    2222from SeismicHandler.basics.analysispar import PhaseList 
     
    218218                    reltimes.append( None ) 
    219219            # Get residuals for each entry. 
    220             resid = self.planeWaveResiduals( reltimes, relx, rely ) 
     220            resid = planeWaveResiduals( reltimes, relx, rely ) 
    221221            #print "dbg: names", names 
    222222            #print "dbg: resid", resid 
     
    263263        return (optcorrinfo,minvecnorm) 
    264264     
    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. 
    268         if None in reltimes: 
    269             rreltimes = [] 
    270             rrelx = [] 
    271             rrely = [] 
    272             for i,relt in enumerate(reltimes): 
    273                 if relt == None: 
    274                     continue 
    275                 rreltimes.append( relt ) 
    276                 rrelx.append( relx[i] ) 
    277                 rrely.append( rely[i] ) 
    278         else: 
    279             rreltimes = reltimes 
    280             rrelx = relx 
    281             rrely = rely 
    282         az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes), 
    283             np.array(rrelx), np.array(rrely) ) 
    284         # Now we have slowness and azimuth, let's compute residuals. 
    285         reltimes_mean = sum(rreltimes)/float(len(rreltimes)) 
    286         # compute beam shifts (as in command 'beam') 
    287         radaz = (90.0 - az) / (180./np.pi) 
    288         kmslo = slo / DEG_TO_KM 
    289         bshifts = [] 
    290         for i,relt in enumerate(reltimes): 
    291             bshifts.append( 
    292                 (relx[i]*np.cos(radaz) + rely[i]*np.sin(radaz)) * kmslo 
    293             ) 
    294         bshifts_mean = sum(bshifts)/float(len(bshifts)) 
    295         resids = [] 
    296         # Residuals consists of sum (difference) of correlation maxima and beam. 
    297         for i,relt in enumerate(reltimes): 
    298             if relt == None: 
    299                 resid = None 
    300             else: 
    301                 resid = abs( (bshifts[i]-bshifts_mean) + (relt-reltimes_mean) ) 
    302             resids.append( resid ) 
    303         return np.array(resids) 
    304  
    305265    def modifyCorrinfo( self, corrinfo, resid, names, maxtol=2., loopcnt=50 ): 
    306266        """ 
  • SHX/trunk/SeismicHandler/commands/locate_tele.py

    r1162 r1163  
    88from SeismicHandler.basics.command import BaseCommand 
    99from SeismicHandler.modules.seismics import locateTeleEvent 
     10from obspy.sh.core import fromUTCDateTime 
    1011 
    1112 
     
    1516    URI:http://www.seismic-handler.org/portal/wiki/ShLocateTele 
    1617    """ 
    17     numberOfParameters = [1,2,4] 
     18    numberOfParameters = [1,2,5] 
    1819    parameterQueries = [ 
    1920        { 
     
    2526            "text": "depth", 
    2627            "type": "float", 
     28            "question": False, 
     29        }, 
     30        { 
     31            "text": "symbol origin", 
     32            "type": "str", 
    2733            "question": False, 
    2834        }, 
     
    5258        if len(self.parameters) > 1: 
    5359            depth = float(self.parameters[1]) 
     60        origsym = None 
     61        if len(self.parameters) > 2: 
     62            origsym = self.parameters[2][1:] 
    5463        latsym = None 
    55         if len(self.parameters) > 2: 
    56             latsym = self.parameters[2][1:] 
    5764        if len(self.parameters) > 3: 
    58             lonsym = self.parameters[3][1:] 
    59         try: 
    60             lat, lon = locateTeleEvent( phase, depth ) 
    61         except: 
    62             return 
    63         print "dbg: location", lat, lon 
     65            latsym = self.parameters[3][1:] 
     66        lonsym = None 
     67        if len(self.parameters) > 4: 
     68            lonsym = self.parameters[4][1:] 
     69        #try: 
     70        origin, lat, lon = locateTeleEvent( phase, depth ) 
     71        #except: 
     72        #    return 
     73        self.symbols.set( origsym, fromUTCDateTime(origin) ) 
    6474        self.symbols.set( latsym, "%g" % lat ) 
    6575        self.symbols.set( lonsym, "%g" % lon ) 
  • SHX/trunk/SeismicHandler/modules/seismics.py

    r1162 r1163  
    1414from SeismicHandler.basics.analysispar import PhaseList 
    1515from SeismicHandler.modules.stations import Stations 
     16from obspy.core.util.geodetics import gps2DistAzimuth 
     17try: 
     18    from geographiclib.geodesic import Geodesic 
     19    geodesiclib_available = True 
     20except: 
     21    geodesiclib_available = False 
     22 
     23 
     24DEG_TO_KM = 111.19 
     25 
     26 
     27#------------------------------------------------------------------------------- 
    1628 
    1729 
    1830def locateTeleEvent( phase, depth=33. ): 
    19     pslo = getPhaseSlownesses( withcenter=True ) 
     31    """ 
     32    Locate teleseismic event using the phases of the given name defined on 
     33    all traces. 
     34    """ 
     35    pslo, presid = getPhaseSlownesses( withcenter=True, withresiduals=True ) 
    2036    if not phase in pslo: 
    2137        return None 
    2238    slowness,azimuth,slowerr,azimerr,center_lat,center_lon = pslo[phase] 
    2339    dist,acc = fit_distance( phase, slowness, depth=depth ) 
    24     return sphereadd( center_lat, center_lon, dist, azimuth ) 
    25  
    26  
    27 def getPhaseSlownesses( withcenter=False ): 
     40    if geodesiclib_available: 
     41        res = Geodesic.WGS84.Direct( center_lat, center_lon, azimuth, 
     42            dist*DEG_TO_KM*1000. ) 
     43        #print "dbg: sphere cmp:", sphereadd(center_lat,center_lon,dist,azimuth) 
     44        elat, elon = (res['lat2'],res['lon2']) 
     45    else: 
     46        # Compute on sphere 
     47        #print "dbg: geodesic lib not available" 
     48        elat, elon = sphereadd( center_lat, center_lon, dist, azimuth ) 
     49     
     50    # origin time 
     51    for s,r in presid[phase]: 
     52        if r != None: 
     53            station = s 
     54            residtime = r 
     55            break 
     56    r = Stations()[(station,None)] 
     57    slat = float( r.latitude ) 
     58    slon = float( r.longitude ) 
     59    dist, az, baz = gps2DistAzimuth( elat, elon, slat, slon ) 
     60    dist /= (DEG_TO_KM*1000.) 
     61    model = TauPyModel(model="iasp91") 
     62    arr = model.get_travel_times( source_depth_in_km=depth, 
     63        distance_in_degree=dist, phase_list=[phase] ) 
     64    ttrav = arr[0].time 
     65    phaselist = PhaseList() 
     66    manual = phaselist.translatePicktype('manual') 
     67    auto = phaselist.translatePicktype('auto') 
     68    xphase = phaselist.findPhase( station, phase, (manual,auto) ) 
     69    if xphase == None: 
     70        origin = None 
     71    else: 
     72        origin = xphase.picktime - ttrav - residtime 
     73 
     74    return (origin,elat,elon) 
     75 
     76 
     77#------------------------------------------------------------------------------- 
     78 
     79 
     80def getPhaseSlownesses( withcenter=False, withresiduals=False ): 
    2881    "Compute slownesses for all phases in phase list (>4 picks)." 
    2982 
    3083    pslo = {} 
     84    presid = {} 
    3185    phaselist = PhaseList() 
    3286    plist = phaselist.getPhaseListsByName() 
     
    3690        lats = [] 
    3791        abstimes = [] 
     92        stalist = [] 
    3893        for station,picktime in plist[phase]: 
    3994            try: 
     
    46101            lats.append( lat ) 
    47102            abstimes.append( picktime ) 
     103            stalist.append( station ) 
    48104        if len(lons) < 3: 
    49105            continue 
     
    62118        else: 
    63119            pslo[phase] = (slowness,azimuth,slowerr,azimerr) 
    64     return pslo 
     120        if withresiduals: 
     121            resid = planeWaveResiduals( times, relx, rely, takeabs=False ) 
     122            presid[phase] = [] 
     123            for i,iresid in enumerate(resid): 
     124                presid[phase].append( (stalist[i],iresid) ) 
     125    if withresiduals: 
     126        return (pslo,presid) 
     127    else: 
     128        return pslo 
     129 
     130 
     131#------------------------------------------------------------------------------- 
    65132 
    66133 
     
    78145        rely.append( y ) 
    79146    return (relx,rely) 
     147 
     148 
     149#------------------------------------------------------------------------------- 
     150 
    80151 
    81152def computeSlownessFromPicks( reltimes, relx, rely, dt=None ): 
     
    132203         
    133204    RAD_TO_DEG = 180.0/np.pi 
    134     DEG_TO_KM = 111.19 
    135205     
    136206    # coordinate shift to center 
     
    209279            norm * norm * RAD_TO_DEG 
    210280    return (azimuth,azimerr,slowness,slowerr) 
     281 
     282 
     283#------------------------------------------------------------------------------- 
     284 
     285 
     286def planeWaveResiduals( reltimes, relx, rely, takeabs=True ): 
     287    "Compute residuals of given shift times to a plane wave." 
     288    # Bit complicated due to possible None values in reltimes. 
     289    if None in reltimes: 
     290        rreltimes = [] 
     291        rrelx = [] 
     292        rrely = [] 
     293        for i,relt in enumerate(reltimes): 
     294            if relt == None: 
     295                continue 
     296            rreltimes.append( relt ) 
     297            rrelx.append( relx[i] ) 
     298            rrely.append( rely[i] ) 
     299    else: 
     300        rreltimes = reltimes 
     301        rrelx = relx 
     302        rrely = rely 
     303    az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes), 
     304        np.array(rrelx), np.array(rrely) ) 
     305    # Now we have slowness and azimuth, let's compute residuals. 
     306    reltimes_mean = sum(rreltimes)/float(len(rreltimes)) 
     307    # compute beam shifts (as in command 'beam') 
     308    radaz = (90.0 - az) / (180./np.pi) 
     309    kmslo = slo / DEG_TO_KM 
     310    bshifts = [] 
     311    for i,relt in enumerate(reltimes): 
     312        bshifts.append( 
     313            (relx[i]*np.cos(radaz) + rely[i]*np.sin(radaz)) * kmslo 
     314        ) 
     315    bshifts_mean = sum(bshifts)/float(len(bshifts)) 
     316    resids = [] 
     317    # Residuals consist of sum (difference) of correlation maxima and beam. 
     318    for i,relt in enumerate(reltimes): 
     319        if relt == None: 
     320            resid = None 
     321        else: 
     322            resid = (bshifts[i]-bshifts_mean) + (relt-reltimes_mean) 
     323            if takeabs: 
     324                resid = abs( resid ) 
     325        resids.append( resid ) 
     326    return np.array(resids) 
     327 
     328 
     329#------------------------------------------------------------------------------- 
    211330 
    212331 
     
    261380 
    262381 
     382#------------------------------------------------------------------------------- 
     383 
     384 
    263385def sphereadd( slat, slon, dist, bazim ): 
    264386 
     
    319441     
    320442    return (elat,elon) 
     443 
     444 
     445#------------------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.