import numpy as np import matplotlib.pyplot as plt from obspy import read # 读取面波数据并画图。 st = read('MASW_DATA/Sample_Data/*.SAC') dt = st[0].stats.delta data = [] scale = 0.05 dx = 2 plt.figure(figsize=(8, 6)) for i, tr in enumerate(st): d = tr.data data.append(d) t = np.arange(len(d)) * dt plt.plot(t, d*scale+(i+1)*dx, lw=1, color='b') plt.xlabel('Time (s)') plt.ylabel('Offset (m)') plt.tight_layout() plt.savefig('Surface_wave.png') plt.show() # 二维FFT。 d = np.array(data) n = len(d[0]) # m为空间方向的采样点数,m增大可以让FK谱光滑一点,以达到插值效果。 m = len(d[:, 0]) * 5 D = np.zeros((m, n)) D[:len(d[:, 0])] = d # 时间采样率。 fs = 1 / dt # 空间采样率 xs = 1 / dx # 频率 (赫兹)。 f = np.arange(-n//2, n//2) * fs / (n-1) # 波数 (每米)。 k = 2 * np.pi * np.arange(-m//2, m//2) * xs / (m-1) # 二维FFT。 fk = np.fft.fft2(D) # 作图。 pmin = -10 P = abs(np.fft.fftshift(fk)); P /= P.max(); P = 10 * np.log10(P) P2 = abs(fk); P2 /= P2.max(); P2 = 10 * np.log10(P2) plt.figure(figsize=(11, 8)) plt.subplot(221) plt.pcolormesh(f, k, P2, cmap='magma', vmin=pmin, vmax=0) plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (2$\pi$m$^{-1}$)') plt.subplot(222) plt.pcolormesh(f, k, P, cmap='magma', vmin=pmin, vmax=0) plt.plot([f[n//2], f[-1], f[-1], f[n//2], f[n//2]], [k[0], k[0], k[m//2], k[m//2], k[0]], lw=2, ls='--', color='r') plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (m$^{-1}$)') plt.subplot(223) plt.pcolormesh(f[n//2:], k[:m//2], P[:m//2, n//2:], cmap='magma', vmin=pmin, vmax=0) plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (m$^{-1}$)') plt.subplot(224) plt.pcolormesh(f[n//2:], abs(k[:m//2][::-1]), P[:m//2, n//2:][::-1], cmap='magma', vmin=pmin, vmax=0) cbar = plt.colorbar() cbar.set_label(r'FK spectra (dB)') plt.xlim(0, 100) plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (m$^{-1}$)') plt.tight_layout() plt.show()
时间: 2024-03-28 11:40:42 浏览: 17
这段代码主要实现了面波数据的读取和处理,包括画出原始数据的波形图和进行二维FFT后绘制FK谱图。具体步骤如下:
1. 读取面波数据:使用obspy中的read函数读取指定路径下的所有.SAC文件,并将各个通道的数据保存到data列表中。
2. 画出原始数据的波形图:将每个通道的数据与对应的时间轴绘制在同一张图上,用不同的颜色区分不同通道。
3. 进行二维FFT:将原始数据进行二维FFT,得到FK谱,用于分析面波的频率和波数分布情况。
4. 绘制FK谱图:将FK谱进行可视化,绘制四个子图。其中,第一个子图展示了整个FK谱的情况,第二个子图在第一个子图的基础上加入了一个红色方框,表示主要分析区域。第三个子图展示了主要分析区域的上半部分,第四个子图展示了主要分析区域的下半部分。
总的来说,这段代码实现了对面波数据的初步处理和分析,为后续的面波分析提供了基础。
相关问题
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d
import numpy是一个用于科学计算的Python库,它提供了大量用于数组操作和数学计算的函数和工具。使用import numpy as np将numpy库导入为np别名,以便在代码中更方便地使用numpy的函数和工具。
而import matplotlib.pyplot as plt是导入matplotlib库的一种常用方式,matplotlib是一个用于绘制图表和数据可视化的库,import语句将matplotlib.pyplot模块导入为plt别名,以便在代码中更方便地使用matplotlib库的绘图函数和工具。
最后,from mpl_toolkits.mplot3d导入模块同样是用于matplotlib库的一种常见操作,它用于导入mpl_toolkits.mplot3d模块,这是一个用于绘制三维图表和可视化三维数据的工具模块。
综合起来,这些import语句在Python代码中的作用是将numpy和matplotlib库以及其相关的模块导入到代码中,以便使用它们提供的各种函数和工具进行数组操作、数学计算和绘制图表等操作。这样可以极大地方便数据科学家、工程师和研究人员们对数据的处理和可视化工作。
import numpy as np import random import matplotlib.pyplot as plt是否有问题
在Python中,通常每个导入语句都应该单独成行,因此在`import numpy as np import random import matplotlib.pyplot as plt`这行代码中,将三个导入语句放在同一行可能会导致语法错误。正确的写法应该是每个导入语句单独成行,例如:
```
import numpy as np
import random
import matplotlib.pyplot as plt
```
这样写可以保证代码的可读性和易于维护。