MATLAB DEM栅格数据平滑处理

本博客适用于数字高程模型的平滑处理，对于其他数据的平滑处理仅作参考
平滑处理公式采用Daly1984
Daly C1994_A statistical topographic model for mapping climatological precipitation over mountainous terrain.
ele(m,n) = 0.5ele(m,n)+0.125(ele(m-1,n)+ele(m+1,n)+ele(m,n-1)+ele(m,n+1))
即中心点的高程由自身和上下左右四个点的高程加权平均得到
（这里就涉及边界位置处的点如何处理的问题，我的做法是给栅格额外再加一圈0值，根据0的个数改变权重，实现边界位置点的高程平滑）

MATLAB代码如下：

%% 此代码是用于平滑处理DEM数据
%计算公式（Daly,1984） 
%ele(m,n) = 0.5ele(m,n)+0.125(ele(m-1,n)+ele(m+1,n)+ele(m,n-1)+ele(m,n+1))
%为了方便处理边界位置处的点，给栅格矩阵额外加一圈
clc
clear

%% 读取数据
[DEM,GeoRef] = geotiffread('Resample.tif');   %GeoRef存储tif文件的地理空间信息
[x,y] = size(DEM);                            %存储tif文件的栅格大小
lat = GeoRef.LatitudeLimits;                  %存储tif文件的维度范围
lon = GeoRef.LongitudeLimits;                 %存储tif文件的经度范围
DEM(DEM==-32768) = 0;                         %将背景值换成0，之后统一处理

%% 给DEM矩阵再加外围一圈，方便之后边界位置的处理
DEM = [zeros(x,1),DEM];    %最左侧加一列
DEM = [DEM,zeros(x,1)];    %最右侧加一列
DEM = [zeros(1,y+2);DEM];  %最上侧加一行
DEM = [DEM;zeros(1,y+2)];  %最下侧加一行
DEM_origine = DEM;         %存储此时的矩阵，方便检查，看有没有出错

%% 平滑处理
for t = 1:40               %平滑处理的重复次数，最多40次
    for i = 2:x+1
        for j = 2:y+1
            if DEM(i,j) == 0                 %0意味着背景值或者最外围一圈，跳过
                continue
            else           
                a = sum([DEM(i-1,j),DEM(i+1,j),DEM(i,j-1),DEM(i,j+1)]==0);   %计算目标栅格周围为0的栅格数
                %计算中心点的高程，4/(4-a)会根据周围0的个数而变化，从而改变权重，默认不出现周围全是0的情况，a在0-3之间
                DEM(i,j) = 0.5*DEM(i,j)+0.125*4/(4-a)*(DEM(i-1,j)+DEM(i+1,j)+DEM(i,j-1)+DEM(i,j+1));
            end
        end
    end
    if t == 8||t == 16||t == 24||t == 32||t == 40     %重复平滑处理8、16等次数后输出图像
        DEM_out = DEM;
        DEM_out(:,16) = [];
        DEM_out(:,1) = [];
        DEM_out(15,:) = [];
        DEM_out(1,:) = [];                   %删掉外围一圈
        DEM_out(DEM_out==0) = -32768;        %将0值变回之前的默认背景值
        R = georasterref('RasterSize', size(DEM_out),'Latlim', [double(min(lat))...
            double(max(lat))], 'Lonlim', [double(min(lon)) double(max(lon))]); %记录空间坐标信息
        R.ColumnsStartFrom = 'north';
        geotiffwrite(['F:\Data of China Flash Flood\DEM_smoothing_test\Resample_tif\smoothed',num2str(t),'times','.tif'],DEM_out,R);     %保存绘制的地理栅格图像
    end
end
disp('finish!')

以上代码的参考文献有：
MATLAB nc文件转tif （可视化范例）
Matlab：Save as GeoTiff Format

效果如下：

后记：

写博客的初衷是分享经验，同时是算是自己对思路和代码的整理，方便日后处理数据，应该可以帮到很多人。
我已免费分享我的心得，如果看官还有其他问题的，那么：知识付费，我的时间和经验正好可以解决你的问题。
咨询问题请添加QQ：819369354
2022年4月20日