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

Objectif¶

L’objectif de cet exercice est de créer le profil longitudinal d’un cours d’eau.

Cette opération nécessite de disposer de points à distance régulière le long du cours d’eau. Une valeur d’altitude va ensuite être récupérée à partir d’un modèle numérique de terrain.

La dernière partie de l’exercice aborde la création d’une nouvelle couche temporaire en mémoire.

Données¶

• Tracé du Canal du Furon Canal_Furon.shp (Ligne - Source : IGN [1]) ➜ Téléchargement

• Modèle Numérique de terrain DEM_N245E395_l93.tif (Raster float 32bits - Source : EAA [2]) ➜ Téléchargement

Les deux jeux de données utilisent la projection Lambert-93 (EPSG:2154).

Procédure¶

Ouvrir QGIS.
Créer un nouveau projet.
Décompresser les deux archives zip qui ont été téléchargées et ajouter les couches Canal_Furon.shp et DEM_N245E395_l93.tif.
Ouvrir la console Python.

Obtenir la référence de chacune de nos couches :

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

Disposer des points à égale distance le long du tracé du cours d’eau :

# Quelle est la longueur du cours d'eau ?
ft = list(line_lyr.getFeatures())[0]
geom = ft.geometry()
length = geom.length() # En unités de la carte (ici en m.)

# Combien de points ?
nb_pts = 50

# Distance entre points successifs
step = length / nb_pts

# Création d'un itérateur approprié et ajout des points créés à une liste :
result = []
for dist in range(0, int(length), int(step)):
    # Le resultat va être de type `QgsGeometry`
    interp = geom.interpolate(dist)
    # Convertit la géométrie en `QgsPointXY` et l'ajoute au résultat
    result.append(interp.asPoint())

Récupérer la valeur de la (ou des) bande(s) raster à une coordonnée donnée :

dem_pr = dem.dataProvider()
pt = QgsPointXY(925189, 6473673)
res = dem_pr.identify(pt)
if res.isValid():
    _values = res.results() # Retourne un dictionnaire
    # Affiche la valeur pour la première bande
    print(_values[1])
    # Mieux formatté... :
    print("Le point ({}, {}) est situé à environ {:.2f}m d'altitude"
        .format(pt.x(), pt.y(), _values[1]))

« Putting it all together »

1. Vous disposez des briques essentielles pour créer un objet python de type list, contenant, pour une cinquantaine de points répartis à égale distance le long du cours d’eau, la distance à la source et l’altitude de chacun d’entre eux, par exemple sous la forme d’un tuple de 2 éléments (les coordonnées précises de chaque point ne nous intéressent ici que pour connaitre son altitude). Ces deux informations, distance à la source et altitude, seront par la suite utilisées, respectivement en abcisses et en ordonnées, pour tracer le profil du cours d’eau.

2. Lisez le tutoriel https://futurestud.io/tutorials/matplotlib-simple-line-plots et essayez de mobiliser ces nouvelles notions pour utiliser la liste obtenue à l’étape (1) afin d’obtenir le résultat suivant :
Déjà fini ? Utilisez maintenant ces points pour créer une couche de point 3D Dans cet exercice et contrairement à l’exemple précédent, les coordonnées de chaque nouveau point nous intéressent : elles vont être enrichie de l’altitude du point. Ces points devront être ajoutés à une nouvelle couche en mémoire qui aura précédemment été créée avec deux champs attributaires, l’un pour l’identifiant du point (de type entier) et l’autre pour la distance à la source (de type double flottant).

Footnotes

💡 Solution de l'exercice 4

Quelques mots supplémentaires sur l'API PyQGIS

Note

L’API exposée par PyQgis n’est pas toujours très « pythonique ». Certaines fonctions relative à la manipulation des données raster ou vecteur 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 (pensez aux exercices précédent lorsque vous avez pu voir True - ou False - affiché dans l’interpréteur lors de certaines opérations telles que l’appel de la méthode updateFeature).

 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.x(), pt.y(), 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)")

Contenu

Sujet précédent

Sujet suivant

Cette page

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

Objectif¶

Données¶

Procédure¶