Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Problem with elevation (not accurate at all) #12

Open
levydanqc opened this issue Jan 22, 2021 · 2 comments
Open

Problem with elevation (not accurate at all) #12

levydanqc opened this issue Jan 22, 2021 · 2 comments

Comments

@levydanqc
Copy link

Hi,
I know it's been some time since this project was publish, but I'm trying to read elevations for these coordinates : 46.75303,-71.40358 and it doesn't work at all.
Any coordinates for Quebec, Qc, CA, return values not even close.

Hoping to find a solution.

@nicholas-fong
Copy link

nicholas-fong commented Feb 4, 2021

The snippet below (tile N46W072.hgt) should work in all four quadrants.
Your POI is (46.75303, -71.40358) = 55 m
Parc naturel du mont Bélair (46.822845, -71.507354) = 483 m
Mont Roland-J.-Auger (46.845798, -71.472764) = 330 m.
on OpenStreetMap:
Parc naturel du mont Bélair = 484 m
Mont Roland-J.-Auger = 335 m.

import math
import numpy as np

def read_elevation(hgt_file, lon, lat):
    with open(hgt_file, 'rb') as hgt_data:
        # Each data point is a 16bit signed integer(i2) big endian(>)
        height = np.fromfile(hgt_data, np.dtype('>i2'), 1201*1201).reshape((1201, 1201))
        if ( lat >= 0.0 ):                    #first and second quadrants
            x = round (1200 - lat%1 * 1200) 
            y = round ( lon%1 * 1200 )
            return height[ x, y ].astype(np.int16)
        if ( lat < 0.0 ):                     #third and fourth quadrants
            x = round ((1 - lat%1) * 1200)
            y = round ( lon%1 * 1200 )
            return height[ x, y ].astype(np.int16)

def find_tile ( longitude, latitude ):
    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" ) 

@nicholas-fong
Copy link

nicholas-fong commented Feb 12, 2021

See this module which is versatile in reading tiles from different sources.

https://github.com/nicholas-fong/SRTM-GeoTIFF

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants