[GIS] Code to read elevation data from HGT file

javasrtm

I'm trying to read elevation data from NASA HGT file. I've started from Extracting elevation from .HGT file? and have changed this code a little to keep it simple and to use 1-ARCSEC (30-meter) data from http://dwtkns.com/srtm30m/ instead the original old data files.

After change the code, I can check it here, but the values are wrong.

I only changed these values but I'm not very sure what I'm doing:

public static final int HGT_RES = 1;  // <<- The new SRTM is 1-ARCSEC
public static final int HGT_ROW_LENGTH = 3601; <<-- New file resolution is 3601x3601

This is the code I'm using ( credits to don-vip from openstreetmap )

//  https://gis.stackexchange.com/questions/43743/extracting-elevation-from-hgt-file
private static final int SECONDS_PER_MINUTE = 60;

// alter these values for different SRTM resolutions
public static final int HGT_RES = 1; // resolution in arc seconds
public static final int HGT_ROW_LENGTH = 3601; // number of elevation values per line
public static final int HGT_VOID = -32768; // magic number which indicates 'void data' in HGT file  


public static void main(String[] args)  throws Exception {
    double coorLat = -22.352951806595854;
    double coorLon = -42.43121349392618;

    String filePath = "C:/Magno/DEMRJ/";
    String fileName = getHgtFileName(coorLat, coorLon);

    String htgFile = filePath + fileName;

    System.out.println( htgFile );

    //File f = new File( htgFile );

    ShortBuffer data = readHgtFile( htgFile );

    double fLat = frac( Math.abs(coorLat) ) * SECONDS_PER_MINUTE;
    double fLon = frac( Math.abs(coorLon) ) * SECONDS_PER_MINUTE;       

    System.out.println("fLat: " + fLat + "   fLon: " + fLon);

    int row = (int) Math.round(fLat * SECONDS_PER_MINUTE / HGT_RES);
    int col = (int) Math.round(fLon * SECONDS_PER_MINUTE / HGT_RES);        


    row = HGT_ROW_LENGTH - row;
    int cell = (HGT_ROW_LENGTH * (row - 1)) + col;

    System.out.println("Read SRTM elevation data from row/col/cell " + row + "," + col + ", " + cell + " == " + data.limit() );

    // valid position in buffer?
    if ( cell < data.limit() ) {

        short ele = data.get(cell);
        System.out.println("==> Read SRTM elevation data from row/col/cell " + row + "," + col + ", " + cell + " = " + ele);

        // check for data voids
        if (ele == HGT_VOID) {
            System.out.println("No Elevation - VOID");
        } else {
            System.out.println("Elevation : " + ele );
        }

    } else {
        System.out.println("No Elevation - Out of Range");
    }        


}


private static ShortBuffer readHgtFile(String file) throws Exception {

    FileChannel fc = null;
    ShortBuffer sb = null;
    try {
        // Eclipse complains here about resource leak on 'fc' - even with 'finally' clause???
        fc = new FileInputStream(file).getChannel();
        // choose the right endianness

        ByteBuffer bb = ByteBuffer.allocateDirect((int) fc.size());
        while (bb.remaining() > 0) fc.read(bb);

        bb.flip();
        //sb = bb.order(ByteOrder.LITTLE_ENDIAN).asShortBuffer();
        sb = bb.order(ByteOrder.BIG_ENDIAN).asShortBuffer();
    } finally {
        if (fc != null) fc.close();
    }

    return sb;
}   




public static String getHgtFileName(double latD, double lonD) {
    int lat = (int) latD;
    int lon = (int) lonD;
    String latPref = "N";
    if (lat < 0) latPref = "S";

    String lonPref = "E";
    if (lon < 0) {
        lonPref = "W";
    }
    String lonT = String.valueOf(lon).replace("-", "");
    if ( lonT.length() == 2 ) {
        lonT = "0" + lonT;
    }
    String ret = latPref + lat + lonPref + lonT + ".hgt";
    return ret.replace("-", "");
}   


public static double frac(double d) {
    long iPart;
    double fPart;
    iPart = (long) d;
    fPart = d - iPart;
    return fPart;
}

What I'm doing wrong?

The output is:

C:/Magno/DEMRJ/S22W042.hgt
fLat: 21.177108395751247   fLon: 25.872809635570775
Read SRTM elevation data from row/col/cell 2330,1552, 8388281 == 12967201
==> Read SRTM elevation data from row/col/cell 2330,1552, 8388281 = 129
Elevation : 129

But if you check the coordinates you can see the elevation value is 1222m or 4011ft ( not 129 – what unit? ).

Best Answer

I think I made some mistake in X,Y order. All seems fine now and I'm able to:

1) Open a HGT file.
2) Get the elevation value at a given Lat,Lon
3) Get the elevation profile from a path (red bars).

Data format

I'll take it as a little exercise in how to program a data reader. Have a look at the documentation:

SRTM data are distributed in two levels: SRTM1 (for the U.S. and its territories and possessions) with data sampled at one arc-second intervals in latitude and longitude, and SRTM3 (for the world) sampled at three arc-seconds.

Data are divided into one by one degree latitude and longitude tiles in "geographic" projection, which is to say a raster presentation with equal intervals of latitude and longitude in no projection at all but easy to manipulate and mosaic.

File names refer to the latitude and longitude of the lower left corner of the tile - e.g. N37W105 has its lower left corner at 37 degrees north latitude and 105 degrees west longitude. To be more exact, these coordinates refer to the geometric center of the lower left pixel, which in the case of SRTM3 data will be about 90 meters in extent.

Height files have the extension .HGT and are signed two byte integers. The bytes are in Motorola "big-endian" order with the most significant byte first, directly readable by systems such as Sun SPARC, Silicon Graphics and Macintosh computers using Power PC processors. DEC Alpha, most PCs and Macintosh computers built after 2006 use Intel ("little-endian") order so some byte-swapping may be necessary. Heights are in meters referenced to the WGS84/EGM96 geoid. Data voids are assigned the value -32768.

How to proceed

For your position, 50°24'58.888"N 14°55'11.377"E, you already found the correct tile, N50E14.hgt. Let's find out which pixel you are interested in. First latitude, 50°24'58.888"N:

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

arc seconds. Divided by three and rounded to the closest integer gives a grid row of 500. The same calculation for longitude results in grid column 1104.

The quickstart documentation lacks information about how rows and columns are organized in the file, but in the full documentation it is stated that

The data are stored in row major order (all the data for row 1, followed by all the data for row 2, etc.)

The first row in the file is very likely the northernmost one, i.e. if we are interested in row 500 from the lower edge, we actually have to look at row

1201 - 500 = 701

from the start if the file. Our grid cell is number

(1201 * 700) + 1104 = 841804

from the start of the file (i.e. skip 700 rows, and in the 701st one take sample 1104). Two bytes per sample means we have to skip the first 1683606 bytes in the file and then read two bytes in order to get our grid cell. The data is big-endian, which means that you need to swap the two bytes on e.g. Intel platforms.

Sample program

A simplistic Python program to retrieve the right data would look like this (see the docs for use of the struct module):

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)

Note that efficient data retrieval would have to look a little more sophisticated (e.g. not opening the file for each and every sample).

Alternatives

You could also use a program which can read the .hgt files out of the box. But that is boring.

QGIS Processing – How to Merge Several .hgt Files into One

(Perhaps grasshopper/elk also accepts IMG and GeoTiff formats, anyway...)

Not a QGIS answer; my choice is SAGA, currently version 7.

Install and Start SAGA

Mosaicking

Drag and drop your hgt files to SAGA window.
In the Tools window, go to Grid > Tools and activate Mosaicking tool.
Mosaicking window has an entry for input files, as >>Grids. Click [...] button and select all your grids. [Edit] Choose Nearest Neighbour for resampling method. Then click on [Apply] button and then [Execute]. You'll have newly merged raster (grid).

Export hgt

Then export the new grid as hgt. In the Tools window, go to Import/Export > GDAL/OGR and activate Export Raster tool.
[Edit] Choose newly created Grid System (usualy 2401x, 2401y for 2x2 mosaic), then select new grid in the >>Grid(s). Give a new filename in File (file extension should be .hgt. In the Format selection, choose ESRI.hdr Labelled as the output format. Then click on [Apply] button and then [Execute] button. Now you have merged hgt file.

In this export process, SAGA is using a GDAL driver for ESRI BIL format. This is sometimes tricky (depends on data), but would work in most cases.

Best Answer

Related Solutions

Extract Elevation from .HGT File – Elevation Data Guide

Data format

How to proceed

Sample program

Alternatives

QGIS Processing – How to Merge Several .hgt Files into One

Related Question