我有一个非常大的(2019x1678双)DEM(数字高程模型)文件放在MATLAB中作为矩阵。它的边包含NaN值。为了在我的代码中考虑边缘效应，我必须在我的DEM周围放置一个1单元缓冲区(与相邻单元相同的值)。在NaN存在的地方，我需要找到NaN值的边缘，以便构建缓冲区。我已经尝试了两种方法:
在第一个中，我获得所有非nan DEM值的行和列坐标，并找到每列的第一行和最后一行编号以获得南北边界，然后找到每行的第一行和最后一行编号以获得东西边界。我在sub2ind()中使用这些来创建我的缓冲区。
[r, c] = find(~isnan(Zb_ext)); %Zb is my DEM matrix
idx = accumarray(c, r, [], @(x) {[min(x) max(x)]});
idx = vertcat(idx{:});
NorthBoundary_row = transpose(idx(:,1)); % the values to fill my buffer with
NorthBoundary_row_ext = transpose(idx(:,1) - 1); % My buffer cells
columnmax = length(NorthBoundary_row);
column1 = min(c);
Boundary_Colu = linspace(column1,column1+columnmax-1,columnmax);
SouthBoundary_row = (transpose(idx(:,2))); % Repeat for south Boundary
SouthBoundary_row_ext = transpose(idx(:,2) + 1);
SouthB_Ind = sub2ind(size(Zb_ext),SouthBoundary_row,Boundary_Colu);
SouthB_Ind_ext = sub2ind(size(Zb_ext),SouthBoundary_row_ext, Boundary_Colu);
NorthB_Ind = sub2ind(size(Zb_ext),NorthBoundary_row, Boundary_Colu);
NorthB_Ind_ext = sub2ind(size(Zb_ext),NorthBoundary_row_ext, Boundary_Colu);
Zb_ext(NorthB_Ind_ext) = Zb_ext(NorthB_Ind);
Zb_ext(SouthB_Ind_ext) = Zb_ext(SouthB_Ind);
% Repeat above for East and West Boundary by reversing the roles of row and
% column
[r, c] = find(~isnan(Zb_ext));
idx = accumarray(r, c, [], @(x) {[min(x) max(x)]});
idx = vertcat(idx{:});
EastBoundary_colu = transpose(idx(:,1)); % Repeat for east Boundary
EastBoundary_colu_ext = transpose(idx(:,1) - 1);
row1 = min(r);
rowmax = length(EastBoundary_colu);
Boundary_row = linspace(row1,row1+rowmax-1,rowmax);
WestBoundary_colu = transpose(idx(:,2)); % Repeat for west Boundary
WestBoundary_colu_ext = transpose(idx(:,2) + 1);
EastB_Ind = sub2ind(size(Zb_ext),Boundary_row, EastBoundary_colu);
EastB_Ind_ext = sub2ind(size(Zb_ext),Boundary_row, EastBoundary_colu_ext);
WestB_Ind = sub2ind(size(Zb_ext),Boundary_row, WestBoundary_colu);
WestB_Ind_ext = sub2ind(size(Zb_ext),Boundary_row, WestBoundary_colu_ext);
Zb_ext(NorthB_Ind_ext) = Zb_ext(NorthB_Ind);
Zb_ext(SouthB_Ind_ext) = Zb_ext(SouthB_Ind);
Zb_ext(EastB_Ind_ext) = Zb_ext(EastB_Ind);
Zb_ext(WestB_Ind_ext) = Zb_ext(WestB_Ind);
这在我的小开发矩阵上工作得很好，但在我的全尺寸DEM上失败了。我不理解我的代码的行为，但看看数据，我的边界有空白。我想知道我是否需要更好地控制最大/最小行/列值的顺序，尽管在我对一个较小的数据集的测试中，所有似乎都是有序的....
第二种方法是我从一个类似的问题中得到的，基本上是使用扩张法。然而，当我转换到完整的数据集时，计算ZbDilated需要几个小时。虽然我的第一个方法不工作，它至少在几秒钟内计算。
[m, n] = size(Zb); %
Zb_ext = nan(size(Zb)+2);
Zb_ext(2:end-1, 2:end-1) = Zb; % pad Zb with zeroes on each side
ZbNANs = ~isnan(Zb_ext);
ZbDilated = zeros(m + 2, n + 2); % this will hold the dilated shape.
for i = 1:(m+2)
if i == 1 %handling boundary situations during dilation
i_f = i;
i_l = i+1;
elseif i == m+2
i_f = i-1;
i_l = i;
else
i_f = i-1;
i_l = i+1;
end
for j = 1:(n+2)
mask = zeros(size(ZbNANs));
if j == 1 %handling boundary situations again
j_f = j;
j_l = j+1;
elseif j == n+2
j_f = j-1;
j_l = j;
else
j_f = j-1;
j_l = j+1;
end
mask(i_f:i_l, j_f:j_l) = 1; % this places a 3x3 square of 1's around (i, j)
ZbDilated(i, j) = max(ZbNANs(logical(mask)));
end
end
Zb_ext(logical(ZbDilated)) = fillmissing(Zb_ext(logical(ZbDilated)),'nearest');
有没有人有任何想法，使这两个可用?
开头:
NaN NaN 2 5 39 55 44 8 NaN NaN
NaN NaN NaN 7 33 48 31 66 17 NaN
NaN NaN NaN 28 NaN 89 NaN NaN NaN NaN
这是在nan极限上缓冲的矩阵:
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
NaN NaN NaN 2 5 39 55 44 8 NaN NaN NaN
NaN NaN NaN NaN 7 33 48 31 66 17 NaN NaN
NaN NaN NaN NaN 28 NaN 89 NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
这是我使用fillmissing后想要得到的(尽管我注意到缓冲区值如何填充…):
NaN NaN 2 2 5 39 55 44 8 17 NaN NaN
NaN NaN 2 2 5 39 55 44 8 17 17 NaN
NaN NaN 2 2 7 33 48 31 66 17 17 NaN
NaN NaN NaN 2 28 33 89 31 66 17 17 NaN
NaN NaN NaN 5 28 55 89 8 NaN NaN NaN NaN
为了试图澄清我正在做的事情的任何混乱，以下是我用于填充缺失
的扩展所得到的逻辑0 0 1 1 1 1 1 1 1 1 0 0
0 0 1 1 1 1 1 1 1 1 1 0
0 0 1 1 1 1 1 1 1 1 1 0
0 0 0 1 1 1 1 1 1 1 1 0
0 0 0 1 1 1 1 1 0 0 0 0
应用3x3膨胀的更快方法如下:这确实涉及到一些大的中间矩阵，这使得它的效率低于使用imdilate。
[m, n] = size(Zb); %
Zb_ext = nan(size(Zb)+2);
Zb_ext(2:end-1, 2:end-1) = Zb; % pad A with zeroes on each side
ZbNANs = ~isnan(Zb_ext);
ZbDilated = ZbNANs; % this will hold the dilated shape.
% up and down neighbors
ZbDilated(2:end, :) = max(ZbDilated(2:end, :), ZbNANs(1:end-1, :));
ZbDilated(1:end-1, :) = max(ZbDilated(1:end-1, :), ZbNANs(2:end, :));
% left and right neighbors
ZbDilated(:, 2:end) = max(ZbDilated(:, 2:end), ZbNANs(:, 1:end-1));
ZbDilated(:, 1:end-1) = max(ZbDilated(:, 1:end-1), ZbNANs(:, 2:end));
% and 4 diagonal neighbors
ZbDilated(2:end, 2:end) = max(ZbDilated(2:end, 2:end), ZbNANs(1:end-1, 1:end-1));
ZbDilated(1:end-1, 2:end) = max(ZbDilated(1:end-1, 2:end), ZbNANs(2:end, 1:end-1));
ZbDilated(2:end, 1:end-1) = max(ZbDilated(2:end, 1:end-1), ZbNANs(1:end-1, 2:end));
ZbDilated(1:end-1, 1:end-1) = max(ZbDilated(1:end-1, 1:end-1), ZbNANs(2:end, 2:end));
这是一个乏味的方式来写它，我相信有一个循环可以写得更短，但我认为这使意图更清晰。
[编辑:因为我们在这里处理的是一个逻辑数组，我们也可以做A | B而不是max(A,B)。我不确定在时间上是否会有差别。
@beaker在评论中说不要使用
mask = zeros(size(ZbNANs));
mask(i_f:i_l, j_f:j_l) = 1; % this places a 3x3 square of 1's around (i, j)
ZbDilated(i, j) = max(ZbNANs(logical(mask)));
而不是
ZbDilated(i, j) = max(ZbNANs(i_f:i_l, j_f:j_l), [], 'all');
[编辑:因为我们在这里处理一个逻辑数组，而不是max(A,[],'all')，我们也可以做any(A,'all')，这应该更快。请看@beaker的其他评论