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軸標註
若有其它需要,請自行更改代碼