Python使用S

2023-04-21 14:31| 来源: 网络整理| 查看: 265

S-G介绍：

Savitzky-Golay（S-G）滤波是一种数字信号处理技术，主要用于平滑和去噪数据。它通过对数据点进行局部多项式拟合来估算平滑值。与其他平滑方法（如移动平均和高斯滤波）相比，S-G滤波在保留信号特征（如峰值、宽度和高度）方面具有优势。

主要参数：

1、窗口大小（window_size）：窗口大小是用于局部拟合的数据点的数量。通常选择奇数，以便在中心点左右有相等数量的数据点。窗口大小越大，平滑效果越明显，但过大的窗口可能会导致信号失真。

2、多项式阶数（polynomial_order）：多项式阶数决定了用于拟合数据的多项式的复杂程度。较低的阶数（如1或2）适用于简单的数据，较高的阶数（如3或4）适用于复杂的数据。阶数过高可能导致过拟合。

优点：

1、保留信号特征：S-G滤波器能够保留信号的峰值、宽度和高度。

2、去噪能力：通过局部多项式拟合，S-G滤波器可以有效地去除噪声。

3、可调参数：通过调整窗口大小和多项式阶数，可以根据具体应用调整滤波器的性能。

1、导入需要的Python模块：

import os import numpy as np import rasterio from scipy.signal import savgol_filter import matplotlib.pyplot as plt

2、定义一个需要处理的文件夹路径。

# 获取文件夹内的所有文件 all_files = os.listdir(ndvi_folder) ndvi_image_files = [os.path.join(ndvi_folder, file) for file in all_files if file.endswith('.tif')]

3、读取第1幅影像ndvi_image_files[0]以获取影像的高度和宽度值，并将影像属性信息赋值给profile变量：

with rasterio.open(ndvi_image_files[0]) as src: height, width = src.shape profile = src.profile

4、初始化一个包含所有影像数据的NumPy数组ndvi_data，并将每个影像中的第1个波段数据读取到此数组中：

ndvi_data = np.empty((len(ndvi_image_files), height, width)) # 将影像文件读取到NumPy数组中 for i, image_file in enumerate(ndvi_image_files): with rasterio.open(image_file) as src: ndvi_data[i] = src.read(1)

5、绘制原始的NDVI数据的趋势图：

mean_ndvi = np.mean(ndvi_data, axis=(1, 2)) plt.plot(mean_ndvi) plt.title('NDVI') plt.xlabel('TIME') plt.ylabel('NDVI') plt.show()

6、应用Savitzky-Golay滤波平滑处理，将平滑后的NDVI数据存储在变量smoothed_ndvi_data中：

indow_size = 7 # 窗口大小，可以根据需要调整 polynomial_order = 3 # 多项式阶数，可以根据需要调整 smoothed_ndvi_data = savgol_filter(ndvi_data, window_size, polynomial_order, axis=0)

7、绘制平滑后的数据

mean_smoothed_ndvi = np.mean(smoothed_ndvi_data, axis=(1, 2)) plt.plot(mean_smoothed_ndvi) plt.title('smoothed_ndvi') plt.xlabel('time') plt.ylabel('NDVI') plt.show()

8、将平滑后的数据另存为新的影像文件

# 将平滑后的NDVI数据保存为新的影像文件 output_folder = os.path.join(ndvi_folder, 'smoothed_ndvi') os.makedirs(output_folder, exist_ok=True) for i, image_file in enumerate(ndvi_image_files): output_file = os.path.join(output_folder, f'smoothed_{os.path.basename(image_file)}') with rasterio.open(output_file, 'w', **profile) as dst: dst.write(smoothed_ndvi_data[i], 1)

效果大概如下：这里参数需要自己修改一下，，我这个就做一个参考。

注：window_size 只能为奇数！！！！

完整代码如下：

import os import numpy as np import rasterio from scipy.signal import savgol_filter import matplotlib.pyplot as plt #z ndvi_folder = 'D:/SY/SJJ/QLSNDVI_Month' # 获取文件夹内的所有文件 all_files = os.listdir(ndvi_folder) ndvi_image_files = [os.path.join(ndvi_folder, file) for file in all_files if file.endswith('.tif')] # 读取第一幅影像以获取维度信息 with rasterio.open(ndvi_image_files[0]) as src: height, width = src.shape profile = src.profile # 初始化一个包含所有影像数据的NumPy数组 ndvi_data = np.empty((len(ndvi_image_files), height, width)) # 将影像文件读取到NumPy数组中 for i, image_file in enumerate(ndvi_image_files): with rasterio.open(image_file) as src: ndvi_data[i] = src.read(1) # 绘制原始NDVI数据的趋势图 mean_ndvi = np.mean(ndvi_data, axis=(1, 2)) plt.plot(mean_ndvi) plt.title('NDVI_Original') plt.xlabel('TIME') plt.ylabel('NDVI') plt.show() # 应用Savitzky-Golay滤波 window_size = 7 # 窗口大小，可以根据需要调整 polynomial_order = 3 # 多项式阶数，可以根据需要调整 smoothed_ndvi_data = savgol_filter(ndvi_data, window_size, polynomial_order, axis=0) # 绘制平滑后的NDVI数据的趋势图 mean_smoothed_ndvi = np.mean(smoothed_ndvi_data, axis=(1, 2)) plt.plot(mean_smoothed_ndvi) plt.title('smoothed_ndvi') plt.xlabel('time') plt.ylabel('NDVI') plt.show() # 将平滑后的NDVI数据保存为新的影像文件 output_folder = os.path.join(ndvi_folder, 'smoothed_ndvi') os.makedirs(output_folder, exist_ok=True) for i, image_file in enumerate(ndvi_image_files): output_file = os.path.join(output_folder, f'smoothed_{os.path.basename(image_file)}') with rasterio.open(output_file, 'w', **profile) as dst: dst.write(smoothed_ndvi_data[i], 1)

【本文地址】

公司简介

联系我们