数据量大的分块处理的代码逻辑[20250310]-优快云博客

#日数据量很大，尤其是处理栅格数据的时候，日数据的计算会很慢，于是我开始琢磨分块的运算。其实分块的逻辑很好想，但是难者不会，会者不难。脑子里要绕好多此才能绕出来；

首先剖一下我想要解决的问题：

如果你的生长季每一个栅格都不一样，那么怎么可以统计生长季内的气候特征？逻辑就是，定位一个栅格的位置，然后在栅格提取生长的开始日和结束日。比如开始日是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