2

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

Vince
  • 20,017
  • 15
  • 45
  • 64
Magno C
  • 2,140
  • 3
  • 24
  • 63
  • This code is incorrect for values below sea level (DTED uses ones-complement encoding). – Vince Apr 25 '18 at 21:59
  • I don't understood very well what you said. Can you tell noobie-friendly? – Magno C Apr 25 '18 at 22:15
  • Consider to check with this working implementation (Python, though, but perhaps good to compare the formula): https://trac.osgeo.org/grass/browser/grass-addons/grass7/raster/r.in.srtm.region/r.in.srtm.region.py (manual: https://grass.osgeo.org/grass74/manuals/r.in.srtm.html) – markusN Jun 30 '18 at 10:03
  • If it helps, here is the Java source code to a tool that reads SRTM data from an hgt file and has been used and tested for many years: https://github.com/jblindsay/whitebox-geospatial-analysis-tools/blob/master/ImportExport/plugins/ImportSRTM.java – WhiteboxDev Oct 16 '18 at 19:32
  • @WhiteboxDev I think it is too coupled for me.... – Magno C Oct 17 '18 at 18:18

1 Answers1

1

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

enter image description here

Magno C
  • 2,140
  • 3
  • 24
  • 63
  • I thik the algoritm to detect the SRTM file is incorrect. The file name refers to the lower-left corner. So 23.0001 SOUTH is located in tile S24xxxx. – Magno C Apr 26 '18 at 10:30