DEM实验|基本地形因子的提取和DEM的可视化
试验目的
了解基本地形因子的分类及提取原理，熟悉 ArcGIS 软件中各种因子的提取方法。
实验环境与准备
ArcGIS软件
实验数据 : 1 : 1 万 DEM 数据
实验内容和步骤
坡度、坡向的提取
坡度
该点的切平面与水平地面的夹角，表示地表面在该点的倾斜程度
【工具】: Spatial Analyst>表面分析>坡度
原始dem：（以下因子分析均基于此dem）
z：高程变换系数
如果 z 单位是英尺而 x,y 单位是米，则应使用 z 因子 0.3048 将 z 单位从英尺转换为米（1 英尺 = 0.3048 米）。
坡向
地表面上一点的切平面的法线矢量在水平面的投影与过该点的正北方向的夹角，表示该点高程值该变量的最大变化方向
【工具】：Spatial Analyst>表面分析>坡向
变率的提取
坡度变化率
Slope of Slope of DEM，即地面坡度的变化
反映地面有没有一个陡坎，有变化的地方才出现大值
反映的是在剖面上的棱角
在相当程度上和剖面曲率是一致的
【工具】：Spatial Analyst>表面分析>坡度
【方法】：二次求导，即求两次坡度
坡向变化率的提取
Slope of Aspect of DEM，即地面坡向的坡度，可以很好地反映等高线的弯曲程度
等高线有转折的地方取得大值，等高线有波动的地方
在山脊线或山谷线SOA取得大值
平面曲率与SOA在数字上相等
误差小的SOA求法:（（[SOA1]+[SOA2]）-Abs（[SOA1]-[SOA2]））/2 （SOA1：原始DEM的坡向变率；SOA2：反地形的坡向变率）
【工具】：
1、Spatial Analyst>邻域分析>焦点统计
2、Spatial Analyst>地图代数>栅格计算器
3、Spatial Analyst>表面分析>坡向
4、Spatial Analyst>表面分析>坡度
【方法】：
1、栅格计算器求DEM反地形（MAX-DEM）
2、对反DEM求得坡度变化率（先坡向再坡度）
3、栅格计算器输入无误差SOA计算公式
曲率的提取
剖面曲率：对地面坡度的沿最大坡降方向地面高程变化率的度量
平面曲率：过地形表面上某点得水平面沿水平方向切地形表面所得的曲线在该点的曲率值，即该点所在的地面等高线的弯曲程度
【工具】：Spatial Analyst>表面分析>曲率
综合曲率结果图
剖面曲率结果图
平面曲率结果图
坡面形态因子的提取
坡形
坡面复杂因子的提取
地形起伏度
在指定区域内，最高海拔点和最低海拔点的差值
【工具】：
1、Spatial Analyst>邻域分析>焦点统计
2、Spatial Analyst>地图代数>栅格计算器
【方法】：
1、焦点统计求DEM的最大值和最小值
2、直接利用焦点统计求像元范围，统计类型设为RANGE
地形粗糙度
使用之前计算坡度数据来提取地形粗糙度，需要注意的是在ArcGIS中，cos使用弧度值作为角度单位，但是使用【坡度】工具提取的坡度值是角度单位，所以在计算时应该提前把角度单位转化为弧度单位，即在角度值后面乘上Π/180。
【工具】：
1、Spatial Analyst>表面分析>坡度
2、Spatial Analyst>地图代数>栅格计算器
【方法】：
1、求DEM的坡度
2、输入地表粗糙度公式求解
地表切割深度
地面某点的邻域范围的平均高程与该邻域范围内的最小高程的差值,反映地表被侵蚀切割的情况
【工具】：
1、Spatial Analyst>邻域分析>焦点统计
2、Spatial Analyst>地图代数>栅格计算器
【方法】：
1、焦点统计求DEM的平均值和最小值
2、求平均值与最小值的
高程变异系数
反映分析区域内地表单元格网各顶点高程变化的指标，它以格网单元顶点的标准差与平均高程的比值来表示
【工具】：
1、Spatial Analyst>邻域分析>焦点统计
2、Spatial Analyst>地图代数>栅格计算器
【方法】：
1、焦点统计求DEM的标准差和平均值
DEM地形渲染
DEM高程分层设色
DEM地形渲染
美术透视技术已经告诉我们怎么样在一个二维纸面上表达一个三维立方体，那就是通过多个侧视的光线汇聚，或者阴影方式来模拟多一个维度。这一原理可以被充分利用在各种制图中。对这个DEM数据进行【山体阴影】的操作，生成新的DEM数据。如下图，为使用【山体阴影】工具处理后的效果。
操作过后非常明显的展示了地形的起伏，而且这种模拟，类似太阳在某一侧照射，产生了另一侧的阴影，使得山脊非常的明显；凹陷的部分是河流，由于河流在地位，所以没有任何的阴影像素产生，所以效果也是比较明显的。人眼对这种做过做出处理的假数据是非常敏感的，比之前的真实DEM辨别起来要容易的多。
DEM地形组合渲染
原始DEM对着色有着非常连贯的效果，也就是说色彩过度非常好；山体阴影数据对起伏非常敏感，高差显示明显。只要将这两个数据充分中和一下就非常好了。
所以可以将色彩连贯的原始数据放在最上层，高地起伏的山体阴影数据放在下层。然后对原始数据使用色彩渲染，同时在【图层属性】—【显示】中，将透明度设置为30%左右，得到最终效果如下图所示。
三维可视化
在ArcScene中将两个数据高度进行两倍拉伸，结果如下
DEM 3D可视化
DEM地形明暗等高线可视化
根据入射光方向将坡向划分为背光面和受光面两个部分假。光源位于区域西北方向，则可将坡向为0-45°、225-360° 的部分划为受光面，坡向为45-225° 的部分分为背光面。可以借助(Spatial Analyst) 工具栏的重分类工具得到受光面和背光面的二值图数据。
其次，将受光面和背光面的二值图数据转为矢量多边形数据。
最后，使用叠置分析工具，将等高线数据与向背光面数据进行相交处理，获得明暗等
高线数据。再将明暗等高线的符号样式更改为，明暗两种颜色。
得到明暗等高线如下
在ArcScene中将明暗等高线叠加地形晕渲图
附加实验
坡度因子提取的不确定性
DEM 数据的格网分辨率、精度、因子提取算法等对各因子的提取结果都有影响, 即不确定性。这里以坡度因子为例来探讨, 首先基1: 1万 DEM数据, 通过resample得到10m、20m、30m、…100m格网分辨率的系列DEM 数据, 分别基于这些数据堤取坡度, 通过对比坡度直方图及各个统计数据的变化特征来分析坡度提取的不确定性。
GDAL简介
Geospatial Data Abstraction Library （GDAL）是使用C/C++编写的用于读写空间数据的一套跨平台开源库。现有的大部分GIS或者遥感软件，不论是商业软件ArcGIS，ENVI还是开源软件GRASS，QGIS，都使用了GDAL作为底层构建库。
GDAL库由OGR和GDAL项目合并而来，OGR主要用于空间要素矢量矢量数据的解析，GDAL主要用于空间栅格数据的读写。此外，空间参考及其投影转换使用开源库 PROJ.4进行。
目前，GDAL主要提供了三大类数据的支持：栅格数据，矢量数据以及空间网络数据（Geographic Network Model）。
实验内容与步骤
1. 批量重采样
因为需要对一张图像重复重采样很多次，为了避免重复劳动，使用Python + GDAL库对dem_grid进行批量重采样，得到从10m至100m不同分辨率的数据。
2. 批量计算坡度
因为需要对多张图像提取坡度，同样为了避免重复劳动，使用Python + GDAL库对不同分辨率的DEM数据提取坡度，得到从10m至100m不同分辨率的坡度数据。
3. 批量绘制直方图并统计
还是为了避免重复劳动，使用Python + GDAL库对不同分辨率的坡度数据，统计他们的方差，均值，最大最小值等数据，并绘制灰度直方图。
观察直方图可以发现，随着分辨率的降低，坡度的灰度分布仍然大体上一致，因此仍依旧能够表达地形的整体趋势，但是地形的细节信息会逐渐消失。从分辨率这方面来分析，分辨率的降低导致DEM表示细部特征的能力降低，细节消失，仅仅表达区域整体高程变化，保留了宏观性、骨架性的地形特征，同时高程变化大的区域因为被概括而趋于平坦化。
随着分辨率的降低，受影响最大的为空间高频组分，也即高程变化最大的部分，受影响最小的是空间低频组分，也即较平坦的部分，造成坡度变化的原因是高频分量的丢失，从直方图中也可以看出，随着分辨率的降低，坡度的最大值越来越小。
4. 绘制不同分辨率数据的对比图
最后，绘制三种统计数据，使用pyecharts绘制对比图，并进行对比分析。其中平均坡度随DEM分辨率降低而降低，且与分辨率大致呈线性关系。
由最大值和最小值可以得到值域，可以看出，随着分辨率的降低，值域在减小，方差也在减小。值域和方差都对 DEM 分辨率敏感，分辨率对值域和方差结果影响很大，在使用时要考虑DEM分辨率的影响。这是因为随着分辨率的降低，平滑程度作用使地貌细部特征被简化，采样间距作用使高程变化大区域被概括而趋于平坦化，所以由 DEM 计算得到坡度值值变得单一，表面也趋于平滑，一些极值被综合掉，所以值域在减小，方差变小，分布趋于集中。
四、附录
实验使用代码
批量重采样部分代码
# -*- coding: utf-8 -*-
"""
PROJECT_NAME: RS_Toolbox
FILE_NAME: batch_resample
AUTHOR: welt
E_MAIL: tjlwelt@foxmail.com
DATE: 2022/10/31
"""
from osgeo import gdal
import os
import numpy as np
path = r"D:\RS_Toolbox"
os.chdir(path)
image_name = 'dem_grid.tif'
times = [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
for time in times:
in_ds = gdal.Open(image_name)
geotrans = list(in_ds.GetGeoTransform())
geotrans[1] *= time # 像元宽度变为原来的两倍
geotrans[5] *= time # 像元高度也变为原来的两倍
in_band = in_ds.GetRasterBand(1)
xsize = in_band.XSize
ysize = in_band.YSize
x_resolution = int(xsize / time) # 影像的行列都变为原来的一半
y_resolution = int(ysize / time)
(filename, extension) = os.path.splitext(image_name)
outname = filename + '_' + str(5 * time) + 'm' + extension
if os.path.exists(outname):
# 如果存在重采样后的影像，则删除之
os.remove(outname)
out_ds = in_ds.GetDriver().Create(outname, x_resolution, y_resolution, 1,
in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection()) # 设置投影信息
out_ds.SetGeoTransform(geotrans) # 设置地理变换信息
data = np.empty((y_resolution, x_resolution), np.int) # 设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值
in_band.ReadAsArray(buf_obj=data)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)
out_band.FlushCache()
out_band.ComputeStatistics(False) # 计算统计信息
out_ds.BuildOverviews('average', [1, 2, 4, 8, 16, 32]) # 构建金字塔
del in_ds
del out_ds
print("This process has succeeded!")
批量求坡度部分代码
# -*- coding: utf-8 -*-
"""
PROJECT_NAME: RS_Toolbox
FILE_NAME: batch_slope
AUTHOR: welt
E_MAIL: tjlwelt@foxmail.com
DATE: 2022/10/31
"""
from osgeo import gdal
import os
import sys
from osgeo import osr
import numpy as np
# 给栅格最外圈加一圈
def assignBCs(elevGrid):
ny, nx = elevGrid.shape
Zbc = np.zeros((ny + 2, nx + 2))
Zbc[1:-1, 1:-1] = elevGrid
Zbc[0, 1:-1] = elevGrid[0, :]
Zbc[-1, 1:-1] = elevGrid[-1, :]
Zbc[1:-1, 0] = elevGrid[:, 0]
Zbc[1:-1, -1] = elevGrid[:, -1]
Zbc[0, 0] = elevGrid[0, 0]
Zbc[0, -1] = elevGrid[0, -1]
Zbc[-1, 0] = elevGrid[-1, 0]
Zbc[-1, -1] = elevGrid[-1, 0]
return Zbc
# 计算dx,dy
def calcFiniteSlopes(elevGrid, dx):
Zbc = assignBCs(elevGrid)
Sx = (Zbc[1:-1, :-2] - Zbc[1:-1, 2:]) / (2 * dx) # WE方向
Sy = (Zbc[2:, 1:-1] - Zbc[:-2, 1:-1]) / (2 * dx) # NS方向
return Sx, Sy
# 投影转换
def convertProjection(data, filename, dx):
landsatData = gdal.Open(filename, gdal.GA_ReadOnly)
oldRef = osr.SpatialReference()
oldRef.ImportFromWkt(data.GetProjectionRef())
newRef = osr.SpatialReference()
newRef.ImportFromWkt(landsatData.GetProjectionRef())
transform = osr.CoordinateTransformation(oldRef, newRef)
tVect = data.GetGeoTransform()
nx, ny = data.RasterXSize, data.RasterYSize
(ulx, uly, ulz) = transform.TransformPoint(tVect[0], tVect[3])
(lrx, lry, lrz) = transform.TransformPoint(tVect[0] + tVect[1] * nx, tVect[3] + tVect[5] * ny)
memDrv = gdal.GetDriverByName('MEM')
dataOut = memDrv.Create('name', int((lrx - ulx) / dx), int((uly - lry) / dx), 1, gdal.GDT_Float32)
newtVect = (ulx, dx, tVect[2], uly, tVect[4], -dx)
dataOut.SetGeoTransform(newtVect)
dataOut.SetProjection(newRef.ExportToWkt())
# Perform the projection/resampling
res = gdal.ReprojectImage(data, dataOut, oldRef.ExportToWkt(), newRef.ExportToWkt(), gdal.GRA_Cubic)
return dataOut
def slope(DEMFilename, slopeFilename):
gdal.AllRegister()
data = gdal.Open(DEMFilename, gdal.GA_ReadOnly)
if data is None:
print('Cannot open this file:' + DEMFilename)
sys.exit(1)
geotrans = list(data.GetGeoTransform())
dx = geotrans[1]
# 投影变换
projData = convertProjection(data, DEMFilename, dx)
projinfo = data.GetProjection()
geotransform = data.GetGeoTransform()
gridNew = projData.ReadAsArray().astype(np.float)
Sx, Sy = calcFiniteSlopes(gridNew, dx)
# 坡度计算
slope = np.arctan(np.sqrt(Sx ** 2 + Sy ** 2)) * 57.29578
# 输出坡度坡向文件
driver = gdal.GetDriverByName('GTiff')
if os.path.exists(slopeFilename):
os.remove(slopeFilename)
ds1 = driver.Create(slopeFilename, slope.shape[1], slope.shape[0], 1, gdal.GDT_Float32)
ds1.SetGeoTransform(geotransform)
ds1.SetProjection(projinfo)
band = ds1.GetRasterBand(1)
band.WriteArray(slope, 0, 0)
del ds1
def main():
dem_path = 'D:/RS_Toolbox/examples'
slope_path = 'D:/RS_Toolbox'
filelist = os.listdir(dem_path)
for file in filelist:
(filename, extension) = os.path.splitext(file)
outname = filename + '_slope' + extension
dem_file = dem_path + '/' + file
slopefile = slope_path + '/' + outname
slope(dem_file, slopefile)
if __name__ == '__main__':
main()
批量绘制直方图部分代码
# -*- coding: utf-8 -*-
"""
PROJECT_NAME: RS_Toolbox
FILE_NAME: batch_his
AUTHOR: welt
E_MAIL: tjlwelt@foxmail.com
DATE: 2022/11/1
"""
import os
import re
import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal
# 通过GDAL读取栅格影像
filepath = 'D:/RS_Toolbox/slope'
filelist = os.listdir(filepath)
for file in filelist:
(name, extension) = os.path.splitext(file)
filename = filepath + '/' + file
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(name)
print("平均值: ", mean_value)
print("方差: ", variance)
print("最小值: ", data.min())
print("最大值: ", data.max())
# 根据影像中最大最小值设定坐标轴
bins = np.linspace(data.min(), data.max(), 100)
resolution = re.findall(r"\d+\.?\d*", name)
# 绘制直方图，设定直方图颜色
plt.hist(data, bins, facecolor="blue")
# 横坐标轴名称
plt.xlabel('像元值')
# 纵坐标轴名称
plt.ylabel('频数')
# 图表头名称
plt.title(resolution[0] + 'm分辨率坡度数据的' + '灰度分布直方图')
# 显示中文字体
plt.rcParams["font.sans-serif"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False
# 导出绘制得到的图片
plt.tight_layout()
plt.savefig(name + '.jpg')
plt.show()
