基于MATLAB的DEM读取显示以及山顶点提取

前言

利用MATLAB强大的数组运算的功能，我们可以利用其进行山顶点的提取。但由于知识比较薄弱，在MATLAB正确读取并显示DEM模型上还未能找到合适的方法，因此在下面我采用的是常规导入并进行拉伸显示。本篇文章只是一个小实验，希望读者能从中获益，若有不足之处，望能指出。

一、山顶点提取原理

所谓山顶点，即是该部分地形区域上的高程最大点，提取山顶点对地形研究有很重要的意义。通过提取山顶点来获取山顶点的几何位置，在GIS分析中可用来作为选址分析或者地形分析。下图是在arcmap中 提取山顶矢量点，可以看出，其提取的结果还是比较满意，但在边缘区域就会出现多余点的情况，而且还会出现多个山顶点重叠的情况。

下面我们就来谈一谈山顶点提取原理，我们利用的源数据是DEM图，其每个像元值代表的是一个高程值。那如何找到一定区域内的最高点呢，按照下述流程，我们可以得到差值矩阵图。

此差值矩阵中的零值就代表的是山顶点，其他值都是低于这个点的其他高程点。我们进行重分类，将零值换为0，其他值赋值为1，通过转换成double数组，就可以得到一个二值化图。

二、MATLAB代码实现（亲测有效）

1.程序以及程序清单

程序清单：
geotiffread：用于读取地理tiff文件
padarray：用于在原矩阵四周扩充零矩阵
deal：用于分身
reshape：改变数组的维度不改变数值
flip：颠倒数组
find：查找数组某值得具体行列位置
numel：计算数组的元素个数

clc;
clear all;
%geotiffread函数可用于读取地理tiff文件格式图片
% data变量存储坐标或是像元数据，info变量用于存储数据的一些详细信息
[data,info]=geotiffread("DEM1.tif");
%拉伸图像，由于dem数据中边缘存在异常值，因此一般要进行拉伸才能够显示，拉伸的大小视实际情况而定
image=data.*30;
%-------------------------%
imshow(image);
colorbar;
title('DEM image')
%-------------------------%
%padarrary函数用于给图像四周扩充零矩阵，扩充的具体范围视窗口大小而定
image1=padarray(image,[7,7]);
%deal函数用于将数组赋给两个变量，其继承了原矩阵的所有属性值，目的是保证后续不会被迭代改变
[x1,x2]=deal(image1);
[rows,cols]=size(image1);
%采用15*15的窗口遍历每一个像元，将窗口内的最大值赋值给该像元
for i=1:rows-14for j=1:cols-14NR=x1(i:i+14,j:j+14);ER=reshape(NR,[1,15*15]);MAX_NR=max(ER);x2(i+6,j+6)=MAX_NR;end
end
diff_image=image1-x2;
[xx1,xx2]=deal(diff_image);
for ii=1:rowsfor jj=1:colsif xx1(ii,jj)<0xx2(ii,jj)=0;elsexx2(ii,jj)=1;endend
end
top=double(xx2);
%------------------------------------------------------------%
%---------------------------------------%
%查找异常值，并对图像进行裁剪
[r,c]=size(data);
first_row_find=0;end_row_find=0;first_col_find=0;end_col_find=0;
first_row=data(1,1:c);
end_row=flip(data(r,1:c));
first_col=flip(data(1:r,1));
end_col=data(1:r,c);
for mm=1:cif first_row(mm)==max(data)first_row_find=first_row_find+1;elsebreakendif end_row(mm)==max(data)end_row_find=end_row_find+1;elsebreakend
end
for nn=1:rif first_col(nn)==max(data)first_col_find=first_col_find+1;elsebreakendif end_col(nn)==max(data)end_col_find=end_col_find+1;elsebreakend
end
%---------------------------------------%
%------------------------------------------------------------%
crop_top=top(7+first_col_find:rows-7-end_col_find,7+first_row_find:cols-7-end_row_find);
imshow(crop_top);
title('the crop top.point image')
%根据find函数查找corp_top图像上山顶点的行列坐标
[rows1,cols1]=find(crop_top==1);
[rows2,cols2]=size(data);
ol=data(first_col_find:rows2-end_col_find,first_row_find:cols2-end_row_find).*45;
imshow(ol);
title('the crop DEM image')
%------------------------------------------------------------%
png=data(first_col_find:rows2-end_col_find,first_row_find:cols2-end_row_find).*45;
png1=deal(png);
number=0;
for iii=1:size(rows1)number=number+1;png1(rows1(number),cols1(number))=1;
end
imshow(png1);
title('the mix top.point image');
%------------------------------------------------------------%
count=numel(find(crop_top==1))

2.运行结果以及分析

下图为源数据DEM图，该图的高程以灰白色条显示，数据类型为int16。由于我的源数据是掩膜提取的，所以该图的四周含有异常值，因此色条的范围很大，另外要在MATLAB中显示，则需要对图像进行拉伸才能显示正常。

针对上图的异常显示，我们需要对图像进行裁剪，去除四周的异常值，结果如下

本文主要采用的窗口大小为15*15，通过对源数据DEM进行遍历之后，再进行相减以及二值化处理，可以得到下面的图，图中的白色点即为山顶点，其高程在一定区域内最高。

下图即按照上图的白色点坐标，在the crop DEM image中显示的结果，不难看出，其山顶点的位置还是较为正确，而且不存在边缘多余山顶点的情况。