从.HGT文件中提取高度？

#1 楼

数据格式
我将作为一个小程序来练习如何对数据读取器进行编程。请看一下文档：

SRTM数据分为两个级别：SRTM1（用于美国及其领土和财产），其数据在纬度和经度上以一弧秒的间隔采样，而SRTM3 （对于世界）以3弧秒为单位进行采样。
在“地理”投影中，数据按纬度和经度图块一分为一，也就是说，在纬度和经度之间没有等距间隔投影，但易于操作和镶嵌。
文件名指的是瓷砖左下角的经度和纬度-例如N37W105的左下角位于北纬37度和西经105度。更确切地说，这些坐标是指左下像素的几何中心，在SRTM3数据的情况下，其范围约为90米。 。
字节以摩托罗拉的“ big-endian”顺序排列，最高有效字节在前，
可由使用Power PC处理器的Sun SPARC，Silicon Graphics和Macintosh计算机等系统直接读取。 DEC Alpha，大多数2006年以后生产的PC和Macintosh计算机都使用Intel（“小端”）顺序，因此可能需要一些字节交换。高度以米为单位，参考WGS84 / EGM96大地水准面。为数据空赋予的值是-32768。 N50E14.hgt。让我们找出您感兴趣的像素。第一纬度，50°24'58.888“ N：

24'58.888" = (24 * 60)" + 58.888" = 1498.888"

快速入门文档缺少有关文件中行和列的组织方式的信息，但是在完整文档中，声明了

数据以行主要顺序存储（第1行的所有数据，后按行2的所有数据等。）

文件中的第一行很可能是最北端的，即，如果我们对下边缘的第500行感兴趣，我们实际上必须查看从文件开始处在行

1201 - 500 = 701

。我们的网格单元格是从文件开头开始的编号

(1201 * 700) + 1104 = 841804

（即跳过700行，在第701行中取样1104）。每个样本两个字节意味着我们必须跳过文件中的前1683606个字节，然后读取两个字节才能获取网格单元。数据是big-endian，这意味着您需要交换两个字节，例如英特尔平台。
示例程序
一个简单的Python程序可以检索正确的数据，如下所示（请参见有关struct模块使用的文档）：

import struct

def get_sample(filename, n, e):
    i = 1201 - int(round(n / 3, 0))
    j = int(round(e / 3, 0))
    with open(filename, "rb") as f:
        f.seek(((i - 1) * 1201 + (j - 1)) * 2)  # go to the right spot,
        buf = f.read(2)  # read two bytes and convert them:
        val = struct.unpack('>h', buf)  # ">h" is a signed two byte integer
        if not val == -32768:  # the not-a-valid-sample value
            return val
        else:
            return None

if __name__ == "__main__":
    n = 24 * 60 + 58.888
    e = 55 * 60 + 11.377
    tile = "N50E14.hgt"  # Or some magic to figure it out from position
    print get_sample(tile, n, e)

请注意，高效数据检索将看起来有些复杂（例如，不为每个样本打开文件）。但这很无聊。

#2 楼

GDAL可以使用SRTMHGT驱动程序读取/写入这些栅格格式。这意味着您可以使用QGIS，ArcGIS查看栅格，或使用GDAL实用程序（如gdallocationinfo）从某个点获取值，例如：

将DMS转换为DD：


纬度：50°24'58.888“ N = 50 +（24/60）+（58.888 / 3600）= 50.4163577778
纬度：14°55'11.377” E = 14 +（55/60）+（ 11.377 / 3600）= 14.9198269444

然后从壳中使用gdallocationinfo file.hgt -wgs84 long lat：

$ gdallocationinfo N50E14.hgt -wgs84 14.9198269444 50.4163577778
Report:
  Location: (1104P,700L)
  Band 1:
    Value: 216

海拔为216 m。

#3 楼

如果您确实使用QGIS，请检查是否安装了python插件“ Point Sampling Tool”。您可以在->增强（Python）->分析中找到它。选择一个新的点形状进行输出。

就这样：-)

  Chris

#4 楼

克里斯（Chris）的回答表明，从QGIS中的图层采样点很简单。查看SRTM文档中的Quickstart PDF。它说明了高程数据的存储方式。总结一下：


SRTM3文件包含一系列的大端整数值。
每个整数值都代表“以WGS84 / EGM96大地水准面为参考的米”的高程。用于表示无数据像素的.hgt的值。
有1201行的1201行样本，因此总共应该有1442401个整数值。坐标和像素，因此获取高程只需从文件中的适当偏移量读取整数值即可。给定相对于场景左上角的像素坐标-32768和x，基本上是y。像素offset = (y * 1201) + x是文件中的第一个整数，像素0,0是文件中的最后一个整数。