【代码】根据经纬度从数字高程模型（DEM）文件获取高度

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

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

应地理位置的高程。

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

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

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

image-20240325090706481

代码解析

以下是我们的Python函数：

def fromDemFileGetHeight(dem, lon, lat):
    try:
        DEMIDataSet = gdal.Open(dem)
    except Exception as e:
        print('Missing DEM file')
        pass
    DEMBand = 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文件中获取指定经纬度位置的高度信息。

这对于遥感应用来说非常有用，这段代码可以用在几何校正、大气校正、水体提取等应用。