气象编程 | 使用Python处理SRTM(.hgt)文件
引言最近在做美赛时，使用了高精度的地形文件（海拔高度），因此在网站上下载了高精度的.hgt文件，并学习了处理。于是记录下来，有了这一篇笔记。
数据类型
航天飞机雷达地形任务（SRTM），顾名思义，是一个研究任务，产生一个通用的免费数字高程模型。SRTM以.hgt为结尾，文件的名字解释了hgt文件的范围。比如，文件名字为N30E11,表示范围时30°N11°E到31°N12°E的正方形网格范围。
分为两种类型,SRTM1,SRTM3,分别对应的网格(1201*1201或3601*3601)。
资料格式
可以简单的理解为hgt文件将每一个经纬度分为了一个1201*1201(3601*3601)的网格，位置与范围如数据类型所示。
因此使用Numpy打开hgt文件的代码下面所示：
import numpy as np
SAMPLES = 1201# Change this to 3601 for SRTM1
def read_hgt(f_name, lat, lon):
with open(f_name, 'rb') as hgt_data:
elevations = np.fromfile(hgt_data, np.dtype('>i2'), SAMPLES * SAMPLES) \
.reshape((SAMPLES, SAMPLES))
lat_range = np.arange(lat, lat + 1 / 1200 + 1, 1 / 1200)
lon_range = np.arange(lon, lon + 1 + 1 / 1200, 1 / 1200)
return lat_range, lon_range, elevations
上面返回经纬网格及对应的海拔高度。
复杂应用：多个.hgt文件绘制地形图
我们在某地形网站下载到了澳大利亚维多利亚洲的地形图，使用循环对多个文件进行处理。代码如下：
import numpy as np
import matplotlib.pyplot as plt
import os
SAMPLES = 1201# Change this to 3601 for SRTM1
def read_return(lat, lon):
f_name = r'./I55/S{}E{}.hgt'.format(lat, lon)
with open(f_name, 'rb') as hgt_data:
elevations = np.fromfile(hgt_data, np.dtype('>i2'), SAMPLES * SAMPLES) \
.reshape((SAMPLES, SAMPLES))
lat_range = np.arange(-lat, -lat + 1 / 1200 + 1, 1 / 1200)[::-1] # 南半球纬度为负数
lon_range = np.arange(lon, lon + 1 + 1 / 1200, 1 / 1200)
return lat_range, lon_range, elevations
def plot_d(ax):
for i in range(33, 41):
for j in range(139, 153):
if os.path.exists(r'./I55/S{}E{}.hgt'.format(i, j)):
lat, lon, ele = read_return(i, j)
lim = np.arange(1, 2401, 200)
ax.contourf(lon, lat, ele, lim, cmap='binary')
else:
print(i, j)
return ax
if __name__ == '__main__':
fig ,ax = plt.subplots()
ax = plot_d(ax)
plt.show()
得到的结果如下图：
加上一些地理信息后得到：
这只是一个简单的应用，还可以使用插值、采样等方法获得别的网格的地形图。
本文作者孟子路，更多请邮件联系： mzll1202@163.com
文章来源于微信公众号：气象学家
页:
[1]