Get the elevation in an ASCII raster DEM format using python

Titiiii — Thu, 30 Jun 2022 09:22:44 GMT

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 🙂

Re: Get the elevation in an ASCII raster DEM format using python

JohannesLindner — Thu, 30 Jun 2022 09:56:10 GMT

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

topic Re: Get the elevation in an ASCII raster DEM format using python in Python Questions

Get the elevation in an ASCII raster DEM format using python

Re: Get the elevation in an ASCII raster DEM format using python