当多个行政单元合并成一个影像时,在行政边界附近可能会产生一些无效值,如下图所示。
示例影像为分类结果图,值在1~5以内,无效值为127。
本文使用最邻近法给无效值赋值(3 * 3区域内类别数最多的值,如果最多的值有多个,则以第一个为准)。
影像的行政单元也不是规则的,即周边也存在大量的无效值。此时要是逐行遍历的话,运算量非常大。
此时的主要思路为:
- 原始影像A,将无效值赋为0后保存为B;
- 由于行政边界所造成的无效值多呈线状分布(可见上图),因此将B进行均值滤波,窗口可以大一些(程序中设置为25);
- 均值滤波后保存为C,C中为0点则为边界的无效值,识别后不参与运算;
- B == 0 && C != 0 的点则为目标点,进行最邻近插值(需要多次最邻近插值,需要外设一个循环)。
- 最邻近插值后,由于位于行政边界外的部分点也被识别为目标点了,因此需要使用行政边界(shp数据)进行裁剪。
第1~4步的参考代码如下。
%%
clc;
clear;
%%
% 获取.m所在的文件夹
fullpath = mfilename('fullpath');
[path, name] = fileparts(fullpath);
%%
filepath = strcat(path, '\oridata.tif');
[A, RA] = readgeoraster(filepath);
% 有效值在1~5之间,将无效值(127)赋为0
A(A > 10) = 0;
info = geotiffinfo(filepath);
B = A;
R = 25;
I_3 = fspecial('average', [51, 51]); %均值滤波
C = imfilter(A, I_3);
B(C == 0) = 127;
C = B;
R = 1;
for iti = 1 : 100
[X, Y] = find(B == 0);
xsi = size(X, 1);
if xsi == 0
continue;
end
for i = 1 : xsi
[iti, i, xsi]
xmin = X(i) - R;
xmax = X(i) + R;
ymin = Y(i) - R;
ymax = Y(i) + R;
if xmin < 1
xmin = 1;
end
if xmax > max(X)
xmax = max(X);
end
if ymin < 1
ymin = 1;
end
if ymax > max(Y)
ymax = max(Y);
end
T = B(xmin : xmax, ymin : ymax);
F = T(T > 0 & T < 10);
if size(F, 1) == 0
continue;
end
C(X(i), Y(i)) = mode(F(:));
end
B = C;
end
outpath = strcat(path, '\fillnan.tif');
geotiffwrite(outpath, B, RA, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);