目的： 現有2001-2020年 第57天—第273天的資料(間隔八天) 想要得到每一天（129 - 274）資料的20年平均 並存儲到新的tif裡面

資料命名方式如下：SIF_YYYYDDD_new.tif （YYYY為年， DDD為天）ex：SIF_2001057_new.tif 

                                                                                                                          SIF_2019185_new.tif

from osgeo import gdal
import numpy as np

file_format = "GTiff"
driver = gdal.GetDriverByName(file_format)

for day in range(57, 274, 8):
    mean = []
    SIF_mean = []
    day = "%03d" % day   # 對於57 自動補為057
    for year in range(2001, 2021):
        input_file = 'F:/paper_graduate/result_tif/new/SIF_{year1}{day1}_new.tif'.format(year1=year, day1=day)
        data = gdal.Open(input_file)
        data_sif = data.GetRasterBand(1)
        data_sif_array = data_sif.ReadAsArray()

        data_sif_array = np.where(data_sif_array >= 0, data_sif_array, np.NAN)
        mean.append(data_sif_array)
    SIF_mean = np.nanmean(mean, axis=0)
    SIF_mean[np.isnan(SIF_mean)] = 255

    mean_filename = "F:/paper_graduate/result_tif/new/MeanSIF_{dd}_new.tif".format(dd=day)
    New_tif = driver.Create(mean_filename, 300, 300, 1, gdal.GDT_Float32)

    # 讀取之前tif資訊 作為新生成tif的座標系統
    mask = gdal.Open('F:/paper_graduate/result_tif/new/SIF_2001057_new.tif')
    # mask_proj = mask.GetProjection()
    # New_tif.SetProjection(mask_proj)
    mask_trans = mask.GetGeoTransform()
    New_tif.SetGeoTransform(mask_trans)

    New_tif.GetRasterBand(1).WriteArray(SIF_mean)
    dataset = None