MODIS系列之NDVI（MOD13Q1）七：时间序列S-G滤波之Python

 时间序列S-G滤波之Python

根据上上篇博文（MODIS系列之NDVI（MOD13Q1）五：NDVI处理流程 ）做出的NDVI。我们求NDVI时间序列图，但该NDVI时序图为地表各土地类型综合的NDVI时序图。（详情同样参考该系列五博文的文底）

建议：大家应该也能发现从网上粘贴的代码，大部分在各自实际运行中会出现报错，不能运行。这其中有代码本身的错误，但也不乏运行环境的欠缺、误操作、电脑自身特点等原因。本博客的所有代码都经过实际运行再上传，哪怕比较熟悉的代码，再上传前都会尽可能实际运行。目的便是给大家减少参考困难及错误信息。因为己所不欲勿施于人么。不足之处，还请见谅，并留下建议。

1. 完整代码如下（实际运行成功）：

import matplotlib.pyplot as pltfrom osgeo import gdalfrom gdalconst import *import matplotlib import numpyimport timeimport osimport os.pathdef Gettiff(getpath):    tiff=[]    os.system("dir "+getpath+"\\"+"*.tif /b /s>tiff.TXT")    tifftxt = open("tiff.TXT").readlines()    for i in tifftxt:        tiff.append(i.strip())    return tiffdef greater(data,dt,r,c):    count=0    for i in range(r):        for j in range(c):            if data[i][j]!=dt:                count=count+1            else:                continue    return countdef getsum(data,c,r,nodata):    sum1=0    for i in range(c):        for j in range(r):            if data[i][j]==nodata:                continue            else:                sum1=sum1+data[i][j]    return sum1def write_img(savepath,filename,im_proj,im_geotrans,im_data):    #gdal数据类型包括    #gdal.GDT_Byte,     #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,    #gdal.GDT_Float32, gdal.GDT_Float64    a = os.path.exists(savepath)    if a== False:        os.mkdir(savepath)    #判断栅格数据的数据类型    if ‘int8‘ in im_data.dtype.name:        datatype = gdal.GDT_Byte    elif ‘int16‘ in im_data.dtype.name:        datatype = gdal.GDT_UInt16    else:        datatype = gdal.GDT_Float32    #判读数组维数    if len(im_data.shape) == 3:        im_bands, im_height, im_width = im_data.shape    else:        im_bands, (im_height, im_width) = 1,im_data.shape     #创建文件    driver = gdal.GetDriverByName("GTiff")            #数据类型必须有，因为要计算需要多大内存空间    dataset = driver.Create(savepath+"\\"+filename, im_width, im_height, im_bands, datatype)    dataset.SetGeoTransform(im_geotrans)              #写入仿射变换参数    dataset.SetProjection(im_proj)                    #写入投影    if im_bands == 1:        dataset.GetRasterBand(1).WriteArray(im_data)  #写入数组数据    else:        for i in range(im_bands):            dataset.GetRasterBand(i+1).WriteArray(im_data[i])    del datasetdef GetValue(TIFF):    x=list()    y=list()    Geo=[]    for  i,b in zip(TIFF,range(len(TIFF))):        filePath=r"D:\ArcMapData\tif\.tif"   #该路径不是读取和保存路径。是时序图X轴数值定义        end_time=os.path.split(i)[-1][5:8]   #截取字符串，截取该文件夹下.tif文件的文件名的部分作为时序图X轴数值        endtime=int(end_time)                 ds = gdal.Open(i,GA_ReadOnly)        rows = ds.RasterYSize        cols = ds.RasterXSize        band = ds.GetRasterBand(1)        Nodata=band.GetNoDataValue()        data = band.ReadAsArray(0, 0, cols, rows)        sum1=getsum(data,rows,cols,Nodata)        count = greater(data,Nodata,rows,cols)        ignore0_pixel = sum1/count        x.append(ignore0_pixel)        y.append(endtime)        del ds    return x,y        # print(rows,cols,im_proj)        # print("\n\n\n")        # print(data)def showtiff(listx,listy,listx1,listy1):    plt.plot(listy,listx,".-b")    plt.plot(listy1,listx1,"*-r")    plt.tick_params(labelsize=10)    plt.xticks(fontsize = 8)    plt.show()    if __name__==‘__main__‘:        #for i in range(2000,2018):    getpath=r"D:\ArcMapData\xiaomaiNDVI"    tiff=Gettiff(getpath)        zhfont1 = matplotlib.font_manager.FontProperties(fname="微软vista黑体.ttf")        x3,y3=GetValue(tiff)    plt.plot(y3,x3,"*-",label=‘NDVI‘,color=‘red‘)    plt.title("2010年河南省3、4、5月冬小麦NDVI",fontproperties=zhfont1)    plt.xlabel("天数", fontproperties=zhfont1)     plt.ylabel("NDVI")           plt.legend(ncol=1, mode="None")    plt.show()

 运行结果如下：

 Figure界面介绍：

 重置原始视图

 返回上一个视图

 前进到下一个视图

 用鼠标左键平移轴，用鼠标右键缩放

 缩放到矩形

 配置子批次

 保存图形（.png格式）

2. 部分代码解读（为解读的代码通常不用更改，若另有改动需求，请便）：

2.1 安装matplotlib module

import matplotlib

本人在运行代码遇到的第一个问题就是这个，相信大家在运行的过程中也可能会遇到。python交互式命令行页面会报出  无 matplotlib module 类型 。那就安装 matplotlib module 。具体安装步骤请参考博文：Python之 module安装

 2.2 该部分为根据各自需要获取X轴数值

filePath=r"D:\ArcMapData\tif\.tif"    #该路径不是读取和保存路径。是时序图X轴数值定义
        end_time=os.path.split(i)[-1][5:8]    #截取字符串，截取该文件夹下.tif文件的文件名的部分作为时序图X轴数值
        endtime=int(end_time)

2.3 数据读取路径（替换自己的数据读取路径）

getpath=r"D:\ArcMapData\xiaomaiNDVI"

2.4 重点介绍 (图形中文显示)

Matplotlib 默认情况不支持中文，我们可以使用以下简单的方法来解决：

首先下载字体（注意系统）：https://www.fontpalace.com/font-details/SimHei/

SimHei.ttf 文件放在当前执行的代码文件中：

2.4.1安装SimHei.ttf 文件

打开上文链接（该系列操作电脑比较慢，如果打不开，请重复操作，国外网站都这样。我大天国太强）

 如果该操作不行。请用网盘链接

链接：https://pan.baidu.com/s/1h_U37kA_dzEIjW4Vy9tX9Q
提取码：7tm1

如图：

 2.5

plt.plot(y3,x3,"*-",label=‘NDVI‘,color=‘red‘)
    plt.title("2010年河南省3、4、5月冬小麦NDVI",fontproperties=zhfont1)
    plt.xlabel("天数", fontproperties=zhfont1) 
    plt.ylabel("NDVI")

2.5.1

*- 为实线样式。

作为线性图的替代，可以通过向 plot() 函数添加格式字符串来显示离散值。 可以使用以下格式化字符。

2.5.2

label=‘NDVI‘        label为标注

2.5.3

color=‘red‘

线条颜色可以用英文或缩写

以下是颜色的缩写：

2.5.4

第二行为标题设置  "2010年河南省3、4、5月冬小麦NDVI" 为本操作标题

fontproperties=zhfont1 为 fontproperties 设置中文显示，fontsize 设置字体大小

2.5.5

第三行为设置X轴标注和中文显示和字体大小

2.5.6

第四行为设置Y轴标注

若有其它需要，请自行更改代码