Skip to content

TD 2 - Manipuler les couches vectorielles et raster

Le coeur de l'API Python pyQGIS réside dans la manipulation des couches de données. On peut créer et modifier les couches et les entités avec plus de contrôle qu'en passant par les algorithmes processing. Dans ce TD nous allons découvrir comment faire une requête de données depuis les couches et comment les modifier grâce aux classes QgsVectorLayer et QgsRasterLayer.

1 - Manipuler les couches vectorielles

Créez un projet QGIS et sauvegardez-le à l'emplacement de votre choix sous le nom « td2.qgz ».

1.1 Charger une couche vectorielle

Avec la classe QgsVectorLayer:

  1. Chargez-la couche « RPG » (fichier RPG.gpkg) dans une variable rpg.
  2. Affichez son nom.
  3. Ajoutez là au projet courant avec la classe QgsProject.

  4. Extrayez les informations suivantes:

    • Le SRC de la couche.
    • Le chemin du fichier de la couche sur le disque.
    • L'identifiant QGIS de la couche.

  5. Changez le nom de la couche en « RPG - parcelles ».

  6. Paramétrez le SRC sur 2134 (SRC équivalent à celui d'origine).
  7. Ajoutez la couche au projet.

Tip

La classe QgsVectorLayer hérite de la classe QgsMapLayer qui défini des méthodes génériques pour tout les types de couches dans QGIS (vecteur et raster). Pour trouver les fonctions qui sont génériques à toutes les couches QGIS il vous faut donc également aller inspecter la documentation de la classe QgsMapLayer.

Solution 
rpg = QgsVectorLayer("/path/to/RPG.gpkg", "RPG", "ogr")
rpg.name()

rpg.crs()
rpg.source()
rpg.id()

rpg.setName("RPG - parcelles")
rpg.setCrs(2134)

QgsProject.instance().addMapLayer(rpg)

1.2 Champs d'une couche vectorielle

Affichez dans la console tout les champs de la couche RPG avec leur noms et leur types.

Tip

Il existe plusieurs manières de récupérer ces informations. Vous pouvez utiliser des méthodes de la classe QgsFields ou bien itérer sur chaque champ (QgsField) individuellement.

Solution 
fields = rpg.fields()

# itérer sur les champs
for f in fields:
    print(f"name: {f.name()}, type: {f.type()}")

# Lister les champs
list(fields)

# lister les noms des champs uniquement
fields.names()

1.3 Modifier les champs

Cours

Il y a deux manières de modifier les couches vectorielles.

1-Utiliser un buffer d'édition

La première est de passer par un buffer d'édition qui permet de conserver en mémoire les modifications et de les annuler si besoin. C'est ce qui se passe quand vous éditez une couche via l'interface QGIS, il faut sauvegarder les modifications pour qu'elles soient effectives. C'est la méthode sûre mais elle a pour inconvénient d'être plus lente. Avec l'API il faut faire les étapes suivantes:

  1. Démarrer l'édition de la couche (startEditing()).
  2. Faire les modifications.
  3. Valider les modifications (commitChanges()).

Tant que l'édition est en cours vous pouvez utiliser rollBack() pour annuler les modifications. Comme c'est plus long à écrire, pyQGIS fourni un outil pour faciliter l'usage de cette méthode, c'est le gestionnaire de contexte edit. Si une erreur est remontée dans le contexte edit, rollBack() sera appliqué automatiquement.

2-Utiliser le fournisseur de données de la couche

Pour accéder aux données les couches vectorielles utilisent des fournisseur de données (dataProvider). En modifiant le fournisseur de données vous appliquez directement les modifications aux données sans possibilité de retour en arrière. C'est plus dangereux mais plus rapide.

Ajoutez un champ « flottant » (ou « double ») "AREA" aux champs de la couche « RPG - parcelles » en utilisant un buffer d'édition sans utiliser edit. Observez ensuite la table d'attributs de la couche « RPG ». Ensuite ajoutez un champ "AREA2" et "AREA3" en utilisant respectivement edit puis le fournisseur de données dataProvider.

Tip

  • Lorsque vous modifier les champs d'une couche vectorielle il faut ensuite les mettre à jour avec updateFields().
  • Dans pyQGIS il existe beaucoup de classes de type Enum, ce sont des entiers qui possèdent un nom associé et qui sont souvent utilisés pour définir des options comme par exemple les types de champs QMetaType. Il faut un peu chercher dans la documentation Python et dans la documentation Cpp pour trouver les correspondances entre les valeurs et les variables. Pour les types de champs vous pouvez suivre ce lien.
Solution 
# Buffer d'edition
rpg.startEditing()
rpg.addAttribute(QgsField("AREA", 6)) # 6 pour QMetaType.Double
rpg.commitChanges()
rpg.updateFields()

# Buffer d'edition avec le contexte edit
with edit(rpg):
    rpg.addAttribute(QgsField("AREA2", 6))
    rpg.updateFields()

# avec le fournisseur de donnees
rpg.dataProvider().addAttributes([QgsField("AREA3", 6)])
rpg.updateFields()

1.4 Index de champs

Retrouvez la position du champ "AREA" à partir de leur noms en utilisant la classe QgsFields. Faites la même chose avec le champ inconnu "XYZ" et observez la position rendue.

Solution 
area_id = rpg.fields().indexFromName("AREA")
xyz_id = rpg.fields().indexFromName("XYZ") # affiche -1 pour les champs inexistants dans la couche

1.5 Supprimer des champs

Supprimez les champs "AREA", "AREA2" et "AREA3" avec la classe QgsVectorLayer en passant par le fournisseur de données ou le buffer d'édition.

Solution 
area_id = rpg.fields().indexFromName("AREA")
area2_id = rpg.fields().indexFromName("AREA2")
area3_id = rpg.fields().indexFromName("AREA3")
rpg.dataProvider().deleteAttributes([area_id, area2_id, area3_id])
rpg.updateFields()

1.6 Métadonnées des entités

Cours

Les entités des couches vectorielles sont gérées par la classe QgsFeature. Cette classe permet de lire les données et géométries associées à chaque entité. Elle dispose des mêmes champs que la couche à laquelle elle est reliée. La classe QgsFeature est centrale dans l'API pyQGIS lorsqu'on travaille sur des couches vectorielles. La méthode id() de QgsFeature permet d'obtenir l'index de l'entité dans la couche, c'est cet index qui va permettre de faire les requêtes sur la couche vectorielles afin de récupérer les entités.

Extrayez le nombre et la liste des entités de la couche « RPG » avec la classe QgsVectorLayer.

Solution 
count = rpg.featureCount()
features = list(rpg.getFeatures())

1.7 Extraire des entités

Récupérez l'entité avec l'index 5. Puis récupérez en une commande les entités 10, 200 et 1500.

Solution 
feat = rpg.getFeature(5)
feat_10, feat_200, feat_1500 = tuple(rpg.getFeatures([10,20,1500]))

Cours

La classe QgsVectorLayer renverra toujours un objet de type QgsFeature même si l'entité correspondant à l'index demandé n'existe pas. Pour vous assurer de la validité d'une entité récupérée vous pouvez appeler la méthode isValid() sur l'entité: 

feat = rpg.getFeature(5)
feat.isValid()
True    
feat = rpg.getFeature(1000000)
feat.isValid()
False

NB: isValid() est une méthode commune à beaucoup de classe pyQGIS, notamment aux classes QgsVectorLayer et QgsRasterLayer.

1.8 Inspecter une entité

Récupérez les informations suivante de l'entité 5:

  • L'identifiant de l'entité (fid).
  • La géométrie de l'entité.
  • La liste des champs (idem que pour les couches).
  • La liste des attributs (les valeurs pour chaque champs).
  • La valeur pour le champ "CODE_CULTU".
  • La géométrie de l'entité.
Solution 
feature = rpg.getFeature(5)
feature.id() 
feature.fields().names()
feature.attributes() 
feature.geometry()
feature["CODE_CULTU"]

1.9 Modifier une entité

Sélectionnez l'entité 5 avec la classe QgsVectorLayer et zoomez dessus avec QgsMapCanvas. Modifiez l'entité en changeant la valeur pour le champ "CODE_CULTU" par "PPH" et mettez à jour la couche avec cette nouvelle entité. Observez le changement de valeur pour le champ "CODE_CULTU".

Tip

On peut modifier les attributs des entités de la même manière qu'on modifie un dictionnaire Python avec: 

feature["field"]
5
feature["field"] = 8
feature["field"] # vaut 8 maintenant

Solution 
rpg.select(5)
iface.mapCanvas().zoomToSelected()
feat = rpg.getFeature(5)
feat["CODE_CULTU"] = "PPH"

with edit(rpg):
    rpg.updateFeature(feat)

1.10 Ajouter et supprimer des entités

Créez une nouvelle entité avec les mêmes champs que ceux de la couche « RPG » puis ajoutez-là à la couche avec la classe QgsVectorLayer. Ensuite supprimez-là. Notez que l'entité ne possède pas de géométrie et donc ne sera pas affichée dans la couche.

Tip

Le champ "fid" a une contrainte d'unicité. Il ne peut pas y avoir deux entités avec le même identifiant donc pour ajoutez une entité à une couche pensez à mettre un "fid" unique.

Solution 
feat = QgsFeature(rpg.fields())
feat.setAttributes([rpg.featureCount() + 1, "156543514", "PPH"])
with edit(rpg):
    rpg.dataProvider().addFeature(feat)

with edit(rpg):
    rpg.deleteFeature(feat.id())

1.11 Inspecter une géométrie

Cours

Les géométries des entités sont gérées par la classe QgsGeometry. De nombreuses méthodes sont disponibles depuis cette classe pour appliquer des transformations aux géométries ou comparer des géométries entre elles. Les géométries stockées dans la classe QgsGeometry peuvent être de différentes natures (points, lignes, polygones, multi polygones) et donc toutes les méthodes de la classe QgsGeometry ne fonctionneront pas toutes selon le type de géométrie réellement stockée. Il faut donc choisir la méthode qui correspond au type de géométrie stockée. Pour récupérer les éléments de la géométrie en objets python manipulables on peut utiliser les méthodes avec la syntaxe as{GeometryType}() telles que asMultiPolygon().

Récupérez la géométrie de l'entité 100 de la couche RPG et inspectez:

  • Son type.
  • Sa représentation WKB.
  • Sa représentation WKT.
Solution 
feature = rpg.getFeature(100)
geometry = feature.geometry()
geometry.type()
geometry.asWkb()
geometry.asWkt()

1.12 Modifier une géométrie

Récupérez l'entité 5 et sa géométrie, appliquez un buffer de 10m sur la géométrie puis remplacer la géométrie de l'entité 5 par cette nouvelle géométrie avec la classe QgsFeature. Ensuite mettez à jour l'entité dans la couche. Répétez l'opération en remplaçant la géométrie par sa boîte englobante orientée la plus petite (orientedMinimumBoundingBox()).

Solution 
feat = rpg.getFeature(5)
geom = feat.geometry()
buffer_geom = geom.buffer(5,5)
feat.setGeometry(buffer_geom)
with edit(rpg):
    rpg.updateFeature(feat)

oriented_box = geometry.orientedMinimumBoundingBox()[0]
feat.setGeometry(oriented_box)

with edit(rpg):
    rpg.updateFeature(feat)

1.13 Comparaison de géométries

Cours

La classe QgsGeometry permet de faire des comparaisons et des calculs entre plusieurs géométries:

  1. Comparaisons:

    • intersects()
    • contains()
    • equal()
    • ...

  2. Calculs:

    • difference()
    • intersection()
    • combine()
    • ...

Récupérez les entités 1310 et 1312 ainsi que leur géométries. Appliquez un buffer de 100m sur la géométrie de l'entité 1312 et assignez cette nouvelle géométrie à la variable geom_buffer. Testez si la géométrie de l'entité 1310 intersecte geom_buffer. Ensuite calculer la géométrie résultante de l'intersection de la géométrie 1310 avec geom_buffer et mettez à jour l'entité 1310 avec cette nouvelle géométrie.

Solution 
feat_1310, feat_1312 = tuple(list(rpg.getFeatures([1310,1312])))
geom_1310, geom_1312 = feat_1310.geometry(), feat_1312.geometry()
geom_buffer = geom_1312.buffer(100,5)
intersection = geom_1310.intersection(geom_buffer)
feat_1310.setGeometry(intersection)
with edit(rpg):
    rpg.updateFeature(feat_1310)

1.14 Créer une géométrie

Récupérez l'entité 50 et l'entité 100. Créez une géométrie à partir de la représentation WKT de la géométrie de l'entité 50. Appliquez cette géométrie à l'entité 100 et mettez à jour la couche.

Solution 
feat_50, feat_100 = tuple(list(rpg.getFeatures([50,100])))
geom_50 = feat_50.geometry()
wkt_string = geom_50.asWkt()
new_geom = QgsGeometry.fromWkt(wkt_string)
feat_100.setGeometry(new_geom)
with edit(rpg):
    rpg.updateFeature(feat_100)

1.15 Modification fine de géométries

Cours

Les structures géométriques en QGIS sont représentées par des listes de points avec différents niveaux pour représenter les composants des géométries. Un polygone par exemple est composée d'une liste de points pour l'anneau extérieur et d'autant de listes de points supplémentaires pour les trous dans le polygone: 

polygon = [
    [QgsPointXY(), ...], # anneau extérieur
    [QgsPointXY, ...], # Trou 1
    [QgsPointXY, ...], # Trou 2
    ...
]
Et par conséquent un multi polygone est une liste de polygones donc contiendra un niveau de liste supplémentaire.

Récupérer la géométrie de l'entité 20 (c'est un MultiPolygone), transformez son polygone (uniquement l'anneau extérieur, vous pouvez ignorer les trous) pour qu'ils ne contienne plus qu'un point sur deux. Modifiez la géométrie de l’entité et mettez à jour la couche.

Tip

Pour recréer une QgsGeometry depuis une structure en liste de points (polygones, etc) vous pouvez utiliser les méthodes avec la syntaxe QgsGeometry.from{StructureType}() et en l’occurrence il s'agit de la méthode QgsGeometry.fromMultiPolygonXY().

Solution 
feat = rpg.getFeature(20)
geometry = feat.geometry()
multipolygons = geometry.asMultiPolygon()
polygon = multipolygons[0][0] # anneau extérieur
selected_points = polygon[::2]
new_geom = QgsGeometry.fromMultiPolygonXY([[selected_points]])
feat.setGeometry(new_geom)
with edit(rpg):
    rpg.updateFeature(feat)

1.16 Requête spatiale

Cours

On peut faire des requêtes d'entités via la classe QgsVectorLayer et sa méthode getFeatures(). Pour cela il faut passer un objet QgsFeatureRequest comme argument de getFeatures().

Récupérez l'entité 300 depuis la couche RPG. Appliquez un buffer de 5km à la géométrie récupérez sont étendue. Extrayez toutes les entités de la couche RPG qui intersectent cette étendue et calculez l'aire totale de ces géométries.

Tip

L'étendue d'une géométrie peut être récupérée avec la méthode extent() de QgsGeometry.

Solution 
feat = rpg.getFeature(300)
extent = feat.geometry().buffer(5000,5).boundingBox()
query = QgsFeatureRequest().setFilterRect(extent)
features = rpg.getFeatures(query)
geometries = [f.geometry() for f in features]
area_tot = sum([g.area() for g in geometries])

1.17 Utiliser un index spatial

Créez un index spatial sur la couche RPG puis répétez l'énoncé de la question 1.15.

NB: L'index spatial améliore significativement les performances de la requête sur de grandes couches de données. La couche « RPG » ne possède pas assez d'entités pour qu'une différence de performance puisse être remarquée.

Tip

C'est la classe SpatialIndex qui permet la création d'index spatiaux et la récupération des identifiants des entités lorsqu'on effectue une requête.

Solution 
feat = rpg.getFeature(537)
geom_buffer = feat.geometry().buffer(5000,5)
spatial_index = QgsSpatialIndex(rpg.getFeatures())
query = spatial_index.intersects(geom_buffer.boundingBox())
features = rpg.getFeatures(query)
geometries = [f.geometry() for f in features]
area_tot = sum([g.area() for g in geometries])

1.18 Requête par expression

Cours

QGIS dispose d'un moteur d'expression permettant de travailler sur l'ensemble des entités d'une couche avec une certaine flexibilité. Elle sont par exemple utilisées pour:

  • Faire des requêtes sur les entités.
  • Calculer de nouveaux champs.
  • Trier len entités.
  • etc.

L'utilisation des expressions en pyQGIS se fait via la classe QgsExpression mais certaines méthodes prennent directement des chaînes de caractères qu'elles transforment en QgsExpression. Dans ce cours nous n'utiliseront pas directement la classe QgsExpression mais si vous souhaitez apprendre comment l'utiliser vous pouvez suivre le tutoriel sur les expressions fourni par la documentation pyQGIS.

Créez une couche qui ne contient que les entités dont le champ "CODE_CULTU" vaut "PPH" en utilisant QgsFeatureRequest et la méthode setFilterExpression().

Tip

La méthode materialize() de QgsVectorLayer (héritée de QgsFeatureSource) permet de créer une couche vectorielle depuis une liste d'index des entités.

Solution 
query = QgsFeatureRequest().setFilterExpression("CODE_CULTU = 'PPH'")
rpg_selected = rpg.materialize(query)

2 - Manipuler les couches raster

Les possibilités de manipulation des couches raster avec pyQGIS sont plus réduites que pour les couches vectorielles. On utilisera majoritairement des algorithmes processing existants pour effectuer des opérations sur les raster.

2.1 Charger une couche raster

Chargez la couche raster « indices » (fichier indices.tiff) avec QgsRasterLayer et ajoutez-là au projet. Ensuite extrayez les informations suivantes:

  • Nombre de bandes.
  • Nom des bandes.
  • Hauteur et largeur du raster.
  • Unités par pixel de hauteur/largeur.
  • Étendue du raster.
Solution 
indices = QgsRasterLayer("/path/to/indices.tiff")
QgsProject.instance().addMapLayer(indices)
indices.bandCount()
indices.bandName(x) # x est l'indice de la bande (de 1 à N bandes)
indices.height()
indices.width()
indices.rasterUnitsPerPixelY()
indices.rasterUnitsPerPixelX()
indices.extent()

2.2 Extraire la valeur d'un pixel

Récupérer le fournisseur de données de la couche avec la méthode dataProvider() puis positionnez votre curseur de souris sur un des pixels du raster et extrayez la valeur de ce pixel pour la bande « NDVI » avec la méthode sample() du fournisseur de données.

Tip

Utilisez iface pour accéder au canvas et à la position de votre curseur.

Solution 
point = iface.mapCanvas().center()
value, _ = indices.dataProvider().sample(point, 1)

2.3 Transformation en array numpy

Cours

Pour une analyse plus fine du raster on peut transformer la couche raster en array numpy avec la méthode as_numpy(). Cependant une array numpy est chargée en mémoire et peut donc surcharger la mémoire si elle trop grande. Il existe des classes pyQGIS pour lire des portions spécifiques du raster mais c'est une partie de l'API pyQGIS qui est assez complexe à utiliser et que nous n'aborderons pas dans ce cours. Si vous souhaitez vous renseigner sur ces options vous pouvez regarder le tutoriel sur les blocks raster.

Extrayez l'array numpy qui correspond à la couche raster et calculez la médiane de chaque bandes.

Solution 
import numpy as np
array = indices.as_numpy()
median = np.nanmedian(array, axis=(1,2))