Zemětřesné eventy a ohniskové koule


In [1]:
%pylab inline
from __future__ import print_function
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
import warnings
warnings.filterwarnings("ignore")

Populating the interactive namespace from numpy and matplotlib
In [3]:
from obspy.core import read, UTCDateTime
from obspy import read_inventory, read_events

Funkce read_events umožňuje načíst eventová metadata uložená v různých formátech do jednotného objektu Catalog.


Podporované formáty: CMTSOLUTION, FNETMT, GSE2, IMS10BULLETIN, MCHEDR, NDK, NLLOC_HYP, NORDIC, QUAKEML, SC3ML, SCARDEC, ZMAP.

Načteme soubor easiaa_p_pick.xml obsahující seznam eventů určených pro pickování teleseismických P vln na stanicích projektu AlpArray-EASI a jeho nejbližšího okolí. Jako zdroj metadat byly použity ISC a NEIC katalogy.

In [4]:
cat = read_events('easiaa_p_pick.xml')
print(cat)

2004 Event(s) in Catalog:
2014-06-16T12:01:08.720000Z | +67.695, -162.146 | 5.7 MW
2014-06-16T13:26:46.510000Z | +11.746,  -81.349 | 5.3 MW
...
2017-10-28T19:10:55.240000Z | +86.891,  +53.301 | 4.4 mb
2017-10-28T19:11:02.760000Z | +86.923,  +55.200 | 5.9 Mww
To see all events call 'print(CatalogObject.__str__(print_all=True))'

Catalog

Objekt Catalog má předdefinováno několik funkcí, základní z nich je plot pro vykreslení zemětřesných lokací.

In [5]:
plt.rcParams['figure.figsize'] = 12, 8
cat.plot(),
Out[5]:
(<matplotlib.figure.Figure at 0xcc92e48>,)

Vybereme jen eventy s magnitudem minimálně 7 a barevně rozlišíme dobu vzniku zemětřesení.

In [6]:
cat.filter('magnitude >= 7').plot(color="date"),
Out[6]:
(<matplotlib.figure.Figure at 0x12b72710>,)

Do obrázku lze začlenit i informaci o poloze stanic, použijeme STATIONXML soubor vytvořený ze stanic ALpArray-EASI. Projekce tentokrát lokální.

In [7]:
cat_eu = cat.filter("latitude >= 30.0", "latitude <= 60.0", "longitude >= -15.0", "longitude <= 30.0")
fig = cat_eu.plot(projection="local", show=False)
inv = read_inventory('XT_all.xml')
inv.plot(fig=fig, label=False, size=50, color_per_network={"XT": "red"}),
Out[7]:
(<matplotlib.figure.Figure at 0x15210128>,)

Zajímá nás azimutální rozložení eventů. Vybereme prostřední stanici z EASI profilu a napočteme back-azimuty a epicentrální vzdálenosti pro všechny eventy. K tomu poslouží funkce calc_vincenty_inverse s eliptickou aproximací tvaru Země.

In [8]:
from obspy.geodetics.base import calc_vincenty_inverse, kilometer2degrees

refstat = inv.get_coordinates(seed_id="XT.AAE28.00.HHZ")  # souřednice z prostřední stanice v profilu

dists, bazs = [], []
for event in cat:
    origin = event.preferred_origin() or event.origins[0]  # finta na výběr nejlépe určeného ohniska
    [dist_m, az, baz] = calc_vincenty_inverse(origin.latitude, origin.longitude, 
                                              refstat['latitude'], refstat['longitude'])
    dist = kilometer2degrees(dist_m/1000.)
    dists.append(dist); bazs.append(baz)

# vykresli histogramy pro BAZ a EP. DIST
fig = plt.figure(figsize=(14.5,4))
ax1 = fig.add_subplot(1,2,1); ax2 = fig.add_subplot(1,2,2)
ax1.hist(bazs, 36, range=(0,360), facecolor='red', alpha=0.75)
ax2.hist(dists, 30, facecolor='green', alpha=0.75)
ax1.set_xlim(0,360); ax1.set_xticks(range(0,361, 60)); ax2.set_xlim(0,max(dists))
ax1.set_xlabel('back-azimuth [$^\circ$]'); ax2.set_xlabel('epic. distance [$^\circ$]')
Out[8]:
Text(0.5,0,u'epic. distance [$^\\circ$]')

Beachballs / ohniskové koule

  • beachball vykreslí ohniskovou kouli
  • beach to samé, navíc ale s umístěním do mapy

Momentový tenzor lze definovat dvěma způsoby:

In [9]:
from obspy.imaging.beachball import beachball

mt1 = [150, 87, 30]    # strike, dip, and rake
mt2 = [0.91, -0.89, -0.02, 1.78, -1.55, 0.47]   # MOMENT TENSOR:
mt3 = [-2.39, 1.04, 1.35, 0.57, -2.94, -0.94]   # M11, M22, M33, M12, M13, M23

beachball(mt1, size=200, linewidth=2, facecolor='b')
beachball(mt2, size=200, linewidth=2, facecolor='r')
beachball(mt3, size=200, linewidth=2, facecolor='g'),
Out[9]:
(<matplotlib.figure.Figure at 0x17cacdd8>,)

Zatím v ObsPy neexistuje jednoduchý způsob, jak dávkově načíst eventy i s ohniskovými mechanismy. Nejjednoduší řešení je stáhnout si např. z https://www.globalcmt.org celý Harvardský katalog (soubor http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/jan76_dec17.ndk.gz pro období 1976 až 2017) a pomocí jednoduchého skriptu např. v shellu z něj vyfiltrovat požadovanou část.

Funkce read_events umí načíst eventová metadata i ve formátu ndk.

In [10]:
cat1 = read_events('GCMT_ALL_30N_50N_000E_030E.ndk')
cat1.plot(projection="local"),
Out[10]:
(<matplotlib.figure.Figure at 0x150eb5f8>,)

Pro hezčí vykreslení mapy s ohniskovými koulemi jsem převzal a upravil původní skript z https://github.com/qingkaikong/blog/tree/master/15_plot_historical_moment_tensors .

In [11]:
from obspy.imaging.beachball import beach
from mpl_toolkits.basemap import Basemap
from matplotlib.lines import Line2D

def plot_mt(cat, latlons, figsize = (16,24), mt_size = 10, mt_aspect = None,\
                 pretty = False, resolution='l'):
    '''
    Plot the historical moment tensor from the query of GCMT
    
    cat - catalog with focal mechanisms (e.g. from GCMT)
    latlons - list or tuple form a map domain: lat1, lat2, long1, long2
    figsize - tuple of the size of the figure
    mt_size - size of the moment tensor on the plot for magnitude 5 (or all, see mt_size_large)
    mt_aspect - magnification of large events, magn=8: size = mt_aspect * mt_size
                    if None, all events have the same size mt_size
    pretty - boolean, whether want to plot nice maps
    resolution - low or high as you want
    '''

    llat = float(latlons[0])
    ulat = float(latlons[1])
    llon = float(latlons[2])
    ulon = float(latlons[3])
    
    lats, lons, depths, mts, mags, mt_sizes = [],[],[],[],[],[]
    for event in cat:
        
        # jak z objektu Catalog dostat lat, long, depth, magnitude, moment tensor
        origin = event.preferred_origin() or event.origins[0]
        magnitude = event.preferred_magnitude() or event.magnitudes[0]
        focmec = event.preferred_focal_mechanism() or event.focal_mechanisms[0]
        tensor = focmec.moment_tensor.tensor
        
        lats.append(origin.latitude)
        lons.append(origin.longitude)
        mts.append([tensor.m_rr, tensor.m_tt, tensor.m_pp, tensor.m_rt, tensor.m_rp, tensor.m_tp])
        depths.append(origin.depth/1000.)
        mags.append(magnitude.mag)
        
        if mt_aspect:
            # velikost koule úměrná magnitudu
            mtsize = max([mt_size * ((mt_aspect-1.)*(magnitude.mag-5.)/3.+1.), mt_size/2.])
            mt_sizes.append(mtsize)
        else:
            mt_sizes.append(mt_size)
    
    plt.figure(figsize=figsize)
    m = Basemap(projection='cyl', resolution=resolution,
        llcrnrlon=llon,llcrnrlat=llat,urcrnrlon=ulon,urcrnrlat=ulat)
    
    m.drawcoastlines()
    m.drawcountries()
    
    if pretty:    
        m.etopo()   # ETOPO podklad
    else:
        m.drawmapboundary()
        m.fillcontinents()
    
    m.drawparallels(np.arange(llat, ulat, (ulat - llat) / 4.0), labels=[1,0,0,0])
    m.drawmeridians(np.arange(llon, ulon, (ulon - llon) / 4.0), labels=[0,0,0,1])   
    
    ax = plt.gca()
    x, y = m(lons, lats)
    
    for i in range(len(mts)):
        
        # barevné rozlišení podle hloubky
        if depths[i] <= 50:
            color = '#FFA500'
        elif depths[i] > 50 and depths[i] <= 100:
            color = 'g'
        elif depths[i] > 100 and depths[i] <= 200:
            color = 'b'
        else:
            color = 'r'
        
        #note here, if the mrr is zero, then you will have an error
        #so, change this to a very small number 
        if mts[i][0] == 0:
            mts[i][0] = 0.001
        
        try:
            b = beach(mts[i], xy=(x[i], y[i]),width=mt_sizes[i], linewidth=1, facecolor=color)
        except:
            pass
            
        b.set_zorder(10)
        ax.add_collection(b)
        
    # add the legend
    circ1 = Line2D([0], [0], linestyle="none", \
            marker="o", alpha=0.6, markersize=10, markerfacecolor="#FFA500")
    circ2 = Line2D([0], [0], linestyle="none", \
            marker="o", alpha=0.6, markersize=10, markerfacecolor="g")
    circ3 = Line2D([0], [0], linestyle="none", \
            marker="o", alpha=0.6, markersize=10, markerfacecolor="b")
    circ4 = Line2D([0], [0], linestyle="none", \
            marker="o", alpha=0.6, markersize=10, markerfacecolor="r")
    plt.legend((circ1, circ2, circ3, circ4), \
                ("depth $\leq$ 50 km", "50 km $<$ depth $\leq$ 100 km", \
                "100 km $<$ depth $\leq$ 200 km", "200 km $<$ depth"), \
                numpoints=1, loc=3)
    plt.show()
    return
In [12]:
plot_mt(cat1, latlons=[30,50,0,30], mt_size=0.5, mt_aspect=2, pretty=True),
Out[12]:
(None,)
In [14]:
cat2 = read_events('GCMT_ALL_10S_10N_090W_070W.ndk')  # výřez pro SZ Jižní Ameriky
plot_mt(cat2, latlons=[-10,10,-90,-70], figsize = (11,11), mt_size=0.4, mt_aspect=2, pretty=True),
Out[14]:
(None,)
In [ ]:
# výřez jen pro AlpArray
#plot_mt(cat1, latlons=[33,42,18,30], mt_size=0.25, mt_aspect=2, pretty=True),