python绘图plt.figure\subplot\add_subplots\Axes3D\contourf

本文详细介绍了Python中matplotlib库的绘图方法,包括plt.figure的参数解释,subplot创建单个和多个子图,面向对象API的add_subplots与add_axes用法,以及Axes3D模块下的3D绘图、填充等高线图和绘制轮廓线的技巧。文章适合对Python绘图感兴趣的读者学习。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

一、plt.figure参数解释

  matplotlib.pyplot.figure() 创建一个新的画布(figure)。

matplotlib.pyplot.figure(num=None, figsize=None, dpi=None, facecolor=None, edgecolor=None, frameon=True, FigureClass=<class 'matplotlib.figure.Figure'>, clear=False, **kwargs)

输入参数:
num:整型或者字符串,可选参数,默认:None。图像编号或名称,数字为编号 ,字符串为名称

  • 如果不提供该参数,一个新的画布(figure)将被创建而且画布数量将会增加。
  • 如果提供该参数,带有id的画布是已经存在的,激活该画布并返回该画布的引用。
  • 如果这个画布不存在,创建并返回画布实例。
  • 如果num是字符串,窗口标题将被设置为该图的数字。
    figsize:整型元组,可选参数 ,默认:None。每英寸的宽度和高度。如果不提供,默认值是figure.figsize。
    dpi:整型,可选参数,默认:None。每英寸像素点。如果不提供,默认是figure.dpi。
    facecolor:背景色。如果不提供,默认值:figure.facecolor。 [c=labels.astype(np.float)]
    edgecolor:边界颜色。如果不提供,默认值:figure.edgecolor。
    framemon:布尔类型,可选参数,默认值:True。如果是False,禁止绘制画图框。
    FigureClass:源于matplotlib.figure.Figure的类。(可选)使用自定义图实例。
    clear:布尔类型,可选参数,默认值:False。如果为True和figure已经存在时,这是清理掉改图。

返回值:
figure:Figure。返回的Figure实例也将被传递给后端的new_figure_manager,这允许将自定义的图类挂接到pylab接口中。附加的kwarg将被传递给图形init函数。

import matplotlib.pyplot as plt
#创建自定义图像
fig=plt.figure(figsize=(4,3),facecolor='blue')
plt.show()

二、subplot创建单个子图

  subplot可以规划figure划分为n个子图,但每条subplot命令只会创建一个子图

import numpy as np  
import matplotlib.pyplot as plt  
x = np.arange(0, 100)  
#作图1
plt.subplot(221)  
plt.plot(x, x)  
#作图2
plt.subplot(222)  
plt.plot(x, -x)  
 #作图3
plt.subplot(223)  
plt.plot(x, x ** 2)  
plt.grid(color='r', linestyle='--', linewidth=1,alpha=0.3)
#作图4
plt.subplot(224)  
plt.plot(x, np.log(x))  
plt.show()  

三、subplots创建多个子图

plt.subplots(
    nrows=1,
    ncols=1,
    sharex=False,
    sharey=False,
    squeeze=True,
    subplot_kw=None,
    gridspec_kw=None,
    **fig_kw,
)

  subplots参数与subplot相似

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
x = np.arange(0, 
import numpy as np import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.animation import FuncAnimation from mpl_toolkits.mplot3d import Axes3D import ipywidgets as widgets from IPython.display import display plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题 plt.rcParams["font.family"] = "SimHei" # 设置中文字体 # 物理常数 c = 3e8 # 光速 (m/s) class RawFDAArraySimulator: def __init__(self, f0=24e8, delta_f=100e3, N=16): """ 去归一化FDA阵列仿真器 f0: 中心频率 (Hz) delta_f: 频率增量 (Hz) - 关键FDA参数 N: 阵元数量 """ self.f0 = f0 self.delta_f = delta_f self.N = N self.d = 0.5 * c / f0 # 阵元间距 (半波长) # 仿真参数 - 无归一化处理 self.theta = np.linspace(-90, 90, 181) # 方位角 (度) self.R = np.linspace(100, 10000, 101) # 距离 (米) self.time_points = np.linspace(0, 1e-3, 100) # 时间序列 (秒) self.theta_rad = np.deg2rad(self.theta) # 转换为弧度 def calculate_fda_beam(self, t): """ 计算FDA阵列因子,体现空间-时间耦合特性 - 无归一化 参数: t: 时间点 (秒) 返回: 原始阵列因子幅度 (距离 × 角度) """ # 初始化阵列因子矩阵 AF = np.zeros((len(self.R), len(self.theta)), dtype=complex) for m in range(self.N): # 遍历每个阵元 # 计算当前阵元的频率 f_m = self.f0 + m * self.delta_f # 空间相位分量 (角度相关) space_phase = 2 * np.pi * f_m * (m * self.d * np.sin(self.theta_rad)) / c # 时间相位分量 (时间相关) time_phase = 2 * np.pi * self.delta_f * m * t # 距离相位分量 (距离相关) range_phase = 2 * np.pi * f_m * self.R[:, None] / c # 综合阵列因子 - 体现空间-时间耦合 AF += np.exp(1j * (space_phase - range_phase + time_phase)) return np.abs(AF) # 直接返回幅度值,无归一化 def plot_spatio_temporal_coupling(self, t=0.5e-3): """ 可视化空间-时间耦合特性 - 无归一化 """ fig = plt.figure(figsize=(15, 10)) # 计算波束方向图 beam_pattern = self.calculate_fda_beam(t) # 1. 距离-角度热力图 - 使用原始幅度值 plt.subplot(2, 2, 1) plt.imshow(beam_pattern, # 直接使用幅度值 extent=[self.theta.min(), self.theta.max(), self.R.min() / 1000, self.R.max() / 1000], aspect='auto', cmap='jet', origin='lower') plt.colorbar(label='阵列因子幅度') plt.xlabel('方位角 (度)') plt.ylabel('距离 (km)') plt.title(f'FDA空间-时间耦合特性 (t={t * 1000:.1f}ms)') # 2. 固定距离的切面 - 使用原始幅度值 plt.subplot(2, 2, 2) fixed_range_idx = 50 # 选择中间距离 plt.plot(self.theta, beam_pattern[fixed_range_idx, :]) plt.xlabel('方位角 (度)') plt.ylabel('阵列因子幅度') plt.title(f'固定距离{self.R[fixed_range_idx] / 1000:.1f}km的波束方向图') plt.grid(True) # 3. 固定角度的切面 - 使用原始幅度值 plt.subplot(2, 2, 3) fixed_angle_idx = 90 # 选择0度方向 plt.plot(self.R / 1000, beam_pattern[:, fixed_angle_idx]) plt.xlabel('距离 (km)') plt.ylabel('阵列因子幅度') plt.title(f'固定方位角{self.theta[fixed_angle_idx]:.0f}°的波束方向图') plt.grid(True) # 4. 空间-时间耦合示意图 - 使用原始幅度值 plt.subplot(2, 2, 4) # 在固定距离和角度上观察幅度随时间变化 fixed_range = 5000 # 固定距离5km fixed_angle = 0 # 固定角度0度 r_idx = np.argmin(np.abs(self.R - fixed_range)) a_idx = np.argmin(np.abs(self.theta - fixed_angle)) # 计算时间序列上的幅度变化 time_amplitude = [] for time_point in self.time_points: beam = self.calculate_fda_beam(time_point) time_amplitude.append(beam[r_idx, a_idx]) plt.plot(self.time_points * 1000, time_amplitude) plt.xlabel('时间 (ms)') plt.ylabel('阵列因子幅度') plt.title(f'固定点(θ={fixed_angle}°, R={fixed_range / 1000:.1f}km)的幅度随时间变化') plt.grid(True) plt.tight_layout() plt.show() def animate_auto_scanning(self): """ 创建自动扫描特性的动画 - 无归一化 """ fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # 创建网格 T, R_grid = np.meshgrid(self.theta, self.R / 1000) # 初始化空图 beam = self.calculate_fda_beam(0) surf = ax.plot_surface(T, R_grid, beam, cmap=cm.jet, alpha=0.8) # 设置坐标轴标签 ax.set_xlabel('方位角 (度)') ax.set_ylabel('距离 (km)') ax.set_zlabel('阵列因子幅度') ax.set_title('FDA自动扫描特性 (t=0 ms)') ax.set_zlim(0, self.N) # 幅度范围0到N # 更新函数用于动画 def update(frame): ax.clear() t = self.time_points[frame] beam = self.calculate_fda_beam(t) surf = ax.plot_surface(T, R_grid, beam, cmap=cm.jet, alpha=0.8) ax.set_xlabel('方位角 (度)') ax.set_ylabel('距离 (km)') ax.set_zlabel('阵列因子幅度') ax.set_title(f'FDA自动扫描特性 (t={t * 1000:.1f} ms)') ax.set_zlim(0, self.N) # 保持幅度范围一致 return surf, # 创建动画 ani = FuncAnimation(fig, update, frames=len(self.time_points), interval=50, blit=False) plt.tight_layout() plt.show() return ani def interactive_parameter_analysis(self): """ 交互式参数分析界面 - 无归一化 """ style = {'description_width': 'initial'} # 创建交互控件 delta_f_slider = widgets.FloatLogSlider( value=self.delta_f, min=3, max=6, step=0.1, description='频率增量Δf (Hz):', readout_format='.0f', style=style ) time_slider = widgets.FloatSlider( value=0, min=0, max=1e-3, step=1e-5, description='时间 t (s):', readout_format='.5f', style=style ) angle_slider = widgets.FloatSlider( value=0, min=-90, max=90, step=1, description='固定角度θ (度):', style=style ) range_slider = widgets.FloatSlider( value=5000, min=100, max=10000, step=100, description='固定距离R (m):', style=style ) # 更新函数 def update_plots(delta_f_val, t_val, fixed_angle, fixed_range): # 临时创建新的仿真器实例 temp_simulator = RawFDAArraySimulator( f0=self.f0, delta_f=delta_f_val, N=self.N ) # 计算波束方向图 beam_pattern = temp_simulator.calculate_fda_beam(t_val) # 创建图形 fig, axs = plt.subplots(1, 2, figsize=(15, 6)) # 1. 距离-角度热力图 - 使用原始幅度值 im = axs[0].imshow(beam_pattern, extent=[self.theta.min(), self.theta.max(), self.R.min() / 1000, self.R.max() / 1000], aspect='auto', cmap='jet', origin='lower') fig.colorbar(im, ax=axs[0], label='阵列因子幅度') axs[0].set_xlabel('方位角 (度)') axs[0].set_ylabel('距离 (km)') axs[0].set_title(f'FDA波束方向图 (Δf={delta_f_val / 1000:.1f}kHz, t={t_val * 1000:.2f}ms)') # 2. 固定角度和距离的幅度变化 r_idx = np.argmin(np.abs(self.R - fixed_range)) a_idx = np.argmin(np.abs(self.theta - fixed_angle)) # 计算时间序列上的幅度变化 time_amplitude = [] for time_point in self.time_points: temp_beam = temp_simulator.calculate_fda_beam(time_point) time_amplitude.append(temp_beam[r_idx, a_idx]) axs[1].plot(self.time_points * 1000, time_amplitude) axs[1].axvline(x=t_val * 1000, color='r', linestyle='--', label=f'当前时间: {t_val * 1000:.2f}ms') axs[1].set_xlabel('时间 (ms)') axs[1].set_ylabel('阵列因子幅度') axs[1].set_title(f'固定点(θ={fixed_angle}°, R={fixed_range / 1000:.1f}km)的幅度随时间变化') axs[1].grid(True) axs[1].legend() plt.tight_layout() plt.show() # 创建交互界面 widgets.interactive( update_plots, delta_f_val=delta_f_slider, t_val=time_slider, fixed_angle=angle_slider, fixed_range=range_slider ) # ========== 主执行函数 ========== if __name__ == "__main__": print("频控阵天线(FDA)波束特性仿真 (去归一化)") # 创建FDA仿真器实例 fda_simulator = RawFDAArraySimulator() print("1. 空间-时间耦合特性 (无归一化)") fda_simulator.plot_spatio_temporal_coupling() print("\n2. 自动扫描特性 (无归一化动画展示)") fda_simulator.animate_auto_scanning() print("\n3. 交互式参数分析 (无归一化)") fda_simulator.interactive_parameter_analysis()
最新发布
07-27
E:\Anaconda\python.exe C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py 生成模拟数据... 开始数据预处理... 数据预处理完成 E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 27969 (\N{CJK UNIFIED IDEOGRAPH-6D41}) missing from font(s) DejaVu Sans. fig.tight_layout() E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans. fig.tight_layout() E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans. fig.tight_layout() E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans. fig.tight_layout() E:\Anaconda\Lib\site-packages\statsmodels\tsa\seasonal.py:360: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans. fig.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 23395 (\N{CJK UNIFIED IDEOGRAPH-5B63}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 33410 (\N{CJK UNIFIED IDEOGRAPH-8282}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 24615 (\N{CJK UNIFIED IDEOGRAPH-6027}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 20998 (\N{CJK UNIFIED IDEOGRAPH-5206}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:149: UserWarning: Glyph 35299 (\N{CJK UNIFIED IDEOGRAPH-89E3}) missing from font(s) DejaVu Sans. plt.tight_layout() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 23395 (\N{CJK UNIFIED IDEOGRAPH-5B63}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 33410 (\N{CJK UNIFIED IDEOGRAPH-8282}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 24615 (\N{CJK UNIFIED IDEOGRAPH-6027}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 20998 (\N{CJK UNIFIED IDEOGRAPH-5206}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:150: UserWarning: Glyph 35299 (\N{CJK UNIFIED IDEOGRAPH-89E3}) missing from font(s) DejaVu Sans. plt.savefig('季节性分解.png', dpi=300) C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 27969 (\N{CJK UNIFIED IDEOGRAPH-6D41}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 37327 (\N{CJK UNIFIED IDEOGRAPH-91CF}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 23395 (\N{CJK UNIFIED IDEOGRAPH-5B63}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 33410 (\N{CJK UNIFIED IDEOGRAPH-8282}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 24615 (\N{CJK UNIFIED IDEOGRAPH-6027}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 20998 (\N{CJK UNIFIED IDEOGRAPH-5206}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 35299 (\N{CJK UNIFIED IDEOGRAPH-89E3}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 21547 (\N{CJK UNIFIED IDEOGRAPH-542B}) missing from font(s) DejaVu Sans. plt.show() C:\Users\cheny\Desktop\PythonProject2\水沙通量周期性模型.py:151: UserWarning: Glyph 27801 (\N{CJK UNIFIED IDEOGRAPH-6C99}) missing from font(s) DejaVu Sans. plt.show() 上面是报错内容 下面是代码 import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.dates as mdates import pywt from statsmodels.tsa.seasonal import seasonal_decompose from scipy import signal class HydrologicalAnalysis: def __init__(self, data_path=None): """ 水文时间序列分析工具 :param data_path: 数据文件路径(CSV格式) """ self.data = None if data_path: self.load_data(data_path) def load_data(self, file_path): """ 加载水文监测数据 :param file_path: CSV文件路径 """ # 读取数据并处理时间格式 self.data = pd.read_csv(file_path, parse_dates=['时间'], index_col='时间') # 检查必要字段 required_columns = ['水位(m)', '流量(m3/s)', '含沙量(kg/m3)'] if not all(col in self.data.columns for col in required_columns): raise ValueError("数据文件缺少必要列:水位(m), 流量(m3/s), 含沙量(kg/m3)") print(f"成功加载数据:{len(self.data)}条记录") return self.data def preprocess_data(self): """ 数据预处理流程 """ if self.data is None: raise ValueError("未加载数据,请先调用load_data方法") print("开始数据预处理...") # 1. 缺失值处理(前向填充+线性插值) self.data = self.data.ffill().interpolate(method='linear') # 2. 异常值处理(3σ原则) for col in ['流量(m3/s)', '含沙量(kg/m3)']: mean = self.data[col].mean() std = self.data[col].std() self.data[col] = np.where( (self.data[col] < mean - 3 * std) | (self.data[col] > mean + 3 * std), mean, self.data[col] ) # 3. 重采样为日数据 daily_data = self.data.resample('D').mean() print("数据预处理完成") return daily_data def seasonal_decomposition(self, series, period=12): """ 时间序列季节性分解 :param series: 时间序列数据 :param period: 季节性周期(月数据默认为12) :return: 分解结果 """ result = seasonal_decompose(series, model='additive', period=period) return result def wavelet_analysis(self, series, title='水文序列'): """ 小波分析主函数 :param series: 时间序列数据 :param title: 分析标题 :return: (小波系数, 主周期) """ # 1. 参数设置 scales = np.arange(1, 365) # 1天到1年尺度 wavelet = 'morl' # Morlet小波 # 2. 小波变换 coef, freqs = pywt.cwt(series, scales, wavelet) # 3. 小波方差计算 variance = np.mean(np.abs(coef) ** 2, axis=1) main_scale = scales[np.argmax(variance)] # 4. 可视化 self.plot_wavelet_results(series, scales, coef, variance, main_scale, title) return coef, main_scale def plot_wavelet_results(self, series, scales, coef, variance, main_scale, title): """ 绘制小波分析结果 """ plt.figure(figsize=(15, 12)) # 1. 原始序列图 plt.subplot(3, 1, 1) plt.plot(series.index, series) plt.title(f'{title} - 原始时间序列') plt.xlabel('日期') plt.ylabel('数值') plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) # 2. 小波系数实部等值线图 plt.subplot(3, 1, 2) plt.contourf(series.index, scales, np.real(coef), cmap='jet', levels=100) plt.colorbar(label='小波系数实部') plt.title(f'{title} - 小波系数实部等值线图') plt.ylabel('时间尺度(天)') plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) # 3. 小波方差图 plt.subplot(3, 1, 3) plt.plot(scales, variance) plt.axvline(main_scale, color='red', linestyle='--', label=f'主周期: {main_scale}天') plt.title(f'{title} - 小波方差分析') plt.xlabel('时间尺度(天)') plt.ylabel('方差') plt.legend() plt.tight_layout() plt.savefig(f'{title}_小波分析.png', dpi=300) plt.show() print(f"{title}主周期: {main_scale}天") def full_analysis(self): """ 完整分析流程 """ # 1. 数据预处理 daily_data = self.preprocess_data() # 2. 季节性分解 plt.figure(figsize=(15, 10)) for i, col in enumerate(['流量(m3/s)', '含沙量(kg/m3)'], 1): plt.subplot(2, 1, i) result = self.seasonal_decomposition(daily_data[col]) result.plot() plt.title(f'{col}季节性分解') plt.tight_layout() plt.savefig('季节性分解.png', dpi=300) plt.show() # 3. 小波分析 flow_coef, flow_period = self.wavelet_analysis( daily_data['流量(m3/s)'], '流量' ) sediment_coef, sediment_period = self.wavelet_analysis( daily_data['含沙量(kg/m3)'], '含沙量' ) # 4. 交叉小波分析(流量与含沙量关系) self.cross_wavelet_analysis( daily_data['流量(m3/s)'], daily_data['含沙量(kg/m3)'], '流量-含沙量' ) return { 'flow_period': flow_period, 'sediment_period': sediment_period } def cross_wavelet_analysis(self, series1, series2, title='交叉分析'): """ 交叉小波分析 :param series1: 第一个时间序列 :param series2: 第二个时间序列 :param title: 分析标题 """ # 1. 计算小波变换 scales = np.arange(1, 365) wavelet = 'morl' coef1, _ = pywt.cwt(series1, scales, wavelet) coef2, _ = pywt.cwt(series2, scales, wavelet) # 2. 计算交叉小波谱 cross_spectrum = coef1 * np.conj(coef2) # 3. 可视化 plt.figure(figsize=(15, 8)) # 交叉小波谱实部 plt.subplot(2, 1, 1) plt.contourf(series1.index, scales, np.real(cross_spectrum), cmap='RdBu_r', levels=100) plt.colorbar(label='实部') plt.title(f'{title} - 交叉小波谱实部') plt.ylabel('时间尺度(天)') plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) # 交叉小波谱相位 plt.subplot(2, 1, 2) phase = np.angle(cross_spectrum) plt.contourf(series1.index, scales, phase, cmap='hsv', levels=100) plt.colorbar(label='相位(弧度)') plt.title(f'{title} - 交叉小波谱相位') plt.xlabel('日期') plt.ylabel('时间尺度(天)') plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m')) plt.tight_layout() plt.savefig(f'{title}_交叉小波分析.png', dpi=300) plt.show() # ====================== # 使用示例 # ====================== if __name__ == "__main__": # 1. 创建分析对象 analyzer = HydrologicalAnalysis() # 2. 加载数据(替换为实际文件路径) # analyzer.load_data('水文监测数据.csv') # 3. 生成模拟数据(实际使用时请注释掉) print("生成模拟数据...") dates = pd.date_range(start='2016-01-01', end='2021-12-31', freq='D') flow = np.sin(2 * np.pi * dates.dayofyear / 365) * 100 + 500 + np.random.normal(0, 50, len(dates)) sediment = np.cos(2 * np.pi * dates.dayofyear / 365) * 2 + 5 + np.random.normal(0, 1, len(dates)) analyzer.data = pd.DataFrame({ '水位(m)': np.random.uniform(40, 45, len(dates)), '流量(m3/s)': flow, '含沙量(kg/m3)': sediment }, index=dates) # 4. 执行完整分析 results = analyzer.full_analysis() print("\n分析结果摘要:") print(f"流量主周期: {results['flow_period']}天") print(f"含沙量主周期: {results['sediment_period']}天")
07-20
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值