### 绘制 NDVI 时间序列谐波分析和平滑前后对比图 在 Python 中绘制 NDVI 时间序列谐波分析结果和平滑前后的对比图，可以通过 `matplotlib` 和自定义或第三方实现的 HANTS 算法完成。以下是一个完整的实现流程，涵盖数据准备、HANTS 平滑处理和可视化展示。 #### 数据准备 假设已经获取了 NDVI 时间序列数据，并存储在一个 NumPy 数组中。时间序列数据可以是按天或按旬采样的 NDVI 值，其中可能包含缺失值（如云覆盖导致的异常点）。数据格式通常为一维数组，时间点均匀或非均匀分布。 #### HANTS 算法实现 以下是一个简化版本的 HANTS 算法实现，适用于 NDVI 时间序列的平滑处理： ```python import numpy as np from scipy.fft import fft, ifft def hants_filter(data, n_frequencies=3, missing_value=-9999): n = len(data) t = np.arange(n) mask = data != missing_value t_valid = t[mask] y_valid = data[mask] # 构建傅里叶基函数 A = np.zeros((len(t_valid), 2 * n_frequencies + 1)) A[:, 0] = 1 for k in range(1, n_frequencies + 1): A[:, 2*k-1] = np.cos(2 * np.pi * k * t_valid / n) A[:, 2*k] = np.sin(2 * np.pi * k * t_valid / n) # 最小二乘拟合 coeff, _, _, _ = np.linalg.lstsq(A, y_valid, rcond=None) # 重构时间序列 A_full = np.zeros((n, 2 * n_frequencies + 1)) A_full[:, 0] = 1 for k in range(1, n_frequencies + 1): A_full[:, 2*k-1] = np.cos(2 * np.pi * k * t / n) A_full[:, 2*k] = np.sin(2 * np.pi * k * t / n) fitted = np.dot(A_full, coeff) return fitted ``` #### 数据可视化 使用 `matplotlib` 可以绘制原始 NDVI 时间序列和平滑后的结果对比图： ```python import matplotlib.pyplot as plt # 假设 ndvi_data 是原始 NDVI 时间序列数据，包含缺失值 ndvi_data = np.array([...]) # 替换为实际数据 smoothed_data = hants_filter(ndvi_data) # 时间轴（例如按旬采样） time_axis = np.arange(len(ndvi_data)) # 绘图 plt.figure(figsize=(12, 6)) plt.plot(time_axis, ndvi_data, label='Original NDVI', color='gray', linestyle='--') plt.plot(time_axis, smoothed_data, label='Smoothed NDVI (HANTS)', color='blue', linewidth=2) plt.xlabel('Time Index') plt.ylabel('NDVI Value') plt.title('NDVI Time Series Before and After HANTS Smoothing') plt.legend() plt.grid(True) plt.show() ``` #### 可视化效果优化 为了增强对比效果，可以在图中加入缺失值标记或异常点标注。此外，如果时间序列数据具有明确的日期信息，可以将时间轴转换为日期格式，使图表更具可读性。 ###