通过DEM图生成坡度图

（1）生成原理以及公式

所谓坡度，即过地面某一点的切平面与水平面的夹角，该夹角就是该点的坡度。而坡度一般有两种表示方法（度数或坡度百分比），本文以度数为例。因此我们只需要知道两点的高程增量以及水平增量，便可以算出这两点所在平面的单一坡度值。如果将高程增量百分比视为高程增量除以水平增量后再乘以 100，就可以更好地理解高程增量百分比。请考虑下面的三角形 B。当角度为 45 度时，高程增量等于水平增量，所以高程增量百分比为 100%。如三角形 C 所示，当坡度角接近直角（90 度）时，高程增量百分比开始接近无穷大。

图1

图2

DEM作为地理常用数据集的4D产品之一，其具有非常广泛的作用，它的数据组织格式即为连续地形高程值集合，在MATLAB中显示为由高程值组成的矩阵N。要利用该DEM生成与DEM同样大小的坡度图，需要满足以下几步：第一步：将DEM数据矩阵外围扩充一圈矩阵（如右下图灰色边缘所示），可以默认其值为零，本文采用其值为所有高程的平均值；第二步：构建一个3x3的窗口，该窗口所对应的中间像元对应DEM每一个位置的像元，除中间像元，其他八个像元称为领域；

计算原理将一个平面与要处理的像元或中心像元周围一个 3 x 3 的像元邻域的 z 值进行拟合。该平面的坡度值通过最大平均值法来计算。该平面的朝向就是待处理像元的坡向。坡度值越小，地势越平坦；坡度值越大，地势越陡峭。通过移动移动窗口来遍历每一个像元，将计算出来的坡度值赋值给每一个像元。计算窗口大小内在x方向上的变化率

在这里插入图片描述

计算窗口大小内在y方向上的变化率

在这里插入图片描述

计算坡度的公式（其中的 57.29578 是对 180/pi 的计算结果进行截断而得到的值）：

在这里插入图片描述

（2）代码段

clc;
clear all;
[data,info]=geotiffread("DEM1.tif");
figure;
[rows,cols]=size(data);
%下面这一步是为了提出原始DEM数据中的异常值，其和裁剪时的边缘有关
%------------------------------------------------------------------------
DEM_COP=data(8:rows-8,6:cols-6);
%------------------------------------------------------------------------
figure;
imagesc(DEM_COP);
axis off;
colorbar;
title('DEM图');

%获取整个DEM的像元平均值，作为边缘扩充的值
DEM_mean=fix(mean(DEM_COP(:)));
MAP=padarray(DEM_COP,[1 1],DEM_mean);

double_map=double(MAP);
[rows1,cols1]=size(double_map);
[result,result2]=deal(double_map);
for i=1:rows1-2
    for j=1:cols1-2
        %------------------------------------------------------------------------
        MR=double_map(i:i+2,j:j+2);
        DZDX=((MR(1,3)+2*MR(2,3)+MR(3,3))-(MR(1,1)+2*MR(2,1)+MR(3,1))) / (8 * 5);
        DZDY=((MR(3,1)+2*MR(3,2)+MR(3,3))-(MR(1,1)+2*MR(1,2)+MR(1,3))) / (8 * 5);
        result(i+1,j+1)=round(atan(sqrt(DZDX^2+DZDY^2))* 57.29578);
        %------------------------------------------------------------------------
        if DZDX ~= 0
            Aspect_rad = atan2(DZDY,-DZDX);
            if Aspect_rad<0
                Aspect_rad = 2*pi + Aspect_rad;
            end
        end
        if DZDX==0
            if DZDY>0
                Aspect_rad = pi/2;
            elseif DZDY<0
                Aspect_rad = 2 * pi - pi / 2;
            end
        else
            Aspect_rad=Aspect_rad;
        end
        angle=round(Aspect_rad*(180/pi));
        result2(i+1,j+1)=angle;
        %------------------------------------------------------------------------
    end
end
%------------------------------------------------------------------------
solpe_map=int16(result(2:rows1-1,2:cols1-1));
figure;
imagesc(solpe_map);
axis off;
colorbar ;
c=colorbar;
set(c,'tickdir','out') 
set(c,'YTick',0:10:100);
set(c,'YTickLabel',{'0\circ','10\circ','20\circ','30\circ',...
    '40\circ','50\circ','60\circ','70\circ','80\circ','90\circ','100\circ'});
title('坡度图');
%进行迭代融合
solpe_mean=fix(mean(solpe_map(:)));
ol=DEM_mean/solpe_mean;
die_map=imadd(DEM_COP,solpe_map.*ol);
figure;
imagesc(die_map);
axis off;
colormap(gray(256));
title('DEM坡度融合图');

DEM生产坡度图（转载）

通过DEM图生成坡度图

（1）生成原理以及公式

（2）代码段

（3）结果