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

Revision 1168, 19.2 KB checked in by klaus, 4 years ago (diff)

fit depth using depth phases

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    print "dbg: depthphases", determineDepthFromDepthPhases( elat, elon )
76
77    return (origin,elat,elon,center_lat,center_lon)
78
79
80#-------------------------------------------------------------------------------
81
82
83def getPhaseSlownesses( withcenter=False, withresiduals=False ):
84    "Compute slownesses for all phases in phase list (>4 picks)."
85
86    pslo = {}
87    presid = {}
88    phaselist = PhaseList()
89    plist = phaselist.getPhaseListsByName()
90    stations = Stations()
91    for phase in plist:
92        lons = []
93        lats = []
94        abstimes = []
95        stalist = []
96        for station,picktime in plist[phase]:
97            try:
98                r = stations[(station,None)]
99                lat = float( r.latitude )
100                lon = float( r.longitude )
101            except:
102                continue
103            lons.append( lon )
104            lats.append( lat )
105            abstimes.append( picktime )
106            stalist.append( station )
107        if len(lons) < 3:
108            continue
109        relx, rely = absToRelPositions( lons, lats )
110        mintime = min( abstimes )
111        times = []
112        for a in abstimes:
113            times.append( a-mintime )
114        azimuth,azimerr,slowness,slowerr = computeSlownessFromPicks(
115            np.array(times), np.array(relx), np.array(rely) )
116        if withcenter:
117            center_lon = sum(lons)/float(len(lons))
118            center_lat = sum(lats)/float(len(lats))
119            pslo[phase] = (slowness,azimuth,slowerr,azimerr,
120                center_lat,center_lon)
121        else:
122            pslo[phase] = (slowness,azimuth,slowerr,azimerr)
123        if withresiduals:
124            resid = planeWaveResiduals( times, relx, rely, takeabs=False )
125            presid[phase] = []
126            for i,iresid in enumerate(resid):
127                presid[phase].append( (stalist[i],iresid) )
128    if withresiduals:
129        return (pslo,presid)
130    else:
131        return pslo
132
133
134#-------------------------------------------------------------------------------
135
136
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
190def absToRelPositions( lons, lats ):
191    "Compute relative distances in km from lon,lat pairs."
192    if len(lons) == 0 or len(lons) != len(lats):
193        return None
194    center_lon = sum(lons)/float(len(lons))
195    center_lat = sum(lats)/float(len(lats))
196    relx = []
197    rely = []
198    for i in range(len(lons)):
199        x, y = utlGeoKm( center_lon, center_lat, lons[i], lats[i] )
200        relx.append( x )
201        rely.append( y )
202    return (relx,rely)
203
204
205#-------------------------------------------------------------------------------
206
207
208def computeSlownessFromPicks( reltimes, relx, rely, dt=None ):
209
210    """
211    /* Computes azimuth and slowness of an event from arrival times
212     * at different locations.  Returns also an average deviation
213     * for both results.
214     *
215     * Formula used:
216     *
217     *    (xi,yi)   relative location for n points
218     *    ti        observed relative travel times
219     *
220     *    axx := sum_i{xi*xi};  axy := sum_i{xi*yi};  ayy := sum_i{yi*yi};
221     *    taux := sum_i{ti*xi};  tauy := sum_i{ti*yi};
222     *    wanted:  sx, sy
223     *
224     *      /  axx    axy   \   / sx \   / taux \
225     *      |               | * |    | = |      | ,   or  A*s = tau
226     *      \  axy    ayy   /   \ sy /   \ tauy /
227     *
228     *    ==>   s = (1/A)*tau
229     *
230     *
231     * Important:
232     *                              _
233     *    Only the geometric center x = (1/n)*(sum_i{xi},sum_i{yi}) as origin is
234     *    stable against a constant time shift dt in observed times:
235     *
236     *    ti' := ti + dt   for all i
237     *    then            _                 _
238     *    tau' = tau + dt*x       = tau for x = 0
239     *    For other origins the slowness is warped by time residuals on the
240     *    reference station.
241     *
242     * parameters of routine
243     * int        listlth;     input; number of locations
244     * REAL       time[];      modify; in: relative time, out: scratch
245     * REAL       x[];         modify; in: x-location (km), out: scratch
246     * REAL       y[];         modify; in: y-location (km), out: scratch
247     * REAL       dt;          input; sample distance in seconds
248     * REAL       *azimuth;    output; computed azimuth (in degrees)
249     * REAL       *azimerr;    output; average deviation of azimuth
250     * REAL       *sloness;    output; computed slowness (in sec/degree)
251     * REAL       *slowerr;    output; average deviation of slowness
252     * int        *status;     output; return status
253     */
254    """
255    if len(reltimes) < 3:
256            return None
257    if len(reltimes) == 3 and dt == None:
258        return None
259       
260    RAD_TO_DEG = 180.0/np.pi
261   
262    # coordinate shift to center
263    times = reltimes - reltimes.mean()
264    xrel = relx - relx.mean()
265    yrel = rely - rely.mean()
266
267    tx = ty = tq = 0.0
268    axx = axy = ayx = ayy = 0.0
269
270    if len(times) == 3:
271        tx = times[1]
272        ty = times[2]
273        axx = xrel[1]
274        axy = yrel[1]
275        ayx = xrel[2]
276        ayy = yrel[2]
277    else:
278        for i in range(len(times)):
279            tq += times[i]*times[i]
280            tx += times[i]*xrel[i]
281            ty += times[i]*yrel[i]
282            axx += xrel[i]*xrel[i]
283            axy += xrel[i]*yrel[i]
284            ayy += yrel[i]*yrel[i]
285        ayx = axy
286
287    # compute determinant, Cramer's rule
288    det = axx*ayy - axy*ayx
289    if abs(det) < 1.e-10:
290        return None
291    sx = (tx*ayy - ty*axy) / det
292    sy = (ty*axx - tx*ayx) / det
293
294    if sy == 0.0:
295        if sx < 0.:
296            azimuth = 90.
297        else:
298            azimuth= 270.
299    else:
300        azimuth = np.arctan( sx/sy ) * RAD_TO_DEG
301        if sy > 0.0:
302            azimuth += 180.0
303        elif  (sx > 0.0):
304            azimuth += 360
305
306    norm = np.sqrt( sx*sx + sy*sy )
307    slowness = DEG_TO_KM * norm
308    if norm < 1.e-10:
309        return None
310    norm = 1.0 / norm
311
312    # deviations
313    det = abs(det)
314    if len(times) == 3:
315        #/* dtc = 2.0 * dt; ??? */
316        r = np.arctan( sx/sy )
317        saz = np.sin( r )
318        caz = np.cos( r )
319        r = norm * dt/det
320        scr1 = saz*axx + caz*axy
321        scr2 = saz*ayx + caz*ayy
322        azimerr = r * np.sqrt( scr1*scr1 + scr2*scr2 ) * RAD_TO_DEG
323        scr1 = caz*axx - saz*axy
324        scr2 = caz*ayx - saz*ayy
325        r *= norm * np.sqrt( scr1*scr1 + scr2*scr2 )
326        slowerr = r * DEG_TO_KM / (norm*norm)
327    else:
328        r = np.sqrt( abs(tq-tx*sx-ty*sy) / float(len(times)-3) )
329        sx2 = sx*sx
330        sy2 = sy*sy
331        sxy = sx*sy
332        slowerr = r * np.sqrt( (ayy*sx2 - 2.0*axy*sxy + axx*sy2) / det ) *\
333            norm * DEG_TO_KM
334        azimerr = r * np.sqrt( (ayy*sy2 + 2.0*axy*sxy + axx*sx2) / det ) *\
335            norm * norm * RAD_TO_DEG
336    return (azimuth,azimerr,slowness,slowerr)
337
338
339#-------------------------------------------------------------------------------
340
341
342def planeWaveResiduals( reltimes, relx, rely, takeabs=True ):
343    "Compute residuals of given shift times to a plane wave."
344    # Bit complicated due to possible None values in reltimes.
345    if None in reltimes:
346        rreltimes = []
347        rrelx = []
348        rrely = []
349        for i,relt in enumerate(reltimes):
350            if relt == None:
351                continue
352            rreltimes.append( relt )
353            rrelx.append( relx[i] )
354            rrely.append( rely[i] )
355    else:
356        rreltimes = reltimes
357        rrelx = relx
358        rrely = rely
359    az, azerr, slo, sloerr = computeSlownessFromPicks( np.array(rreltimes),
360        np.array(rrelx), np.array(rrely) )
361    # Now we have slowness and azimuth, let's compute residuals.
362    reltimes_mean = sum(rreltimes)/float(len(rreltimes))
363    # compute beam shifts (as in command 'beam')
364    radaz = (90.0 - az) / (180./np.pi)
365    kmslo = slo / DEG_TO_KM
366    bshifts = []
367    for i,relt in enumerate(reltimes):
368        bshifts.append(
369            (relx[i]*np.cos(radaz) + rely[i]*np.sin(radaz)) * kmslo
370        )
371    bshifts_mean = sum(bshifts)/float(len(bshifts))
372    resids = []
373    # Residuals consist of sum (difference) of correlation maxima and beam.
374    for i,relt in enumerate(reltimes):
375        if relt == None:
376            resid = None
377        else:
378            resid = (bshifts[i]-bshifts_mean) + (relt-reltimes_mean)
379            if takeabs:
380                resid = abs( resid )
381        resids.append( resid )
382    return np.array(resids)
383
384
385#-------------------------------------------------------------------------------
386
387
388def fit_distance( phase, slo, depth=33., minacc=0.01 ):
389    dist_range = {
390        'P'   : (10.0,97.0),
391    }
392    if phase not in dist_range:
393        return None
394    model = TauPyModel(model="iasp91")
395    lo, hi = dist_range[phase]
396    arr = model.get_travel_times( source_depth_in_km=depth,
397        distance_in_degree=lo, phase_list=[phase] )
398    loslo = arr[0].ray_param_sec_degree
399    arr = model.get_travel_times( source_depth_in_km=depth,
400        distance_in_degree=hi, phase_list=[phase] )
401    hislo = arr[0].ray_param_sec_degree
402    if loslo < hislo:
403        if slo < loslo:
404            return None
405        if slo > hislo:
406            return None
407    else:
408        if slo < hislo:
409            return None
410        if slo > loslo:
411            return None
412    for i in range(50):
413        mid = (lo+hi)/2.
414        arr = model.get_travel_times( source_depth_in_km=depth,
415            distance_in_degree=mid, phase_list=[phase] )
416        midslo = arr[0].ray_param_sec_degree
417        #print loslo, midslo, hislo, '|', lo, hi
418        if loslo < hislo:
419            if slo < midslo:
420                hi = mid
421                hislo = midslo
422            else:
423                lo = mid
424                loslo = midslo
425        else:
426            if slo < midslo:
427                lo = mid
428                loslo = midslo
429            else:
430                hi = mid
431                hislo = midslo
432        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
470        if abs(acc) < minacc:
471            return (mid,acc)
472    return (mid,acc)
473
474
475#-------------------------------------------------------------------------------
476
477
478def sphereadd( slat, slon, dist, bazim ):
479
480    """
481    /* computes new location from input location (slat,slon) and distance
482     * "dist" and azimuth "bazim"
483     *
484     * parameters of routine
485     * float      slat, slon;    input; (station) input location in degrees
486     * float      dist;          input; distance in degrees
487     * float      bazim;         input; (back-)azimuth in degrees
488     * float      *elat, *elon;  output; (epicentre) location in degrees
489     */
490    #define DEG2RAD(x) ((x)*BC_PI/180.0)
491    #define RAD2DEG(x) ((x)/BC_PI*180.0)
492    """
493    # transformation to mathematical system
494    slat = 90.0 - slat
495    if  (slon < 0.0):
496        slon = 360 + slon
497
498    invert = (bazim > 180.0)
499    if  invert:
500        bazim = 360 - bazim
501
502    # degrees to radian
503    slat = slat*np.pi/180.
504    dist = dist*np.pi/180.
505    bazim = bazim*np.pi/180.
506
507    arg = np.cos(dist)*np.cos(slat) + np.sin(dist)*np.sin(slat)*np.cos(bazim)
508    if arg > 1.0:
509        arg = 1.0
510    if arg < -1.0:
511        arg = -1.0
512    elat = np.arccos( arg );
513    arg = (np.cos(dist)-np.cos(elat)*np.cos(slat))/(np.sin(elat)*np.sin(slat))
514    if arg > 1.0:
515        arg = 1.0
516    if arg < -1.0:
517        arg = -1.0
518    gamma = np.arccos( arg )
519
520    # back to degrees
521    elat = elat*180./np.pi
522    gamma = gamma*180./np.pi
523
524    if invert:
525        gamma = 360.0 - gamma
526    elon = slon + gamma
527    if elon > 360.0:
528        elon -= 360.0
529
530    # transformation back to geo-system
531    elat = 90.0 - elat
532    if elon > 180.0:
533        elon -= 360.0
534   
535    return (elat,elon)
536
537
538#-------------------------------------------------------------------------------
539
540
541def compareLocation( origin, lat, lon, depth, cmplist,
542    centerlat=None, centerlon=None ):
543
544    """
545    Compare location with results from other sources using FDSN webservice.
546    """
547
548    toltime = 120.
549    evpar = {}
550    evpar['latitude'] = lat
551    evpar['longitude'] = lon
552    evpar['maxradius'] = 5.
553    evpar['minmagnitude'] = 4.0
554    evpar['starttime'] = origin - toltime
555    evpar['endtime'] = origin + toltime
556    evpar['includeallorigins'] = False
557    evpar['includeallmagnitudes'] = False
558    resultstr = ""
559    for src in cmplist.split(','):
560        client = Client( src )
561        evlist = client.get_events( **evpar )
562        for ev in evlist:
563            if len(ev.origins) < 0:
564                continue
565            orig = ev.origins[0]
566            magn = ev.magnitudes[0]
567            origdiff = origin - orig.time
568            depdiff = depth - orig.depth/1000.
569            dist, az, baz = gps2DistAzimuth( lat, lon, orig.latitude,
570                orig.longitude )
571            kmdist = dist / 1000.
572            degdist = kmdist / DEG_TO_KM
573            resultstr += "%s: locdiff %3.1f km or " % (src,kmdist)\
574                +"%4.2f deg, timediff %3.1f s, " % (degdist,origdiff)\
575                +"depdiff %3.1f\n" % depdiff
576            if centerlat == None or centerlon == None:
577                continue
578            dist, az, baz = gps2DistAzimuth( centerlat, centerlon,
579                orig.latitude, orig.longitude )
580            evdist = dist / (1000.*DEG_TO_KM)
581            pslo = getPhaseSlownesses()
582            model = TauPyModel(model="iasp91")
583            arr = model.get_travel_times( source_depth_in_km=depth,
584                distance_in_degree=evdist, phase_list=pslo.keys() )
585            for phase in pslo.keys():
586                ownslowness, ownaz, sloerr, azerr = pslo[phase]
587                print "dbg:", phase, ownslowness
588                cmpval = None
589                for a in arr:
590                    print "dbg: arr dump", a.phase.name, a.ray_param_sec_degree
591                    if a.phase.name == phase:
592                        cmpval = a.ray_param_sec_degree
593                        break
594                resultstr += "   slowness %s %3.1f" % (phase,ownslowness)
595                if cmpval != None:
596                    print "dbg: slownesses", phase, ownslowness, cmpval
597                    resultstr += ", diff %3.1f" % (ownslowness-cmpval)
598                resultstr += '\n'
599    print "dbg:", resultstr
600    return resultstr
601
602
603#-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.