MODIS系列之NDVI(MOD13Q1)七:時間序列S-G濾波之Python

 時間序列S-G濾波之Python

根據上上篇博文(MODIS系列之NDVI(MOD13Q1)五:NDVI處理流程 )做出的NDVI。我們求NDVI時間序列圖,但該NDVI時序圖為地表各土地類型綜合的NDVI時序圖。(詳情同樣參考該系列五博文的文底)

建議:大家應該也能發現從網上粘貼的代碼,大部分在各自實際運行中會出現報錯,不能運行。這其中有代碼本身的錯誤,但也不乏運行環境的欠缺、誤操作、電腦自身特點等原因。本博客的所有代碼都經過實際運行再上傳,哪怕比較熟悉的代碼,再上傳前都會儘可能實際運行。目的便是給大家減少參考困難及錯誤信息。因為己所不欲勿施於人么。不足之處,還請見諒,並留下建議。

1. 完整代碼如下(實際運行成功):

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




 運行結果如下:

 

 

 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 默認情況不支持中文,我們可以使用以下簡單的方法來解決:

首先下載字體(注意系統)://www.fontpalace.com/font-details/SimHei/

SimHei.ttf 文件放在當前執行的代碼文件中:

2.4.1安裝SimHei.ttf 文件

打開上文鏈接(該系列操作電腦比較慢,如果打不開,請重複操作,國外網站都這樣。我大天國太強)

 

 

 

 

 

 如果該操作不行。請用網盤鏈接

鏈接://pan.baidu.com/s/1h_U37kA_dzEIjW4Vy9tX9Q
提取碼:7tm1

如圖:

 

 

 2.5

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") 

2.5.1

*- 為實線樣式。

作為線性圖的替代,可以通過向 plot() 函數添加格式字符串來顯示離散值。 可以使用以下格式化字符。

字符 描述
'-' 實線樣式
'--' 短橫線樣式
'-.' 點劃線樣式
':' 虛線樣式
'.' 點標記
',' 像素標記
'o' 圓標記
'v' 倒三角標記
'^' 正三角標記
'<' 左三角標記
'>' 右三角標記
'1' 下箭頭標記
'2' 上箭頭標記
'3' 左箭頭標記
'4' 右箭頭標記
's' 正方形標記
'p' 五邊形標記
'*' 星形標記
'h' 六邊形標記 1
'H' 六邊形標記 2
'+' 加號標記
'x' X 標記
'D' 菱形標記
'd' 窄菱形標記
'|' 豎直線標記
'_' 水平線標記

2.5.2

label=’NDVI’        label為標註

2.5.3

color=’red’

線條顏色可以用英文或縮寫

以下是顏色的縮寫:

字符 顏色
'b' 藍色
'g' 綠色
'r' 紅色
'c' 青色
'm' 品紅色
'y' 黃色
'k' 黑色
'w' 白色

2.5.4

第二行為標題設置  “2010年河南省3、4、5月冬小麥NDVI” 為本操作標題

fontproperties=zhfont1 為 fontproperties 設置中文顯示,fontsize 設置字體大小

2.5.5

第三行為設置X軸標註和中文顯示和字體大小

2.5.6

第四行為設置Y軸標註

 

若有其它需要,請自行更改代碼

 

Tags: