DEM实验|DEM的建立

试验目的

通过DTM的创建与插值、DTM的模型转化实验深入掌握数字地面模型的概念与建模方法。

实验环境与准备

PC、ArcGIS软件

实验内容和步骤

基于不规则、 规则分布采样点的DEM建立

导入数据

导入散点数据（feapt）

插值分析

打开【ArcToolBox】, 选择【3D Analyst Tools】, 选择【栅格插值】，这里使用反距离权重法对散点数据进行插值

image-20221013153625991

结果生成

使用反距离权重法对散点数据进行内揷计算生成 DEM如下：

DEM6

基于等高线数据的 DEM建立

数据导入

添加等高线数据（terlk)

线要素转点要素

使用Feature Vertices To Points（要素折点转点）工具，将线要素转为点要素

image-20221017212320869

生成的点要素如下：

image-20221017212916044

插值分析

使用反距离权重法对线要素转成的点要素进行插值分析：

image-20221017212531258

结果生成

使用反距离权重法对等高线数据进行内揷计算生成 DEM如下：

DEM7

TIN的建立及TIN与GRID的转换

TIN的建立

数据导入

添加等高线数据（terlk) 与散点数据（feapt）

创建TIN

打开【ArcToolBox】, 选择【3D Analyst Tools】, 选择【TIN Creation】, 双击【Create TIN】。在弹出的窗口中,【Output TIN】中输人新建 TIN 的名称及保存路径; 【Spatial Reference 】输人 TIN 的空间坐标, 没有则不填。【Input Feature Class】输人构建 TIN 的矢量数据层, 这里选择了等高线（terlk），在属性栏选项中, height_field 属性选择高程 (ELVE)，SF_type属性选择软断线 (softline)，其他都为缺省值。

image-20221013162130498
由矢量数据生成TIN

单独使用等高线数据生成TIN，同时使用等高线数据和散点数据组合生成TIN，比较二者的区别：

由等高线数据生成的TIN如下：

DEM8

由等高线数据和散点数据生成的TIN如下：

DEM11

TIN转换为GRID

数据导入

导入之前由等高线数据数据生成的TIN

TIN转换为GRID

打开【ArcToolBox】, 选择【3D Analyst Tools】, 选 择【Conversion】, 选择【From TIN】, 双击【TIN to Raster】
在弹出的窗口中,【Input TIN】中输人需要被转换的 TIN 数据; 【Output Raster】 中输人计算结果数据的名称和路径; 【Output Data Type (optional) 】选择转换得到棚格数据的高程数据属性, 有浮点型 (FLOAT) 和整型 (INT) 两种可选;【Method (optional)】选择内揷方法, 有线性内揷 (LINEAR) 和邻域内揷 (NATURAL_NEIGHBORS) 两种可选; 【Sampling Distance (optional)】输人栅格的分辨率; 【Z Factor (optional) 】输人垂直放 大因子, 通常选 1, 表示维持原始地形。

image-20221013165819107
结果生成

由TIN生成的GRID如下图所示。

DEM5

GRID转换为TIN

数据导入

导入栅格DEM数据（dem）

GRID转换为TIN

打开【ArcToolBox】, 选择【3D Analyst Tools】, 选 择【Conversion】, 选择【From Raster】, 双击【Raster to TIN】。在弹出的窗口中, 【Input Raster】中输人需要被转换的 DEM 数据; 【Output TIN】 中输人计算结果数据的名称和路径; 【Z Tolerance (optional) 】 为高程容限值; 【 Maximum Number of Points (optional)】为可选值; 【Z Factor (optional)】输人垂直放大因子, 通常选 1 , 表示维持原始地形。

image-20221017123412156
结果生成

由GRID转换而来的TIN如下所示：DEM11

提取等高线

数据导入

导入栅格DEM数据（dem）

从GRID中提取等高线

打开【ArcToolBox】, 选择【3D Analyst Tools】, 选择【栅格表面】, 选择【等值线】。在弹出的窗口中, 【输入栅格】中输入需要提取等高线的DEM数据; 【等值线间距】中输入等值线间的间距；【起始等值线】中输入等值线的起始值。

image-20221018003715237

结果生成

使用不同的等值线间距提取等高线，并比较它们之间的区别

由DEM数据，等值线间的间距为10m提取的等高线如下：

DEM3

由DEM数据，等值线间的间距为5m提取的等高线如下：

DEM12

实验结论

对于DEM的建立，定性的来看，由等高线和TIN建立的DEM都比较精细，对于地表的凸凹情况还原的比较好。而由散点建立的DEM则比较粗糙，丢失了许多细节。

对于TIN的建立，单独使用由等高线数据生成的TIN和使用等高线数据和散点数据生成的TIN在大体上没有什么变化。但在细微的地方，使用等高线数据和散点数据生成的TIN，细节更加丰富，表面更加粗糙，不像单独使用等高线一样光滑。因此可以认为使用等高线数据和散点数据生成的TIN细节更加丰富。

对于等高线的提取，使用较大的等值线间距时，提取出的等高线条数较少，有部分细节缺失，但二者大体上的形状均是一致的。使用较小的等值线间距时，提取出等高线的条数较多，细节更丰富，特别是对于谷底和山顶的细节。

从以上的对比可明显看出不管时对DEM的建立还是TIN的建立，还是等高线的绘制，不同的方法图像呈现的效果不同，精度也不同。

思考与总结

不同来源数据DEM建立的误差分析

以现有 1：1 万 DEM 数据为真值, 利用同一区域的等高线数据、 散点数据、TIN 各生成一幅该区域 DEM, 将不同数据源生成的 DEM 与现有 DEM 数据相减, 统计结果的最大值、最小值、平均值、均方差等统计值, 使用Python+GDAL库实现了栅格图像的相减与相减后图像最大值、最小值、平均值、均方差等统计值的获取，并绘制了DEM图像的灰度直方图。

等高线数据生成的该区域DEM误差分析

等高线数据生成的该区域DEM与现有DEM数据相减后得到的图像
DEM2

图中白色或者黑色的区域为生成的DEM与现有DEM相比差距较大的地方，可以发现在这张图中，大部分区域均为灰色，也即生成的DEM与现有的DEM相差较小，因此等高线数据生成的该区域DEM误差较小。可以发现图中白色和黑色的区域主要分布于山谷的谷底、山峰顶部和山脊等部位。个人推测出现这种现象是由于在山顶和山谷谷底缺少相应的等高线，导致在谷底生成的DEM较标准DEM的值偏高，而在山顶生成的DEM较标准DEM的值偏低。区域中等高线的分布情况如下图所示：

DEM10
等高线数据生成的该区域DEM与现有DEM数据相减后的统计值

平均值: 1.0985761
方差: 0.81780744
最小值: 0.00012207031
最大值: 8.347046

等高线数据生成的该区域DEM灰度直方图

可以发现，等高线数据生成的该区域DEM灰度直方图中，中等的DEM值这一范围内的点是最多的，大致呈现出近似正态分布的特征。图像中存在频数特别大的DEM值，这里个人推测是由于插值所导致的。

image-20221017221117610

散点数据生成的该区域DEM误差分析

散点数据生成的该区域DEM与现有DEM数据相减后得到的图像
DEM1

图中白色或者黑色的区域为生成的DEM与现有DEM相比差距较大的地方，可以发现在这张图中，特征线被比较清晰的刻画出来，个人推测出现这种现象是因为采样点主要分布在特征线上，因此特征线附近生成的DEM比较准确。但特征线之外的区域，由于采样点的数量较少，导致生成的DEM与标准DEM之间存在较大的差距，在图像的上部的中间区域与图像的右下角，这种情况尤为明显。区域中采样点的分布情况如下图所示：

DEM9
散点数据生成的该区域DEM与现有DEM数据相减后的统计值

与平均值: 4.058907
方差: 20.889515
最小值: -85.730225
最大值: 95.529175

散点数据生成的该区域DEM灰度直方图

可以发现，散点数据生成的该区域DEM灰度直方图中，分布较不均匀，DEM值在1150附近的点的数量明显偏少。

image-20221017124414095

TIN生成的该区域DEM误差分析

TIN生成的该区域DEM与现有DEM数据相减后得到的图像
DEM

图中白色或者黑色的区域为生成的DEM与现有DEM相比差距较大的地方，可以发现在这张图中，大部分区域均为灰色，也即生成的DEM与现有的DEM相差较小，因此TIN生成的该区域DEM误差较小。与由等高线生成的DEM类似，可以发现图中白色和黑色的区域主要分布于山谷的谷底、山峰顶部和山脊等部位。个人推测出现这种现象是由于使用的TIN是由等高线生成的，而在山顶和山谷谷底缺少相应的等高线，导致在谷底生成的DEM较标准DEM的值偏高，而在山顶生成的DEM较标准DEM的值偏低。同时可以发现图像中出现了明显的分块的现象，导致这种现象推测是因为在由TIN生成DEM的过程中使用的是区域插值算法导致的。

TIN生成的该区域DEM与现有DEM数据相减后的统计值

平均值: 0.007681254
方差: 1.7352749
最小值: -10.950562
最大值: 11.114014

TIN生成的该区域DEM灰度直方图

可以发现，TIN生成的该区域DEM灰度直方图中，中等的DEM值这一范围内的点是最多的，大致呈现出近似正态分布的特征，且不像等高线数据生成的该区域DEM灰度直方图一样有很多极其突出的部分。推测是因为在TIN转GRID时，使用区域插值的缘故。

image-20221017161739668

不同来源数的DEM与现有DEM数据相减后的各项统计指标对比

Bar_DEM

可以发现，在方差和平均值上，散点数据得到的数据均高于其他两组，同时其数据的最小值和最大值之间的差值也更大，因此由散点数据生成的DEM是误差最大的，其次是由TIN生成的DEM，误差最小的是由等高线数据生成的DEM。可以得知，在本次实验中，使用等高线转点后插值生成DEM的方法最好。

附录

实验代码

图像相减部分代码

# -*- coding: utf-8 -*-
"""
PROJECT_NAME: RS_Toolbox 
FILE_NAME: Raster_subtraction 
AUTHOR: welt 
E_MAIL: tjlwelt@foxmail.com
DATE: 2022/10/17 
"""

from osgeo import gdal
import numpy as np
from tqdm import tqdm

def subtract_raster(DEM_standard_path, DEM_path, out_path):
   dataset = gdal.Open(DEM_standard_path)
   projinfo = dataset.GetProjection()
   geotransform = dataset.GetGeoTransform()

   format = "GTiff"
   driver = gdal.GetDriverByName(format)
   depth = dataset.RasterCount
   rows = int(dataset.RasterYSize)
   colmns = int(dataset.RasterXSize)

   DEM_standard = dataset.ReadAsArray()
   src2 = gdal.Open(DEM_path)
   DEM = src2.ReadAsArray()
   Result = np.zeros((rows, colmns))
   Result_DEM = driver.Create(out_path, colmns, rows, depth, gdal.GDT_Float32)
   Result_DEM.SetGeoTransform(geotransform)
   Result_DEM.SetProjection(projinfo)


   for row in tqdm(range(rows)):
      for col in range(colmns):
         Result[row, col] = DEM_standard[row, col] - DEM[row, col]

   Result_DEM.GetRasterBand(1).WriteArray(Result)


if __name__ == '__main__':
   # 修改路径
   Standard_DEM_Path = r"dem_standard.tif"  # 标准的DEM
   MY_DEM_Path = r"Extract_Idw_fleap.tif"  # 由点，线，TIN得到的DEM
   OutTif = r"out.tif"
   subtract_raster(Standard_DEM_Path, MY_DEM_Path, OutTif)

计算均值，方差，最值，绘制灰度直方图部分的代码

# -*- coding: utf-8 -*-
"""
PROJECT_NAME: RS_Toolbox 
FILE_NAME: histogram_statistics 
AUTHOR: welt 
E_MAIL: tjlwelt@foxmail.com
DATE: 2022/10/17 
"""

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

# 通过GDAL读取栅格影像
filename = "out_terlk.tif"
dataset = gdal.Open(filename)
im_width = dataset.RasterXSize
im_height = dataset.RasterYSize
im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
print(im_data.shape)

# 显示灰度直方图
# 遍历影像中的每一个像元的像元值
data = []
for i in range(im_data.shape[0]):
	for j in range(im_data.shape[1]):
		if im_data[i][j] > 0:
			data.append(im_data[i][j])
data.sort()

# 统计最大最小值
data = np.array(data)
mean_value = np.mean(data)
variance = np.std(data)
print("平均值: ", mean_value)
print("方差: ", variance)
print("最小值: ", data.min())
print("最大值: ", data.max())

# 根据影像中最大最小值设定坐标轴
bins = np.linspace(data.min(), data.max(), 100)
# 绘制直方图，设定直方图颜色
plt.hist(data, bins, facecolor="blue")
# 横坐标轴名称
plt.xlabel('像元值')
# 纵坐标轴名称
plt.ylabel('频数')
# 图表头名称
plt.title('灰度分布直方图')
# 显示中文字体
plt.rcParams["font.sans-serif"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False
# 导出绘制得到的图片
plt.savefig('./test.jpg')
plt.show()

绘制不同来源数据DEM的统计指标对比图部分的代码

# -*- coding: utf-8 -*-
"""
PROJECT_NAME: PLOT 
FILE_NAME: Bar_DEM 
AUTHOR: welt 
E_MAIL: tjlwelt@foxmail.com
DATE: 2022/10/17 
"""

from pyecharts import options as opts
from pyecharts.charts import Bar

if __name__ == '__main__':
   bar_base = (
      Bar(init_opts=opts.InitOpts(bg_color='white'))
         .add_xaxis(['平均值', '方差', '最小值', '最大值'])
         .add_yaxis("等高线", [1.0986, 0.8178, 0.0001, 8.3471], gap="0%")
         .add_yaxis("散点", [4.0589, 20.8895, -85.7302, 95.5292], gap="0%")
         .add_yaxis("TIN", [0.0077, 1.7353, -10.9506, 11.1140], gap="0%")
         .set_global_opts(
         title_opts=opts.TitleOpts(title="不同来源数据DEM的统计指标分析图", subtitle=""),
      )
         .render("bar_DEM.html")
   )
posted @ 2022-10-18 22:07  Weltㅤ  阅读(1196)  评论(0)  编辑  收藏  举报