Exercice 4 : Ouverture d'un raster et enrichissement d'une couche vecteur ========================================================================= .. figure:: img/ex5_mnt_color.png :width: 40% :align: center :figwidth: 98% 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**. .. figure:: img/ex5_matplotlib.png :width: 40% :align: center :figwidth: 98% Données ------- | • Tracé du Canal du Furon ``Canal_Furon.shp`` (Ligne - Source : IGN [#f1]_) ➜ `Téléchargement `_ | • Modèle Numérique de terrain ``DEM_N245E395_l93.tif`` (Raster float 32bits - Source : EAA [#f2]_) ➜ `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 : .. figure:: img/ex5_matplotlib.png :width: 40% :align: center :figwidth: 98% - **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*). .. rubric:: Footnotes .. [#f1] Institut national de l'information géographique et forestière : http://professionnels.ign.fr/bdtopo. Données téléchargées le 28/02/2019. .. [#f2] European Environment Agency : https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem. Données téléchargées le 28/02/2019. .. raw:: html

💡 Solution de l'exercice 4

.. raw:: html


Quelques mots supplémentaires sur l'API PyQGIS

.. container:: spoiler blq .. 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``). .. code-block:: python :emphasize-lines: 3,4,10,11 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)")