Handle Negative Latitude and Longitude Correctly

[JohnBolander]

srtm-python/srtm.py

Line 36 in ea07b6d

return elevations[SAMPLES - 1 - lat_row, lon_row].astype(int)

To handle negative values, shouldn't you be taking the absolute value of lon_row and lat_row prior to passing them into elevations[SAMPLES - 1 - lat_row, lon_row].astype(int)? Otherwise I think lon_row will be off and lat_row will not work. For example: latitude -62.1935316, longitude -58.9801776 yields lat_row:-232,lon_row:-1176. lat_row of -232 will break since you are looking up elevations[1432, -1176]. Am I missing something?

[JohnBolander]

After further review, to handle negative latitude and longitude correctly you need to get correct file. get_file_name should look the following:
def get_file_name(self,lon, lat):
"""
Returns filename such as N27E086.hgt, concatenated
with HGTDIR where these 'hgt' files are kept
"""
if lat >= 0:
ns = 'N'
lat_offset = 0
elif lat < 0:
ns = 'S'
lat_offset = 1

    if lon >= 0:
        ew = 'E'
        lon_offset = 0
    elif lon < 0:
        ew = 'W'
        lon_offset = 1

    hgt_file = "%(ns)s%(lat)02d%(ew)s%(lon)03d.hgt" % {'lat': abs(lat) + lat_offset, 
                 'lon': abs(lon) + lon_offset, 'ns': ns, 'ew': ew}

and you'll then need to adjust the longitude to account for the fact that the file you referenced has a corner that is more negative
if lat_row < 0:
lat_row = abs(lat_row)
else:
lat_row = self.samples - 1 - lat_row
return self.elevations[lat_row, lon_row].astype(int)

[aatishnn]

@JohnBolander Thanks. This is indeed a bug. I was working on this yesterday and this is the test I wrote to accommodate all possible latititude, longitude scenarios. If you have already fixed these, maybe you can try testing it with this code:

import re
import unittest
import subprocess
from srtm import get_elevation, get_file_name


'''
These elevation values were taken using gdallocationinfo command which is
a part of gdal-bin package. You can install it in Ubuntu (or derivatives)
using:
    sudo apt-get install gdal-bin
'''

TEST_DATA = [
    {
        'name': 'Mt. Everest (NE)',
        'lat': 27.988056,
        'lon': 86.925278,
        'filename': 'hgt/N27E086.hgt',
        # gdallocationinfo N27E086.hgt -wgs84 86.925278 27.988056
    },
    {
        'name': 'Mt. Kanchanjunga (NE)',
        'lat': 27.7025,
        'lon': 88.146667,
        'filename': 'hgt/N27E088.hgt',
        # gdallocationinfo N27E088.hgt -wgs84 88.146667 27.7025
    },
    {
        'name': 'A point in NW corner',
        'lat': 52.25,
        'lon': -125.22,
        'filename': 'hgt/N52W126.hgt',
        # gdallocationinfo hgt/N52W126.hgt -wgs84 -125.22 52.25
    },
    {
        'name': 'Another point in NW corner',
        'lat': 52.68,
        'lon': -125.67,
        'filename': 'hgt/N52W126.hgt',
        # gdallocationinfo hgt/N52W126.hgt -wgs84 -125.22 52.25
    },
    {
        'name': 'A point in SE corner',
        'lat': -27.25,
        'lon': 133.22,
        'filename': 'hgt/S28E133.hgt',
    },
    {
        'name': 'Another point in SE corner',
        'lat': -27.65,
        'lon': 133.69,
        'filename': 'hgt/S28E133.hgt',
    },
    {
        'name': 'Extreme point in SE corner',
        'lat': -28,
        'lon': 133,
        'filename': 'hgt/S28E133.hgt',
    },
    {
        'name': 'A point in SW corner',
        'lat': -50.12,
        'lon': -74.65,
        'filename': 'hgt/S51W075.hgt',
    },

]


def get_elevation_from_gdallocationinfo(filename, lat, lon):
    output = subprocess.check_output([
        'gdallocationinfo', filename, '-wgs84', str(lon), str(lat)
    ])
    return int(re.search('Value: (\d+)', str(output)).group(1))


class TestSRTMMethods(unittest.TestCase):

    def test_get_file_name(self):
        for mountain in TEST_DATA:
            filename = get_file_name(mountain['lat'], mountain['lon'])
            self.assertEqual(filename, mountain['filename'])


def create_get_elevation_test(mountain):
    def do_test(self):
        elevation = get_elevation(mountain['lat'], mountain['lon'])
        gdal_elevation = get_elevation_from_gdallocationinfo(
            mountain['filename'], mountain['lat'], mountain['lon']
        )
        self.assertEqual(elevation, gdal_elevation)
    return do_test


for mountain in TEST_DATA:
    test_method = create_get_elevation_test(mountain)
    test_method.__name__ = 'test_elevation_for_%s(%d, %d)' \
                % (mountain['name'], mountain['lat'], mountain['lon'])
    setattr(TestSRTMMethods, test_method.__name__, test_method)


if __name__ == '__main__':
    unittest.main()

I will probably work on this next week.

Thanks!

[andrewcooke]

you can use floor() rather than int() in a bunch of places to get things to work. it gives simpler code.

thanks for this. i tried to comment on the website with the article but couldnt get disqus to work, so came here and found this issue.

[jayveesea]

was this bug ever fixed? i'm currently getting errors with south america data. i'll try messing with what @JohnBolander has above.

[andrewcooke]

if it's any help https://github.com/andrewcooke/choochoo/tree/master/py/ch2/srtm includes code that works, with bilinear and spline interpolation.

[jayveesea]

if i'm reading the above correctly i've got the following revisions:

def read_elevation_from_file(hgt_file, lat, lon):
    with open(hgt_file, 'rb') as hgt_data:
        # HGT is 16bit signed integer(i2) - big endian(>)
        elevations = np.fromfile(
            hgt_data,  # binary data
            np.dtype('>i2'),  # data type
            SAMPLES * SAMPLES  # length
        ).reshape((SAMPLES, SAMPLES))

        lat_row = int(round((lat - int(lat)) * (SAMPLES - 1), 0))
        lon_row = int(round((lon - int(lon)) * (SAMPLES - 1), 0))

        if lat_row < 0:
            lat_row = abs(lat_row)
        else:
            lat_row = SAMPLES - 1 - lat_row
        return elevations[lat_row, lon_row].astype(int)

and

def get_file_name(lon, lat):
    """
    Returns filename such as N27E086.hgt, concatenated
    with HGTDIR where these 'hgt' files are kept
    """
    if lat >= 0:
        ns = 'N'
        lat_offset = 0
    elif lat < 0:
        ns = 'S'
        lat_offset = 1

    if lon >= 0:
        ew = 'E'
        lon_offset = 0
    elif lon < 0:
        ew = 'W'
        lon_offset = 1

    hgt_file = "%(ns)s%(lat)02d%(ew)s%(lon)03d.hgt" % {'lat': abs(lat) + lat_offset, 
                 'lon': abs(lon) + lon_offset, 'ns': ns, 'ew': ew}
    hgt_file_path = os.path.join(HGTDIR, hgt_file)
    print(hgt_file)
    if os.path.isfile(hgt_file_path):
        return hgt_file_path
    else:
        return None

...it looks like get_file_name returns the correct file but read_elevation_from_file does not return the correct elevation. With S33W071.hgt and using get_elevation(-70.011111, -32.653611) it returns 3117, but elevation should be closer to 6961.

[jayveesea]

if it's any help https://github.com/andrewcooke/choochoo/tree/master/py/ch2/srtm includes code that works, with bilinear and spline interpolation.

thanks! this might help if I get stuck.

[JohnBolander]

was this bug ever fixed? i'm currently getting errors with south america data. i'll try messing with what @JohnBolander has above.

Sorry no, I don't think it ever got fixed in this repository. This was a long time ago. I no longer have access to the solution that I wrote. There might have been another flip I had to do for the Southern Hemisphere.

[nicholas-fong]

This snippet solves the negative latitude and negative longitude problems. It works in all four quadrants. Referring to above lat=-32.653611, lon=-70.011111 (S33W071.hgt), it returns an elevation of 6928. The core of the snippet uses gdal's GeoTransform to deal with the translation of lat/lon to pixel. It can handle other file formats and tile sizes, providing a migration path to modern file formats. https://github.com/nicholas-fong/SRTM-GeoTIFF

Save this file as issues10.py

from osgeo import gdal
import math

def  read(my_file, lat, lon):
    data = gdal.Open(my_file)
    band1 = data.GetRasterBand(1)
    GT = data.GetGeoTransform()
    # call gdal's Affine Transformation (GetGeoTransform method)
    # GetGeoTransform translates latitude, longitude to pixel indices
    x_pixel_size = GT[1]
    y_pixel_size = GT[5]
    xP = int((lon - GT[0]) / x_pixel_size )
    yL = int((lat - GT[3]) / y_pixel_size )
    return ( int( band1.ReadAsArray(xP,yL,1,1) ) )

def   find( latitude, longitude ):
    if  ( latitude >= 0.0 and longitude >= 0.0 ):
        hemi, meri = "N", "E"
        t1 = f"{math.floor(latitude):02d}"
        t2 = f"{math.floor(longitude):03d}"
    elif ( latitude >= 0.0 and longitude < 0.0 ):
        hemi, meri = "N", "W"
        t1 = f"{math.floor(latitude):02d}"
        t2 = f"{math.ceil(abs(longitude)):03d}"
    elif ( latitude < 0.0 and longitude < 0.0 ):
        hemi, meri = "S", "W"
        t1 = f"{math.ceil(abs(latitude)):02d}"
        t2 = f"{math.ceil(abs(longitude)):03d}"
    elif ( latitude < 0.0 and longitude >= 0.0 ):
        hemi, meri = "S", "E"
        t1 = f"{math.ceil(abs(latitude)):02d}"
        t2 = f"{math.floor(longitude):03d}"
    return( f"{hemi}{t1}{meri}{t2}.hgt" ) 

Test snippet with -32.653611, -70.01111 (highest mountain in South America)

>>> import issues10
>>> issues10.find(-32.653611, -70.011111)
'S33W071.hgt'
>>> issues10.read( 'S33W071.hgt', -32.653611, -70.011111 )
6928
>>>

[TFehlh]

I'm having problems with one of the test points: lat 27.988056, lon: 27.988056.

I think it's a problem with the hgt file, but can someone else confirm? What are YOU getting?

I downloaded the hgt file from https://dds.cr.usgs.gov/srtm/version2_1/SRTM3. Is there another source? I tried the version1 data with the same result.

Here's the gdallocationinfo that shows the same answer that I am getting.

gdallocationinfo N27E086.hgt -wgs84 86.925278 27.988056
Report:
Location: (1110P,14L)
Band 1:
Value: -32768

[LaurentBerder]

I've been using the code as is without a worry for some months now.
However, I'm now stumbling on a set of particular cases with points just west of the Greenwich meridian (with -1 <= longitude <= 0).

We can take as an example the point {45.837883, -0.882169} in Western France.

For these points, the function get_file_name returns a name with W000 (N45W000.hgt for our point in example), which does not exist.
The elevations for these points are actually located in the E000 files (N45E000.hgt for the example point).

So I corrected the code as such:

def get_file_name(lat, lon):
    """
    Returns filename such as N27E086.hgt, concatenated
    with HGTDIR where these 'hgt' files are kept
    """
    if pd.isna(lat) or pd.isna(lon):
        return None

    if lat >= 0:
        ns = 'N'
    elif lat < 0:
        ns = 'S'

    # The East/West distinction is not at 0, but at -1
    if lon > -1:
        ew = 'E'
    elif lon <= -1:
        ew = 'W'

    hgt_file = "%(ns)s%(lat)02d%(ew)s%(lon)03d.hgt" % \
               {'lat': abs(lat), 'lon': abs(lon), 'ns': ns, 'ew': ew}
    hgt_file_path = os.path.join(HGTDIR, hgt_file)
    if os.path.isfile(hgt_file_path):
        return hgt_file_path
    else:
        return None

A similar change in the latitude would probably be necessary for points close to the equator, but I haven't tested it as I don't have any data concerned by that issue.

[Pinnh]

I solved this using math.floor() method ,codes available at here