如何在DTED中查询感兴趣矩形区域内的最高海拔？

一、使用Python的dted库实现 #

dted库是专门处理DTED格式的轻量库，适合单文件小范围处理，直接按以下步骤操作：

1. 加载目标DTED文件，获取其地理覆盖范围和高程数据矩阵
2. 根据你的矩形ROI经纬度，计算出在DTED数据中的对应行列范围
3. 提取该行列范围内的所有高程值，取最大值

代码示例 #

import dted

# 替换为你的DTED文件路径
dted_path = "your_file.dt2"
dted_data = dted.DTED(dted_path)

# 定义你的矩形ROI经纬度（示例为西南角和东北角坐标）
roi_lon_min = 116.0
roi_lon_max = 116.1
roi_lat_min = 39.0
roi_lat_max = 39.1

# 验证ROI是否在当前DTED文件覆盖范围内
if not (dted_data.lon_min <= roi_lon_min and dted_data.lon_max >= roi_lon_max and
        dted_data.lat_min <= roi_lat_min and dted_data.lat_max >= roi_lat_max):
    raise ValueError("ROI超出当前DTED文件的覆盖范围")

# DTED纬度从北到南递减，北边界对应更小的行号
row_start = dted_data.lat_to_index(roi_lat_max)
row_end = dted_data.lat_to_index(roi_lat_min)
col_start = dted_data.lon_to_index(roi_lon_min)
col_end = dted_data.lon_to_index(roi_lon_max)

# 提取ROI内的高程数据并计算最大值
roi_elevations = dted_data.elevations[row_start:row_end+1, col_start:col_end+1]
max_elevation = roi_elevations.max()

print(f"ROI内最高海拔：{max_elevation} 米")

二、使用GDAL实现 #

GDAL适合处理多DTED文件拼接、大范围区域的场景，兼容性更强，操作步骤如下：

1. 用GDAL打开DTED数据集（支持单文件或虚拟数据集VRT）
2. 通过地理变换参数，将经纬度坐标转换为像素行列号
3. 读取对应像素区域的高程数据，排除NoData值后取最大值

代码示例 #

from osgeo import gdal
import numpy as np

# 替换为你的DTED文件路径
dted_path = "your_file.dt2"
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
dataset = gdal.Open(dted_path, gdal.GA_ReadOnly)

if not dataset:
    raise FileNotFoundError(f"无法打开DTED文件：{dted_path}")

# 获取地理变换参数：(左上角经度, 经度分辨率, 0, 左上角纬度, 0, 纬度分辨率)
geotransform = dataset.GetGeoTransform()
lon_origin, lon_res, _, lat_origin, _, lat_res = geotransform

# 获取高程波段和NoData值
band = dataset.GetRasterBand(1)
nodata_value = band.GetNoDataValue()

# 定义你的矩形ROI经纬度
roi_lon_min = 116.0
roi_lon_max = 116.1
roi_lat_min = 39.0
roi_lat_max = 39.1

# 计算对应像素行列号
col_start = int((roi_lon_min - lon_origin) / lon_res)
col_end = int((roi_lon_max - lon_origin) / lon_res)
# 纬度分辨率为负数，需注意计算方向
row_start = int((roi_lat_max - lat_origin) / lat_res)
row_end = int((roi_lat_min - lat_origin) / lat_res)

# 确保行列号在数据集范围内
col_start = max(0, col_start)
col_end = min(dataset.RasterXSize - 1, col_end)
row_start = max(0, row_start)
row_end = min(dataset.RasterYSize - 1, row_end)

# 读取ROI区域高程数据
width = col_end - col_start + 1
height = row_end - row_start + 1
elev_data = band.ReadAsArray(col_start, row_start, width, height)

# 排除NoData值后计算最大值
valid_data = elev_data[elev_data != nodata_value]
max_elevation = valid_data.max() if valid_data.size > 0 else np.nan

print(f"ROI内最高海拔：{max_elevation} 米")

# 关闭数据集
dataset = None

注意事项 #

• DTED文件的纬度分辨率为负数，计算行号时要注意方向，避免行列颠倒
• 如果ROI跨多个DTED文件，GDAL可以通过创建VRT虚拟数据集合并后再处理，无需手动拼接
• 若遇到NoData值，需先过滤再计算最大值，避免结果异常

内容的提问来源于stack exchange，提问作者carraro