Matlab 将三维矩阵分块计算然后组合计算结果

---

前言

如果矩阵过大，有三维时间序列，则需要将矩阵分块成为无数子个数然后计算，这里计算的是一个文件夹下格式为tif和nc的大矩阵的第三维的分位数。注意是一个文件夹下有若干tif文件或者nc文件。比如1990-2020年每一天有一副影像。这里是针对计算机一次性无法读取那么多的文件，内存不够的情况，但是分块读取计算是可以的，请注意，计算在栅格尺度上，相互之间不受影响。如果计算机存储性能足够，请忽略。

代码

1.直接上代码

%如果矩阵过大，有三维时间序列，则需要将矩阵分块成为无数子个数然后计算
%这里计算的是一个文件夹下大矩阵的第三维的分位数
clc
clear all
tic  % 开始计时
folder = 'H:\MSWEP_V280\PastDaily\';
tifFiles = dir(fullfile(folder, '*.nc'));

firstImage = ncread(fullfile(folder, tifFiles(1).name),'precipitation');
[Rows, Cols] = size(firstImage);

% 设置块的大小
blockSizeRows = 200; % 选择适当的行大小，最好可以除尽，如果不能整除，采取加0的行或列，最终再还回去。
blockSizeCols = 200; % 选择适当的列大小

%prc5=nan(Rows, Cols);
prc10=nan(Rows, Cols);
%prc15=nan(Rows, Cols);
%prc20=nan(Rows, Cols);
prc90=nan(Rows, Cols);

for i=1:Rows/blockSizeRows
   for j=1:Cols/blockSizeCols
        stack = nan(blockSizeRows, blockSizeCols, numel(tifFiles), 'double');
        startRow = (i - 1) * blockSizeRows + 1;
        endRow = min(startRow + blockSizeRows - 1, Rows);
        startCol = (j - 1) * blockSizeCols + 1;
        endCol =  min(startCol + blockSizeCols - 1, Cols);

        for k = 1:numel(tifFiles)
    % 读取块
        %[tifdata,R]=geotiffread(fullfile(folder, tifFiles(k).name));
        %dataBlock = imread(fullfile(folder, tifFiles(k).name), 'PixelRegion', {
   
   [startRow, endRow], [startCol, endCol]});
        tifdata=ncread(fullfile(folder, tifFiles(k).name),'precipitation');
        dataBlock = tifdata(startRow:endRow, startCol:endCol);
    % 填充 stack
        stack(:, :, k) = dataBlock;
        end
        stack(stack == 0)=nan; % 选择大于0的值

        %prc5(startRow:endRow, startCol:endCol) = prctile(stack, 5,3);
        prc10(startRow:endRow, startCol:endCol) = prctile(stack, 10,3);
        %prc15(startRow:endRow, startCol:endCol) = prctile(stack, 15,3);
        %prc20(startRow:endRow, startCol:endCol) = prctile(stack, 20,3);
        prc90(startRow:endRow, startCol:endCol) = prctile(stack, 90,3);
        toc
   end
end

  disp('well done')  

% clc
% clear all
% tic  % 开始计时
% folder = 'E:\BaiduSyncdisk\precipitationCal\';
% tifFiles = dir(fullfile(folder, '*.tif'));
% 
% firstImage = imread(fullfile(folder, tifFiles(1).name));
% [Rows, Cols] = size(firstImage);
% 
% % 设置块的大小
% blockSizeRows = 900; % 选择适当的行大小
% blockSizeCols = 900; % 选择适当的列大小
% 
% gpmprc5=nan(Rows, Cols);
% gpmprc10=nan(Rows, Cols);
% gpmprc15=nan(Rows, Cols);
% gpmprc20=nan(Rows, Cols);
% gpmprc30=nan(Rows, Cols);
% 
% for i=1:Rows/blockSizeRows
%         startRow = (i - 1) * blockSizeRows + 1
%         endRow = min(startRow + blockSizeRows - 1, Rows);
%    for j=1:Cols/blockSizeCols
%         stack = nan(blockSizeRows, blockSizeCols, numel(tifFiles), 'double');
%         %startRow = (i - 1) * blockSizeRows + 1
%         %endRow = min(startRow + blockSizeRows - 1, Rows);
%         startCol = (j - 1) * blockSizeCols + 1
%         endCol =  min(startCol + blockSizeCols - 1, Cols);
% 
%         for k = 1:numel(tifFiles)
%     % 读取块
%         %[tifdata,R]=geotiffread(fullfile(folder, tifFiles(k).name));
%        dataBlock = imread(fullfile(folder, tifFiles(k).name), 'PixelRegion', {
   
   [startRow, endRow], [startCol, endCol]});
%         %dataBlock = tifdata(startRow:endRow, startCol:endCol);
%     % 填充 stack
%         stack(:, :, k) = dataBlock;
%         end
%         stack(stack == 0)=nan; % 选择大于0的值
% 
%         gpmprc5(startRow:endRow