Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

how to read NASA .hgt binary files

Tags:

file

binary

A tested numpy example:

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))

Since the records are fixed length (16-bit signed integers) and you know the grid size (1201 x 1201 or 3601x3601), Python's struct module seems ideally suited (untested code):

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

Describing it in more generic terms, basically you want to read the file in 2 bytes at a time, parse the bytes read as a 16-bit signed integer and process it. Since you know the grid size already you can read it in row by row or in any other manner that is convenient to your application. It also means you can randomly seek to specific coordinates inside the data file.


If you want a little more speed than you get from millions of calls to struct.unpack, have a look at array.array. While the "struct-and-for-loop" implementation takes several seconds on my admittedly slow laptop, the following is near instantaneous:

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()

https://gdal.org/drivers/raster/srtmhgt.html

    Input_HGT = 'N30E120.hgt'
    import gdal
    Raster = gdal.Open(Input_HGT) 

All functions available with GDAL on raster files could be applied on this 'Raster' like Functions available with the variable, 'Raster'