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

import matplotlib.pyplot as plt
from osgeo import gdal
from gdalconst import *
import matplotlib 
import numpy
import time
import os
import os.path



def 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 tiff

def 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 count
def 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 sum1


def 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 dataset

def 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()

import matplotlib 

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

getpath=r"D:\ArcMapData\xiaomaiNDVI"

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