DEM生产坡度图、坡向图、山体阴影图、地形图、等高线图原理以及MATLAB实现
2024-05-03 18:49:03
文章目录
- 前言
- 一、通过DEM图生成坡度图
- （1）生成原理以及公式
- （2）代码段
- （3）结果
- 二、生成坡向图
- （1）生成原理以及公式
- （2）代码段
- （3）结果
- 三、生成山体阴影图
- （1）生成原理以及公式
- （2）代码段
- （3）结果
- 四、通过DEM数据生成三维地形图、伪彩色图以及等高线图
- 代码段
- 结果
前言
一、通过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-2for 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 ~= 0Aspect_rad = atan2(DZDY,-DZDX);if Aspect_rad<0Aspect_rad = 2*pi + Aspect_rad;endendif DZDX==0if DZDY>0Aspect_rad = pi/2;elseif DZDY<0Aspect_rad = 2 * pi - pi / 2;endelseAspect_rad=Aspect_rad;endangle=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:(请参考上面的坡度代码段，里面囊括以下内容)
if DZDX ~= 0Aspect_rad = atan2(DZDY,-DZDX);if Aspect_rad<0Aspect_rad = 2*pi + Aspect_rad;endendif DZDX==0if DZDY>0Aspect_rad = pi/2;elseif DZDY<0Aspect_rad = 2 * pi - pi / 2;endelseAspect_rad=Aspect_rad;endangle=round(Aspect_rad*(180/pi));result2(i+1,j+1)=angle;
- 段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:rows2for jj=1:cols2result3(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)<0result4(ii,jj)=0;elseresult4(ii,jj)=1;endif result3(ii,jj)<0result5(ii,jj)=0;elseresult5(ii,jj)=result3(ii,jj);endend
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:rows3X(iii,:)=mx(iii);
end
for jjj=1:cols3Y(:,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('等高线图');
结果
|
|
|
DEM生产坡度图、坡向图、山体阴影图、地形图、等高线图原理以及MATLAB实现相关推荐
- GEE学习笔记：在Google Earth Engine(GEE)中计算坡度、坡向、山体阴影
本次实验使用的 SRTM 数字高程数据是 30 米分辨率数据,对某地区的坡度.坡向和山体阴影信息进行提取. 目录 1.获取SRTMGL1_003 数据 2.计算地形特征 3.分别提取各地形因子 4.完 ...
- 【GDAL】python读取DEM计算坡度与坡向
利用GDAL读入DEM与Landsat影像,由于DEM是WG84坐标系,Landsat是WGS84坐标系UTM投影,因此处理在实际应用中需要将DEM进行投影转换. 大概分为以下几个步骤: 读取DEM, ...
- DEM中坡度和坡向的计算
1.坡度的计算 地表单元的坡度就是其切平面的法线方向与Z轴的夹角.若需求格网点上的坡度时,可取3×3的格网单元进行计算.也可求出该格网点八个方向上的坡度,再取其平均值. 2.坡向的计算 坡向是地表单元 ...
- python计算坡度_基于python实现利用DEM数据计算坡度、坡向
1.Python的地形三维可视化--简介Matplotlib和gdal https://blog.csdn.net/allenlu2008/article/details/51880333 2.Pyc ...
- 基于python实现利用DEM数据计算坡度、坡向
基本概念 DEM数据 DataMark:CNSDTF-DEM Version:1.0 Unit:M Alpha:0.000000 Compress:0.000000 X0:258000.000 Y0: ...
- GIS操作之高程、坡度、坡向
高程数据的处理: 获取研究区域高程数据(地理空间数据云) 将获取数据解压并且导入ArcGis当中,开始处理高程数据 打开ArcToolbox当中的数据管理工具中下的栅格,栅格当中的栅格数据集下的镶嵌或 ...
- CentOS 7.2下ELK分析Nginx日志生产实战(高清多图)
注:本文系原创投稿 本文以api.mingongge.com.cn域名为测试对象进行统计,日志为crm.mingongge.com.cn和risk.mingongge.com.cn请求之和(此二者域名 ...
- 【GlobalMapper精品教程】036：基于DEM的流域计算生成流域图
Globalmapper基于DEM的流域计算生成流域图教程. 文章目录 一.加载DEM 二.流域分析 一.加载DEM 加载配套实验数据. 二.流域分析 GM中的流域分析工具位于分析→生成流域,如下所示 ...
- Tableau绘制漏斗图、甘特图、瀑布图、镶边面积图、阴影坡度图
Tableau绘制漏斗图.甘特图.瀑布图.镶边面积图.阴影坡度图 本文首发于博客冰山一树Sankey,去博客浏览效果更好.直接右上角搜索该标题即可 一. 漏斗图 数据源 1.1 分色直条漏斗图 (1) ...
最新文章
- Highcharts Pie 饼图提示标签IE下重叠解决方法，及json数据绑定方法
- extjs window显示在顶层
- 无人驾驶图像数据集_自动驾驶数据集
- Java学习笔记十五
- 回归_英国酒精和香烟关系
- linux rpm 校验软件包中的文件
- linux课堂笔记（7）
- MYSQL的索引类型：PRIMARY, INDEX,UNIQUE,FULLTEXT,SPAIAL 有什么区别？各适用于什么场合？
- 如何巧妙的申请换部门_如何设置户外广告？市城管局局长体验户外广告审批流程...
- linux c 获取网卡状态（UP or DOWN）
- 开课吧python小课学了有用吗-好消息！今天，审计、会计、税务、财务主管彻底沸腾了……...
- InnoDB缓存相关优化
- 这家初创公司用端到端安全保护物联网设备
- 用泰勒展开式计算sin（x）的值
- 零代码与低代码快速开发平台的区别
- 学生专用计算机游戏怎么按,学生计算器怎么玩
- 测试人员的基本技能要求 - 快速掌握业务知识的能力
- 关闭iOS上京东app不停询问“京东想从MF839粘贴”，您允许这样做吗？ - 允许iphone应用访问剪切板
- 基于 NCNN 的 Chinese-Lite 模型测试
- U盘EFI分区删不掉怎么办