Script Python pour obtenir la différence d'altitude entre deux points [fermé]

J'ai plusieurs segments de cours d'eau de 1000 km de long. J'ai besoin de trouver la différence d'altitude entre deux points de distance consécutifs de 1 km à partir de l'amont vers l'aval. Comment puis-je obtenir la différence d'altitude à partir de DEM? J'ai des segments de flux au format raster et également au format vectoriel. Ce serait mieux si j'avais une idée du script Python.

En tant que géologue, j'utilise souvent cette technique pour réaliser des coupes géologiques en Python pur. J'ai présenté une solution complète en Python: en utilisant le vecteur et les couches raster dans un point de vue géologique, sans logiciel SIG (en français)

Je présente ici un résumé en anglais:

• pour vous montrer comment extraire les valeurs d'élévation d'un DEM
• comment traiter ces valeurs

Si vous ouvrez un DEM avec le module Python GDAL / OGR:

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)

Par conséquent, vous avez le nombre de bandes et les paramètres de géotransformation. Si vous souhaitez extraire la valeur du raster sous un point xy:

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)
array([[222]]) 

Comme il s'agit d'un DEM, vous obtenez la valeur d'élévation sous le point. Avec 3 bandes raster avec le même point xy, vous obtenez 3 valeurs (R, G, B). Vous pouvez donc créer une fonction qui permet d'obtenir les valeurs de plusieurs rasters sous un point xy:

def Val_raster(x,y,layer,bands,gt):
    col=[]
    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)
        col.append(data[0][0])
  return col

application

# 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  

Après cela, vous traitez le profil de ligne (qui peut avoir des 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)

Pour générer des points équidistants sur la ligne, vous pouvez utiliser le module Shapely avec interpolation (plus facile que l'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
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)[0]
     dista.append(currentdistance)

et les résultats (avec aussi les valeurs RVB d'une carte géologique) avec les valeurs x, y, z, distance des listes En 3D avec matplotlib et Visvis ( valeurs x, y, z)

Coupes (x, élévation de la distance actuelle (liste dista )) avec matplotlib :