GeophyAI

Python与地球物理数据处理

0%

地震数据波形变密度绘图

发表于 分类于 Valine：
地震数据的波形变密度显示。

地震数据变密度显示

以下代码节选自fatiando包的vis.mpl.seismic_wiggle模块

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
from matplotlib import pyplot
import numpy
def seismic_wiggle(section, dt=0.004, ranges=None, scale=1.,
                   color='k', normalize=False):
    """
    Plot a seismic section (numpy 2D array matrix) as wiggles.

    Parameters:

    * section :  二维数组，第一个维度为时间，第二个维度为导数
        matrix of traces (first dimension time, second dimension traces)
    * dt : float
        时间采样率（默认四毫秒）
    * ranges : (x1, x2)
        水平最大最小值（默认为道号）
    * scale : float
        比例因子
    * color : tuple of strings
        Color for filling the wiggle, positive  and negative lobes.
    * normalize :
        True: 使用全局最大和最小值将所有道的数据归一化到[-0.5, 0.5]
    .. warning::
        Slow for more than 200 traces, in this case decimate your
        data or use ``seismic_image``.

    """
    npts, ntraces = section.shape  # time/traces
    if ntraces < 1:
        raise IndexError("Nothing to plot")
    if npts < 1:
        raise IndexError("Nothing to plot")
    t = numpy.linspace(0, dt*npts, npts)
    amp = 1.  # normalization factor
    gmin = 0.  # global minimum
    toffset = 0.  # offset in time to make 0 centered
    if normalize:
        gmax = section.max()
        gmin = section.min()
        amp = (gmax-gmin)
        toffset = 0.5
    pyplot.ylim(max(t), 0)
    if ranges is None:
        ranges = (0, ntraces)
    x0, x1 = ranges
    # horizontal increment
    dx = float((x1-x0)/ntraces)
    pyplot.xlim(x0, x1)
    for i, trace in enumerate(section.transpose()):
        tr = (((trace-gmin)/amp)-toffset)*scale*dx
        x = x0+i*dx  # x positon for this trace
        pyplot.plot(x+tr, t, 'k')
        pyplot.fill_betweenx(t, x+tr, x, tr > 0, color=color)