MODIS系列之NDVI(MOD13Q1)七:时间序列S-G滤波之Python

Posted 青灯常伴古佛

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了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轴标注

 

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

 

以上是关于MODIS系列之NDVI(MOD13Q1)七:时间序列S-G滤波之Python的主要内容,如果未能解决你的问题,请参考以下文章

MODIS系列之NDVI(MOD13Q1)九:单一土地类型NDVI及时序图(以耕地类型为例)

MODIS系列之NDVI(MOD13Q1)九:单一土地类型NDVI及时序图(以耕地类型为例)

ArcGIS Python实现Modis NDVI批量化月最大合成

R:MODIS 原数据下载

Python应用实战案例-Python使用MODIS数据实现温度植被干旱指数TVDI的计算

Python应用实战案例-Python使用MODIS数据实现温度植被干旱指数TVDI的计算