Timesat提取物候信息并绘图

---

前言

Timesat 操作过程记录

---

一、准备数据

Timesat接受：

• img、dat格式的栅格数据；
• 单个像素点的text时序数据；

1. 将 Tiff 数据转成 dat 数据并生成 Filelist

先使用python根据shp裁剪原始tiff文件，然后使用IDL编程，将tiff转成dat文件，并生成Filelist.tex；用ENVI查看dat文件行列数，右键菜单，选View Metadata,可以得到rows=16，columns=34；
每一幅TIFF图像是以8字节的 IFH 开始的，IFH 的构成：

• Byte 0-1: 字节顺序标志位， 值为II（ASCII：4949）或者MM（ASCII：4D4D）。II表示小字节在前， 又称为little-endian。MM表示大字节在前，又称为big-endian。

使用UltraEdit打开一个tiff文件，看到“4949”，即“little-endian”

2. 使用python生成单个像素点text时序数据

Timesat 可以处理单个ASCII文件中的多个时序数据。该ASCII文件的第一行需要给出时序数据的年数nyear，每年的期数nptperyear，文件中时序数据的总数nts。第一行下面是时间序列y1;y2;…,yN,$N=nyear\ast nptperyear$,每个时序数据占一行，共有nts行时序数据，文件结构见下图。

def Tiff2text(lon,lat):
    column,row = Lonlat2Rowcol(lon, lat)#函数定义见前博文
    sif = []
    file_list = os.listdir(tiff_path)
    for file in file_list:
        if file.endswith('.tif'):
            dataset = gdal.Open(os.path.join(tiff_path, file))
            img = dataset.GetRasterBand(1).ReadAsArray()
            sif.append(img[row,column])
    sif = str(sif).replace(',',' ')
    fw = open(textfile_path + "haibeisif.txt", 'w+')
    fw.write(str(sif))  # 转化为str
    fw.close()

运行结果见下图，手工修改成 Timesat 所需格式即可

二、Timesat 打开时序数据

参考文献¹²³

1. 处理dat文件，提取物候信息

选好处理区域后，Load data，选择处理参数，有很多文献可以参考：

• 选好参数后，点“Settings”，保存成set文件，退出该窗口
• 点击“TSF_process”,选择刚才保存的set文件，系统会自动处理，完成后注意看下生成文件的路径，退出该窗口
• 点击“TSM_printseason”，在matlab窗口按提示要求输入信息，可得到物候信息的text文件

2. 处理text文件，提取物候信息

处理过程类似上面dat的处理过程，生成物候信息的text文件

三、Python 提取 text 中物候信息并绘图

1. Python 提取 text 中物候信息

from osgeo import gdal, osr, ogr
import numpy as np
import os
import shapefile
import csv
import matplotlib.pyplot as plt

def PhenoTextDraw():

    # 准备装入返青期、枯黄期、生长季长度
    Begtv = np.zeros(shape=(20, 17, 35), dtype=np.float32)
    Endtv = np.zeros(shape=(20, 17, 35), dtype=np.float32)
    Lengthv = np.zeros(shape=(20, 17, 35), dtype=np.float32)

    with open(textfile, 'r') as f:
        f_csv = csv.reader(f)
        for row in f_csv:
            for str in row:
                str = ','.join(str.split())
                str = str.split(',')
                # print(str)
                if str[0].isdigit():
                    season, row, column = int(str[2]), int(str[0]), int(str[1])
                    Begtv[season][row][column] = (float(str[3]) - ((season - 1) * 46)) * 8
                    Endtv[season][row][column] = (float(str[4]) - ((season - 1) * 46)) * 8
                    Lengthv[season][row][column] = float(str[5])*8

2. 绘图

# 进行contourf绘图
        fig, ax = plt.subplots()
        levels = np.arange(120, 200, 5)  # 对颜色渐进细致程度进行设置，其中50与200是色条显示的数据范围，5(也可0.05)是颜色显示的细致程度
        cs = ax.contourf(lon, lat,  Begtv[9][:][:], levels, cmap=plt.get_cmap('Spectral'))#Spectral,nipy_spectral,gist_rainbow
        # 添加colorbar
        cbar = fig.colorbar(cs, fraction=0.1, pad=0.15, shrink=0.9, anchor=(0.0, 0.3))  # 对colorbar的大小进行设置
        # 设置颜色条的刻度
        tick_locator = ticker.MaxNLocator(nbins=10)  # colorbar上的刻度值个数
        cbar.locator = tick_locator
        cbar.ax.tick_params(labelsize=12.5)
        # 设置颜色条的title
        cbar.ax.set_title('doy', fontsize=12.5)
        cbar.update_ticks()  # 显示colorbar的刻度值
        # 设置坐标刻度及间隔
        ax.set_xlim(np.min(lon[:]), np.max(lon[:]))
        ax.set_ylim(np.min(lat[:]), np.max(lat[:]))
        x_major_locator = ticker.MultipleLocator(0.5)  # 刻度间隔
        y_major_locator = ticker.MultipleLocator(0.1)
        ax = plt.gca()
        ax.xaxis.set_major_locator(x_major_locator)
        ax.yaxis.set_major_locator(y_major_locator)
        # 设置坐标轴标签
        font3 = {'family': 'Arial',
                 'weight': 'normal',
                 'size': 14,
                 }
        ax.set_xlabel('$\t{longitude}$ (°)', font3)
        ax.set_ylabel('$\t{latitude}$ (°)', font3)
        # 设置坐标刻度字体大小
        ax.tick_params(labelsize=14)
        labels = ax.get_xticklabels() + ax.get_yticklabels()
        [label.set_fontname('Arial') for label in labels]
        # 设置图像像素及大小
        plt.rcParams['figure.figsize'] = (6.0, 4.0)
        plt.rcParams['savefig.dpi'] = 300  # 图片像素
        plt.rcParams['figure.dpi'] = 300  # 分辨率
        plt.show()

---

总结

要提取区域物候信息时，需要准备该区域的dat或img数据
要提取单个点的物候信息时，需要准备text文件
使用Timesat提取出物候信息后，使用python读入并绘图分析

---

1. 参考http://the7.net/news/show-27332.html ↩︎

2. 参考https://developer.aliyun.com/article/751569 ↩︎

3. 参考https://blog.csdn.net/veraveraveravera/article/details/104043496 ↩︎