地形因子提取（MATLAB）

最新推荐文章于 2023-09-29 10:36:45 发布

xm00

最新推荐文章于 2023-09-29 10:36:45 发布

阅读量4k

点赞数 3

文章标签：数学建模 matlab 开发语言

本文链接：https://blog.youkuaiyun.com/weixin_56305254/article/details/131384311

版权

以某区域的DEM数据（.tif）为基础，通过MATLAB实现地形因子的提取，并进行可视化表达。

本次提取的地形因子主要包括：坡度、坡向、坡度变率、坡向变率、曲率、粗糙度、凹凸系数（面积、体积）。

一准备

（1）要提取区域的.tif文件。

（2）MATLAB 2022b（其他版本也可以）。

（3）新建脚本文件。

二文件读取

MATLAB在读取文件时，根目录为文件所在位置，所以.tif数据可以直接放在与,m文件相同的位置，可以避免路径写错的问题。

使用readgeoraster函数读取DEM数据，readgeoraster函数可以读取Z和R两个信息，Z表示地形高度数据；R表示空间参考对象。（注：要进行提取的tif数据要具有这两项）

代码如下：

% 读取 DEM 数据
[Z, R] = readgeoraster('你的DEM数据.tif');

三进行计算

（1）坡度

使用gradient函数计算DEM的坡度，再使用atan函数进行坡度的计算。

坡度的计算原理为：i = h / l x 100%

其中：h-高度差，l-水平距离，i-坡度。即两点的高程差与其水平距离的百分比。

% 计算 DEM 数据的梯度
[fx, fy] = gradient(Z, R.CellExtentInWorldX, R.CellExtentInWorldY); 
% 计算坡度
slope = atan(sqrt(fx.^2 + fy.^2));

（2）坡向

坡向以度为单位按顺时针方向进行测量，角度范围介于 0 度(正北)到 360 度(仍是正北，循环一周)之间。坡向格网中各像元的值均表示该像元的坡度所面对的方向。

注意：平坡没有方向，平坡的值被指定为 -1。

 % 计算 DEM 数据的梯度
[fx, fy] = gradient(Z, R.CellExtentInWorldX, R.CellExtentInWorldY);
aspect = atan2(-fx, -fy) * 180 / pi; % 计算坡向（以角度为单位）

（3）坡度变率

坡度变率是地面坡度在微分空间的变化率，是根据坡度的求算原理，在坡度的基础上再求一次坡度。代码如下：

% 计算 DEM 的坡度变率
[slopeX, slopeY] = gradient(slope, R.CellExtentInWorldX, R.CellExtentInWorldY);

（4）坡向变率

与坡度变率的计算原理相同，坡向变率就是地面坡向在微分空间的变化率，即在坡向的基础上再求一次坡向，代码如下：

% 计算 DEM 的坡向变率
[aspectX, aspectY] = gradient(aspect, R.CellExtentInWorldX, R.CellExtentInWorldY);

（5）曲率

地面曲率是描述地形起伏变化的一个重要指标，可以用于地形分析、水文模拟等。使用MATLAB进行地面曲率的计算，要对DEM进行两次梯度运算，代码如下：

% 计算 DEM 的地面曲率
[dx, dy] = gradient(Z, R.CellExtentInWorldX, R.CellExtentInWorldY);
[dxx, ~] = gradient(dx, R.CellExtentInWorldX, R.CellExtentInWorldY);
[~, dyy] = gradient(dy, R.CellExtentInWorldX, R.CellExtentInWorldY);
H = -sqrt(dxx.^2 + dyy.^2);

（6）粗糙度

地表粗糙度也是描述地表起伏变化的一个重要指标，可以用于土地利用评价、水文模拟等。计算过程为：首先计算DEM的梯度矩阵，再计算得到幅值矩阵（d），然后使用std2函数计算梯度矩阵的标准差（粗糙度指标），代码如下：

% 计算 DEM 的粗糙度
[dx, dy] = gradient(Z, R.CellExtentInWorldX, R.CellExtentInWorldY);
d = sqrt(dx.^2 + dy.^2);
r = std2(d);

（7）凹凸系数

地表曲率、地表粗糙度、凹凸系数都是描述地表起伏变化的重要指标，我们使用grandient函数计算DEM的梯度，得到梯度矩阵（H），然后代入公式计算DEM的凹凸系数，代码如下：

% 计算 DEM 的凸凹系数
[dx, dy] = gradient(Z, R.CellExtentInWorldX, R.CellExtentInWorldY);
[dxx, ~] = gradient(dx, R.CellExtentInWorldX, R.CellExtentInWorldY);
[~, dyy] = gradient(dy, R.CellExtentInWorldX, R.CellExtentInWorldY);
H = dxx + dyy;
C = abs(H) ./ sqrt(1 + dx.^2 + dy.^2).^3;

四结果可视化

使用figure函数对结果进行可视化，并导出可视化结果。

有多种表达形式，使用多子图（subplot），使用绘图函数（imagesc、plot），或者直接设置路径保存到电脑的指定位置，或者通过设置figure将读取的DEM数据和提取的地形因子分别显示。

其次，设置title、axis、legend、grid等属性对可视化结果的标题、坐标轴、图例、格网等进行显示，以坡度的可视化为例，代码如下：

% 可视化结果
figure;
imagesc(slope);
colorbar;
title('Slope');

五实验结果

（1）坡度提取结果

（2）坡向提取结果

（3）坡度变率提取结果

（4）坡向变率提取结果

（5）曲率提取结果

（6）粗糙度提取结果

（7）凹凸系数提取结果

六说明

最后，在本次实验中对DEM提取面积和体积主要用到： geotiffread函数读取DEM 数据； abs 函数和R.cellExtentInlworldx、R.cellExtentInWorldY属性计算每个像元的面积；使用ones函数将每个像元的面积乘以高程赋值给一个与 DEM 数据相同大小的矩阵；使用sum函数将所有像元的体积相加，得到DEM的体积。

本次实验设计的算法能够有效地提取出DEM数据中的地形特征，并能够对地形特征进行可视化和分析。该算法在地形分析、地质勘探、环境监测等领域具有广泛的应用前景。

（需要完整代码及运行结果请留言~