问题:CSIF数据,月尺度上只有0.5°空间分辨率的
目标:将0.05°空间分辨率、4d时间分辨率的数据转成0.05° monthly
思路:
- 首先,读取路径下的文件
- 其次,将每个月的tif文件先归拢在一起(涉及到将DOY转为datetime
- 然后,基于每个月的tif文件求均值(注意使用np.nanmean,剔除nan值
- 最后,将均值输出为tif,每月一幅
代码如下:
# -*- coding: utf-8 -*-
"""
- 将每个月的tif归拢到一起,求平均值
- 每个月生成一幅tif输出
20240316, 15:45, YMJ
"""
import os
from datetime import datetime, timedelta
from osgeo import gdal
import numpy as np
# 文件路径
path = 'D:\\CSIF\\SIF_daily'
# 获取文件列表
files = os.listdir(path)
# 创建一个字典来存储每个月份的数据
monthly_data = {}
# 读取每个文件
for file in files:
if file.endswith('.tif'):
# 解析文件名以获取年份和DOY
year = int(file[0:4])
doy = int(file[4:7])
# 使用datetime模块计算具体日期
date = datetime(year=year, month=1, day=1) + timedelta(days=doy - 1)
# 获取月份
month = date.month
# 使用GDAL读取.tif文件
dataset = gdal.Open(os.path.join(path, file))
array = dataset.ReadAsArray()
# 将数据添加到月份字典中
key = (year, month)
if key not in monthly_data:
monthly_data[key] = [array]
else:
monthly_data[key].append(array)
# 计算每个月份的平均值
for key, data in monthly_data.items():
year, month = key
monthly_mean = np.nanmean(data, axis=0)
# print(monthly_mean[200,450])
monthly_mean[monthly_mean < -10000] = np.nan
# 创建一个新的.tif文件来保存平均值
driver = gdal.GetDriverByName('GTiff')
if month < 10:
month = str("0") + str(month)
else:
month = str(month)
output_path = f'D:\\CSIF\\SIF_monthly\\{year}{month}.tif'
print(output_path)
output_dataset = driver.Create(output_path, dataset.RasterXSize, dataset.RasterYSize, 1, gdal.GDT_Float32)
output_dataset.GetRasterBand(1).WriteArray(monthly_mean)
output_dataset.SetProjection(dataset.GetProjection())
output_dataset.SetGeoTransform(dataset.GetGeoTransform())
output_dataset.FlushCache()
output_dataset = None