Exercice 3 : Géométries et calcul de champ

_images/ex3_intro.png

Objectif

L’objectif de cet exercice est double : manipuler les géométries relatives à chaque entité et éditer une couche pour l’enrichir.

À partir de la couche vectorielle contenant les pays du monde, nous cherchons à trouver tous les polygones limitrophes à chaque polygone (tous les voisins de chaque pays).

La couche d’entrée va être enrichie de deux nouveaux champs, l’un contenant le noms des voisins et l’autre contenant la population totale de ces voisins.

Données

Les données proviennent de Natural Earth [1] :

ne_10m_admin_0_countries.shp (Polygones - Pays du monde) ➜ Téléchargement

Procédure

La méthode à utiliser peut être décrite de la manière suivante :

Pour chaque pays P:

  • trouver les polygones qui intersectent sa bounding box

  • parmis ces polygones :

    • trouver ceux qui sont réelement voisins du pays P

  • enregistrer le nom des voisins de P dans le champ NEIGHBORS

  • enregistrer la somme de la population de ces voisins dans le champ SUM

Cet exerice mobilise plusieurs notions nouvelles :

  • édition d’une couche

# Au début de l'édition :
layer.startEditing()

# À la fin de l'édition pour valider :
layer.commitChanges()

# Pour abandonner les changements en cours :
layer.rollBack()
  • création de nouveaux champs attributaires

from PyQt5.QtCore import QVariant

# La couche `layer` doit être en cours d'édition !
# Création d'un champ de chacun des 3 principaux types
# qui seront utilisés :
layer.dataProvider().addAttributes(
    [QgsField("Champ_1", QVariant.String),
     QgsField("Champ_2", QVariant.Int),
     QgsField("Champ_3", QVariant.Double)])
layer.updateFields()
  • mise à jour d’une entité après l’avoir modifiée

# La couche `layer` doit être en cours d'édition !
# `ft` est une variable de type QgsFeature
layer.updateFeature(ft)
  • utilisation de prédicats spatiaux (intersects, disjoint, etc.)

# Obtenir la géométrie d'une entité :
geom1 = feature1.geometry()
geom2 = feature2.geometry()

# Prédicats spatiaux :
geom1.intersects(geom2) # -> True/False
geom1.disjoint(geom2) # -> True/False
geom1.touches(geom2) # -> True/False
# Création de l'index :
index = QgsSpatialIndex()
for f in layer.getFeatures():
    index.insertFeature(f)

# Interrogation de l'index
# (l'argument de la méthode intersects de l'index
#  spatial doit être de type `QgsRectangle`)
intersecting_ids = index.intersects(geom.boundingBox())

# Une des façon d'obtenir les entités qui correspondent à ces ids :
request = QgsFeatureRequest().setFilterFids(intersecting_ids)
for feature in vector_layer.getFeatures(request):
    # ...

À vous de jouer ! Essayer d’écrire le script permettant d’atteindre l’objectif de cet exercice

La méthode à utiliser peut être décrite de la manière suivante :

Création de deux champs NEIGHBORS et SUM (attention au type !)

Création de l’index spatial

Parcours des entités de la couche, pour chaque pays P:

  • trouver les polygones qui intersectent sa bounding box (grace à l’index spatial)

  • parmis ces polygones :

    • trouver ceux qui sont réelement voisins du pays P

    • enregistrer le nom des voisins de P dans le champ NEIGHBORS

    • enregistrer la somme de la population de ces voisins dans le champ SUM

    • mettre à jour l’entité en question

Validation des changements effectués sur la couche.

➜ Quelle est la population des pays voisins de la Hongrie ?

➜ Combien de pays n’ont aucun voisin ?


💡 Solution de l'exercice 3


Footnotes