使用Python求多副tif影像平均值 并输出到新的tif中

本文主要是介绍使用Python求多副tif影像平均值 并输出到新的tif中，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

目的： 现有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

这篇关于使用Python求多副tif影像平均值 并输出到新的tif中的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！