Changeset 1162 for SHX


Ignore:
Timestamp:
06.02.2016 01:55:22 (4 years ago)
Author:
klaus
Message:

locate tele command (to be checked)

Location:
SHX/trunk/SeismicHandler
Files:
2 added
2 edited

Legend:

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

    r1158 r1162  
    1111import numpy as np 
    1212from obspy.signal.util import utlGeoKm 
     13from obspy.taup import TauPyModel 
    1314from SeismicHandler.basics.analysispar import PhaseList 
    1415from SeismicHandler.modules.stations import Stations 
    1516 
    1617 
    17 def getPhaseSlownesses(): 
     18def locateTeleEvent( phase, depth=33. ): 
     19    pslo = getPhaseSlownesses( withcenter=True ) 
     20    if not phase in pslo: 
     21        return None 
     22    slowness,azimuth,slowerr,azimerr,center_lat,center_lon = pslo[phase] 
     23    dist,acc = fit_distance( phase, slowness, depth=depth ) 
     24    return sphereadd( center_lat, center_lon, dist, azimuth ) 
     25 
     26 
     27def getPhaseSlownesses( withcenter=False ): 
    1828    "Compute slownesses for all phases in phase list (>4 picks)." 
    1929 
     
    4555        azimuth,azimerr,slowness,slowerr = computeSlownessFromPicks( 
    4656            np.array(times), np.array(relx), np.array(rely) ) 
    47         pslo[phase] = (slowness,azimuth,slowerr,azimerr) 
     57        if withcenter: 
     58            center_lon = sum(lons)/float(len(lons)) 
     59            center_lat = sum(lats)/float(len(lats)) 
     60            pslo[phase] = (slowness,azimuth,slowerr,azimerr, 
     61                center_lat,center_lon) 
     62        else: 
     63            pslo[phase] = (slowness,azimuth,slowerr,azimerr) 
    4864    return pslo 
    4965 
     
    193209            norm * norm * RAD_TO_DEG 
    194210    return (azimuth,azimerr,slowness,slowerr) 
     211 
     212 
     213def fit_distance( phase, slo, depth=33., minacc=0.01 ): 
     214    dist_range = { 
     215        'P'   : (10.0,97.0), 
     216    } 
     217    if phase not in dist_range: 
     218        return None 
     219    model = TauPyModel(model="iasp91") 
     220    lo, hi = dist_range[phase] 
     221    arr = model.get_travel_times( source_depth_in_km=depth, 
     222        distance_in_degree=lo, phase_list=[phase] ) 
     223    loslo = arr[0].ray_param_sec_degree 
     224    arr = model.get_travel_times( source_depth_in_km=depth, 
     225        distance_in_degree=hi, phase_list=[phase] ) 
     226    hislo = arr[0].ray_param_sec_degree 
     227    if loslo < hislo: 
     228        if slo < loslo: 
     229            return None 
     230        if slo > hislo: 
     231            return None 
     232    else: 
     233        if slo < hislo: 
     234            return None 
     235        if slo > loslo: 
     236            return None 
     237    for i in range(50): 
     238        mid = (lo+hi)/2. 
     239        arr = model.get_travel_times( source_depth_in_km=depth, 
     240            distance_in_degree=mid, phase_list=[phase] ) 
     241        midslo = arr[0].ray_param_sec_degree 
     242        #print loslo, midslo, hislo, '|', lo, hi 
     243        if loslo < hislo: 
     244            if slo < midslo: 
     245                hi = mid 
     246                hislo = midslo 
     247            else: 
     248                lo = mid 
     249                loslo = midslo 
     250        else: 
     251            if slo < midslo: 
     252                lo = mid 
     253                loslo = midslo 
     254            else: 
     255                hi = mid 
     256                hislo = midslo 
     257        acc = midslo-slo 
     258        if abs(acc) < minacc: 
     259            return (mid,acc) 
     260    return (mid,acc) 
     261 
     262 
     263def sphereadd( slat, slon, dist, bazim ): 
     264 
     265    """ 
     266    /* computes new location from input location (slat,slon) and distance 
     267     * "dist" and azimuth "bazim" 
     268     * 
     269     * parameters of routine 
     270     * float      slat, slon;    input; (station) input location in degrees 
     271     * float      dist;          input; distance in degrees 
     272     * float      bazim;         input; (back-)azimuth in degrees 
     273     * float      *elat, *elon;  output; (epicentre) location in degrees 
     274     */ 
     275    #define DEG2RAD(x) ((x)*BC_PI/180.0) 
     276    #define RAD2DEG(x) ((x)/BC_PI*180.0) 
     277    """ 
     278    # transformation to mathematical system 
     279    slat = 90.0 - slat 
     280    if  (slon < 0.0): 
     281        slon = 360 + slon 
     282 
     283    invert = (bazim > 180.0) 
     284    if  invert: 
     285        bazim = 360 - bazim 
     286 
     287    # degrees to radian 
     288    slat = slat*np.pi/180. 
     289    dist = dist*np.pi/180. 
     290    bazim = bazim*np.pi/180. 
     291 
     292    arg = np.cos(dist)*np.cos(slat) + np.sin(dist)*np.sin(slat)*np.cos(bazim) 
     293    if arg > 1.0: 
     294        arg = 1.0 
     295    if arg < -1.0: 
     296        arg = -1.0 
     297    elat = np.arccos( arg ); 
     298    arg = (np.cos(dist)-np.cos(elat)*np.cos(slat))/(np.sin(elat)*np.sin(slat)) 
     299    if arg > 1.0: 
     300        arg = 1.0 
     301    if arg < -1.0: 
     302        arg = -1.0 
     303    gamma = np.arccos( arg ) 
     304 
     305    # back to degrees 
     306    elat = elat*180./np.pi 
     307    gamma = gamma*180./np.pi 
     308 
     309    if invert: 
     310        gamma = 360.0 - gamma 
     311    elon = slon + gamma 
     312    if elon > 360.0: 
     313        elon -= 360.0 
     314 
     315    # transformation back to geo-system 
     316    elat = 90.0 - elat 
     317    if elon > 180.0: 
     318        elon -= 360.0 
     319     
     320    return (elat,elon) 
  • SHX/trunk/SeismicHandler/modules/wx_.py

    r1160 r1162  
    208208                return 
    209209        shortcut = { 
     210            'a': plotter.OnSimulateWoodAnderson, 
    210211            'b': plotter.OnBeam, 
    211212            'd': plotter.OnDelTimeWindow, 
     
    217218            's': plotter.OnSetTimeWindow, 
    218219            'u': plotter.OnSimulateUndo, 
    219             'w': plotter.OnSimulateWoodAnderson, 
     220            'w': plotter.OnPlaneWave, 
    220221            '>': plotter.OnMoveZoomRight, 
    221222            '<': plotter.OnMoveZoomLeft, 
     
    19911992     
    19921993    Shortcuts 
    1993            b(Beam) d(DelTimeWdw) g(PhasePg) l(locate) m(mlSingle) n(PhasePn) 
    1994            p(PhaseP) s(SetTimeWdw) u(Undo Filter) w(FilterWA) 0(PrevEvent) 
    1995            1(NextEvent) 2(PhaseSg) 8(Filter 1-8Hz) >(zoom->right) <(zoom->left) 
     1994           a(FilterWA) b(Beam) d(DelTimeWdw) g(PhasePg) l(locate) m(mlSingle) 
     1995           n(PhasePn) p(PhaseP) s(SetTimeWdw) u(Undo Filter) w(planewave) 
     1996           0(PrevEvent) 1(NextEvent) 2(PhaseSg) 8(Filter 1-8Hz) >(zoom->right) 
     1997           <(zoom->left) 
    19961998    Ctrl-  f(FK) l(FilterSROLP) n(Norm) o(Overlapping) q(Quit) s(SaveT&P) 
    19971999           w(FilterWWSSNSP) y(SortByDistance) 
     
    21532155        self.addEntry( arrayMenu, 'Beam / Del Beam\tB', 
    21542156            'Compute/delete beam for all traces', self.OnBeam ) 
    2155         self.addEntry( arrayMenu, 'Plane Wave', 
     2157        self.addEntry( arrayMenu, 'Plane Wave\tW', 
    21562158            'Compute slowness and azimuth using picks', self.OnPlaneWave ) 
    21572159        self.addEntry( arrayMenu, 'Correlation Picker', 
     
    21602162        # simulate menu 
    21612163        simulateMenu = wx.Menu() 
    2162         self.addEntry( simulateMenu, 'Simulate &Wood-Anderson\tW', 
     2164        self.addEntry( simulateMenu, 'Simulate &Wood-Anderson\tA', 
    21632165            'Simulate Wood-Anderson instrument on all traces', 
    21642166            self.OnSimulateWoodAnderson ) 
     
    21942196        eventsMenu.AppendSeparator() 
    21952197        self.addEntry( eventsMenu, 'Event Info Parser ...', 
    2196             'Read event from text line', self.EventInfoParser ) 
     2198            'Read event from text line', self.OnEventInfoParser ) 
    21972199        eventsMenu.AppendSeparator() 
    2198         self.addEntry( eventsMenu, 'Locate Event\tl', 
     2200        self.addEntry( eventsMenu, 'Locate Local Event\tl', 
    21992201            'Use phase list for locating event', self.OnLocateEvent ) 
     2202        self.addEntry( eventsMenu, 'Locate Tele Event with P', 
     2203            'Use P arrivals for teleseismic location', self.OnLocateTeleEvent ) 
    22002204        self.addEntry( eventsMenu, 'Plot Location', 
    22012205            'Simple map with stations and epicenter', self.OnPlotLocation ) 
     
    26842688        _sendShCommand( 'shx_menu_next_event -1' ) 
    26852689     
    2686     def EventInfoParser( self, e ): 
     2690    def OnEventInfoParser( self, e ): 
     2691        "plotterWindow: events menu entry." 
    26872692        dlg = EventParseDialog( self ) 
    26882693        dlg.ShowModal() 
     
    26922697        "plotterWindow: events menu entry." 
    26932698        _sendShCommand( 'locate shxvar /showresult' ) 
     2699     
     2700    def OnLocateTeleEvent( self, e ): 
     2701        "plotterWindow: events menu entry." 
     2702        _sendShCommand( "shx_menu_locate_tele" ) 
    26942703         
    26952704    def OnPlotLocation( self, e ): 
Note: See TracChangeset for help on using the changeset viewer.