1

From How to extract elevation from .HGT file? I know how to get the elevation from SRTM1 and SRTM3 files. But how to get elevation from SRTM30 Plus ftp://topex.ucsd.edu/pub/srtm30_plus/srtm30/data/ ?
I tried to use the following code, but it returns the incorrect values:

def read_elevation_from_file_srtm30_plus(hgt_file, lon, lat):
with open('e060n40.Bathymetry.srtm', 'rb') as hgt_data:
    # HGT is 16bit signed integer(i2) - big endian(>)
    elevations = np.fromfile(hgt_data, np.dtype('>i2'), 6000*4800)\
                            .reshape((6000, 4800))

    lat_y = int(round((lat - int(lat)) * 6000, 0))
    lon_x = int(round((lon - int(lon)) * 4800, 0))

    return elevations[6000 - lat_y, lon_x].astype(int)

What am I doing wrong?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
LuFang
  • 21
  • 3
  • What is the code supposed to do exactly? Which value do you expect and which do you get? – bugmenot123 Jun 21 '17 at 12:21
  • @bugmenot123 For example, if I execute the code print(get_elevation_srtm30_plus(98.06, 23.06)), then I get -4056 instead of the correct 853 (get this value with the help of Landserf program). – LuFang Jun 21 '17 at 12:37
  • Well, explain step by step what your script does and why. What are those constants, did you make sure that x/y -> lon/lat mapping is correct? etc. – bugmenot123 Jun 21 '17 at 14:45
  • 1
    SRTM defines data as 16bit signed integer in Big-Endian and has NrOfCellsPerLine = 4800 and NrOfLines = 6000 (4800*6000 elevations). The script opens 'e060n40.Bathymetry.srtm' file and converts it to 2d array. Next, the script converts a given latitude-longitude to data position in array. That's all. And yes, I'm not sure that I'm doing it right. – LuFang Jun 22 '17 at 05:23
  • Double check your lat_y and lon_x calculations, shouldn't the data's minimum lat and lon play a part? – bugmenot123 Jun 22 '17 at 11:36
  • Well, I changed some code to get the correct array index. But that did not solve my problem...
        lat_y = int(round((lat - int(lat)) * 5999, 0))
    lon_x = int(round((lon - int(lon)) * 4799, 0))

    return elevations[5999 - lat_y, lon_x].astype(int)</pre>
    – LuFang Jun 22 '17 at 12:43

0 Answers0