DEM生产坡度图、坡向图、山体阴影图、地形图、等高线图原理以及MATLAB实现

前言

一、通过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坡度融合图');

（3）结果

二、生成坡向图

（1）生成原理以及公式

• 坡向就是地表面一点的切平面的法线矢量n在水平面的投影xoy与过该点正北方向X的夹角Aspect （如图一所示），用于识别出从每个像元到其相邻像元方向上值的变化率最大的下坡方向。坡向可以被视为坡度方向。输出栅格中各像元的值可指示出各像元位置处表面的朝向的罗盘方向 （如图三所示）。将按照顺时针方向进行测量，角度范围介于 0（正北）到 360（仍是正北）之间，即完整的圆 （如图四所示）。不具有下坡方向的平坦区域将赋值为 -1 （如图五所示）。
• 坡向的计算类似于坡度的计算方法，公式如下

• 但是为了其罗盘方向能准确，所以坡向的计算牵扯到循环 （如图二所示），在编程时一定要注意循环。

图一

图二

图三

图四

图五

坡向图一般有以下几种作用：

• 在寻找最适合滑雪的山坡的过程中，查找某座山所有朝北的坡。
• 在统计各地区生物多样性的研究中，计算某区域中各个位置的日照强度。
• 作为判断最先遭受洪流袭击的居住区位置研究的一部分，在某山区中查找所有朝南的山坡，从而判断出雪最先融化的位置。
• 识别出地势平坦的区域，以便从中挑选出可供飞机紧急着陆的一块区域。

（2）代码段

1. 段1:(请参考上面的坡度代码段，里面囊括以下内容)

 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;

2. 段2：

angle_map=result2(2:rows1-1,2:cols1-1);
figure;
imagesc(angle_map);
axis off;
title('坡向图');
%------------------------------------------------------------------------
%进行颜色重分类
%------------------------------------------------------------------------
c_map=zeros(360,3);
c_map(1:22,:)=repmat([0,1,0],22,1);
c_map(338:360,:)=repmat([0,1,0],23,1);

c_map(23:67,:)=repmat([0,1,1],45,1);
c_map(68:112,:)=repmat([1,1,0],45,1);
c_map(113:157,:)=repmat([1,0,0],45,1);
c_map(158:202,:)=repmat([1,1,1],45,1);
c_map(203:247,:)=repmat([0,0,0],45,1);
c_map(248:292,:)=repmat([0,0,1],45,1);
c_map(293:337,:)=repmat([1,0,1],45,1);
colormap(c_map);
c=colorbar;
set(c,'tickdir','out') 
set(c,'YTick',0:45:360); 
set(c,'YTickLabel',{'北','东北','东','东南','南','西南','西','西北','北'})
%------------------------------------------------------------------------
%------------------------------------------------------------------------
angle_mean=fix(mean(angle_map(:)));
ol2=DEM_mean/angle_mean;
die_map2=(imadd(DEM_COP,int16(angle_map.*ol)))/2;
figure;
imagesc(die_map2);
axis off;
colormap(gray(256));
title('DEM坡向融合图');

（3）结果

三、生成山体阴影图

（1）生成原理以及公式

山体阴影图就是已知高度角和方向角的光照源的照射下，去假定来获取表面的假定照明度。通过设置假定光源的位置和计算与相邻像元相关的每个像元的照明度值，即可得出假定照明度。由于坡度和坡向已经解决，山体阴影还需要以下两种定义：

太阳方位角： 是太阳在方位上的角度，它通常被定义为从北方沿着地平线顺时针量度的角。可以理解为太阳与中心地物的连线在水平面的投影线与正北方向的夹角，在本文中默认为315°；（如图一所示）

太阳高度角： 是指某地太阳光线与通过该地与地心相连的地表切面的夹角。当太阳高度角为90°时，太阳辐射强度最大；太阳斜射地面程度越大（即太阳高度角越小），太阳辐射强度就越小。本文默认60°。（如图二所示）

图一

图二

山体阴影计算公式：

注：式中Hillshade表示山体阴影（-255<Hillshade<255）

通过计算出山体阴影我们可以生成以下两种图：

1）山体阴影渲染图（山体阴影渲染可更好的表示地形）

2）山体阴影二值图（山体阴影二值图可应用于通视分析）

（2）代码段

%------------------------------------------------------------------------
%------------------------------------------------------------------------
%计算入射角度
Altitude=30;%默认高度角
Zenith_deg=90-Altitude;%将高度角转换为天顶角
Zenith_rad=Zenith_deg*(pi/180);%将角度转换为弧度
%计算入射方向
Azimuth=315;%默认方位角
Azimuth_math=360.0-Azimuth+90;%将方位角转换为天顶角
Azimuth_rad=Azimuth_math*(pi/180);%将角度转换为弧度
[rows2,cols2]=size(DEM_COP);
[result3,result4,result5]=deal(zeros(rows2,cols2));
%将坡向坡度数据转换为弧度数据
solpe_map_rad=double(solpe_map).*(pi/180);
angle_map_rad=double(angle_map).*(pi/180);
for ii=1:rows2
    for jj=1:cols2
        result3(ii,jj)=255*((cos(Zenith_rad)*cos(solpe_map_rad(ii,jj)))...
            +(sin(Zenith_rad)*sin(solpe_map_rad(ii,jj))*...
            cos(Azimuth_rad-angle_map_rad(ii,jj))));
        if result3(ii,jj)<0
            result4(ii,jj)=0;
        else
            result4(ii,jj)=1;
        end
        if result3(ii,jj)<0
            result5(ii,jj)=0;
        else
            result5(ii,jj)=result3(ii,jj);
        end
    end
end
figure;
imagesc(result5);
axis off;
colormap(gray(256));
colorbar;
title('315°太阳方位角下的山体阴影渲染图');
figure;
imshow(result4);
title('315°太阳方位角下的山体阴影二值图');

（3）结果

四、通过DEM数据生成三维地形图、伪彩色图以及等高线图

我们知道，DEM格网数据是由高程值集合组成的规则格网数据集，而每个格网对应一个行列号，因此我们
需要利用循环构建行矩阵A和列矩阵B，在已知DEM矩阵的情况下，利用surf(X,Y,Z)函数直接绘制，也可以采用[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(y)),linspace(min(y),max(y)),‘v4’)函数生成格网数据，再使用surf函数。

代码段

%------------------------------------------------------------------------
%------------------------------------------------------------------------
DEM_COP=data(8:rows-8,6:cols-6);
[X,Y]=deal(DEM_COP);
[rows3,cols3]=size(DEM_COP);
my=(linspace(1,cols3,cols3))';
mx=(linspace(1,rows3,rows3))';
for iii=1:rows3
    X(iii,:)=mx(iii);
end
for jjj=1:cols3
    Y(:,jjj)=my(jjj);
end
figure;
surf(X,Y,DEM_COP);
% el=85;  %设置仰角为75度。
% for az=0:1:360  %让方位角从0变到360，绕z轴一周
%     view(az,el);
%     drawnow;
% end
% title('三维地形图');

figure;
A=pcolor(X,Y,DEM_COP);
shading interp;
view(90,90);
axis off;
colormap hot;
title('伪彩色图');
figure;
B=contour(X,Y,DEM_COP);
view(90,90);
axis off;
title('等高线图');

结果