source: SHX/trunk/SeismicHandler/modules/seismics.py @ 1167

Revision 1167, 16.1 KB checked in by klaus, 4 years ago (diff)

bugfix in amplitude display; little progress in compare_location

Line 
1# -*- coding: utf-8 -*-
2
3#    This file is part of Seismic Handler eXtended (SHX). For terms of use and
4#    license information please see license.txt and visit
5#    http://www.seismic-handler.org/portal/wiki/Shx/LicenseTerms
6
7"""
8Seismic algorithms.
9"""
10
11import numpy as np
12from obspy.signal.util import utlGeoKm
13from obspy.taup import TauPyModel
14from SeismicHandler.basics.analysispar import PhaseList
15from SeismicHandler.modules.stations import Stations
16from obspy.core.util.geodetics import gps2DistAzimuth
17from obspy.fdsn.client import Client
18try:
19    from geographiclib.geodesic import Geodesic
20    geodesiclib_available = True
21except:
22    geodesiclib_available = False
23
24
25DEG_TO_KM = 111.19
26
27
28#-------------------------------------------------------------------------------
29
30
31def locateTeleEvent( phase, depth=33. ):
32    """
33    Locate teleseismic event using the phases of the given name defined on
34    all traces.
35    """
36    pslo, presid = getPhaseSlownesses( withcenter=True, withresiduals=True )
37    if not phase in pslo:
38        return None
39    slowness,azimuth,slowerr,azimerr,center_lat,center_lon = pslo[phase]
40    dist,acc = fit_distance( phase, slowness, depth=depth )
41    if geodesiclib_available:
42        res = Geodesic.WGS84.Direct( center_lat, center_lon, azimuth,
43            dist*DEG_TO_KM*1000. )
44        #print "dbg: sphere cmp:", sphereadd(center_lat,center_lon,dist,azimuth)
45        elat, elon = (res['lat2'],res['lon2'])
46    else:
47        # Compute on sphere
48        #print "dbg: geodesic lib not available"
49        elat, elon = sphereadd( center_lat, center_lon, dist, azimuth )
50   
51    # origin time
52    for s,r in presid[phase]:
53        if r != None:
54            station = s
55            residtime = r
56            break
57    r = Stations()[(station,None)]
58    slat = float( r.latitude )
59    slon = float( r.longitude )
60    dist, az, baz = gps2DistAzimuth( elat, elon, slat, slon )
61    dist /= (DEG_TO_KM*1000.)
62    model = TauPyModel(model="iasp91")
63    arr = model.get_travel_times( source_depth_in_km=depth,
64        distance_in_degree=dist, phase_list=[phase] )
65    ttrav = arr[0].time
66    phaselist = PhaseList()
67    manual = phaselist.translatePicktype('manual')
68    auto = phaselist.translatePicktype('auto')
69    xphase = phaselist.findPhase( station, phase, (manual,auto) )
70    if xphase == None:
71        origin = None
72    else:
73        origin = xphase.picktime - ttrav - residtime
74
75    return (origin,elat,elon,center_lat,center_lon)
76
77
78#-------------------------------------------------------------------------------
79
80
81def getPhaseSlownesses( withcenter=False, withresiduals=False ):
82    "Compute slownesses for all phases in phase list (>4 picks)."
83
84    pslo = {}
85    presid = {}
86    phaselist = PhaseList()
87    plist = phaselist.getPhaseListsByName()
88    stations = Stations()
89    for phase in plist:
90        lons = []
91        lats = []
92        abstimes = []
93        stalist = []
94        for station,picktime in plist[phase]:
95            try:
96                r = stations[(station,None)]
97                lat = float( r.latitude )
98                lon = float( r.longitude )
99            except:
100                continue
101            lons.append( lon )
102            lats.append( lat )
103            abstimes.append( picktime )
104            stalist.append( station )
105        if len(lons) < 3:
106            continue
107        relx, rely = absToRelPositions( lons, lats )
108        mintime = min( abstimes )
109        times = []
110        for a in abstimes:
111            times.append( a-mintime )
112        azimuth,azimerr,slowness,slowerr = computeSlownessFromPicks(
113            np.array(times), np.array(relx), np.array(rely) )
114        if withcenter:
115            center_lon = sum(lons)/float(len(lons))
116            center_lat = sum(lats)/float(len(lats))
117            pslo[phase] = (slowness,azimuth,slowerr,azimerr,
118                center_lat,center_lon)
119        else:
120            pslo[phase] = (slowness,azimuth,slowerr,azimerr)
121        if withresiduals:
122            resid = planeWaveResiduals( times, relx, rely, takeabs=False )
123            presid[phase] = []
124            for i,iresid in enumerate(resid):
125                presid[phase].append( (stalist[i],iresid) )
126    if withresiduals:
127        return (pslo,presid)
128    else:
129        return pslo
130
131
132#-------------------------------------------------------------------------------
133
134
135def absToRelPositions( lons, lats ):
136    "Compute relative distances in km from lon,lat pairs."
137    if len(lons) == 0 or len(lons) != len(lats):
138        return None
139    center_lon = sum(lons)/float(len(lons))
140    center_lat = sum(lats)/float(len(lats))
141    relx = []
142    rely = []
143    for i in range(len(lons)):
144        x, y = utlGeoKm( center_lon, center_lat, lons[i], lats[i] )
145        relx.append( x )
146        rely.append( y )
147    return (relx,rely)
148
149
150#-------------------------------------------------------------------------------
151
152
153def computeSlownessFromPicks( reltimes, relx, rely, dt=None ):
154
155    """
156    /* Computes azimuth and slowness of an event from arrival times
157     * at different locations.  Returns also an average deviation
158     * for both results.
159     *
160     * Formula used:
161     *
162     *    (xi,yi)   relative location for n points
163     *    ti        observed relative travel times
164     *
165     *    axx := sum_i{xi*xi};  axy := sum_i{xi*yi};  ayy := sum_i{yi*yi};
166     *    taux := sum_i{ti*xi};  tauy := sum_i{ti*yi};
167     *    wanted:  sx, sy
168     *
169     *      /  axx    axy   \   / sx \   / taux \
170     *      |               | * |    | = |      | ,   or  A*s = tau
171     *      \  axy    ayy   /   \ sy /   \ tauy /
172     *
173     *    ==>   s = (1/A)*tau
174     *
175     *
176     * Important:
177     *                              _
178     *    Only the geometric center x = (1/n)*(sum_i{xi},sum_i{yi}) as origin is
179     *    stable against a constant time shift dt in observed times:
180     *
181     *    ti' := ti + dt   for all i
182     *    then            _                 _
183     *    tau' = tau + dt*x       = tau for x = 0
184     *    For other origins the slowness is warped by time residuals on the
185     *    reference station.
186     *
187     * parameters of routine
188     * int        listlth;     input; number of locations
189     * REAL       time[];      modify; in: relative time, out: scratch
190     * REAL       x[];         modify; in: x-location (km), out: scratch
191     * REAL       y[];         modify; in: y-location (km), out: scratch
192     * REAL       dt;          input; sample distance in seconds
193     * REAL       *azimuth;    output; computed azimuth (in degrees)
194     * REAL       *azimerr;    output; average deviation of azimuth
195     * REAL       *sloness;    output; computed slowness (in sec/degree)
196     * REAL       *slowerr;    output; average deviation of slowness
197     * int        *status;     output; return status
198     */
199    """
200    if len(reltimes) < 3:
201            return None
202    if len(reltimes) == 3 and dt == None:
203        return None
204       
205    RAD_TO_DEG = 180.0/np.pi
206   
207    # coordinate shift to center
208    times = reltimes - reltimes.mean()
209    xrel = relx - relx.mean()
210    yrel = rely - rely.mean()
211
212    tx = ty = tq = 0.0
213    axx = axy = ayx = ayy = 0.0
214
215    if len(times) == 3:
216        tx = times[1]
217        ty = times[2]
218        axx = xrel[1]
219        axy = yrel[1]
220        ayx = xrel[2]
221        ayy = yrel[2]
222    else:
223        for i in range(len(times)):
224            tq += times[i]*times[i]
225            tx += times[i]*xrel[i]
226            ty += times[i]*yrel[i]
227            axx += xrel[i]*xrel[i]
228            axy += xrel[i]*yrel[i]
229            ayy += yrel[i]*yrel[i]
230        ayx = axy
231
232    # compute determinant, Cramer's rule
233    det = axx*ayy - axy*ayx
234    if abs(det) < 1.e-10:
235        return None
236    sx = (tx*ayy - ty*axy) / det
237    sy = (ty*axx - tx*ayx) / det
238
239    if sy == 0.0:
240        if sx < 0.:
241            azimuth = 90.
242        else:
243            azimuth= 270.
244    else:
245        azimuth = np.arctan( sx/sy ) * RAD_TO_DEG
246        if sy > 0.0:
247            azimuth += 180.0
248        elif  (sx > 0.0):
249            azimuth += 360
250
251    norm = np.sqrt( sx*sx + sy*sy )
252    slowness = DEG_TO_KM * norm
253    if norm < 1.e-10:
254        return None
255    norm = 1.0 / norm
256
257    # deviations
258    det = abs(det)
259    if len(times) == 3:
260        #/* dtc = 2.0 * dt; ??? */
261        r = np.arctan( sx/sy )
262        saz = np.sin( r )
263        caz = np.cos( r )
264        r = norm * dt/det
265        scr1 = saz*axx + caz*axy
266        scr2 = saz*ayx + caz*ayy
267        azimerr = r * np.sqrt( scr1*scr1 + scr2*scr2 ) * RAD_TO_DEG
268        scr1 = caz*axx - saz*axy
269        scr2 = caz*ayx - saz*ayy
270        r *= norm * np.sqrt( scr1*scr1 + scr2*scr2 )
271        slowerr = r * DEG_TO_KM / (norm*norm)
272    else:
273        r = np.sqrt( abs(tq-tx*sx-ty*sy) / float(len(times)-3) )
274        sx2 = sx*sx
275        sy2 = sy*sy
276        sxy = sx*sy
277        slowerr = r * np.sqrt( (ayy*sx2 - 2.0*axy*sxy + axx*sy2) / det ) *\
278            norm * DEG_TO_KM
279        azimerr = r * np.sqrt( (ayy*sy2 + 2.0*axy*sxy + axx*sx2) / det ) *\
280            norm * norm * RAD_TO_DEG
281    return (azimuth,azimerr,slowness,slowerr)
282
283
284#-------------------------------------------------------------------------------
285
286
287def planeWaveResiduals( reltimes, relx, rely, takeabs=True ):
288    "Compute residuals of given shift times to a plane wave."
289    # Bit complicated due to possible None values in reltimes.
290    if None in reltimes:
291        rreltimes = []
292        rrelx = []
293        rrely = []
294        for i,relt in enumerate(reltimes):
295            if relt == None:
296                continue
297            rreltimes.append( relt )
298            rrelx.append( relx[i] )
299            rrely.append( rely[i] )
300    else:
301        rreltimes = reltimes
302        rrelx = relx
303        rrely = rely
304    az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes),
305        np.array(rrelx), np.array(rrely) )
306    # Now we have slowness and azimuth, let's compute residuals.
307    reltimes_mean = sum(rreltimes)/float(len(rreltimes))
308    # compute beam shifts (as in command 'beam')
309    radaz = (90.0 - az) / (180./np.pi)
310    kmslo = slo / DEG_TO_KM
311    bshifts = []
312    for i,relt in enumerate(reltimes):
313        bshifts.append(
314            (relx[i]*np.cos(radaz) + rely[i]*np.sin(radaz)) * kmslo
315        )
316    bshifts_mean = sum(bshifts)/float(len(bshifts))
317    resids = []
318    # Residuals consist of sum (difference) of correlation maxima and beam.
319    for i,relt in enumerate(reltimes):
320        if relt == None:
321            resid = None
322        else:
323            resid = (bshifts[i]-bshifts_mean) + (relt-reltimes_mean)
324            if takeabs:
325                resid = abs( resid )
326        resids.append( resid )
327    return np.array(resids)
328
329
330#-------------------------------------------------------------------------------
331
332
333def fit_distance( phase, slo, depth=33., minacc=0.01 ):
334    dist_range = {
335        'P'   : (10.0,97.0),
336    }
337    if phase not in dist_range:
338        return None
339    model = TauPyModel(model="iasp91")
340    lo, hi = dist_range[phase]
341    arr = model.get_travel_times( source_depth_in_km=depth,
342        distance_in_degree=lo, phase_list=[phase] )
343    loslo = arr[0].ray_param_sec_degree
344    arr = model.get_travel_times( source_depth_in_km=depth,
345        distance_in_degree=hi, phase_list=[phase] )
346    hislo = arr[0].ray_param_sec_degree
347    if loslo < hislo:
348        if slo < loslo:
349            return None
350        if slo > hislo:
351            return None
352    else:
353        if slo < hislo:
354            return None
355        if slo > loslo:
356            return None
357    for i in range(50):
358        mid = (lo+hi)/2.
359        arr = model.get_travel_times( source_depth_in_km=depth,
360            distance_in_degree=mid, phase_list=[phase] )
361        midslo = arr[0].ray_param_sec_degree
362        #print loslo, midslo, hislo, '|', lo, hi
363        if loslo < hislo:
364            if slo < midslo:
365                hi = mid
366                hislo = midslo
367            else:
368                lo = mid
369                loslo = midslo
370        else:
371            if slo < midslo:
372                lo = mid
373                loslo = midslo
374            else:
375                hi = mid
376                hislo = midslo
377        acc = midslo-slo
378        if abs(acc) < minacc:
379            return (mid,acc)
380    return (mid,acc)
381
382
383#-------------------------------------------------------------------------------
384
385
386def sphereadd( slat, slon, dist, bazim ):
387
388    """
389    /* computes new location from input location (slat,slon) and distance
390     * "dist" and azimuth "bazim"
391     *
392     * parameters of routine
393     * float      slat, slon;    input; (station) input location in degrees
394     * float      dist;          input; distance in degrees
395     * float      bazim;         input; (back-)azimuth in degrees
396     * float      *elat, *elon;  output; (epicentre) location in degrees
397     */
398    #define DEG2RAD(x) ((x)*BC_PI/180.0)
399    #define RAD2DEG(x) ((x)/BC_PI*180.0)
400    """
401    # transformation to mathematical system
402    slat = 90.0 - slat
403    if  (slon < 0.0):
404        slon = 360 + slon
405
406    invert = (bazim > 180.0)
407    if  invert:
408        bazim = 360 - bazim
409
410    # degrees to radian
411    slat = slat*np.pi/180.
412    dist = dist*np.pi/180.
413    bazim = bazim*np.pi/180.
414
415    arg = np.cos(dist)*np.cos(slat) + np.sin(dist)*np.sin(slat)*np.cos(bazim)
416    if arg > 1.0:
417        arg = 1.0
418    if arg < -1.0:
419        arg = -1.0
420    elat = np.arccos( arg );
421    arg = (np.cos(dist)-np.cos(elat)*np.cos(slat))/(np.sin(elat)*np.sin(slat))
422    if arg > 1.0:
423        arg = 1.0
424    if arg < -1.0:
425        arg = -1.0
426    gamma = np.arccos( arg )
427
428    # back to degrees
429    elat = elat*180./np.pi
430    gamma = gamma*180./np.pi
431
432    if invert:
433        gamma = 360.0 - gamma
434    elon = slon + gamma
435    if elon > 360.0:
436        elon -= 360.0
437
438    # transformation back to geo-system
439    elat = 90.0 - elat
440    if elon > 180.0:
441        elon -= 360.0
442   
443    return (elat,elon)
444
445
446#-------------------------------------------------------------------------------
447
448
449def compareLocation( origin, lat, lon, depth, cmplist,
450    centerlat=None, centerlon=None ):
451
452    """
453    Compare location with results from other sources using FDSN webservice.
454    """
455
456    toltime = 120.
457    evpar = {}
458    evpar['latitude'] = lat
459    evpar['longitude'] = lon
460    evpar['maxradius'] = 5.
461    evpar['minmagnitude'] = 4.0
462    evpar['starttime'] = origin - toltime
463    evpar['endtime'] = origin + toltime
464    evpar['includeallorigins'] = False
465    evpar['includeallmagnitudes'] = False
466    resultstr = ""
467    for src in cmplist.split(','):
468        client = Client( src )
469        evlist = client.get_events( **evpar )
470        for ev in evlist:
471            if len(ev.origins) < 0:
472                continue
473            orig = ev.origins[0]
474            magn = ev.magnitudes[0]
475            origdiff = origin - orig.time
476            depdiff = depth - orig.depth/1000.
477            dist, az, baz = gps2DistAzimuth( lat, lon, orig.latitude,
478                orig.longitude )
479            kmdist = dist / 1000.
480            degdist = kmdist / DEG_TO_KM
481            resultstr += "%s: locdiff %3.1f km or " % (src,kmdist)\
482                +"%4.2f deg, timediff %3.1f s, " % (degdist,origdiff)\
483                +"depdiff %3.1f\n" % depdiff
484            if centerlat == None or centerlon == None:
485                continue
486            dist, az, baz = gps2DistAzimuth( centerlat, centerlon,
487                orig.latitude, orig.longitude )
488            evdist = dist / (1000.*DEG_TO_KM)
489            pslo = getPhaseSlownesses()
490            model = TauPyModel(model="iasp91")
491            arr = model.get_travel_times( source_depth_in_km=depth,
492                distance_in_degree=evdist, phase_list=pslo.keys() )
493            for phase in pslo.keys():
494                ownslowness, ownaz, sloerr, azerr = pslo[phase]
495                print "dbg:", phase, ownslowness
496                cmpval = None
497                for a in arr:
498                    print "dbg: arr dump", a.phase.name, a.ray_param_sec_degree
499                    if a.phase.name == phase:
500                        cmpval = a.ray_param_sec_degree
501                        break
502                resultstr += "   slowness %s %3.1f" % (phase,ownslowness)
503                if cmpval != None:
504                    print "dbg: slownesses", phase, ownslowness, cmpval
505                    resultstr += ", diff %3.1f" % (ownslowness-cmpval)
506                resultstr += '\n'
507    print "dbg:", resultstr
508    return resultstr
509
510
511#-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.