matlab/arcpy（貌似会死机）对ERA5（tif文件）每小时数据求年平均值

       写在前头

        这是本人的第一篇CSDN博客，如有什么做的不太好的地方，敬请批评指正

        这代码应该挺基础的，发出来也只是：

        ①记录一下自己的学习情况进展，希望大佬嘴下留情，不喜勿喷。

        ②翻阅了csdn，发现似乎并没有一个针对地理信息tif数据进行批量均值合成的（arcgis手动点击好麻烦）matlab代码，只找到了一篇用matlab进行最大值合成NDVI的csdn文章

1.缘由

        帮一个同学处理数据，他有8760个小时的era5数据，类型是tif，刚好是一年的数据，将8760个tif文件求均值，得到一年的均值tif数据。翻阅csdn后，发现好像没我想要的代码，于是迷迷糊糊写了以下的代码

2.matlab代码（求均值，可行）

%2023.5.7晚间与实验室，一开始没有弄清楚矩阵，现在想想大概是明白了
% 心得
% 1.matlab处理tif文件，必须得有tif的sample例子，所以[a,R]，info必不可少，否则丢失地理信息
% 2.关于矩阵，单波段的tif文件，是m*n的二维矩阵，如果要求tif均值，实际上是吧所有二维矩阵，对应的位置相加
% 然后再让每个位置除以总数，则得到的均值二维矩阵
% 3.将得到的均值二维矩阵写入地理信息，就可以输出为带有地理信息的tif了

% 1.设置指定文件夹，打印出来你要均值处理的tif文件的个数
folder_path = 'D:\000-wangsiqi\tcwv_avg\2019'; % 替换为实际的文件夹路径
tif_files = dir(fullfile(folder_path, '*.tif')); % 获取所有tif文件
num_tif_files = length(tif_files); % 获取tif文件个数
fprintf('There are %d tif files in the folder %s.\n', num_tif_files, folder_path); % 打印结果

% 2.导入sample文件的地理信息
[a,R]=geotiffread('D:\000-wangsiqi\tcwv_avg\2019\era5_1.tif');%先投影信息
info=geotiffinfo('D:\000-wangsiqi\tcwv_avg\2019\era5_1.tif');

% 3.创建一个二维零矩阵，行列数=sample的行列数
datasum=zeros(size(a,1),size(a,2)); %size(a,1)和size(a,2)分别是进行合成的图像的行列号
% 4.循环读取tif，并且求sum
for i=1:num_tif_files
    data=importdata(strcat('D:\000-wangsiqi\tcwv_avg\2019\era5_',int2str(i),'.tif'));
    %data=reshape(data,size(a,1)*size(a,2),1);
    datasum=datasum+data
    disp(i)%用于提醒进行到第几次循环了
end
% 5.求均值
datamean=datasum/num_tif_files
% 6.设置输出路径
filenameet=strcat('D:\000-wangsiqi\tcwv_avg\2019_shuchu\2019_average_matlab_v02.tif');
% 7.写入文件夹
geotiffwrite(filenameet,datamean,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
disp('all done')

3.arcpy尝试

本人尝试了两串arcpy代码，但是都有一定问题，代码如下。

arcpy代码1（求均值，最多支持1000个tif文件以内的均值计算）

# 这个代码写于2023.5.7，是用arcpy批量mean合成每小时的era5数据，获取每年的平均值
# 笑死，文件多了的话arcpy处理不了，arcpy最多一次处理1000个文件的平均值
import arcpy
import os
from  arcpy import env
arcpy.CheckOutExtension("spatial")
arcpy.env.workspace=r"D:\000-wangsiqi\tcwv_avg\2019"# 工作空间
tif_list = arcpy.ListRasters("*","tif")# 读取文件夹下的所有提防文件
#东西太多了，别一个一个打印，打印最后一项看看对不对
print(tif_list[-1])
out = r"D:\000-wangsiqi\tcwv_avg\2019_shuchu\2019_arcpy_average.tif"#输出路径和名称
arcpy.gp.CellStatistics_sa(tif_list, out, "MEAN", "DATA")#取平均值，想求最大值可以将MEAN换位MAXIMUM
print("All done,please close!")

报错如图

arcpy代码2（求均值，我在arcgis中运行以下代码，直接死机了哈哈哈，应该是8760个文件太多了，运行不过来，少点应该可以）

import arcpy
from arcpy.sa import *
from arcpy import env


# 设置工作空间
arcpy.env.workspace = r"D:\000-wangsiqi\tcwv_avg\2019"

# 获取所有tif文件的文件名
tif_list = arcpy.ListRasters("*", "TIF")
count=1
# 如果存在tif文件，则进行平均值合成
if tif_list:
    # 创建一个空的Raster对象，用于存储平均值合成后的结果
    output_raster = arcpy.Raster(tif_list[0])
    output_raster = output_raster * 0

    # 遍历所有tif文件，计算平均值并合成
    for tif_file in tif_list:
        raster = arcpy.Raster(tif_file)
        output_raster = output_raster + raster
        print(count)
        count=count+1

    output_raster = output_raster / len(tif_list)

    # 保存结果
    output_raster.save(r"D:\000-wangsiqi\tcwv_avg\2019_shuchu\2019_arcpy_average.tif")
else:
    print("No TIFF files found in the workspace.")