如何阅读NASA .hgt二进制文件

Question

如果你对二进制文件有所了解,我确信这很简单,但我在这个分数上是新手.

我如何从NASA .hgt文件中提取数据？以下是来自www2.jpl.nasa.gov/srtm/faq.html的说明:

SRTM数据文件的名称类似于"N34W119.hgt".字母和数字是指什么,什么是".hgt"格式？

每个数据文件覆盖地球表面的一度经度块的一度纬度.前七个字符表示街区的西南角,N,S,E和W指北,南,东,西.因此,"N34W119.hgt"文件涵盖纬度34至35北纬和西经118-119(此文件包括加利福尼亚州洛杉矶市中心).文件扩展名".hgt"只代表单词"height",意思是提升.它不是格式类型.这些文件采用"原始"格式(无标题且未压缩),16位有符号整数,以海拔高度为单位测量的高程,在"地理"(纬度和经度阵列)投影中,数据空白由-32768表示.国际3弧秒文件具有1201列和1201行数据,总文件大小为2,884,802字节(= 1201 x 1201 x 2).美国1弧秒文件有3601列和3601行数据,总文件大小为25,934,402字节(= 3601 x 3601 x 2).有关更多信息,请阅读http://edcftp.cr.usgs.gov/pub/data/srtm/Readme.html上的文本文件"SRTM_Topo.txt".

谢谢你的帮助!我将在python脚本中使用这些数据,所以如果你不能使用任何其他语言的任何语言特定的技巧,那将是非常棒的.

Answer 1

经过测试的numpy示例:

import os
import math
import numpy

fn = 'DMV/N51E000.hgt'

siz = os.path.getsize(fn)
dim = int(math.sqrt(siz/2))

assert dim*dim*2 == siz, 'Invalid file size'

data = numpy.fromfile(fn, numpy.dtype('>i2'), dim*dim).reshape((dim, dim))

Answer 2

由于记录是固定长度(16位有符号整数)并且您知道网格大小(1201 x 1201或3601x3601),因此Python的struct模块似乎非常适合(未经测试的代码):

from struct import unpack,calcsize

# 'row_length' being 1201 or 3601 and 'row' being the raw data for one row
def read_row( row, row_length ):
    format = 'h'  # h stands for signed short

    for i in range(0, row_length):
        offset = i * calcsize(format)
        (height,) = unpack(format, row[offset : offset+calcsize(format))
        # do something with the height

用更通用的术语来描述它,基本上你想一次读取2个字节的文件,将读取的字节解析为16位有符号整数并进行处理.由于您已经知道网格大小,因此您可以逐行或以适合您的应用程序的任何其他方式读取它.它还意味着您可以随机搜索数据文件中的特定坐标.

Answer 3

如果您想要比从数百万次调用 struct.unpack 中获得的速度更快一点，请查看 array.array。虽然“struct-and-for-loop”实现在我公认的慢速笔记本电脑上需要几秒钟，但以下内容几乎是即时的：

from array import array

f = open(filename, 'rb')
format = 'h'
row_length = 1201
data = array(format)
data.fromfile(f, row_length*row_length)
data.byteswap()
f.close()