HGT 文件解析：从二进制数据到高程网格的实现方案

原创于 2025-09-25 11:24:27 发布 · 465 阅读

9 ·

CC 4.0 BY-SA版权

文章标签：

#python #高程数据 #hgt #GIS

python操作专栏收录该内容

4 篇文章

订阅专栏

一、问题起源

在地理信息系统 (GIS) 开发和地形分析中，经常需要处理高程数据。HGT 格式是美国地质调查局 (USGS) 提供的一种常用高程数据格式，广泛应用于各种地形分析场景。然而，HGT 文件采用二进制存储，且二进制文件内部不带有坐标信息，无法直接读取和使用，这给许多开发者带来了困扰。本文将详细解析 HGT 文件的结构，并提供一个完整的 Python 函数，将 HGT 文件解析为可直接用于可视化和分析的高程网格数据。

二、解决路径

1、HGT 文件格式分析

HGT 文件命名遵循特定规则，例如 "N36E117.hgt" 表示该文件包含北纬 36-37 度、东经 117-118 度范围内的高程数据。文件采用二进制存储，每个高程值使用 2 字节的有符号整数表示。HGT 文件有多种分辨率：1 角秒分辨率：1201×1201 像素，3 角秒分辨率：401×401 像素，1/3 角秒分辨率：3601×3601 或 3600×3600 像素。

三、解析函数实现

下面是解析 HGT 文件的完整函数实现：

import struct
import numpy as np

def parse_hgt_to_grid(filename):
    # 从文件名提取经纬度范围（N36E117对应北纬36-37度，东经117-118度）
    lat_str = filename[1:3]
    lon_str = filename[4:7]
    start_lat = int(lat_str)  # 36°N
    start_lon = int(lon_str)  # 117°E

    # 读取文件
    try:
        with open(filename, 'rb') as f:
            data = f.read()
    except FileNotFoundError:
        raise FileNotFoundError(f"找不到文件: {filename}")
    
    file_size = len(data)
    
    # 确定分辨率和像素数量
    if file_size == 2 * 1201 * 1201:
        pixels = 1201  # 1角秒分辨率
        step = 1 / 3600
    elif file_size == 2 * 401 * 401:
        pixels = 401   # 3角秒分辨率
        step = 3 / 3600
    elif file_size in [2 * 3601 * 3601, 2 * 3600 * 3600]:
        pixels = 3601 if file_size == 2 * 3601 * 3601 else 3600  # 1/3角秒分辨率
        step = 1 / (3 * 3600)
    else:
        raise ValueError(f"不支持的HGT文件大小：{file_size}字节")

    # 解析高程数据
    fmt = f">{pixels * pixels}h"
    elevations_int = np.array(struct.unpack(fmt, data), dtype=np.int16)
    elevations = elevations_int.astype(np.float32)
    elevations[elevations_int == -32768] = np.nan  # 将无效值设为NaN
    elevation_grid = elevations.reshape((pixels, pixels))

    # 计算完整经纬度范围（1×1度）
    lat_range = np.linspace(start_lat + 1, start_lat, pixels)  
    lon_range = np.linspace(start_lon, start_lon + 1, pixels)  

    return elevation_grid, lat_range, lon_range, step

四、函数解析

1、文件名解析：

从文件名中提取起始经纬度信息，确定该 HGT 文件覆盖的地理范围。

lat_str = filename[1:3]
lon_str = filename[4:7]
start_lat = int(lat_str)  # 36°N
start_lon = int(lon_str)  # 117°E

2、文件读取：

以二进制模式读取文件内容，并处理文件不存在的异常。

try:
    with open(filename, 'rb') as f:
        data = f.read()
except FileNotFoundError:
    raise FileNotFoundError(f"找不到文件: {filename}")

3、分辨率判断：

根据文件大小判断 HGT 文件的分辨率，不同分辨率对应不同的像素数量和地理步长。

file_size = len(data)

if file_size == 2 * 1201 * 1201:
    pixels = 1201  # 1角秒分辨率
    step = 1 / 3600
# 其他分辨率判断...

4、数据解析：

使用 struct 模块解析二进制数据，将其转换为 numpy 数组，并将无效值 (-32768) 转换为 NaN，最后重塑为二维高程网格。

fmt = f">{pixels * pixels}h"
elevations_int = np.array(struct.unpack(fmt, data), dtype=np.int16)
elevations = elevations_int.astype(np.float32)
elevations[elevations_int == -32768] = np.nan
elevation_grid = elevations.reshape((pixels, pixels))

5、经纬度范围计算：

lat_range = np.linspace(start_lat + 1, start_lat, pixels)  
lon_range = np.linspace(start_lon, start_lon + 1, pixels)

五、总结

上述通过一个完整的 Python 函数实现了从二进制 HGT 文件到高程网格数据的转换。该函数能够自动识别不同分辨率的 HGT 文件，并处理无效值，最终返回可直接用于可视化和分析的高程数据及对应的经纬度范围。使用该函数解析 HGT 文件后，可以方便地进行地形可视化、坡度分析、等高线绘制等后续操作。例如，结合 Matplotlib 可以快速绘制高程图，并在地图上标记特定点位，为地理信息分析提供有力支持。在实际应用中，还可以根据需要扩展该函数，例如添加数据裁剪、重采样等功能，以满足不同场景的需求。