Script de Python para obtener la diferencia de altitud entre dos puntos

Question

Tengo varios segmentos de arroyos de 1000 km de longitud. Necesito encontrar la diferencia de elevación entre dos puntos consecutivos de distancia 1 Km comenzando desde el río arriba hasta el río abajo. ¿Cómo puedo obtener la diferencia de elevación a partir del MDE? Tengo segmentos de arroyos en formato raster y también en formato vectorial. Sería mejor si tengo alguna idea sobre el script de Python.

Answer 1

Como geólogo, a menudo utilizo esta técnica para hacer secciones transversales geológicas en Python puro. He presentado una solución completa en Pitón: Uso de capas vectoriales y rasterizadas en una perspectiva geológica, sin software SIG (en francés)

Aquí presento un resumen en inglés:

• para mostrarle cómo extraer los valores de elevación de un MDE
• cómo tratar estos valores

Si abre un DEM con el módulo GDAL/OGR Python:

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)

Como resultado, tiene el número de bandas y los parámetros de geotransformación. Si quieres extraer el valor del raster bajo un punto 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]]) 

Al tratarse de un MDE, se obtiene el valor de elevación bajo el punto. Con 3 bandas raster con el mismo punto xy se obtienen 3 valores (R,G,B). Así que se podría hacer una función que permita obtener los valores de múltiples rásteres bajo un punto 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

aplicación

# 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  

Después, se procesa el perfil de la línea (que puede tener segmentos):

# 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)

Para generar puntos equidistantes en la línea, se puede utilizar la función Shapely módulo con interpolación (más fácil que 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)

y los resultados (con también los valores RGB de un mapa geológico) con los valores x,y,z, de distancia de las listas En 3D con matplotlib y Visvis (valores x,y,z)

Secciones transversales (x, elevación de la distancia actual (lista dista)) con matplotlib :