《15节课Python绘图从入门到精通》8——时间序列与频率分析实操

import numpy as npimport pandas as pdfrom datetime import datetime, timedeltanp.random.seed(2024)# ----------时间戳数据（多日温度监测）----------# 生成2024年3月1日至3月7日，每10分钟一条的温度记录start_date = datetime(2024, 3, 1, 0, 0, 0)n_days = 7n_points = n_days * 24 * 6# 每10分钟timestamps = [start_date + timedelta(minutes=10*i) for i in range(n_points)]# 模拟温度：日周期波动 + 整体趋势 + 噪声time_hours = np.arange(n_points) / 6# 以小时计temp = ( 15 + 5 * np.sin(2*np.pi*time_hours/24 - np.pi/2) +  # 日变化0.5 * np.sin(2*np.pi*time_hours/(24*7)) +           # 周趋势0.02 * time_hours +                                 # 轻微升温         np.random.normal(0, 0.8, n_points) )df_temp_time = pd.DataFrame({'时间戳': timestamps, '温度': temp})df_temp_time.to_csv('时间序列温度数据.csv', index=False, encoding='utf-8-sig')# ---------- 8.2 模拟EEG信号（多通道脑电）----------# 生成10秒的模拟EEG，采样率500 Hzfs = 500# Hzduration = 10# secondst_eeg = np.linspace(0, duration, int(fs * duration))# 模拟多个频段成分：delta (0.5-4 Hz), theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz), gamma (30-50 Hz)eeg_signal = ( 0.5 * np.sin(2*np.pi*2*t_eeg) +                         # delta0.8 * np.sin(2*np.pi*6*t_eeg) +                         # theta1.2 * np.sin(2*np.pi*10*t_eeg) * (1 + 0.3*t_eeg) +      # alpha（幅度变化）0.6 * np.sin(2*np.pi*20*t_eeg) +                        # beta0.3 * np.sin(2*np.pi*40*t_eeg) * np.exp(-t_eeg/5) )     # gamma（衰减）# 添加瞬态伪迹（眨眼）eeg_signal += 2.0 * np.exp(-((t_eeg-2.5)**2)/0.01)  # 眨眼1eeg_signal += 1.5 * np.exp(-((t_eeg-7.0)**2)/0.02)  # 眨眼2# 背景噪声eeg_signal += np.random.normal(0, 0.15, len(t_eeg))np.savez('eeg_signal.npz', t=t_eeg, signal=eeg_signal, fs=fs)# ---------- 8.3 风向风速数据（一年逐时数据）----------# 模拟某地一年的风向（0-360°）和风速（m/s）n_hours_year = 365 * 24hours = np.arange(n_hours_year)# 风速：Weibull分布特征，季节变化wind_speed = ( 5 + 2 * np.sin(2*np.pi*hours/(365*24)*2 - 1) +  # 半年周期1.5 * np.sin(2*np.pi*hours/(24)) +             # 日变化               np.random.weibull(2, size=n_hours_year) * 3 )wind_speed = np.clip(wind_speed, 0.5, 25)# 风向：主导风向为西南(225°)和东北(45°)，加随机扰动wind_dir_base = np.where(np.random.rand(n_hours_year) < 0.6, 225, 45)wind_dir = (wind_dir_base + 30 * np.random.randn(n_hours_year)) % 360df_wind = pd.DataFrame({'小时序': hours, '风向': wind_dir, '风速': wind_speed})df_wind.to_csv('风向风速数据.csv', index=False, encoding='utf-8-sig')# ---------- 8.4 二维流场数据（圆柱绕流模拟简化）----------# 在矩形区域内生成速度场 (u, v)x_flow = np.linspace(-2, 4, 30)y_flow = np.linspace(-2, 2, 25)Xf, Yf = np.meshgrid(x_flow, y_flow)# 均匀来流 + 圆柱扰动（简化势流叠加）U_inf = 1.0radius = 0.5# 极坐标r = np.sqrt(Xf**2 + Yf**2)theta = np.arctan2(Yf, Xf)# 势流解：u = U_inf * (1 - (R^2/r^2)*cos(2θ)), v = -U_inf * (R^2/r^2)*sin(2θ)mask = r > radiusu_flow = np.where(mask, U_inf * (1 - (radius**2 / r**2) * np.cos(2*theta)), 0)v_flow = np.where(mask, -U_inf * (radius**2 / r**2) * np.sin(2*theta), 0)# 在圆柱内部速度设为零u_flow[~mask] = np.nanv_flow[~mask] = np.nan# 添加一些随机扰动模拟湍流np.random.seed(42)u_flow += np.random.normal(0, 0.05, u_flow.shape)v_flow += np.random.normal(0, 0.05, v_flow.shape)np.savez('flow_field.npz', x=x_flow, y=y_flow, u=u_flow, v=v_flow, X=Xf, Y=Yf)print("数据生成完成：")print("- 时间序列温度数据.csv（8.1节用）")print("- eeg_signal.npz（8.2节用）")print("- 风向风速数据.csv（8.3节用）")print("- flow_field.npz（8.4节用）")

import pandas as pdimport matplotlib.pyplot as plt# 学术样式设置defset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('时间序列温度数据.csv', parse_dates=['时间戳'])df.set_index('时间戳', inplace=True)fig, ax = plt.subplots(figsize=(5.0, 3.2))# 旋转标签，避免重叠（学术图常用45°/90°，这里用45°更美观）plt.xticks(rotation=45, ha='right')  # ha='right'让标签右对齐，更整齐ax.plot(df.index, df['温度'], color='#4472C4', linewidth=1.0)ax.set_xlabel('时间')ax.set_ylabel('温度 (°C)')ax.grid(alpha=0.3, linestyle='--', linewidth=0.5)plt.tight_layout()plt.show()fig.savefig('时间序列_基础.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时间序列_基础.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

pandas自动将时间戳列解析为datetime64类型，plot时横轴显示完整日期时间，但可能过于密集。

温度曲线呈现清晰的昼夜波动，叠加了长期趋势。默认横轴标签可能重叠，上例采用了plt.xticks(rotation=45, ha='right')进行格式化。

进阶1：自定义日期刻度格式（主/次刻度与日期格式化）

使用mdates设置主刻度为每天，次刻度为每6小时，并格式化日期标签。

import pandas as pdimport matplotlib.pyplot as pltimport matplotlib.dates as mdates# 学术样式设置defset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()# 读取数据df = pd.read_csv('时间序列温度数据.csv', parse_dates=['时间戳'])df.set_index('时间戳', inplace=True)# 创建画布fig, ax = plt.subplots(figsize=(6.5, 3.8))ax.plot(df.index, df['温度'], color='#4472C4', linewidth=1.0)# ===================== 核心：双行标签（日期 + 时间 自动换行）=====================# 每天只显示 0点、12点 两个标签（避免太多）ax.xaxis.set_major_locator(mdates.HourLocator(byhour=[0,  12], interval=1))# 关键：\n 表示换行！ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d\n%H:%M'))# 不旋转！居中！完全不重叠！plt.xticks(rotation=0, ha='center', fontsize=8)# ==============================================================================# 次刻度（只显示刻度线，不显示文字）ax.xaxis.set_minor_locator(mdates.HourLocator(byhour=[0,6,12,18]))ax.tick_params(axis='x', which='minor', labelbottom=False)# 坐标轴与网格ax.set_xlabel('时间', labelpad=8)ax.set_ylabel('温度 (°C)', labelpad=8)ax.grid(which='major', alpha=0.4, linestyle='-', linewidth=0.5)ax.grid(which='minor', alpha=0.2, linestyle='--', linewidth=0.3)plt.tight_layout()plt.show()# 保存fig.savefig('时间序列_双行标签.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时间序列_双行标签.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

刻度显示“月-日”和“时:分”，清晰展示时间跨度。

用分行来避免重叠，网格线区分主次刻度。

适用于多日连续监测数据的展示。

进阶2：双纵轴时间序列

在相同时间轴上绘制两个不同量纲的变量，使用双纵轴。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport matplotlib.dates as mdatesdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('时间序列温度数据.csv', parse_dates=['时间戳'])df.set_index('时间戳', inplace=True)# 模拟相对湿度（与温度负相关）np.random.seed(2024)humidity = 70 - 0.8 * (df['温度'] - 15) + np.random.normal(0, 5, len(df))humidity = np.clip(humidity, 30, 95)fig, ax1 = plt.subplots(figsize=(5.0, 3.5))ax1.plot(df.index, df['温度'], color='#C00000', linewidth=1.2, label='温度')ax1.set_xlabel('时间')ax1.set_ylabel('温度 (°C)', color='#C00000')ax1.tick_params(axis='y', labelcolor='#C00000')ax2 = ax1.twinx()ax2.plot(df.index, humidity, color='#0066CC', linewidth=1.2, label='相对湿度')ax2.set_ylabel('相对湿度 (%)', color='#0066CC')ax2.tick_params(axis='y', labelcolor='#0066CC')# 日期格式化ax1.xaxis.set_major_locator(mdates.DayLocator())ax1.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))ax1.tick_params(axis='x', rotation=30)# 图例lines1, labels1 = ax1.get_legend_handles_labels()lines2, labels2 = ax2.get_legend_handles_labels()ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper left', fontsize=8)ax1.grid(alpha=0.3, linestyle='--', linewidth=0.5)plt.tight_layout()plt.show()fig.savefig('时间序列_双轴.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时间序列_双轴.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

左轴温度（红色）右轴湿度（蓝色），颜色与坐标轴标签对应，避免混淆。

两条曲线反向变化，直观展示温湿负相关关系。

双轴图需注意纵轴范围匹配，避免误导视觉。

进阶3：时间区间高亮标注

使用axvspan在时间轴上高亮特定时段（如每日18:00至次日06:00），突出昼夜差异。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport matplotlib.dates as mdatesfrom datetime import datetime, timedeltadefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('时间序列温度数据.csv', parse_dates=['时间戳'])df.set_index('时间戳', inplace=True)fig, ax = plt.subplots(figsize=(5.0, 3.5))ax.plot(df.index, df['温度'], color='#4472C4', linewidth=1.0)# 高亮每日夜间时段（18:00至次日06:00）start_date = df.index.min().date()end_date = df.index.max().date()current = start_datewhile current <= end_date:    night_start = datetime.combine(current, datetime.strptime('18:00', '%H:%M').time())    night_end = datetime.combine(current + timedelta(days=1), datetime.strptime('06:00', '%H:%M').time())# 转换为pandas Timestamp以便索引    night_start_ts = pd.Timestamp(night_start)    night_end_ts = pd.Timestamp(night_end)    ax.axvspan(night_start_ts, night_end_ts, facecolor='gray', alpha=0.15, edgecolor='none')    current += timedelta(days=1)ax.xaxis.set_major_locator(mdates.DayLocator())ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d'))ax.set_xlabel('日期')ax.set_ylabel('温度 (°C)')ax.tick_params(axis='x', rotation=30)ax.grid(alpha=0.3, linestyle='--', linewidth=0.5)plt.tight_layout()plt.show()fig.savefig('时间序列_区间高亮.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时间序列_区间高亮.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

灰色半透明区域标记了每日夜间，温度明显低于白天，周期规律一目了然。

适用于展示周期性事件（如夜班、休息时段）对变量的影响。

可结合文本标注“夜间”于图例或图中。

进阶4：多时区时间转换与显示（UTC与本地时间）

若原始数据为UTC时间，需转换为本地时间并正确标注时区。

import pandas as pdimport matplotlib.pyplot as pltimport matplotlib.dates as mdatesdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()# 模拟UTC时间数据（假设原始数据是UTC）df = pd.read_csv('时间序列温度数据.csv', parse_dates=['时间戳'])# 假设读取时未指定时区，现将其设为UTCdf['时间戳'] = df['时间戳'].dt.tz_localize('UTC')# 转换为北京时间 (UTC+8)df['北京时间'] = df['时间戳'].dt.tz_convert('Asia/Shanghai')df.set_index('北京时间', inplace=True)fig, ax = plt.subplots(figsize=(5.0, 3.5))ax.plot(df.index, df['温度'], color='#4472C4', linewidth=1.0)ax.xaxis.set_major_locator(mdates.DayLocator())ax.xaxis.set_major_formatter(mdates.DateFormatter('%m-%d %H:%M'))ax.set_xlabel('北京时间 (CST)')ax.set_ylabel('温度 (°C)')ax.tick_params(axis='x', rotation=30)ax.grid(alpha=0.3, linestyle='--', linewidth=0.5)ax.set_title('UTC+8 时间轴', fontsize=10)plt.tight_layout()plt.show()fig.savefig('时间序列_时区转换.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时间序列_时区转换.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

使用tz_localize和tz_convert完成时区转换，索引变为带时区信息的DatetimeIndex。

横轴标签显示北京时间，标题注明时区，避免歧义。

适用于国际合作项目或跨时区观测数据。

实操任务2：绘制脑电信号时频谱图

本任务从基础的specgram绘制时频谱入手，逐步进阶到对数频率轴、小波变换、多通道对比及自定义归一化。

基础版：matplotlib.specgram绘制时频谱

使用specgram函数快速生成时频谱图。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()data = np.load('eeg_signal.npz')t = data['t']signal = data['signal']fs = int(data['fs'])fig, ax = plt.subplots(figsize=(5.0, 3.5))# 绘制时频谱，返回频谱矩阵用于后续自定义Pxx, freqs, bins, im = ax.specgram(signal, NFFT=256, Fs=fs, noverlap=200,                                    cmap='viridis', scale='dB')ax.set_xlabel('时间 (s)')ax.set_ylabel('频率 (Hz)')ax.set_ylim(0, 60)  # 关注0-60Hzcbar = fig.colorbar(im, ax=ax)cbar.set_label('功率谱密度 (dB)', fontsize=9)plt.tight_layout()plt.show()fig.savefig('时频谱_基础.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时频谱_基础.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

时频谱清晰显示alpha频段（10Hz）能量随时间增强，delta（2Hz）和theta（6Hz）持续存在，gamma（40Hz）衰减。

两个眨眼伪迹在低频段表现为宽带高能量。

线性频率轴下，低频区域被压缩，细节不易观察。

进阶1：对数频率轴与自定义频率刻度

将Y轴设置为对数刻度，并手动指定刻度位置以符合EEG频段习惯（如1, 4, 8, 13, 30, 50 Hz）。

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.ticker import ScalarFormatterdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()data = np.load('eeg_signal.npz')t, signal, fs = data['t'], data['signal'], int(data['fs'])fig, ax = plt.subplots(figsize=(5.0, 3.8))Pxx, freqs, bins, im = ax.specgram(signal, NFFT=256, Fs=fs, noverlap=200,                                    cmap='plasma', scale='dB')ax.set_yscale('log')ax.set_ylim(1, 60)  # 避免log(0)# 设置自定义刻度（EEG常用频段边界）ax.set_yticks([1, 4, 8, 13, 30, 50])ax.get_yaxis().set_major_formatter(ScalarFormatter())ax.set_xlabel('时间 (s)')ax.set_ylabel('频率 (Hz)')cbar = fig.colorbar(im, ax=ax)cbar.set_label('功率谱密度 (dB)', fontsize=9)# 添加频段标注for y, label in zip([2.5, 6, 10.5, 21.5, 40], ['δ', 'θ', 'α', 'β', 'γ']):    ax.text(0.01, y, label, transform=ax.get_yaxis_transform(), fontsize=8,             va='center', color='white', fontweight='bold')plt.tight_layout()plt.show()fig.savefig('时频谱_对数频率轴.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时频谱_对数频率轴.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

对数轴拉伸了低频区域，使delta、theta频段变化更清晰。

右侧标注希腊字母对应EEG经典频段，提升专业感。

颜色条使用plasma，比viridis对比度更高，适合功率谱展示。

进阶2：使用连续小波变换（CWT）绘制时频谱

相比短时傅里叶变换，小波变换在时频分辨率上更灵活，尤其适合非平稳信号。

import numpy as npimport matplotlib.pyplot as pltimport pywt #pip install PyWaveletsdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()data = np.load('eeg_signal.npz')t, signal, fs = data['t'], data['signal'], int(data['fs'])# 使用复Morlet小波进行连续小波变换scales = np.arange(1, 128)  # 尺度wavelet = 'cmor1.5-1.0'sampling_period = 1/fscoeffs, freqs_cwt = pywt.cwt(signal, scales, wavelet, sampling_period=sampling_period)# 计算功率power = np.abs(coeffs) ** 2fig, ax = plt.subplots(figsize=(5.0, 4.0))im = ax.contourf(t, freqs_cwt, power, levels=50, cmap='inferno', extend='both')ax.set_yscale('log')ax.set_ylim(1, 60)ax.set_yticks([1, 4, 8, 13, 30, 50])ax.set_xlabel('时间 (s)')ax.set_ylabel('频率 (Hz)')cbar = fig.colorbar(im, ax=ax)cbar.set_label('小波功率', fontsize=9)plt.tight_layout()plt.show()fig.savefig('时频谱_小波变换.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时频谱_小波变换.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

小波变换对瞬态事件（如眨眼）的时域定位更精确，高频成分随时间变化更平滑。

计算量大于STFT，但时频分辨率更优。

可调整小波基函数以适应不同信号特性。

进阶3：多通道时频谱子图对比（如Fz, Cz, Pz电极）

模拟三个通道的EEG，并排绘制时频谱，共享颜色条与坐标轴。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()# 模拟三个通道（振幅略有差异）np.random.seed(2024)data = np.load('eeg_signal.npz')t, signal_base, fs = data['t'], data['signal'], int(data['fs'])channels = {'Fz': signal_base * 0.8, 'Cz': signal_base, 'Pz': signal_base * 1.2 + 0.5*np.sin(2*np.pi*8*t)}fig, axes = plt.subplots(3, 1, figsize=(5.5, 7.0), sharex=True)for ax, (ch_name, sig) in zip(axes, channels.items()):    Pxx, freqs, bins, im = ax.specgram(sig, NFFT=256, Fs=fs, noverlap=200,                                        cmap='viridis', scale='dB', vmin=-20, vmax=20)    ax.set_yscale('log')    ax.set_ylim(1, 60)    ax.set_yticks([1, 4, 8, 13, 30, 50])    ax.set_ylabel('频率 (Hz)')    ax.set_title(f'通道 {ch_name}', fontsize=9)axes[-1].set_xlabel('时间 (s)')# 共享颜色条cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7])fig.colorbar(im, cax=cbar_ax, label='功率谱密度 (dB)')plt.tight_layout(rect=[0, 0, 0.91, 1])plt.show()fig.savefig('时频谱_多通道对比.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时频谱_多通道对比.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

三个通道时频谱垂直排列，便于比较空间分布差异。

统一颜色映射范围（vmin/vmax）确保颜色可跨子图比较。

适用于EEG地形图配套的时频分析展示。

进阶4：基线归一化与事件相关谱扰动（ERSP）显示

以刺激前基线为参考，绘制相对功率变化的时频谱。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()data = np.load('eeg_signal.npz')t, signal, fs = data['t'], data['signal'], int(data['fs'])# 计算时频谱Pxx, freqs, bins, im = plt.specgram(signal, NFFT=256, Fs=fs, noverlap=200,                                     scale='dB', mode='psd')plt.close()  # 不显示临时图# 定义基线区间（前2秒为基线）baseline_idx = bins < 2.0baseline_mean = np.mean(Pxx[:, baseline_idx], axis=1, keepdims=True)# 转换为dB变化 (10*log10(ratio))Pxx_norm = 10 * np.log10(Pxx / baseline_mean)fig, ax = plt.subplots(figsize=(5.0, 4.0))im = ax.contourf(bins, freqs, Pxx_norm, levels=50, cmap='RdBu_r', vmin=-5, vmax=5)ax.set_yscale('log')ax.set_ylim(1, 60)ax.set_yticks([1, 4, 8, 13, 30, 50])ax.set_xlabel('时间 (s)')ax.set_ylabel('频率 (Hz)')cbar = fig.colorbar(im, ax=ax)cbar.set_label('相对功率变化 (dB)', fontsize=9)ax.axvline(2.0, color='k', linestyle='--', linewidth=1, label='事件起点')ax.legend(fontsize=8)plt.tight_layout()plt.show()fig.savefig('时频谱_基线归一化.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('时频谱_基线归一化.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

使用发散色盘RdBu_r，蓝色表示功率降低，红色表示增加。

基线归一化消除了背景活动差异，突出事件相关变化。

垂直虚线标记事件起点（如刺激呈现时刻），是事件相关电位/场分析的标准可视化。

实操任务3：处理风向与风速数据

风向玫瑰图（Wind Rose）是气象学中展示风向频率与风速关系的标准工具。本任务从基础柱状玫瑰图进阶到风速分级、概率密度、极坐标KDE及多子图对比。

基础版：极坐标柱状图（风向频率玫瑰图）

将风向划分为16个方位，统计各方位出现频率，绘制极坐标柱状图。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('风向风速数据.csv')# 将风向划分为16个方位（每个22.5°）bins = np.arange(0, 361, 22.5)labels = ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE','S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW']df['风向方位'] = pd.cut(df['风向'], bins=bins, labels=labels, right=False)# 计算各方位频率freq = df['风向方位'].value_counts(normalize=True).reindex(labels) * 100# 极坐标柱状图theta = np.linspace(0, 2*np.pi, 16, endpoint=False)width = 2*np.pi / 16fig = plt.figure(figsize=(5.0, 5.0))ax = fig.add_subplot(111, projection='polar')bars = ax.bar(theta, freq.values, width=width, bottom=0,               color='#4472C4', alpha=0.7, edgecolor='black', linewidth=0.5)# 设置角度标签ax.set_xticks(theta)ax.set_xticklabels(labels, fontsize=8)ax.set_theta_zero_location('N')  # 0度指向北ax.set_theta_direction(-1)       # 顺时针ax.set_ylabel('频率 (%)', fontsize=9)ax.set_title('风向频率玫瑰图', fontsize=10, pad=20)plt.tight_layout()plt.show()fig.savefig('风向玫瑰图_基础.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('风向玫瑰图_基础.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

极坐标柱状图清晰显示主导风向为西南（SW）和东北（NE），符合模拟设定。

0度指向北，顺时针旋转，符合气象学惯例。

各方位频率百分比标注于径向轴。

进阶1：分级风速玫瑰图（不同颜色代表风速区间）

将每个风向方位按风速大小进一步细分，用堆叠柱状图或不同颜色展示。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('风向风速数据.csv')bins_dir = np.arange(0, 361, 22.5)labels_dir = ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE','S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW']df['风向方位'] = pd.cut(df['风向'], bins=bins_dir, labels=labels_dir, right=False)# 风速分级（Beaufort风力等级简化）speed_bins = [0, 2, 4, 6, 10, 15, 30]speed_labels = ['0-2', '2-4', '4-6', '6-10', '10-15', '>15']df['风速级'] = pd.cut(df['风速'], bins=speed_bins, labels=speed_labels, right=False)# 计算交叉频率（百分比）cross_tab = pd.crosstab(df['风向方位'], df['风速级'], normalize='index') * 100# 绘制分级玫瑰图theta = np.linspace(0, 2*np.pi, 16, endpoint=False)width = 2*np.pi / 16colors = ['#2E75B6', '#70AD47', '#FFC000', '#ED7D31', '#C00000', '#7030A0']fig = plt.figure(figsize=(6.0, 6.0))ax = fig.add_subplot(111, projection='polar')bottom = np.zeros(16)for i, (speed_cat, color) in enumerate(zip(speed_labels, colors)):    values = cross_tab[speed_cat].values    ax.bar(theta, values, width=width, bottom=bottom, label=speed_cat,           color=color, alpha=0.8, edgecolor='white', linewidth=0.3)    bottom += valuesax.set_xticks(theta)ax.set_xticklabels(labels_dir, fontsize=8)ax.set_theta_zero_location('N')ax.set_theta_direction(-1)ax.set_ylabel('频率 (%)', fontsize=9)ax.legend(title='风速 (m/s)', bbox_to_anchor=(1.1, 0.5), loc='center left', fontsize=8)ax.set_title('分级风速玫瑰图', fontsize=10, pad=20)plt.tight_layout()plt.show()fig.savefig('风向玫瑰图_风速分级.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('风向玫瑰图_风速分级.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

堆叠柱状图展示了每个风向方位内不同风速的频率占比。

西南风不仅频率高，且高风速（红色/紫色）比例也大，说明该方向风强劲。

图例置于右侧，颜色区分风速区间，符合气象玫瑰图标准。

进阶2：概率密度玫瑰图（极坐标KDE）

使用核密度估计在极坐标下绘制风向的概率密度分布，比离散方位更平滑。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import gaussian_kdedefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('风向风速数据.csv')angles = np.radians(df['风向'].values)# 在圆周上进行核密度估计（需处理周期性）# 方法：在[0, 2π]外复制数据以处理边界angles_ext = np.concatenate([angles, angles + 2*np.pi, angles - 2*np.pi])kde = gaussian_kde(angles_ext, bw_method=0.1)theta_grid = np.linspace(0, 2*np.pi, 360)density = kde.evaluate(theta_grid)density = density * len(angles) / density.sum()  # 缩放为近似频率fig = plt.figure(figsize=(5.5, 5.5))ax = fig.add_subplot(111, projection='polar')ax.plot(theta_grid, density, color='#C00000', linewidth=2.5)ax.fill(theta_grid, density, color='#C00000', alpha=0.3)ax.set_xticks(np.radians(np.arange(0, 360, 45)))ax.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'], fontsize=9)ax.set_theta_zero_location('N')ax.set_theta_direction(-1)ax.set_ylabel('概率密度', fontsize=9)ax.set_title('风向概率密度玫瑰图 (KDE)', fontsize=10, pad=20)plt.tight_layout()plt.show()fig.savefig('风向玫瑰图_KDE.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('风向玫瑰图_KDE.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

KDE曲线平滑展示了风向分布的双峰特征（东北和西南），峰宽反映风向变异程度。

相比离散柱状图，KDE不受分箱数量影响，更适合展示连续角度分布。

需注意边界连续性处理，否则在0°/360°处可能出现断裂。

进阶3：联合风速-风向玫瑰图（颜色映射风速大小）

在极坐标散点图中，用颜色表示风速，展示所有样本点。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()df = pd.read_csv('风向风速数据.csv')# 为了可视化清晰，抽取一部分点df_sample = df.sample(n=2000, random_state=42)angles = np.radians(df_sample['风向'].values)speeds = df_sample['风速'].valuesfig = plt.figure(figsize=(6.0, 6.0))ax = fig.add_subplot(111, projection='polar')sc = ax.scatter(angles, speeds, c=speeds, cmap='viridis', s=3, alpha=0.6)cbar = plt.colorbar(sc, ax=ax, pad=0.1)cbar.set_label('风速 (m/s)', fontsize=9)ax.set_xticks(np.radians(np.arange(0, 360, 45)))ax.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'], fontsize=9)ax.set_theta_zero_location('N')ax.set_theta_direction(-1)ax.set_ylabel('风速 (m/s)', fontsize=9)ax.set_title('风速-风向联合分布散点图', fontsize=10, pad=20)plt.tight_layout()plt.show()

fig.savefig('风向玫瑰图_联合散点.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('风向玫瑰图_联合散点.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

散点图展示了每个观测的风向和风速，颜色越亮风速越大。

可观察到西南风不仅频率高，且高风速点集中。

径向轴直接表示风速大小，无需额外变换。

进阶4：多季节对比玫瑰图（子图排列）

将一年数据按季节划分，绘制四个子玫瑰图，对比风向季节性变化。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()# ===================== 1. 读取数据 =====================df = pd.read_csv('风向风速数据.csv')# 季节划分（修复：先生成全列，不产生切片副本）df['季节'] = pd.cut(    df['小时序'] % (365 * 24),    bins=[0, 90 * 24, 180 * 24, 270 * 24, 365 * 24],    labels=['春季', '夏季', '秋季', '冬季'],    right=False)# ===================== 2. 风向划分配置 =====================bins_dir = np.arange(0, 361, 22.5)labels_dir = ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE','S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW']theta = np.linspace(0, 2 * np.pi, 16, endpoint=False)width = 2 * np.pi / 16# ===================== 3. 绘制极坐标子图 =====================fig, axes = plt.subplots(2, 2, figsize=(7.0, 7.0), subplot_kw={'projection': 'polar'})axes = axes.flatten()seasons = ['春季', '夏季', '秋季', '冬季']for ax, season in zip(axes, seasons):# 筛选季节数据（使用 .copy() 彻底消除警告）    sub_df = df[df['季节'] == season].copy()# 安全新增列，无警告    sub_df['风向方位'] = pd.cut(        sub_df['风向'],        bins=bins_dir,        labels=labels_dir,        right=False    )# 计算频率    freq = sub_df['风向方位'].value_counts(normalize=True).reindex(labels_dir).fillna(0) * 100# 绘图    ax.bar(        theta, freq.values,        width=width,        bottom=0,        color='#4472C4',        alpha=0.7,        edgecolor='black',        linewidth=0.3    )# 极坐标设置：正北为0，顺时针    ax.set_xticks(theta)    ax.set_xticklabels(labels_dir, fontsize=7)    ax.set_theta_zero_location('N')    ax.set_theta_direction(-1)# 标题与标签    ax.set_title(season, fontsize=10, pad=15)    ax.set_ylabel('')  # 极坐标图不建议显示ylabel，会重叠# 调整布局，防止标签重叠plt.tight_layout(pad=2.0)plt.show()# 保存高清图fig.savefig('风向玫瑰图_四季对比.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('风向玫瑰图_四季对比.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

四季玫瑰图对比可揭示风向的季节性变化（本模拟数据未设季节差异，故各图相似）。

实际应用中可观察到季风转换或局地环流特征。

子图共享相同角度刻度，便于横向比较。

实操任务4：用quiver绘制二维流场速度矢量箭头图

矢量场可视化在流体力学、电磁学中至关重要。本任务从基础quiver绘制进阶到颜色映射、流线叠加、自适应缩放及多子图对比。

基础版：quiver绘制速度矢量场

读取流场数据，使用quiver绘制速度矢量。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()data = np.load('flow_field.npz')X, Y = data['X'], data['Y']U, V = data['u'], data['v']fig, ax = plt.subplots(figsize=(6.0, 4.0))# 绘制矢量，scale控制箭头大小q = ax.quiver(X, Y, U, V, scale=30, width=0.002, color='#4472C4')ax.set_xlabel('x')ax.set_ylabel('y')ax.set_aspect('equal')ax.set_title('圆柱绕流速度场', fontsize=10)# 绘制圆柱轮廓theta = np.linspace(0, 2*np.pi, 100)ax.plot(0.5*np.cos(theta), 0.5*np.sin(theta), 'k-', linewidth=1.5, label='圆柱')ax.legend(fontsize=8)plt.tight_layout()plt.show()fig.savefig('流场矢量_基础.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('流场矢量_基础.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

quiver在网格点绘制箭头，箭头长度和方向表示速度矢量。

scale=30控制箭头整体缩放，值越大箭头越短；width调整箭头杆粗细。

圆柱绕流特征清晰：上游均匀来流，绕过圆柱后形成尾流区（低速）。

进阶1：箭头颜色映射速度大小（标量场）

将箭头颜色与速度幅值关联，并添加颜色条。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 300set_academic_style()data = np.load('flow_field.npz')X, Y = data['X'], data['Y']U, V = data['u'], data['v']speed = np.sqrt(U**2 + V**2)fig, ax = plt.subplots(figsize=(6.0, 4.0))# 使用c参数映射颜色，cmap选择q = ax.quiver(X, Y, U, V, speed, scale=30, width=0.002, cmap='hot')cbar = fig.colorbar(q, ax=ax)cbar.set_label('速度大小 (m/s)', fontsize=9)ax.set_xlabel('x')ax.set_ylabel('y')ax.set_aspect('equal')ax.set_title('圆柱绕流速度场（颜色表示速度大小）', fontsize=10)theta = np.linspace(0, 2*np.pi, 100)ax.plot(0.5*np.cos(theta), 0.5*np.sin(theta), 'w-', linewidth=1.5)plt.tight_layout()plt.show()fig.savefig('流场矢量_颜色映射.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('流场矢量_颜色映射.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

箭头颜色从深红（高速）到亮黄（低速）变化，直观展示圆柱前缘驻点速度为零，侧缘加速。

颜色条标明数值范围，便于定量解读。

白色圆柱轮廓与深色背景对比清晰。

进阶2：quiver与streamplot叠加（流线+箭头）

同时绘制流线和箭头，流线展示轨迹，箭头指示方向与相对大小。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 1200set_academic_style()data = np.load('flow_field.npz')X, Y = data['X'], data['Y']U, V = data['u'], data['v']speed = np.sqrt(U**2 + V**2)fig, ax = plt.subplots(figsize=(6.0, 4.0))# 背景流线（灰色）strm = ax.streamplot(X, Y, U, V, color='gray', linewidth=0.8, density=1.5, arrowsize=0)# 叠加箭头（半透明）q = ax.quiver(X[::2, ::2], Y[::2, ::2], U[::2, ::2], V[::2, ::2],               speed[::2, ::2], scale=30, width=0.003, cmap='plasma', alpha=0.9)cbar = fig.colorbar(q, ax=ax)cbar.set_label('速度大小 (m/s)', fontsize=9)ax.set_xlabel('x')ax.set_ylabel('y')ax.set_aspect('equal')ax.set_title('流线与速度矢量叠加', fontsize=10)theta = np.linspace(0, 2*np.pi, 100)ax.plot(0.5*np.cos(theta), 0.5*np.sin(theta), 'k-', linewidth=1.5)plt.tight_layout()plt.show()fig.savefig('流场矢量_流线叠加.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('流场矢量_流线叠加.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

灰色流线提供了连续的运动轨迹，箭头在稀疏网格上指示方向与强度。

流线密度通过density参数控制，避免过度拥挤。

结合两者优势，既看到整体流动模式，又能读取局部速度值。

进阶3：自适应箭头缩放与缺失值处理

当速度场动态范围很大时，统一缩放会导致低速区箭头过小或高速区箭头过大。可对箭头长度进行非线性变换（如对数或幂律）或设置最小阈值。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 300set_academic_style()data = np.load('flow_field.npz')X, Y = data['X'], data['Y']U, V = data['u'], data['v']speed = np.sqrt(U**2 + V**2)# 对速度大小进行幂律压缩（0.5次方），使箭头长度差异减小speed_norm = speed ** 0.5scale_factor = 5# 调整整体大小U_scaled = U / (speed + 1e-6) * speed_norm * scale_factorV_scaled = V / (speed + 1e-6) * speed_norm * scale_factorfig, ax = plt.subplots(figsize=(6.0, 4.0))q = ax.quiver(X, Y, U_scaled, V_scaled, speed, scale=1, width=0.003, cmap='viridis')# 注意：scale=1因为我们已手动缩放，quiver不再进一步缩放cbar = fig.colorbar(q, ax=ax)cbar.set_label('原始速度大小 (m/s)', fontsize=9)ax.set_xlabel('x')ax.set_ylabel('y')ax.set_aspect('equal')ax.set_title('自适应箭头缩放（幂律压缩）', fontsize=10)theta = np.linspace(0, 2*np.pi, 100)ax.plot(0.5*np.cos(theta), 0.5*np.sin(theta), 'k-', linewidth=1.5)plt.tight_layout()plt.show()fig.savefig('流场矢量_自适应缩放.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('流场矢量_自适应缩放.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

通过幂律压缩（指数0.5），高速区箭头长度被抑制，低速区箭头仍可见。

手动缩放后设置scale=1，避免二次缩放干扰。

适用于速度跨度几个数量级的情况（如边界层流动）。

进阶4：多子图矢量场对比（不同雷诺数或时刻）

并排绘制多个流场状态，对比参数变化对流动的影响。

import numpy as npimport matplotlib.pyplot as pltdefset_academic_style():    plt.rcParams['font.family'] = ['Times New Roman', 'SimSun']    plt.rcParams['font.size'] = 9    plt.rcParams['axes.unicode_minus'] = False    plt.rcParams['axes.linewidth'] = 1.0    plt.rcParams['xtick.major.width'] = 1.0    plt.rcParams['ytick.major.width'] = 1.0    plt.rcParams['xtick.major.size'] = 3.5    plt.rcParams['ytick.major.size'] = 3.5    plt.rcParams['xtick.direction'] = 'in'    plt.rcParams['ytick.direction'] = 'in'    plt.rcParams['legend.frameon'] = False    plt.rcParams['pdf.fonttype'] = 42    plt.rcParams['savefig.dpi'] = 300set_academic_style()data = np.load('flow_field.npz')X, Y = data['X'], data['Y']U_base, V_base = data['u'], data['v']# 模拟不同雷诺数下的流场：通过放大扰动模拟湍流增强np.random.seed(42)U1 = U_base + np.random.normal(0, 0.03, U_base.shape)V1 = V_base + np.random.normal(0, 0.03, V_base.shape)U2 = U_base + np.random.normal(0, 0.1, U_base.shape)V2 = V_base + np.random.normal(0, 0.1, V_base.shape)U3 = U_base + np.random.normal(0, 0.2, U_base.shape)V3 = V_base + np.random.normal(0, 0.2, V_base.shape)fields = [(U_base, V_base, 'Re = 100 (层流)'),          (U1, V1, 'Re = 1000'),          (U2, V2, 'Re = 5000'),          (U3, V3, 'Re = 20000 (湍流)')]fig, axes = plt.subplots(2, 2, figsize=(7.0, 6.0))axes = axes.flatten()for ax, (U, V, title) in zip(axes, fields):    speed = np.sqrt(U**2 + V**2)    q = ax.quiver(X[::2, ::2], Y[::2, ::2], U[::2, ::2], V[::2, ::2],                   speed[::2, ::2], scale=25, width=0.003, cmap='coolwarm')    ax.set_title(title, fontsize=10)    ax.set_xlabel('x')    ax.set_ylabel('y')    ax.set_aspect('equal')    theta = np.linspace(0, 2*np.pi, 100)    ax.plot(0.5*np.cos(theta), 0.5*np.sin(theta), 'k-', linewidth=1)# 共享颜色条cbar_ax = fig.add_axes([0.92, 0.15, 0.02, 0.7])norm = plt.Normalize(vmin=0, vmax=2.0)sm = plt.cm.ScalarMappable(cmap='coolwarm', norm=norm)sm.set_array([])fig.colorbar(sm, cax=cbar_ax, label='速度大小 (m/s)')plt.tight_layout(rect=[0, 0, 0.91, 1])plt.show()fig.savefig('流场矢量_多状态对比.pdf', bbox_inches='tight', pad_inches=0.05)fig.savefig('流场矢量_多状态对比.png', bbox_inches='tight', pad_inches=0.05)

执行结果分析：

随着雷诺数增加，随机扰动增强，流场从有序变为紊乱，尾流区涡结构复杂化。

子图共享颜色条，保证色彩映射一致。

适用于展示参数化研究或时间序列中的流场演化。