Python地理数据处理 十七：植被物候提取和分析（Savitzky-Golay）

Savitzky-Golay滤波

1. 引子

import numpy as npimport matplotlib.pyplot as pltfrom scipy.signal import savgol_filternp.set_printoptions(precision=0, suppress=True)hist_equal = np.array([ 0, 0, 0, 0, 2, 7, 13, 23, 35, 46, 56, 64, 72, 78, 82, 86, 89, 92, 95, 99, 102, 105, 108, 111, 113, 115, 117, 119, 121, 122, 124, 126,128, 129, 131, 133, 134, 136, 137, 139, 140, 142, 143, 145, 146, 148, 149, 151, 153, 155, 157, 159, 161, 163, 166, 168, 170, 172, 174, 175, 177, 178, 179, 180,181, 182, 182, 183, 184, 185, 185, 186, 187, 188, 189, 189, 190, 191, 192, 193, 193, 194, 195, 196, 197, 197, 198, 199, 200, 200, 201, 202, 202, 203, 204, 204,205, 205, 206, 207, 207, 208, 208, 209, 209, 210, 210, 211, 211, 212, 213, 213, 214, 214, 215, 215, 216, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222,222, 223, 224, 225, 225, 226, 227, 227, 228, 229, 230, 230, 231, 232, 233, 234, 235, 236, 236, 237, 238, 239, 239, 240, 241, 241, 242, 242, 243, 243, 244, 244,245, 245, 246, 246, 247, 247, 248, 248, 249, 249, 250, 250, 250, 251, 251, 251, 251, 251, 251, 251, 251, 251, 251, 251, 251, 251, 252, 252, 252, 252, 252, 252,252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 253,253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 255])# window_length即窗口长度：取值为奇数且不能超过len(x)。它越大，则平滑效果越明显；越小，则更贴近原始曲线。# polyorder为多项式拟合的阶数：它越小，则平滑效果越明显；越大，则更贴近原始曲线hist_equal_savg = savgol_filter(hist_equal, 90, 4)x = np.arange(256)fig, axes = plt.subplots(1, 1)plt.plot(x, hist_equal, label = "hist")plt.plot(x, hist_equal_savg, label = "Savitzky Golay")plt.legend()plt.show()plt.show()

结果显示：

2. Savitzky-Golay滤波提取物候信息

import osimport mathimport numpy as npimport pandas as pdfrom osgeo import gdalfrom datetime import datetimefrom scipy.signal import savgol_filter

def jd_to_time(time):dt = datetime.strptime(time,'%Y%j').date()fmt = '%Y%m%d'return dt.strftime(fmt)def readTif(filename):dataset = gdal.Open(filename)if dataset == None:print(filename + "文件无法打开")return datasetdef writeTiff(img, img_geotrans, img_proj, path):if 'int8' in img.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in img.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(img.shape) == 3:img_bands, img_height, img_width = img.shapeelif len(img.shape) == 2:img = np.array([img])img_bands, img_height, img_width = img.shape# 创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(path, int(img_width), int(img_height), int(img_bands))if(dataset != None):dataset.SetGeoTransform(img_geotrans)# 写入仿射变换参数dataset.SetProjection(img_proj)#写入投影for i in range(img_bands):dataset.GetRasterBand(i+1).WriteArray(img[i])del datasetdef SG_filter(tifFolder, suffix, efolder):'''tifFolder:tif所在文件夹suffix：生成结果文件后缀'''# 获取文件夹内的文件名tifNameList = os.listdir(tifFolder)phe = []for p in tifNameList:pt = os.path.splitext(p)[0] # 获取tiff图像的文件名pt = int(pt[10:])phe.append(pt)# 获取第一张TIFF图像的信息：宽度、高度、投影、坐标系tifPath = tifFolder + "/" + tifNameList[0]dataset = readTif(tifPath) # 获取第一个tiff图像文件信息width = dataset.RasterXSize # 栅格矩阵的列数height = dataset.RasterYSize # 栅格矩阵的行数Tif_geotrans = dataset.GetGeoTransform()Tif_proj = dataset.GetProjection()Tif_data = dataset.ReadAsArray(0, 0, width, height)Tif_datas = np.zeros((len(tifNameList),Tif_data.shape[0], Tif_data.shape[1])) # 多张TIFF数据保存到 Tif_datas中for i in range(len(tifNameList)):tifPath = tifFolder + "/" +tifNameList[i]dataset = readTif(tifPath)Tif_data = dataset.ReadAsArray(0, 0, width, height)Tif_datas[i] = Tif_data# 切换维度:宽、高、图像个数# [0,1,2]->[1,0,2]->[1,2,0]类似于汉诺塔[宽，高，46]Tif_datas = Tif_datas.swapaxes(1, 0)Tif_datas = Tif_datas.swapaxes(1, 2)#定义空集合SGfilter = np.zeros(Tif_datas.shape)sos = np.zeros(Tif_data.shape)eos = np.zeros(Tif_data.shape)# 读取每张图像的信息,所有图像相同位置的计算for i in range(Tif_datas.shape[0]): # 行for j in range(Tif_datas.shape[1]): # 列value = Tif_datas[i][j]m = 0for k in value:if math.isnan(k) == True: # 空值判断，若有空值，则不计算该像素（全部图像的该像素）m = m+1if m == 0:sif = Tif_datas[i][j]sif = pd.DataFrame(sif)sif[sif<0] = 0sif[sif>32766] = 0sif = sif * 0.0001sif = sif.values.flatten()SGfilter[i][j] = savgol_filter(sif, window_length = 9, polyorder = 3) # 拟合阶数为3ysg = SGfilter[i][j]_y2 = list(ysg)maxy = max(_y2)miny = min(_y2)th = 0.2 # 设置阈值amplitude = maxy - miny # 振幅thresh = amplitude * th + minynewnums = list(filter(lambda x: x >= thresh, _y2))r = newnums[0]r2 = newnums[-1]index1 = _y2.index(r)index2 = _y2.index(r2)sos[i][j] = phe[index1]eos[i][j] = phe[index2]sos[np.where(sos == 1)] = np.NANsos[np.where(eos == 1)] = np.NANsavePaths = efolder + "sos.tif"savePathe = efolder + "eos.tif"writeTiff(sos, Tif_geotrans, Tif_proj, savePaths)writeTiff(eos, Tif_geotrans, Tif_proj, savePathe)# 维度还原为原来的维度SGfilter = SGfilter.swapaxes(1,0)SGfilter = SGfilter.swapaxes(0,2)for i in range(SGfilter.shape[0]):savePath = os.path.splitext(tifNameList[i])[0] + suffix + ".tif"savePath = efolder + savePathwriteTiff(SGfilter[i], Tif_geotrans, Tif_proj, savePath)

数据导入：

np.random.seed(10)import timestart = time.time()SG_filter(r"E:/dataset/sif/sif_SG", "_SGFilter", "E:/dataset/sif/sifphe/")end = time.time()print(end - start)

结果显示：