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 默认情况不支持中文,我们可以使用以下简单的方法来解决:
首先下载字体(注意系统):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
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轴标注
若有其它需要,请自行更改代码