(Python)根据经纬度从数字高程模型（DEM）文件获取高度

基本介绍

在地理信息系统（GIS）和遥感中，数字高程模型（Digital Elevation Model，简称DEM）是一种表示

地表或地形高程信息的重要数据。DEM数据通常以栅格（raster）形式存在，其中每个像素的值表示相

应地理位置的高程。

之前介绍了如何进行6S大气校正，其中有一段代码是自动获取DEM的高度。

原理及过程是，输入参数分别是经纬度、DEM文件。输入参数是高度。

先把DEM文件读取为栅格，再把经纬度转换为行列号，根据行列号定位到DEM栅格，读取该栅格的数值，返回数据。

代码

以下是我们的Python函数：

def fromDemFileGetHeight(dem, lon, lat):try:DEMIDataSet = gdal.Open(dem)except Exception as e:print('Missing DEM file')passDEMBand = DEMIDataSet.GetRasterBand(1)geotransform = DEMIDataSet.GetGeoTransform()# DEM分辨率pixelWidth = geotransform[1]pixelHight = geotransform[5]# DEM起始点：左上角，X：经度，Y：纬度originX = geotransform[0]originY = geotransform[3]# 研究区左上角在DEM矩阵中的位置yoffset1 = int((originY - lat) / pixelWidth)xoffset1 = int((lon - originX) / (-pixelHight))DEMRasterData = DEMBand.ReadAsArray(xoffset1, yoffset1, 1, 1)DEMRasterData = np.mean(DEMRasterData)return DEMRasterData

代码解析

这个函数接受三个参数：DEM文件路径、经度和纬度。它返回指定经纬度位置的高度值。

首先，函数尝试打开DEM文件。如果文件不存在或无法打开，它将打印一条错误消息并退出。

然后，函数获取DEM数据的第一波段（在大多数DEM数据中，高程数据都存储在第一波段）。它还获取DEM数据的地理转换参数，这些参数描述了像素大小和左上角像素的地理位置。

接下来，函数计算指定的经纬度位置在DEM数据矩阵中的像素坐标。这是通过将经纬度位置与DEM数据的起始点进行比较并除以像素大小来实现的。

最后，函数读取DEM数据中相应像素的值，并返回这个值作为高度。

小结

这个函数提供了一种快速方便的方法，可以从DEM文件中获取指定经纬度位置的高度信息。

这对于遥感研究、应用来说非常有用，例如几何校正、大气校正、水体提取等。