以OMI NO2 L3 2013-2019年的数据为例，为了于2020年的NO2排放进行对比，因此需要求出7年每天的平均值作为基线，代码如下：
"""
计算2013-2019年每天的OMI NO2 L3数据的平均值
"""
import glob
import os
import rasterio
import numpy as np
out_dir = r'G:\DATA\OMINO2_L3\2013_2019_daily_mean'
file_list = glob.glob(r'G:\DATA\OMINO2_L3\tif\*.tif')
# obtain raste meta data
with rasterio.open(file_list[0]) as src:
meta = src.meta
meta.update(dtype=rasterio.float32)
# make a key word list from input files
key = []
for file in file_list:
filename = file.split('\\')[-1]
date = filename.split('_')[1]
date = date.split('.')[0]
key.append(date)
# key_list = set(key)
key_list = key[0:365]
print('The keys are: ')
print(key_list)
# read raster data function
def read_raster(raster):
with rasterio.open(raster) as src:
return(src.read(1))
# select the files match each key and calculate mean and write to disk
file_dict = {}
for key in key_list:
print('The key word is ------> ' + key)
file_dict[key] = []
for file in file_list:
if key in file:
print('Found files ... ' + file)
file_dict[key].append(file)
# calculate mean
for k, v in file_dict.items():
array_list = [read_raster(x) for x in v]
array_out = np.nanmean(array_list, axis=0) # ignoring the NA values
# write to disk
out_file = os.path.join(out_dir, k + '.tif')
with rasterio.open(out_file, 'w', **meta) as dst:
dst.write(array_out.astype(rasterio.float32), 1)
print('Writing files to ... ' + out_file)
print('... ... ... COMPLETED ... ... ...')
