#plot 2 harmonic lines
import os
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from obspy import read
g = [14]
# 定义bin的角度间隔和范围
bin_angle_interval = 5
bin_angle_range = np.arange(0, 360, bin_angle_interval)
csv_rfani = "/media/zekun/Software/fenglin_niu/rfani.csv"
ani = pd.read_csv(csv_rfani)
csv_rfcan = "/media/zekun/Software/pictures/anis_yn_final/rfcan.csv"
rfcan = pd.read_csv(csv_rfcan)
# 创建绘图对象
#%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
# 创建一个包含两个子图的图形窗口 (1行2列布局)
fig, axs = plt.subplots(1, 2)
# 创建绘图对象
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
fig, axs[0] = plt.subplots(figsize=(9, 12))
fig, axs[1] = plt.subplots(figsize=(9, 12))
for i in g:
if len(str(i))==1:
n='00'+str(i)
else:
n='0'+str(i)
# 创建存储叠加波形和事件数目的字典
stacked_waveforms = {angle: np.zeros(10004) for angle in bin_angle_range}
event_counts = {angle: 0 for angle in bin_angle_range}
# SAC文件所在的文件夹路径
sac_folder = '/media/zekun/Software/fenglin_niu/RFcan/examples/YN.'+str(n)+'/RFcan'
png_folder='/media/zekun/Software/fenglin_niu/yn-pic/YN.'+str(n)+'.png'
# 遍历文件夹中的所有SAC文件
sac_files = glob.glob(os.path.join(sac_folder, '*BHR'))
for sac_file in sac_files:
# 读取SAC文件
st = read(sac_file)
tr = st[0]
# 获取反方位角信息
back_azimuth = tr.stats.sac['baz']
# 找到对应的bin
bin_angle = int((back_azimuth-2.5) // bin_angle_interval) * bin_angle_interval+bin_angle_interval
# 对波形进行叠加
stacked_waveforms[bin_angle] += tr.data
event_counts[bin_angle] += 1
# 对叠加后的波形进行归一化处理
for angle in bin_angle_range:
stacked_waveforms[angle] /= event_counts[angle]
stacked_waveforms[angle]*=30
# 绘制每个bin的波形
for angle in bin_angle_range:
# 获取波形数据
waveform = stacked_waveforms[angle]
#print(waveform)
# 设置绘图的纵轴刻度位置
y = ((angle) // bin_angle_interval) * 5
# 绘制波形
axs[0].plot(tr.times()-10, waveform + y, color='black',linewidth=0.1)
# 填充颜色
axs[0].fill_between(tr.times()-10, waveform + y, y, where=(waveform >= 0), color='lightcoral')
axs[0].fill_between(tr.times()-10, waveform + y, y, where=(waveform < 0), color='lightblue')
#next
for index, row in ani.iterrows():
# 检查sta列的值是否为7
if row['sta'] == i:
# 获取对应的k列的值
t = row['t']
for index, row in rfcan.iterrows():
# 检查sta列的值是否为7
if row['sta'] == i:
# 获取对应的k列的值
phi_ra = row['phi(ra)']
phi_jof = row['phi(jof)']
dt_ra = row['dt(ra)']
dt_jof = row['dt(jof)']
# 绘制竖线
axs[0].axvline(x=t, linestyle='dashed', color='gray')
#cos
y = np.linspace(0, 361, 1000)
#x = t - dt_ra * np.cos(2*np.deg2rad(phi_ra - y))/2
axs[0].plot(x, y, linewidth=2, label='ra', zorder=1)
x = t - dt_jof * np.cos(2*np.deg2rad(phi_jof - y))/2
axs[0].plot(x, y, color='gold', linewidth=2, label='jof', alpha=0.9, zorder=1)
# 绘制每个bin的波形
k=0
for angle in bin_angle_range:
# 获取波形数据
waveform = stacked_waveforms[angle]
#print(waveform)
# 设置绘图的纵轴刻度位置
y = ((angle) // bin_angle_interval) * 5
# 绘制波形
#ax.plot(tr.times()-10, waveform + y, color='black',linewidth=0.1)
# 填充颜色
#ax.fill_between(tr.times()-10, waveform + y, y, where=(waveform >= 0), color='lightcoral')
#ax.fill_between(tr.times()-10, waveform + y, y, where=(waveform < 0), color='lightblue')
if str(waveform[0])!='nan':
k+=1
#print(waveform)
#peak value
max_values = []
trace_data = waveform
trace_times = tr.times()
for j in range(0,57):
if float(ani['sta'][j])==float(str(i)):
start_time=float(ani['t'][j])-float(ani['n5'][j])+10
end_time=float(ani['t'][j])+float(ani['n5'][j])+10
indices = [i for i, t in enumerate(trace_times) if start_time <= t <= end_time]
#print(trace_data[indices])
if indices:
max_value = max(trace_data[indices])
max_indices = [i for i, value in enumerate(trace_data) if value == max_value]
max_times = [(trace_times[i]-10) for i in max_indices]
max_values.extend(max_times)
if max_values[-1]==max_values[0]+(len(max_values)-1)*0.005:
max_time=(max_values[0]+max_values[-1])/2
else:
max_time=max_values[0]
if k==1:
axs[0].scatter(max_values[0], y, s=30, color='green', alpha=0.9, label='peak', zorder=2)
else:
axs[0].scatter(max_values[0], y, s=30, color='green', alpha=0.9, zorder=2)
# 设置图像范围和标签
axs[0].set_xlim(-2, 20)
axs[0].set_ylim(-5, len(bin_angle_range) * 5 + 5)
axs[0].set_xlabel('Time (s)', fontsize=7)
axs[0].set_ylabel('Back Azimuth (°)', fontsize=7)
axs[0].set_title('R RFs, before correction', fontsize=8)
# SAC文件所在的文件夹路径
sac_folder = '/media/zekun/Software/fenglin_niu/RFcan/examples/YN.'+str(n)+'/RFcan'
#png_folder='/media/zekun/Software/fenglin_niu/yn-pic/BHR.F/YN.'+str(n)+'.png'
# 创建存储叠加波形和事件数目的字典
stacked_waveforms = {angle: np.zeros(6002) for angle in bin_angle_range}
event_counts = {angle: 0 for angle in bin_angle_range}
# 遍历文件夹中的所有SAC文件
sac_files = glob.glob(os.path.join(sac_folder, '*BHR.F'))
for sac_file in sac_files:
# 读取SAC文件
st = read(sac_file)
tr = st[0]
# 获取反方位角信息
back_azimuth = tr.stats.sac['baz']
# 找到对应的bin
bin_angle = int(back_azimuth // bin_angle_interval) * bin_angle_interval
# 对波形进行叠加
stacked_waveforms[bin_angle] += tr.data
event_counts[bin_angle] += 1
# 对叠加后的波形进行归一化处理
for angle in bin_angle_range:
stacked_waveforms[angle] /= event_counts[angle]
stacked_waveforms[angle]*=30
# 创建绘图对象
#%matplotlib inline
#%config InlineBackend.figure_format = 'svg'
# 绘制每个bin的波形
for angle in bin_angle_range:
# 获取波形数据
waveform = stacked_waveforms[angle]
# 设置绘图的纵轴刻度位置
y = (angle // bin_angle_interval + 1) * 5
# 绘制波形
axs[1].plot(tr.times()-10, waveform + y, color='black',linewidth=0.1)
# 填充颜色
axs[1].fill_between(tr.times()-10, waveform + y, y, where=(waveform >= 0), color='lightcoral')
axs[1].fill_between(tr.times()-10, waveform + y, y, where=(waveform < 0), color='lightblue')
for index, row in ani.iterrows():
# 检查sta列的值是否为7
if row['sta'] == i:
# 获取对应的k列的值
t = row['t']
# 绘制竖线
axs[1].axvline(x=float(t), linestyle='dashed', color='gray')
# 设置图像范围和标签
axs[1].set_xlim(-2, 20)
axs[1].set_ylim(-5, len(bin_angle_range) * 5 + 5)
axs[1].set_xlabel('Time (s)', fontsize=7)
#axs[1].set_ylabel('Amplitude')
axs[1].set_title('R RFs, after correction', fontsize=8)
for ax in axs.flat: # 遍历所有的子图 (如果是二维数组的话需要flatten化)
# 设置横轴的大刻度和小刻度
ax.yaxis.set_major_locator(plt.MultipleLocator(50))
ax.yaxis.set_minor_locator(plt.MultipleLocator(10))
#ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, pos: f'{x:.1f}' if pos % 2 != 0 else ''))
plt.tick_params(axis='y', top=True, labeltop=True)
ax.tick_params(axis='y', which='minor', top=True)
ax.tick_params(axis='both', which='major', labelsize=6)
#plt.subplots_adjust(wspace=0.1)
# 显示图像
plt.savefig(png_folder,dpi=1000)
plt.show()
我想用这段代码绘制一行二列的两个子图,为什么在jupyter notebook中得到的是二行一列的图,甚至生成的png图像中只有第二个子图