# Python Script for getting elevation difference between two points

I have several stream segments 1000 Km long. I need to find the elevation difference between two consecutive points of distance 1 Km starting from the upstream to downstream. How can I get the elevation difference from DEM. I have stream segments in raster format and also in vector format. Hoping for your kind suggestion. It would be better if i got some idea on python script.

As a geologist, I often use this technique to make geological cross section in pure Python. I presented a complete solution in Python: Using vector and raster layers in a geological perspective, without GIS software (in French)

I present here a summary in English:

• to show you how to extract the elevation values ​​of a DEM
• how to treat these values

If you open a DEM with GDAL/OGR Python module:

``````from osgeo import gdal
# raster dem10m
file = 'dem10m.asc'
layer = gdal.Open(file)
gt =layer.GetGeoTransform()
bands = layer.RasterCount
print bands
1
print gt
(263104.72544800001, 10.002079999999999, 0.0, 155223.647811, 0.0, -10.002079999999999)
``````

As a result, you have the number of bands and the geotransform parameters. If you want to extract the value of the raster under a xy point:

``````x,y  = (263220.5,155110.6)
# transform to raster point coordinates
rasterx = int((x - gt) / gt)
rastery = int((y - gt) / gt)
# only one band here
array([])
``````

As it is a DEM, you get the elevation value under the point. With 3 raster bands with the same xy point you get 3 values (R,G,B). So you could make a function that allows to get the values of multiple rasters under a xy point:

``````def Val_raster(x,y,layer,bands,gt):
col=[]
px = int((x - gt) / gt)
py =int((y - gt) / gt)
for j in range(bands):
band = layer.GetRasterBand(j+1)
col.append(data)
return col
``````

application

``````# with a DEM (1 band)
px1 = int((x - gt1) / gt1)
py1 = int((y - gt1) / gt1)
print Val_raster(x,y,layer, band,gt)
 # elevation
# with a geological map (3 bands)
px2 = int((x - gt2) / gt2)
py2 = int((y - gt2) / gt2)
print Val_raster(x,y,couche2, bandes2,gt2)
[253, 215, 118] # RGB color
``````

After that, you process the line profile (which may have segments):

``````# creation of an empty ogr linestring to handle all possible segments of a line with  Union (combining the segements)
profilogr = ogr.Geometry(ogr.wkbLineString)
# open the profile shapefile
source = ogr.Open('profilline.shp')
cshp = source.GetLayer()
# union the segments of the line
for element in cshp:
geom =element.GetGeometryRef()
profilogr = profilogr.Union(geom)
``````

To generate equidistant points on the line, you can use the Shapely module with interpolate (easier than ogr)

``````from shapely.wkb import loads
# transformation in Shapely geometry
# creation the equidistant points on the line with a step of 20m
lenght=profilshp.length
x = []
y = []
z = []
# distance of the topographic profile
dista = []
for currentdistance  in range(0,lenght,20):
# creation of the point on the line
point = profilshp.interpolate(currentdistance)
xp,yp=point.x, point.y
x.append(xp)
y.append(yp)
# extraction of the elevation value from the MNT
z.append(Val_raster(xp,yp,layer, bands,gt)
dista.append(currentdistance)
``````

and the results (with also the RGB values of a geologic map) with the x,y,z, distance values of the lists In 3D with matplotlib and Visvis (x,y,z values) Cross sections (x, elevation from currentdistance (dista list)) with matplotlib:  • +1 for creating a great pythonic solution and for using matplotlib to create great looking figures. – Fezter Apr 27 '13 at 21:09
• Isn't this possible with arcpy ? – Ja Geo Apr 28 '13 at 9:01
• I don't know, I do not use ArcPy – gene Apr 28 '13 at 9:53
• You typed rasterx twice, it should be rasterx,rastery – Metiu Dec 29 '16 at 13:29