Hi guys !
I have some problems to get the elevation, based on the latitude and longitude of one point, from an ASCII DEM file.
The .asc file looks like this :
ncols 1000
nrows 1000
xllcorner 1213999.500000000000
yllcorner 6153000.500000000000
cellsize 1.000000000000
NODATA_value -99999.00
486.41 486.57 486.47 486.27 486.17 485.98 485.74 485.66
485.57 485.51 485.43 485.34 485.18 485.06 485.07 485.05
484.85 484.84 485.04 485.37 485.73 485.98 486.15 486.51
486.85 487.12 487.53 487.93 488.19 488.56 488.89 489.23
489.65 490.02 490.37 490.72 491.05 491.31 491.69 492.07
492.41 492.82 493.31 493.76 494.10 494.51 494.90 495.27
495.45 495.83 496.35 496.72 497.06 497.17 497.26 497.42
497.68 497.92 498.04 498.15 498.32 498.56 498.77 499.02
First I import the file in a numpy array, having 1000 rows and 1000 columns.
Based on the latitude and longitude of one point I calculate the index on x and y to find the altitude in the numpy array.
The code to do this is as follow :
#! usr/bin/env python3
#! -*- coding : utf-8 -*-
import numpy as np
Tile = 'T5821_1213999.5_6153000.5.asc'
rep = 'PathToFile' + Tile
par = open(rep)
dat = np.loadtxt(par, skiprows=6)
nc, nr = 1000, 1000
CellsizeX, CellsizeY = 1, 1
Ref_x = 1213999.5
Ref_y = 6153000.5
Max_x = 1213999.5 + nc*CellsizeX
Max_y = 6153000.5 + nr*CellsizeY
def Get_Altitude(dat, nr, GPS_x, GPS_y, Ref_x, Ref_y, CellsizeX, CellsizeY, index_x=None, index_y=None):
if ((index_x==None) and (index_y==None)):
# Get the distance between tile reference and point
dx = np.abs(GPS_x - Ref_x)
dy = np.abs(GPS_y - Ref_y)
# Get the index of the current point
index_x = int(np.round((dx/CellsizeX),0))
index_y = int(np.round((dy/CellsizeY),0))
# Get the altitude in the DEM
alt = dat[np.abs(-nr+index_y)][index_x]
else:
# Get the altitude in the DEM
alt = dat[np.abs(-nr+index_y)][index_x]
if ((alt == -99999.00) or (alt == -9999)):
alt = 0
return alt
lt = Get_Altitude(dat, nr, Max_x, Max_y, Ref_x, Ref_y, CellsizeX, CellsizeY, index_x=None, index_y=None)
print(lt)
But when I want to get the altitude of the reference of the tile (bottom left corner), or the maximum (top right corner), I have an error saying that the index is out of the array :
alt = dat[np.abs(-nr+index_y)][index_x]
IndexError: index 1000 is out of bounds for axis 0 with size 1000
The index of my numpy array goes from to 0 to 999, so 1000 columns (and rows).
I do not understand why the index calculated for the reference or the maximum is 1000... It should be 999.
Does anybody have an idea about how to fix this problem ? Thanks 🙂
