数据量大的分块处理的代码逻辑[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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值