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

Savitzky-Golay滤波

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

1. 引子

import numpy as np
import matplotlib.pyplot as plt
from 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 os
import math
import numpy as np
import pandas as pd
from osgeo import gdal
from datetime import datetime
from 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 time
start = time.time()SG_filter(r"E:/dataset/sif/2008sif_SG", "_SGFilter", "E:/dataset/sif/2008sifphe/")end = time.time()
print(end - start)

结果显示：

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

长时间序列遥感数据植被物候提取/遥感数据产品分析暨MODIS NDVILAI多年产品数据批处理分析/Python长时间序列遥感数据处理及在全球变化、物候提取、植被变绿与固碳分析、生物量估算与趋势分析
基于MATLAB长时间序列遥感数据植被物候提取与分析 1.本课程基于matlab语言 2.提供所有代码 3.以实践案例为课程内容主线,原理与操作相结合 4.根据讲解内容,布置作业,巩固所学内容及拓展在 ...
Python读写矢量数据（2）矢量数据写入（属性数据）——Python地理数据处理学习分享
这一节主要介绍矢量数据的写入(只有属性数据,无几何),如果有读者没有读取的基础建议先看一下上一篇文章,需要对矢量数据读取有一定的了解才能继续学习本节.在这里我们用到的数据仍为goble文件夹下的数据, ...
《Python 地理数据处理》by Chris Garrard
科研小白从头开始学习用Python处理栅格数据.矢量数据等. 下面就跟我一起开启<Python 地理数据处理>的学习之旅吧!!! 本书主要利用Python + GDAL 等相关库进行地理空 ...
《Python地理数据处理》——导读
前言本书可以帮助你学习使用地理空间数据的基础知识,主要是使用GDAL / OGR.当然,还有其他选择,但其中一些都是建立在GDAL的基础之上,所以如果你理解了本书中的内容,就可以很轻松地学习其他知识 ...
python地理数据处理下载_python-doc/将Python用于地理空间数据处理.md at master · zhuxinyizhizun/python-doc · GitHub...
毫无疑问,Python是当今最流行,最通用的编程语言之一.这有很多种强有力的原因,但在我看来,最重要的是:开源定义,语法简单,包括电池的理念(batteries included philosophy ...
Python地理数据处理二：Python基础知识
目录 1.编写执行代码 2.脚本结构 3.变量 4.数据类型 4.1 布尔型 4.2 数值型 4.3 字符串 4.3.1 连接字符串 4.3.2 转义字符 4.4 列表和元组 4.5 集合 4.6 字 ...
Python地理数据处理十一：空间参照系统(SRS)
目录 1. 空间参照 2. OSR空间参考 2.1 空间参考对象 2.2 创建控件参考对象 2.3 分配SRS 2.4 重投影 2.5 更改基准 2.6 重投影整个图层 3. 使用pyproj空间参考 ...
python地理数据处理下载_【实用书】Python地理数据处理，362页pdf，Geoprocessing with Python...
关于本书我编写了<Geoprocessing for Python> 来帮助您学习处理地理空间数据的基础知识,主要使用GDAL/OGR.当然,还有其他的选择,但是其中一些是在GDAL之上 ...
Python地理数据处理库shapely支持函数总结
Shapely是一个Python库,用于操作和分析笛卡尔坐标系中的几何对象. 本文通过部分示例介绍了空间处理库Shape的部分概念与操作函数. 官方文档:https://shapely.readthe ...

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

Savitzky-Golay滤波

1. 引子

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

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

最新文章

热门文章

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

Savitzky-Golay滤波

1. 引子

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

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

最新文章

热门文章

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

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