使用TIMESAT软件遇到的各种问题

在做时间序列重建工作中，由于实验中必须用到A-G滤波和D-L滤波，这两种经典方法在python中实现很难，资料也很少，所以只能用TIMESAT软件封装好的函数，在此过程中遇到了非常多的坑和问题…记录一下。


正常的操作流程可以参考文章TIMSAT拟合LST数据，NDVI处理方式类似。


1.首先是TIMESAT的下载与安装，基于MATLAB2018b环境，整体流程都参考了青灯常伴古佛大佬的文章时间序列滤波之软件TIMESAT。
由于MATLAB的破解安装使用iso镜像文件，两个iso文件的安装逻辑跟插入光盘一样，使用disk1安装会导致安装到一半提示需要插入disk2。
最后使用虚拟光驱软件DAEMON-TOOLS才正常安装，参考了文章Windows版iso格式MATLAB安装文件使用方法说明


2.成功打开TIMESAT软件后界面如下所示：
其中TSM_ imageview可以加载单张或多张图像预览其效果，No of rows in image代表图像的行数，columns代表列数，且只能输入单波段图像，正确输入后点击Draw可以加载出单张图像。

但是如果行列数错误或者图像本身的格式、编码方式有问题都会导致图像加载错误，由于TIMESAT对于图像格式的严苛要求，导致我很多Python处理后的TIF格式图像在ENVI5.3中能正确加载显示，而在TIMESAT中却显示为乱码或无法对齐。


3.Python中批处理输出TIF格式图像和ENVI中的格式转换：
由于我的图像都在python中做过相关处理，需要以数组的形式存储到TIF格式文件中，如果使用Opencv之类的常规图像处理库来读取和存储，图像将会保存为RGB三通道而非单通道灰度图像，投影信息也会丢失，完全无法被TIMESAT识别。
所以此处只能使用GIS相关库gdal来读取原TIF文件并写入相关波段，参考文章Python下遥感图像文件的读写，此处将波段数限制为1，批量输出到目标文件夹。gdal库的安装和使用也存在相当多的问题，需要与系统适配的版本才能正常import。

#保存单波段tif文件，下一步还需要在ENVI中做格式转换，其他相关库的引用此处未全部列举
from osgeo import gdal
def FileSave():
    
    lists = list(range(230))
    def read(i):
        despath = 'test/NDVI_%d.tif'%i
        outpath = 'NDVI/NDVI_%d.tif'%i
        dataset = gdal.Open(despath)
        if dataset == None:
            print(despath + "文件无法打开")
        im_width = dataset.RasterXSize # 栅格矩阵的列数
        im_height = dataset.RasterYSize # 栅格矩阵的行数
        im_bands = dataset.RasterCount # 波段数
        im_geotrans = dataset.GetGeoTransform() # 获取仿射矩阵信息
        im_pro = dataset.GetProjection()
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height) # 获取数据
        
        driver = gdal.GetDriverByName("GTiff")
        datatype = gdal.GDT_Int16#TIMESAT需要规定usigned8位或signed16位图
        dataset = driver.Create(outpath , im_width, im_height, 1, datatype,
        options=["INTERLEAVE=BAND"])
        if (dataset != None and im_geotrans != None and im_pro != None):
            dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
            dataset.SetProjection(im_pro) # 写入投影
        
        dataset.GetRasterBand(1).WriteArray(im_data[0])
        del dataset

    """多进程读取写法"""  
    start = time.time()
    pool = ThreadPool()
    results = pool.map(read, lists)
    pool.close()
    pool.join()
    end = time.time()
    print("总共用时{}秒".format((end - start)))

输出后的单波段TIF图像能被ENVI读取，但仍然无法被TIMESAT识别，需要在ENVI中使用批处理插件ENVI Raster Processing Batch Tools统一转换为dat格式。具体操作参考文章TIMSAT拟合LST数据，插件下载地址ENVI扩展工具。
在ENVI5.3 64位版本中打开所有需要处理的文件，使用Convert Interleave Batch批量格式转换成.dat文件，最后在TIMESAT中终于可以识别了。

4.时间序列图像的输入与处理：
在TIMESAT中处理时间序列主要使用TSM_GUI功能，需要输入一个list文件代表所有时间序列图像。此文件为txt格式，第一行写入总的文件数量，其余行为所有文件的地址。
在输入设置中，No of years代表一年有多少幅图像，后一个代表一共多少年。其余参数与前面提到的图像行列数、位数一致，但需要在Rows to process中限定需要处理的区域，此处只处理800-805这一块小区域，如果不做限定可能会导致计算量过大程序崩溃。

数据加载完成后，时间序列会显示在界面上，在滤波方法一栏就可以选择相关滤波函数了，用Output可以输出某像素的时间序列到txt文件中，最后可以在python中使用numpy读取。

log = np.loadtxt("logistic.txt")
log = log/10000
raw = ndvi[800][802]

pic = {
    'raw':raw,
    'log':log
}
pic = pd.DataFrame(pic)
pic.plot()