#日数据量很大,尤其是处理栅格数据的时候,日数据的计算会很慢,于是我开始琢磨分块的运算。其实分块的逻辑很好想,但是难者不会,会者不难。脑子里要绕好多此才能绕出来;
首先剖一下我想要解决的问题:
如果你的生长季每一个栅格都不一样,那么怎么可以统计生长季内的气候特征?逻辑就是,定位一个栅格的位置,然后在栅格提取生长的开始日和结束日。比如开始日是X1a月X2a日,结束日是X2a月X2b日,每一个栅格都有这样两个日期,那么一年365或366天,这两个日期换成日序DOY;那么一年365天的数据被你储存在一个mat 中,我需要处理20年的数据,那么每一个栅格都要调用年数据,那么大,多次调用且要运算快且不卡,那么就需要分块处理。
我的运行主代码如下:
%提取所有频率类别的点坐标及对应频率值
all_points = [];
%生成每个点的idx
n_points = size(all_points, 1);
%% 加入累积降水量的数据
for year_idx = 1:n_years
year = years(year_idx);
last_day = all_last_day{year_idx};%加载每个栅格点每年的生长季结束日期
maxday=dynamic_data1(:,year_idx);%选择所有栅格最大值的点
% 处理单一年份
pre_data_final(:,year_idx) = process_single_year_mean(...
year, last_day, maxday, precip_paths, all_points);
fprintf('完成年份: %d\n', year);
save('somewhere you like.mat', 'pre_data_final', '-v7.3'); % 定期保存
end
我处理一年的数据的代码如下,这部分代码的逻辑主要是:分块提取降水数据的运行代码,
function SM_mean = process_single_year_mean(year, last_day, maxday, precip_paths, all_points)
n_points = size(all_points, 1);
% 确定降水文件索引
path_idx = year - 2000; % 例如year=2001对应precip_paths{2}
% 使用 imfinfo 获取波段数(直接读取文件元数据)
info_current = imfinfo(precip_paths{path_idx+1});
total_bands_current = info_current.SamplesPerPixel;%不同的数据不一样
% 计算实际需要的波段范围
max_fd = max(maxday); % 这年所有栅格的最大值
% 分块读取降水数据(每次最多读10天)
band_block_size = 50; % 可根据内存调整
% 处理当前年降水(返回累积值和有效天数)
[pre_f, valid_days_f] = process_precipitation_f(... % 修改返回值
precip_paths{path_idx+1},...
1,...
max_fd,...
maxday,...
all_points,...
band_block_size,...
n_points...
);
% 合并结果(使用真实有效天数)
total_valid_days = valid_days_f; % 关键修改点
SM_mean = pre_f ./ total_valid_days; % 用真实有效天数做分母
% 处理除零情况
SM_mean(total_valid_days == 0) = NaN;
end
第三步,则是处理单个栅格的加和或平均,因为累积加和的比较简单,所以我们来看求平均的情况:求平均的情况是这样的~
第三部分代码有问题,先保存草稿,去调整了!!!
% 定义通用处理函数
function [sum_precip, valid_days] = process_precipitation_f(precip_path, band_start_init, band_end_limit, day_array, all_points, band_block_size, n_points)
sum_precip = zeros(n_points, 1, 'single');
valid_days = zeros(n_points, 1, 'single'); % 统计有效天数
for band_start = band_start_init:band_block_size:band_end_limit
band_end = min(band_start + band_block_size - 1, band_end_limit);% band_end为最大值,因为分块,所以最大值分大小
bands = band_start:band_end; %开始时间到分块的结束部分
PRE_block = single(readgeoraster(precip_path, 'Bands', bands));
for i = 1:n_points
day_i = day_array(i);
if isnan (day_i)
continue; % 统一跳过无效值
end
start_band_in_block = 1;
end_band_in_block = band_end - band_start + 1;
% 提取数据段
data_segment = PRE_block(all_points(i,1), all_points(i,2), start_band_in_block:end_band_in_block);
% 统计非NaN值
valid_mask = ~isnan(data_segment);
sum_precip(i) = sum_precip(i) + sum(data_segment(valid_mask), 'omitnan');
valid_days(i) = valid_days(i) + sum(data_segment(:)); % 新增:累加有效天数
end
clear PRE_block; % 释放内存
end
end