Correction - Exercice 4 : Ouverture d’un raster et enrichissement d’une couche vecteur

Afficher le profil en long du cours d’eau :

import matplotlib.pyplot as plt

# Nos 2 couches de référence à utiliser :
line_lyr = QgsProject.instance().mapLayersByName("Canal_Furon")[0]
dem = QgsProject.instance().mapLayersByName("DEM_N245E395_l93")[0]

# Récupération de la géométrie du cours d'eau
# (la couche ne contient qu'une entité, on peut donc prendre directement
# la première)
ft = list(line_lyr.getFeatures())[0]
geom = ft.geometry()
length = geom.length()

# Nombre de points / distance entre chaque point
nb_pts = 50
step = length / nb_pts

# Le fournisseur de données de la couche raster
dem_pr = dem.dataProvider()

# La liste qui va acceuillir les résultats
result = []

# Construction du résultat
for dist in range(0, int(length), int(step)):
    interp = geom.interpolate(dist) # Le resultat est de type `QgsGeometry`
    res = dem_pr.identify(interp.asPoint(), 1)
    if res.isValid():
        _values = res.results() # Retourne un dictionnaire
        # Distance à la source et la valeur de la 1ere bande (altitude)
        result.append((dist, _values[1]))

# Affichage avec matplotlib
plt.plot([i[0] for i in result], [i[1] for i in result])
plt.xlabel('Distance a la source')
plt.ylabel('Altitude')
plt.show()

Créer une couche de point 3D :

from PyQt5.QtCore import QVariant

line_lyr = QgsProject.instance().mapLayersByName("Canal_Furon")[0]
dem = QgsProject.instance().mapLayersByName("DEM_N245E395_l93")[0]

# Récupération de la géométrie du cours d'eau
ft = list(line_lyr.getFeatures())[0]
geom = ft.geometry()
length = geom.length()

# Nombre de points / distance entre chaque point
nb_pts = 50
step = length / nb_pts

# Le fournisseur de données de la couche raster
dem_pr = dem.dataProvider()

# Nouvelle couche (en mémoire) de points et son fournisseur de données
layer = QgsVectorLayer('Point?crs=EPSG:2154', 'PointsWithAltitude' , 'memory')
pr = layer.dataProvider()

# Spécification de deux champs
# (un entier pour l'ID et un entier pour la distance à la source)
pr.addAttributes([
    QgsField("ID",  QVariant.Int),
    QgsField("dist_source", QVariant.Double),
])
layer.updateFields()

for ix, dist in enumerate(range(0, int(length), int(step))):
    pt = geom.interpolate(dist).asPoint()
    # Avec la methode sample on recupere direction la
    # valeur a l'endroit du point 'pt' (1er argument)
    # pour la premiere bande (2eme argument)
    # Cette fonction renvoi un 'tuple' de deux valeurs
    # la premiere est la valeur a l'endroit du point
    # et la seconde une valeur booleene qui indique
    # que l'operation s'est bien deroulee
    alt, valid = dem_pr.sample(pt, 1)
    if valid:
        new_pt = QgsPoint(pt.x(), pt.y(), alt)
        feat = QgsFeature()
        feat.setGeometry(QgsGeometry(new_pt))
        feat.setAttributes([ix, dist])
        ok = pr.addFeature(feat)
        if not ok:
            # Peut se produire si le nombre ou le type d'attribut ne correspondent
            # pas à ceux attendus pour cette couche
            print("Erreur lors de la création de l'entité {}".format(ix))

# Mise à jour de l'étendue de la couche
layer.updateExtents()

# Ajout de la couche
QgsProject.instance().addMapLayer(layer)

Note

L’API exposée par PyQgis n’est pas toujours très « pythonique ». Certaines fonctions vues au cours de cet exercice ne vont pas générer une Exception comme c’est la convention en Python mais vont retourner une valeur booléenne qu’il faudra lire pour s’assurer de la bonne exécution de la fonction en question.

 for ix, dist in enumerate(range(0, int(length), int(step))):
     pt = geom.interpolate(dist).asPoint()
     alt, valid = data_pr.sample(pt, 1)
     if valid:
         new_pt = QgsPoint(pt)
         new_pt.setZ(alt)
         feat = QgsFeature()
         feat.setGeometry(QgsGeometry(new_pt))
         feat.setAttributes([ix, dist])
         ok = pr.addFeature(feat)
         if not ok:
             print("Erreur lors de la création de l'entité {}".format(ix))
     else:
         print(
           "Erreur lors de la récupération de l'altitude pour l'entité {}"
           .format(ix)")