数字地面模型|DEM的建立

试验目的

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

实验环境与准备

PC、ArcGIS软件

实验内容和步骤

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

导入数据

导入散点数据（feapt）

插值分析

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

结果生成

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

基于等高线数据的 DEM建立

数据导入

添加等高线数据（terlk)

线要素转点要素

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

生成的点要素如下：

插值分析

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

结果生成

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

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)，其他都为缺省值。

由矢量数据生成TIN

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

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

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

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, 表示维持原始地形。

结果生成

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

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 , 表示维持原始地形。

结果生成

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

提取等高线

数据导入

导入栅格DEM数据（dem）

从GRID中提取等高线

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

结果生成

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

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

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

实验结论

对于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数据相减后得到的图像

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

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

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

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

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

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

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

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

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

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

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

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

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

TIN生成的该区域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时，使用区域插值的缘故。

不同来源数的DEM与现有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")
   )

原文地址：http://www.cnblogs.com/tangjielin/p/16804337.html