Aller au contenu

Fonctions sur une couche vecteur#

Utilisation des expressions QGIS#

  • Les expressions sont très présentes dans QGIS, tant dans l'interface graphique que dans l'utilisation en Python.
  • Nous partons de la couche des COMMUNES uniquement chargé dans QGIS.

Tip

Adapter le numéro des codes INSEE ou des départements selon votre BDTOPO 😉

Sélection d'entité#

Nous souhaitons sélectionner les entités dont le code INSEE commence par 77. Commençons par faire cela graphiquement dans QGIS Bureautique. À l'aide d'une expression QGIS, sélectionner les codes INSEE qui commencent par 77 (à choisir un code INSEE propre au jeu de données).

Sélectionner par expression

Solution en mode graphique :

1
"INSEE_COM" LIKE '77%'

Nous allons faire la même chose, mais en utilisant Python. Pensez à désélectionner les entités.

Il va falloir "échapper" un caractère à l'aide de \. Voir la page Wikipédia sur l'échappement ou ce meme pour les devs 🫢

1
2
3
4
5
6
7
8
from qgis.utils import iface

layer = iface.activeLayer()
layer.removeSelection()
layer.selectByExpression(f"\"INSEE_COM\" LIKE '77%'")
# layer.selectByExpression(f'"INSEE_COM" LIKE \'77%\'')  # Résultat identique
layer.invertSelection()
layer.removeSelection()

Le raccourci iface.activeLayer() est très pratique, mais de temps en temps, on a besoin de plusieurs couches qui sont déjà dans la légende. Il existe dans QgsProject plusieurs méthodes pour récupérer des couches dans la légende :

1
2
3
4
5
from qgis.core import QgsProject

projet = QgsProject.instance()
communes = projet.mapLayersByName('communes')[0]
insee = projet.mapLayersByName('tableau INSEE')

Notons le s dans mapLayersByName. Il peut y avoir plusieurs couches avec ce même nom de couche. La fonction retourne donc une liste de couches. Il convient alors de regarder si la liste est vide ou si elle contient plusieurs couches avec len(communes) par exemple.

Warning

mapLayersByName fait uniquement une recherche stricte, sensible à la casse. Il faut passer par du code Python "pure" en itérant sur l'ensemble des couches, indépendamment de leur nom si l'on souhaite faire une recherche plus fine. Si vraiment, on a besoin, on peut utiliser le module re (lien du Docteur Python).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
from qgis.core import QgsProject

projet = QgsProject.instance()
communes = projet.mapLayersByName('communes')

if len(communes) == 0:
    print("Pas de couches dans la légende qui se nomme 'communes'")
    layer = None
elif len(communes) >= 1:
    # TODO FIX ME, pas forcément la bonne couche 'communes'
    layer = communes[0]

Exemple d'une sélection avec un export#

On souhaite pouvoir exporter les communes par département. On peut créer une variable depts = ('34', '30') puis boucler dessus pour exporter les entités sélectionnées dans un nouveau fichier.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from pathlib import Path
from qgis.core import QgsProject, QgsVectorFileWriter
from qgis.utils import iface

layer = iface.activeLayer()

options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = 'ESRI Shapefile'
options.fileEncoding = 'UTF-8'
options.onlySelectedFeatures = True  # Nouvelle option pour la sélection

depts = ('34', '30')
for dept in depts:
    print(f"Dept {dept}")
    layer.selectByExpression(f"\"INSEE_DEP\"  =  '{dept}'")
    result = QgsVectorFileWriter.writeAsVectorFormatV3(
        layer,
        str(Path(QgsProject.instance().homePath()).joinpath(f'{dept}.shp')),
        QgsProject.instance().transformContext(),
        options
    )
    print(result)
    if result[0] == QgsVectorFileWriter.WriterError.NoError:
        print(" → OK")

Bonus

Si l'on souhaite parcourir automatiquement les départements existants, on peut récupérer les valeurs uniques. Pour cela, il faut modifier deux lignes :

1
2
index = layer.fields().indexFromName("INSEE_DEP")
for dept in layer.uniqueValues(index):

Boucler sur les entités d'une couche sans expression#

Si besoin, pour que la suite de l'exercice soit plus rapide, on peut utiliser une couche ARRONDISSEMENT par exemple.

On peut parcourir les entités d'une couche QgsVectorLayer à l'aide de getFeatures().

Info

Avec PyQGIS, on peut accéder aux attributs d'une QgsFeature simplement avec l'opérateur [] sur l'objet courant comme s'il s'agissait d'un dictionnaire Python :

1
2
# Pour accéder au champ "NOM" de l'entité "feature" :
print(feature['NOM'])

On peut le voir dans les exemples attribute de QgsFeature : https://qgis.org/pyqgis/3.34/core/QgsFeature.html#qgis.core.QgsFeature.attribute

1
2
3
4
5
6
7
from qgis.utils import iface

layer = iface.activeLayer()
for feature in layer.getFeatures():
    print(feature)
    print(feature['NOM'])
    print(feature.attribute('NOM'))

Boucler sur les entités à l'aide d'une expression#

L'objectif est d'afficher dans la console le nom des communes dont la code département INSEE_DEP correspond uniquement à un seul département arbitraire.

L'exemple à ne pas faire, même si cela fonctionne (car on peut l'optimiser très facilement) :

1
2
3
4
5
6
from qgis.utils import iface

layer = iface.activeLayer()
for feature in layer.getFeatures():
    if feature['INSEE_DEP'] == '84':
        print(f'{feature['NOM']} : département {feature['INSEE_DEP']}')
  1. Imaginons qu'il s'agisse d'une couche PostgreSQL, sur un serveur distant
  2. On demande à QGIS de récupérer l'ensemble de la table distante, équivalent à SELECT * FROM ma_table
  3. Puis, on filtre dans QGIS (toute la donnée est présente dans QGIS Bureautique désormais)

Tip

Ce qui prend du temps lors de l'exécution, c'est surtout le print en lui-même. Si vous n'utilisez pas print, mais un autre traitement, cela sera plus rapide. Un simple print ralenti l'exécution d'un script.

Optimisation de la requête#

Dans la documentation, observez bien la signature de la fonction getFeatures. Que remarquez-vous ? Utilisons donc une expression pour limiter les résultats.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
from qgis.utils import iface
from qgis.core import QgsFeatureRequest

layer = iface.activeLayer()

request = QgsFeatureRequest()
# Équivalent à SELECT * FROM ma_table WHERE "INSEE_DEP" = '84'
request.setFilterExpression('"INSEE_DEP" = \'84\'')

for feature in layer.getFeatures(request):
    print(f'{feature['NOM']} : département {feature['INSEE_DEP']}')

Nous pouvons accessoirement ordonner les résultats et surtout encore optimiser la requête en :

  • Ne demandant pas de charger la géométrie
  • Ne demandant pas de charger tous les attributs, par exemple, on souhaite afficher que le nom de la commune et sa population.
La solution pour les experts 
1
2
3
4
5
6
7
8
9
request = QgsFeatureRequest()
request.setFilterExpression('"INSEE_DEP" = \'84\'')
request.addOrderBy('NOM')
request.setFlags(QgsFeatureRequest.NoGeometry)
# request.setSubsetOfAttributes([1, 4]) autre manière moins pratique, historique
request.setSubsetOfAttributes(['NOM', 'POPULATION'], layer.fields())
# # Équivalent à SELECT NOM, POPULATION FROM ma_table WHERE "INSEE_DEP" = '84' ORDER BY NOM
for feature in layer.getFeatures(request):
    print('{commune} : {nombre} habitants'.format(commune=feature['NOM'], nombre=feature['POPULATION']))
  • Faire le test en affichant un champ qui n'est pas dans la requête.

Enregistrement d'une requête dans une couche en mémoire#

Si l'on souhaite "enregistrer" le résultat de cette expression QGIS, on peut la matérialiser dans une nouvelle couche :

1
2
memory_layer = layer.materialize(request)
QgsProject.instance().addMapLayer(memory_layer)

Warning

Attention à la ligne iface.activeLayer() qui peut changer lors de l'ajout d'une nouvelle couche dans la légende.

Regardons le résultat et corrigeons ce problème d'export afin d'obtenir les géométries et les attributs, il faut supprimer la ligne NoGeometry si vous l'avez.

Valeur NULL#

En PyQGIS, il existe la valeur NULL qui peut être présente dans la table attributaire d'une couche vecteur.

1
2
3
4
5
6
7
8
from qgis.PyQt.QtCore import NULL

if feature['nom_attribut'] == NULL:
    # Traiter la valeur NULL
    pass
else:
    # Continuer
    pass

Calculer un champ "densite"#

Nous souhaitons avoir une colonne densite dans notre table attributaire, avec la densité de population.

Mais regardons avant la gestion des erreurs lors d'un traitement. En effet, nous allons vouloir "caster" (transformer le type) de la variable population en entier, mais attention, il y a des valeurs NC dans les valeurs.

_Note, il n'y a désormais plus de valeur NC dans le champ POPULATION dans la donnée, mais imaginons. Il peut s'agir d'une autre couche dont on ne connait pas la provenance et le contenu.

Les exceptions en Python#

Avant de traiter cet exercice, nous devons voir ce qu'est une exception en Python.

À plusieurs reprises depuis le début de la formation, il est fort à parier que nous ayons des messages en rouges dans la console de temps en temps. Ce sont des exceptions. C'est une notion de programmation qui existe dans beaucoup de languages.

Dans le langage informatique, une exception peut-être :

  • levée ("raise" en anglais) pour déclencher une erreur
  • attrapée ("catch" en anglais, ou plutôt "except" en Python) pour traiter l'erreur

Essayons dans la console de faire une opération 10 / 2 :

1
10 / 2

Essayons cette fois-ci 10 / 0, ce qui est mathématiquement impossible : 

1
10 / 0

Passons cette fois-ci dans un script pour que cela soit plus simple, et voir que le script s'arrête brutalement 😉

1
2
3
print('Début')
print(10 / 2)
print('Fin')

On peut "attraper" cette erreur Python à l'aide d'un try ... except... :

1
2
3
4
5
6
print('Début')
try:
    print(10 / 2)
except ZeroDivisionError:
    print('Ceci est une division par zéro !')
print('Fin')

Le try permet d'essayer le code qui suit. Le except permet d'attraper en filtrant s'il y a des exceptions et de traiter l'erreur si besoin.

Tip

On peut avoir une ou plusieurs lignes de code dans chacun de ces blocs. On peut appeler des fonctions, etc.

Une exception remonte le fil d'exécution du programme#

Important, une exception remonte tant qu'elle n'est pas attrapée : 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
def function_3():
    print("Début fonction 3")
    a = 10
    b = 2
    print(f"→ {a} / {b} = {a/b}")
    print("Fin fonction 3")

def function_2():
    print("Début fonction 2")
    function_3()
    print("Fin fonction 2")

def function_1():
    print("Début fonction 1")
    function_2()
    print("Fin fonction 1")

function_1()

Testons désormais d'attraper l'erreur dans la fonction 1 :

1
2
3
4
try:
    function_2()
except ZeroDivisionError:
    print("Fin de l'exception")

On voit que Python, quand il peut, nous indique la "stacktrace" ou encore "traceback", c'est-à-dire une sorte de fil d'ariane.

Héritage des exceptions#

Toutes les exceptions héritent de Exception donc le code ci-dessous fonctionne, mais n'est pas recommandé, car il masque d'autres erreurs :

1
2
3
4
5
6
7
8
try:
    a = 10 / 5

    mois = ['janvier', 'février']
    b = mois[0]

except Exception:
    print('Erreur inconnue')

On peut par contre "enchaîner" les exceptions, afin de filtrer progressivement les exceptions.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
try:
    a = 10 / 5
    # a = 10 / 0
    # a = 10 / int("NC")
    # a = 10 / "NC"
except ZeroDivisionError:
    print('Erreur, division par 0')
except ValueError:
    print("Erreur, il n'y a avait pas que des chiffres.")
except Exception:
    print('Erreur inconnue')

Il existe d'autres mots-clés en Python pour les exceptions comme finally: et else:. Voir un autre tutoriel.

Évidement, on peut vérifier la valeur de b en amont si c'est égal à 0. Mais ceci est pour présenter le concept des exceptions en Python.

Retour à l'exercice#

On souhaite donc savoir si un nombre est transformable en entier, dans le cas de la population (s'il y a NC par exemple) :

1
2
int('10')
int('NC')

Correction possible de l'exercice :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
from qgis.utils import iface
from qgis.core import QgsFeatureRequest

layer = iface.activeLayer()
request = QgsFeatureRequest()
# request.setLimit(5)  # Pour aller plus vite si-besoin
request.addOrderBy('NOM')
request.setSubsetOfAttributes(['NOM', 'POPULATION'], layer.fields())
for feature in layer.getFeatures(request):
    area = feature.geometry().area() / 1000000
    try:
        population = int(feature['POPULATION'])
        # Par exemple, pour nettoyer une chaîne de caractère en Python des espaces avant/après : "  Bonjour  ".strip()
        # "rstrip"/"lstrip" existent également
    except ValueError:
        population = 0

    densite = population/area

    print(f"{feature['NOM']} : {densite} habitants/km²")

Nous souhaitons enregistrer ces informations dans une vraie table avec un nouveau champ densite_population.

Solution possible :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
from qgis.utils import iface
from qgis.core import QgsFeatureRequest, QgsField, edit

from qgis.PyQt.QtCore import QVariant

layer = iface.activeLayer()

if 'densite' not in layer.fields().names():
    with edit(layer):
        field = QgsField('densite', QVariant.Double, prec=2, len=2)
        layer.addAttribute(field)

index = layer.fields().indexFromName('densite')
layer.startEditing()
request = QgsFeatureRequest()
# request.setLimit(5)  # Pour aller plus vite si-besoin
request.addOrderBy('NOM')
request.setSubsetOfAttributes(['NOM', 'POPULATION'], layer.fields())
# SELECT NOM, POPULATION FROM communes
for feature in layer.getFeatures(request):
    area = feature.geometry().area() / 1000000
    try:
        population = int(feature['POPULATION'])
    except ValueError:
        population = 0

    densite = population/area

    # Cette ligne n'aura aucun effet, contrairement à l'exercice d'export au format CSV précédent.
    # https://docs.3liz.org/formation-pyqgis/fonctions-scripts/#solution-possible
    # La variable "feature" est une copie, comme un peu le résultat du SELECT * FROM ma_table LIMIT 5
    # Un SELECT est en lecture seule. Ce n'est pas comme ça que cela se passe, c'est pour imager.
    feature['densite'] = densite

    # Uniquement l'appel à "changeAttributeValue" fonctionne
    # Pour information, il existe "changeGeometry" pour la même raison
    # Un peu comme la commande SQL UPDATE, sur une entité existante, bien qu'il ne faille pas oublier la session d'édition.
    layer.changeAttributeValue(feature.id(), index, densite)
    # print('{commune} : {densite} habitants/km²'.format(commune=feature['NOM'], densite=round(population/area,2)))

layer.commitChanges()

Calculer deux champs en utilisant la géométrie et une reprojection à la volée#

Manipulons désormais la géométrie en ajoutant le centroïde de la commune dans une colonne latitude et longitude en degrées.

Warning

TODO, en cours de correction, suppression de la variable petite_communes

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
from qgis.utils import iface
from qgis.core import QgsFeatureRequest, QgsField, edit

layer = iface.activeLayer()

request = QgsFeatureRequest()
request.setFilterExpression('to_int( "POPULATION" ) < 1000')
petites_communes = layer.materialize(request)

with edit(petites_communes):
    petites_communes.addAttribute(QgsField('densite_population', QVariant.Double))

    # /!\ Ajouter les 2 lignes ci-dessous
    petites_communes.addAttribute(QgsField('longitude', QVariant.Double))
    petites_communes.addAttribute(QgsField('latitude', QVariant.Double))

request = QgsFeatureRequest()
request.setSubsetOfAttributes([4])

# /!\ Ajouter les 2 lignes ci-dessous à propos de la transformation
transform = QgsCoordinateTransform(
    QgsCoordinateReferenceSystem("EPSG:2154"), QgsCoordinateReferenceSystem("EPSG:4326"), QgsProject.instance())

with edit(petites_communes):
    for feature in petites_communes.getFeatures(request):
        area = feature.geometry().area() / 1000000
        population = int(feature['POPULATION'])
        densite=population/area
        petites_communes.changeAttributeValue(feature.id(), 5, densite)

        # /!\ Ajouter les lignes ci-dessous
        geom = feature.geometry()
        # La transformation affecte directement l'objet Python en cours, mais pas l'entité dans la couche
        geom.transform(transform)
        centroid = geom.centroid().asPoint()
        petites_communes.changeAttributeValue(feature.id(), 6, centroid.x())
        petites_communes.changeAttributeValue(feature.id(), 7, centroid.y())

QgsProject.instance().addMapLayer(petites_communes)