Get the elevation in an ASCII raster DEM format using python

136
1
06-30-2022 02:22 AM
Titiiii
New Contributor

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 🙂

0 Kudos
1 Reply
JohannesLindner
MVP Regular Contributor

I do not understand why the index calculated for the reference or the maximum is 1000... It should be 999.

When you use Max_X as GPS_X input, dx will be equal to nc*CellsizeX, which would be 1000. index_x = dx/CellsizeX will still be 1000. See lines 16, 22, and 25 in your code, analogous for dy.

 

Does anybody have an idea about how to fix this problem ?

The root of the problem are lines 16 and 17. Here, you're telling Python that your raster has 1001 rows/columns: the start row + 1000 rows after that.

This should fix it:

Max_x = Ref_X + (nc - 1) * CellsizeX

 


Have a great day!
Johannes
0 Kudos