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. It would be better if I got some idea on Python script.


1 Answer 1


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
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[0]) / gt[1])
rastery = int((y - gt[3]) / gt[5])
# only one band here
print layer.GetRasterBand(1).ReadAsArray(rasterx,rastery, 1, 1)

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):
    px = int((x - gt[0]) / gt[1])
    py =int((y - gt[3]) / gt[5])
    for j in range(bands):
        band = layer.GetRasterBand(j+1)
        data = band.ReadAsArray(px,py, 1, 1)
  return col


# with a DEM (1 band)
px1 = int((x - gt1[0]) / gt1[1])
py1 = int((y - gt1[3]) / gt1[5])
print Val_raster(x,y,layer, band,gt)
[222] # elevation
# with a geological map (3 bands)
px2 = int((x - gt2[0]) / gt2[1])
py2 = int((y - gt2[3]) / gt2[5])
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
profilshp = loads(profilogr.ExportToWkb())
# creation the equidistant points on the line with a step of 20m
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
     # extraction of the elevation value from the MNT
     z.append(Val_raster(xp,yp,layer, bands,gt)[0]

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)

enter image description here

Cross sections (x, elevation from currentdistance (dista list)) with matplotlib: enter image description here

enter image description here

  • 3
    +1 for creating a great pythonic solution and for using matplotlib to create great looking figures.
    – Fezter
    Apr 27, 2013 at 21:09
  • Isn't this possible with arcpy ?
    – Ja Geo
    Apr 28, 2013 at 9:01
  • 3
    I don't know, I do not use ArcPy
    – gene
    Apr 28, 2013 at 9:53
  • You typed rasterx twice, it should be rasterx,rastery
    – Metiu
    Dec 29, 2016 at 13:29

Not the answer you're looking for? Browse other questions tagged or ask your own question.