Tutoriel avancé QGIS et Python (PyQGIS)

Détails: 19 octobre 2024; 90005

Index de l'article

Mots-clefs : interopérabilité Python/QGIS, variable, boucle, fonction, condition, liste, argument, chaîne formatée, fonctions géométriques, création de fichiers, génération de cartes, fonctions cartographiques, fonctions standalone, traitement de données cartographiques, workflow, milieu alpin

Enseignant responsable : Georges Hinot

Pré-requis : une bonne connaissance d'un logiciel SIG supportant Python (QGIS ou ArcGIS par exemple) et des problématiques/besoins/enjeux de la géomatique.

Compétences visées : mettre en place un protocole de développement autour d'une problématique géomatique, de la « re-contextualisation du besoin » jusqu'à l'écriture d'un script Python répondant à ce besoin.

Enjeux du cours : rendre les étudiants autonomes dans la création et l'entretien de scripts « automatisants ». Dans ce type de tâche, l’autonomie consiste à savoir interpréter un besoin humain, ré-interpréter ses propres problématiques de programmation et à écrire/chercher/modifier les blocs de code nécessaires.

Programme pédagogique : initiation aux fondamentaux du langage puis mise en situation (générer des cartes à partir de données géographiques). En plus de la simple automatisation des cartographies, le cours nous donnera l'occasion, toujours en Python, d'utiliser des API, de traiter des données et de mettre en place des modes d'affichages customisés sur QGIS. Selon l'avancée du cours, nous pourrons aussi aborder le mode « standalone » permis par le duo QGIS-Python, la création de plugins QGIS et les possibilités de webmapping/analyse/ingénierie de données permises par Python.

Exemple de carte produite par Python

Versions :

Tutoriel mis à jour sur QGIS 3.28.13 et QGIS 3.34.11.
Tutoriel écrit à l'origine sur QGIS 3.14, QGIS 3.16 et QGIS 3.2, mais j'ai constaté qu'il fonctionne de façon assez similaire sur des versions précédentes (sans remonter trop loin non plus bien sûr).

Sources et liens divers :

Advanced PyQGIS tutorial
Tips on PyQGIS in standalone mode
Full code to generate maps
PyQGIS Developer Cookbook
Tutoriel d'initiation à Python dans QGIS3 (débutant)
QGIS StackExchange
Blog spécialisé d'Underdark, aka Anita Graser
OSMDownloader
QuickOSM
Wikidata
Landscape Archaeology
AppDividend
OpenSourceOptions
webgeodatavore

Découverte de la console Python dans QGIS

Ce 1^er chapitre est inspiré de l'excellent site qgistutorials.

Ouvrez les couches peaks et eau (jointes en bas de cet article) dans un projet QGIS, ainsi que la console Python du logiciel (un bouton à cet effet est disponible dans le menu QGIS).

La couche de points en premier plan, le linéaire en second, avec une symbologie évidente. Ce sont les principaux sommets des Écrins et le réseau hydrographique français*.

* Les sommets sont issus d'une rapide extraction depuis OSM, et le réseau hydrographique est une couche IGN peu à jour. Les sommets sont encodés en UTF-8, pour bien lire les accents, précisez-le éventuellement dans QGIS. Pour augmenter les difficultés, ces 2 couches sont dans des projections différentes.

Sélectionnez la couche de points dans le panneau calque de QGIS, puis dans la console Python entrez :

layer = iface.activeLayer()

Vous venez de créer une variable arbitrairement nommée layer et contenant la couche active (celle que vous avez sélectionnée dans QGIS).

iface est une variable globale, qui n'est autre que le canvas courant de votre projet QGIS, les données visibles dans votre projet. Vous n'avez pas créé cette variable vous-même, mais elle s'est automatiquement créée à l'ouverture de QGIS.

activeLayer est une méthode de Python, une action à réaliser, qui ici va chercher la couche active de votre projet. Nous reviendrons bientôt sur le point (.) une fois que les conditions de notre compréhension seront plus explicites.

Votre 1^ère variable existe donc déjà, et vous pouvez en connaître le contenu simplement puisque :

layer

Vous renverra quelques informations sur le contenu de cette variable :

<QgsMapLayer: 'peaks' (ogr)>

Vous venez de lire (d'un point de vue système) le contenu de la variable. Si vous sélectionnez maintenant la couche linéaire, et lancez les mêmes commandes, vous lirez bien-entendu un contenu différent pour la même variable. En effet elle contiendra autre chose, la méthode activeLayer allant chercher la couche sélectionnée dans QGIS. Bienvenue dans Python pour QGIS !

Astuce
Depuis la console, nul besoin de retaper les mêmes commandes déjà utilisées. Utilisez la flèche ↑ (flèche du haut) de votre clavier pour remonter la liste des commandes exécutées. Vous pouvez également tester de lancer vos commandes depuis l'éditeur Python de QGIS (survolez les boutons disponibles en haut de la console Python de QGIS, l'éditeur s'ouvrira à droite de la console). Une interprêtation plus rigoureuse du code est à l'œuvre depuis l'éditeur Python, avec un comportement parfois un peu différent. Dans le cas de nos 2 commandes ci-dessus par exemple, il faudra ajouter un print pour lire la couche visée.

Vous pouvez aussi interroger les quelques 500 méthodes disponibles pour cet objet :

dir(layer)

Ci-dessus, la variable layer est soumise en tant que paramètre (ce mot n'est pas anodin en programmation) de la fonction dir, entre parenthèse. La fonction sait donc sur quoi s'appliquer.

L'une des méthodes disponibles ets la méthode fields(), nous pouvons consulter toutes ces options sur la documentation officielle pyqgis :

https://qgis.org/pyqgis/3.0/core/other/QgsFields.html

Exemples :

Compter le nombre de champs

layer.fields().count()

Afficher les noms des champs

layer.fields().names()

Les listes

La méthode names() nous renvoie les champs entre crochets, c'est donc une liste (ce mot n'est pas anodin en programmation).

Nous pouvons accéder aux élements indicés des listes en ajoutant [NUMÉRO DE L'INDICE]. Pour afficher le 1^er champ par exemple :

layer.fields().names()[0]

Astuce
Les éléments d'une liste sont ordonnés de 0 à X, le 1^er élément ayant toujours le numéro 0. Comme en histoire, quand nous sommes en 50 après JC, nous sommes au 1^er siècle après JC. Et quand nous avons 20 ans, nous sommes dans notre 21^ème année.

Boucle et itération

L'une des méthodes disponibles pour nos couches s'appelle getFeatures. Nous allons nous en servir pour boucler sur les 341 entités contenues dans la couche de points :

layer = iface.activeLayer()
 
for f in layer.getFeatures():
    print(f['name'], f['osm_id'])

Le mot-clé for signifie que nous démarrons une boucle. Chaque contenu de cette boucle sera stocké dans une variable f (arbitrairement nommée ainsi, mais correspondant à une convention Python, entre autre). Le in regarde où lancer la boucle : sur les entités contenues dans notre variable.

Et qu'est-ce que nous y lancerons ? Des instructions itératives.

Astuce
Si à ce stade la boucle ne s'est pas exécutée, c'est que vous devez sortir de la boucle, en tapant sur la touche ↵ (Entrée).

Cheminement

Dans la 2^ème ligne notamment, remarquez l'usage du point (.), qui concatène variables et instructions Python (là ou PHP concatène les chaînes).

Ce n'est donc plus une simple concaténation, mais plutôt le signe d'un cheminement du code : Python voit votre variable layer, mais étant suivie d'un . il sait qu'il y a quelque chose à faire dessus. Quoi donc ? Y lancer la méthode qui suit.

Mais au fait, si le point sert à Python pour concaténer des instructions, comment concatène-t-il les chaînes de texte ? Avec le signe +.

Indentation

Observez l'indentation devant la 3^ème ligne. L'indentation est partie intégrante du langage Python, elle permet de structurer son code et plus encore, de signifier des choses dans le code.

Ici l'indentation précédant le print signifie qu'il est physiquement contenu dans la boucle.

L'indentation a donc son rôle à jouer. Avec Python une mauvaise indentation soulèvera une erreur d'exécution (cela contrairement à PHP, où l'indentation n'est là que pour aider à la lecture, au mieux).

Après saisie complète du code ci-dessus dans la console, pour exécuter la boucle, il vous faudra taper une dernière fois sur ↵ (Entrée) , pour sortir de la boucle. Les noms des principaux sommets de la zone de glaciers des Écrins, et leur identifiant OSM, défilent dans la fenêtre de résultats.

Astuce
Une erreur courante est parfois soulevée, suite à des copier-coller ou des modifications du code depuis un autre logiciel :

<span style="color: #ff0000; font-size: 10pt;">Inconsistent use of tabs and spaces in indentation.</span>

Le log est explicite : il s'agit d'un mélange de tabulations et d'espaces dans votre indentation. Visuellement votre code est propre, mais pas pour la machine, qui voit TOUT. En effet si vous avez commencé par utiliser des espaces ou des tabulations pour indenter, alors vous devez vous y tenir. De même si vous utilisez des indentation à 2 espaces, alors 2 indentations feront 4 espaces ; si vous faîtes des indentations à 4 espaces, alors 2 indentations feront 8 espaces. Souvenez-vous en pour la suite de ce tutoriel !

Corrigez donc la ligne concernée et prenez de bonnes habitudes : Notepad possède un paramètre pour signifier quel caractère systématiquement utiliser dans les indentations de la touche Tab du clavier.

1^ère fonction géométrique

De même que nous affichons les noms des sommets, nous pouvons les géocoder, en listant leurs coordonnées géographiques (xy) :

for f in layer.getFeatures():
    geom = f.geometry()
    print(geom.asPoint())

Astuce
Si à ce stade la boucle ne s'est pas exécutée, c'est que vous devez sortir de la boucle, en tapant sur la touche ↵ (Entrée).

Nous avons cette fois créé une variable directement dans notre boucle (f).

La méthode geometry récupère certaines des informations spatiales de nos entités, toutes stockées et bien organisées dans notre variable geom, même si nous ne les voyons pas.

Puis le print interroge notre variable sur l'un de ses contenus, asPoint, méthode qui récupére les latitudes et longitudes.

Note personnelle : extraire un réseau hydrographique plus complet depuis OSM, ainsi que les sommets tagués différemment (les pointes par exemple).

Écriture de fichiers

Écriture de fichier Text

Nous allons maintenant exporter les noms des sommets, leurs identifiants OSM ainsi que leur coordonnées géographiques dans un joli fichier texte.

Avant cela, récupérez l'emplacement de votre machine où vous souhaitez écrire votre fichier, car nous mentionnerons son chemin complet dans la 1^ère ligne ci-dessous, suivi du nom de fichier voulu.

Vous pouvez aussi écrire directement dans le répertoire courant (si vous le connaissez ! Voir plus bas), dans ce cas le seul nom du fichier suffira :

Astuce
Remarquez que nous utilisons de simples slashs pour définir nos chemins, contrairement à Windows, qui utilise des slahs inversés. SI vous récupérez vos chemins via les Propriétés de vos fichiers Windows, vous devez donc remplacer vos slahs inversés par de simples slahs, ou alors doubler vos slahs inversés (\\).

with open('C:/Users/georg/Downloads/SommetsEcrins.txt', 'w') as file:
    for f in layer.getFeatures():
        geom = f.geometry()
        line = '{};{};{:.4f};{:.4f}\n'.format(f['name'], f['osm_id'], geom.asPoint().y(), geom.asPoint().x())
        file.write(line)

Astuce
Si à ce stade la boucle ne s'est pas exécutée, c'est parce que ...

Chaînes formatées

Observez la syntaxe de la 4^ème ligne ci-dessus, un peu particulière.

Dans la variable line nous indiquons d'abord des zones vides ({}), ce sont les zones dans lesquelles nous souhaitons écrire les valeurs des champs. Nous en profitons parfois pour y indiquer un 1^er formatage. La notation :.4f signifie qu'ici les valeurs doivent être arrondies à 4 chiffres décimaux.

Nous utilisons également le caractère spécial \n issu des expressions régulières, pour faire en sorte qu'après chaque écriture d'entité, une ligne soit sautée dans le fichier avant d'écrire la suivante.

Nous y retrouvons donc 4 blocs avant le .format().

Ensuite seulement, après le .format(), nous y indiquons ce que les blocs doivent contenir, en respectant le même ordre (le 1^er champ mentionné sera contenu dans la 1^er bloc, et ainsi de suite).

Enfin, nous lançons l'écriture de notre variable line.

Il y a plusieurs façons de formater des chaînes avec Python. Ici nous utilisons la méthode {} mais nous en verrons d'autres.

Écriture de CSV

Pour écrire un beau fichier CSV nous allons procéder différemment, sachant par ailleurs qu'il existe plusieurs façons de faire.

out_file = open('C:/Users/georg/Downloads/SommetsEcrins.csv', 'w', encoding='cp1252')
out_file.write("name" + "," + "osm_id" + "," + "y" + "," + "x" + "\n")
 
for f in layer.getFeatures():
    geom = f.geometry()
    line = '{},{},{:.4f},{:.4f}\n'.format(f['name'], f['osm_id'], geom.asPoint().y(), geom.asPoint().x())
    out_file.write(line)
 
out_file.close()

Astuce
Avec le jeu des copier-coller depuis ce site web ou un éditeur de texte gérant différemment les indentations, séparateur Tab, whitespace et saut de ligne, il peut y avoir des erreurs d'interprêtation de Python.

Dans ce cas utilisez plutôt l'éditeur Python de QGIS (survolez les boutons disponibles en haut de la console Python de QGIS, l'éditeur s'ouvrira à droite de la console) et exécutez-y le code ci-dessus. Une meilleure interprêtation du code est à l'œuvre depuis l'éditeur Python.

La 1^ère ligne crée puis ouvre (d'un point de vue système) un fichier CSV vide. Nous ajoutons un encodage cp1252 sans doute nécessaire depuis les machines Windows.

La 2^nd ligne ajoute les noms de champs. Cela de façon manuelle, à ce niveau vous pouvez les renommer à votre guise.

Les suite des opérations est plus classique, hormis la dernière ligne qui vient fermer le fichier.

Écriture de Shape

Nous allons maintenant créer un shapefile. Pour que ce soit drôle, et ne pas créer un shape de shape, nous allons plutôt exporter une sélection des aiguilles de notre shape des sommets.

Cette fois inutile de sélectionner la couche dans le panneau des calques, nous utiliserons la méthode mapLayersByName qui ira chercher la couche concernée par son nom (1^ère ligne).

layer = QgsProject.instance().mapLayersByName('peaks')[0]
layer.selectByExpression('"NAME" LIKE \'Aiguille%\'')
 
root = r'C:/Users/georg/Downloads/Aiguilles/'
if not os.path.exists(root):
    os.makedirs(root)
 
file = str(root)+'Aiguilles.shp'
 
writer = QgsVectorFileWriter.writeAsVectorFormat\
(layer, file, 'utf-8', driverName='ESRI Shapefile', onlySelected=True)

La 2^nd ligne sélectionne nos aiguilles à la manière du langage SQL.

Des lignes 4 à 6 nous créons le dossier qui va contenir le shape. Avant cela nous vérifions qu'il n'existe pas déjà. De plus nous stockons son chemin complet dans une variable dédiée, car elle nous servira juste après à la génération du shape lui-même. Nous aurions pu procéder plus simplement, mais cette partie du code est optimisée pour éviter la redondance. Nous aurions même pu faire plus court en utilisant un chemin relatif, nous y reviendrons plus tard.

Retenez également le paramètre onlySelected=True dans la dernière ligne, qui vous rappelle sans doute des choses.

Observez le slash inversé (\) qui clôt la 10^ème ligne. Puisque les sauts de ligne sont interdits dans le code Python, certaines lignes peuvent être atrocement longue, en particulier les chaînes de texte par exemple. Le slah inversé permet de signifier à Python un saut de ligne sans impact dans le code.

Export d'images

Vous pouvez très rapidement exporter l'image du canvas courant, en format png par exemple :

iface.mapCanvas().saveAsImage("C:/Users/georg/Downloads/SommetsEcrins.png")

Remarquez que Pyhon-QGIS nous exporte également un fichier World File, correspondant notamment à l'étendue géographique de l'image exportée. Si cela ne vous intéresse pas, supprimez-le à la suite de votre export. Comme ici en format jpg cette fois :

iface.mapCanvas().saveAsImage("C:/Users/georg/Downloads/SommetsEcrins.jpg")
os.remove("../SommetsEcrins.jgw")

Remarquez l'usage d'un chemin relatif dans la suppression du fichier World File ci-dessus. Mon répertoire courant étant situé dans /Download, pour le remonter une fois j'utilise ../ qui me permet de viser directement le répertoire Download sans mentionner le chemin complet.

Une bonne maîtrise de la syntaxe de vos chemins vous permettra de raccourcir votre code.

Répertoire courant

En cas de doute vous pouvez connaître facilement votre répertoire courant en lançant le code ci-dessous.

path = os.getcwd()
print("Le répertoire courant est : " + path)

Changer de répertoire courant

path = os.getcwd()
print("Le répertoire courant est :\n" + path)
 
os.chdir('C:/Users/georg/Downloads/')
 
path = os.getcwd()
print("Le répertoire courant a changé et c'est maintenant :\n" + path)

Export de PDF

Pour exporter un PDF il faut d'abord avoir créé un layout, une carte. Créez-en un via le Layouts Manager de QGIS. Dans le code ci-dessous mon layout s'appelle Map1.

Pour connaître les layouts disponible dans un projet QGIS, une petite boucle fera l'affaire :

myLayouts = QgsProject.instance().layoutManager()
for layout in myLayouts.printLayouts():
    print(layout.name())

Pour exporter une carte en PDF :

manager = QgsProject.instance().layoutManager()
layout = manager.layoutByName("Map1")
exporter = QgsLayoutExporter(layout)
exporter.exportToPdf("C:/Users/georg/Downloads/MaJolieCarte1.pdf",QgsLayoutExporter.PdfExportSettings())

Sélection et fonction

Vous verrez que le mot fonction n'est pas anodin en programmation Python ! En revanche quand nous parlons de sélection, nous parlons d'une simple sélection géographique/attributaire.

Sélection sans fonction

Ouvrez les couches data_BDTOPO_V3_Dep05_adresse et chemin_de_fer (jointes en bas de cet article) dans un projet QGIS, ainsi que la console Python du logiciel. Pour plus de clarté, mettez la couche de points en 1^er plan et la couche linéaire en 2^nd plan.

Exécutez le code suivant, que nous allons ensuite décrypter :

mylayer = QgsProject.instance().mapLayersByName("ADRESSE_05")[0]
myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"CODE_POST" = \'05000\'' ) )
mylayer.selectByIds( [ f.id() for f in myselect ] )
iface.mapCanvas().zoomToSelected(mylayer)

La 1^ère ligne crée une variable mylayer, qui regarde le projet QGIS courant, notre instance, et y identifie la couche des points ADRESSE_05.

Le 2^nd ligne créé une variable myselect, qui à l'intérieur de la 1^ère variable contenant la couche spatiale mylayer, identifie des entités de cette couche en fonction d'une expression (CODE_POST = 05000, pour simplifier. Nous y reviendrons). Ce sont les adresses de la ville de Gap et des communes alentours (code postal 05000).

La 3^ème ligne sélectionne les entités concernées avec une boucle for. La 4^ème, pour notre confort personnel, zoome sur ces entités.

Nous reviendrons plus tard sur la notion de paramètre, mais remarquez que la méthode zoomToSelected prend déjà un paramètre (mylayer), qui est lui même notre 1^ère variable. Situé juste après la méthode et entre parenthèse, il mentionne que la méthode doit s'appliquer sur cette couche spatiale.

Les paramètres des méthodes, ou des fonctions, ne sont pas toujours obligatoires (essayez de le supprimer par exemple, en laissant les parenthèses vides, ré-exécutez le code). Mais dans la mesure où nous avons ici plusieurs couches d'ouvertes dans QGIS, il est nécessaire de préciser sur quelle couche le logiciel doit zoomer. Vous avez ici un 1^er exemple du caractère dynamique du langage Python (un autre paramètre, dans cette même fonction, aurait agit ailleurs. Nous reviendrons bientôt sur ce point).

Sélection avec fonction mais sans paramètre

Désélectionnez puis dézoomez. Nous allons maintenant faire la même chose mais en passant par une fonction. Comme son nom l'indique, cela devrait nous permettre d'avoir à notre disposition une vraie fonction, un outil en somme, que nous pourrons appeler et utiliser à notre guise. Exécutez le code suivant :

def SelectAndZoom():
    mylayer = QgsProject.instance().mapLayersByName("ADRESSE_05")[0]
    myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"CODE_POST" = \'05000\'' ) )
    mylayer.selectByIds( [ f.id() for f in myselect ] )
    iface.mapCanvas().zoomToSelected(mylayer)

Absolument rien ne s'est produit, c'est tout-à-fait normal.

Observez l'indentation. De toute évidence quelque chose est contenu à l'intérieur de ce bout de code.

Nous avons ici créé une fonction arbitrairement nommée SelectAndZoom (vous commencez à comprendre l'intérêt des fonctions...). Nous savons qu'il s'agit d'une fonction car elle est précédée d'un def, qui définit la présence d'une fonction. Nous l'avons créée, nous l'avons même enregistrée dans le cache de notre instance QGIS, mais nous ne l'avons pas exécutée.

Elle est tout de même là, prête à bondir à chaque instant !

Maintenant appelez-là :

SelectAndZoom()

Astuce
Avec les jeu des copier-coller depuis ce site web ou un éditeur de texte gérant différemment les indentations - séparateur tab, whitespace et saut de ligne - il peut y avoir des erreurs d'interprêtation/compilation de Python à l'appel de vos fonctions.

Dans ce cas utilisez plutôt l'éditeur Python de QGIS (survolez les boutons disponibles en haut de la console Python de QGIS, l'éditeur s'ouvrira à droite de la console).

Hop, nous voici de retour à Gap ! Amusez-vous à dézoomer, déselectionner, vous déplacer dans votre projet QGIS ou même ouvrir des tables attributaires... Puis ré-exécutez la fonction en l'appelant simplement (SelectAndZoom()). Elle est toujours là, et opérationnelle, agissant tant dans votre canvas que dans vos données attributaires.

Sélection avec fonction et un paramètre

Maintenant nous n'allons même plus prendre la peine de mentionner la couche cible dans la fonction.

À la place, nous mentionnerons la nécessité d'un paramètre dans le nom même de la fonction, entre parenthèse et arbitrairement nommé layername, et nous remplaçerons le nom de la couche cible par ce même paramètre :

def SelectAndZoom(layername):
    mylayer = QgsProject.instance().mapLayersByName(layername)[0]
    myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"CODE_POST" = \'05000\'' ) )
    mylayer.selectByIds( [ f.id() for f in myselect ] )
    iface.mapCanvas().zoomToSelected(mylayer)

La nouvelle fonction du même nom replace la 1^ère, et l'appel se fait maintenant avec un argument, le nom de la couche cible :

SelectAndZoom('ADRESSE_05')

Sélection avec fonction et plusieurs paramètres

Maintenant nous n'allons pas nous contenter d'appeler la fonction et la couche cible, mais carrément la couche cible, le champ et la valeur voulue. Attention nous allons devoir jouer avec la notion de chaîne formatée :

def SelectAndZoom(layername, field, value):
    mylayer = QgsProject.instance().mapLayersByName(layername)[0]
    myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"%s" = \'%s\'' % (field, value)) )
    mylayer.selectByIds( [ f.id() for f in myselect ] )
    iface.mapCanvas().zoomToSelected(mylayer)

L'appel requiert maintenant 3 arguments :

SelectAndZoom('ADRESSE_05', 'CODE_POST', '05100')

Chaînes formatées

Les chaînes formatées sont pénibles, hein ? Ainsi que les caractères d'échappement. C'est vrai.

u'"%s" = \'%s\'' % (field, value)

Mais qu'est-ce que c'est que ça ???

Pénible mais très performant, cela nous permet d'inclure du code hors le code. Ici Python attendait une chaîne. Mais nous avons voulu lui fournir une expression à la place, une expression qui aille chercher nos paramètres avant de les remplacer par des arguments (les chaînes finales). Il faut donc lui dire attention, ici ce n'est pas du simple texte, regarde d'abord les paramètres.

Il y a plusieurs façons de formater des chaînes avec Python. Ici nous utilisons la méthode %.

Le u précise que notre chaine sera en caractère Unicode.

Langage dynamique

Bien, vous utilisez maintenant pleinement le côté dynamique de Python : en appelant la même fonction mais avec des argument différents, vous opérerez des actions différentes. Essayez par exemple d'appeler une autre couche, avec un autre champ et une autre valeur :

SelectAndZoom('chemin_de_fer', 'NATURE', 'TGV')

Vous avez compris : un formulaire, autrement dit une interface graphique, pourrait fournir ces arguments suite à une saisie humaine. C'est comme ça que les logiciels et applications web fonctionnent !

Récupération d'informations d'entités déjà sélectionnées

Récupérez, dézippez et ouvrez la couche decathlon_france.

Sélectionnez quelques magasins au hasard, puis exécutez la code suivant.

layer_selected = QgsProject.instance().mapLayersByName("decathlon_france")[0]
selected_features = layer_selected.selectedFeatures()
 
for i in selected_features:
    attrs = i.__geo_interface__
    id_mag = i['id']
    print (attrs)
    print (id_mag)

Hop ! La console nous affiche certaines informations isues de notre sélection.

Remarquez la méthode attrs, qui récupère tous les attributs.

Sélection intelligente

Téléchargez les couches decathlon_france et iso_iris. Ce sont des points et leurs zones isochrones constituées de portions d'IRIS à moins 15 minutes en voiture de ces points.

Mettez une transparence complète sur les zones isochrones.

Maintenant que vous savez faire des sélections et enregistrer des fonctions dans la mémoire Python, voyons un exemple concret d'usage de fonction dans un SIG.

Imaginez que vous ayez une couche contenant des entités se superposant. Il est malaisé de les examiner sans faire de sélection attributaire, jouer sur la symbologie, les styles ou l'ordre d'affichage... Créons donc une fonction afin de les faire resortir suite à un événement (ici, un clic sur une entité, ou une simple sélection manuelle de quelques entités).

Mais certaines de ces zones se superposent, et les IRIS, tout ou partie, existent autant de fois qu'il appartiennent à une zone. Le code suivant va sélectionner les IRIS d'une zone isochrone après sélection d'un magasin (et donc faire apparaître la zone isochrone concernée).

layer_selected = QgsProject.instance().mapLayersByName("decathlon_france")[0]
layer_to_select = QgsProject.instance().mapLayersByName("iso_iris")[0]
 
def SelectionAuto():
    selected_features = layer_selected.selectedFeatures()
    for i in selected_features:
        attrs = i.__geo_interface__
        id_mag = i['id']
        #print (id_mag)
 
        myselect = layer_to_select.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"id_mag" = \'%s\'' % id_mag ) )
        layer_to_select.selectByIds( [ f.id() for f in myselect ] )
        #iface.mapCanvas().zoomToSelected(layer_to_select)
 
layer_selected.selectionChanged.connect(SelectionAuto)

Dans la suite de cet article, vous apprendrez à faire des symbologies en Python, vous pourrez donc enrichir ce type de fonction.

Note personnelle pour plus tard : enrichir la fonction d'une symbologie et d'un ordre d'affichage des entités sélectionnées. L'objectif étant de ne plus avoir à mettre de la transparence sur la couche pour rendre cette fonction utile.

Action enregistrée

Mais à ce stade votre sélection intelligente ne fonctionnera que tant que votre projet QGIS est ouvert. Fermez-le, et le code Python ne sera pas pris en compte à sa prochaine ouverture.

Pour enregistrer ce code dans votre projet, il faudra l'enregistrer au niveau des Actions d'une couche. Allez dans les propriétés d'une de vos couches (n'importe laquelle pour ce qui nous concerne, mais celle des magasins est la plus à propos), onglet Actions, choisissez le type Python (😄), nommez votre action, décrivez-la brièvement et collez-y le précédent code. Cochez également l'option Canvas Scope afin de la faire apparaître dans le menu de QGIS.

Votre code est désormais activable dans votre projet QGIS sous l'icône Actions, par un clic sur celle-ci, suivi d'un autre sur la couche stockant l'action (avec la petite croix qui apparaîtra). Sélectionnez maintenant un magasin, toute sa zone isochrone sera aussi sélectionnée.

Les 3 lois d'Asimov

Cependant si les machines n'ont a priori pas de sentiment, elles sont tout-de-même capricieuses, et des événements inopportuns peuvent se produire.

En effet votre fonction n'a de sens qu'après sélection d'un seul magasin, d'ailleurs elle ne sélectionne finalement que les entités liées à un seul magasin. Alors que se passe-t-il en cas de sélection de plusieurs magasins ? De plus si jamais vous sélectionnez tous vos magasins en même temps (Ctrl+A), QGIS ne lancera-t-il pas un important calcul, avec risque de plantage ? Faites quelques tests pour comprendre le problème.

Il serait donc bel et bon d'ajouter un comptage, une éventuelle désélection, ainsi que des conditions (nos 3 lois robotiques ! 😄) pour que la sélection ne se déclenche que si un seul magasin a été sélectionné, et désélectionne d'éventuelles entités précédentes. Ceci afin de ne pas induire l'utilisateur en erreur.

Utilisons les méthodes selectedFeatureCount() et removeSelection(), dans des conditions if.

layer_selected = QgsProject.instance().mapLayersByName("decathlon_france")[0]
layer_to_select = QgsProject.instance().mapLayersByName("iso_iris")[0]
 
def SelectionAuto():
    selected_features = layer_selected.selectedFeatures()
    for i in selected_features:
        if layer_selected.selectedFeatureCount() == 1:
            id_mag = i['id']
            myselect = layer_to_select.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"id_mag" = \'%s\'' % id_mag ) )
            layer_to_select.selectByIds( [ f.id() for f in myselect ] )
        if layer_selected.selectedFeatureCount() > 1:
            layer_to_select.removeSelection()
 
layer_selected.selectionChanged.connect(SelectionAuto)

Dans la ligne 7, une 1^ère condition if vérifie que le nombre de magasins sélectionnés est égal à 1, cela grâce à la fonction selectedFeatureCount(), qui renvoie le nombre d'entités sélectionnées. Si c'est bien le cas, le reste du contenu de cette condition s'exécute.

À la ligne 11, un autre if se contente de désélectionner (removeSelection()) la zone éventuellement sélectionnée auparavant, à condition qu'il y ait plus d'1 magasin sélectionné.

Note personnelle pour plus tard : faire la même chose mais sans fonction Def, dans une action activable au clic dans la table attributaire ou via l'icône action. Ceci afin de ne pas être confronté à une sélection des zones en permanence (impossibilité de désactiver cette action Def).

Variante pour afficher plusieurs zones isochrones après sélection de plusieurs magasins

layer_selected = QgsProject.instance().mapLayersByName("decathlon_france")[0]
layer_to_select = QgsProject.instance().mapLayersByName("iso_iris")[0]
 
def SelectionAuto():
    myList = []
    selected_features = layer_selected.selectedFeatures()
    for i in selected_features:
        id_mag = i['id']
        myList.append(id_mag)
 
    myList = str(myList).replace("[", "").replace("]", "")
    print(myList)
 
    myselect = layer_to_select.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"id_mag" IN (%s)' % myList ) )
    layer_to_select.selectByIds( [ f.id() for f in myselect ] )
 
layer_selected.selectionChanged.connect(SelectionAuto)

Générer une carte

Nous allons maintenant générer une carte PDF à partir du canvas QGIS. À chaque déplacement dans le canvas, une nouvelle exécution du script final mettra à jour la carte.

Environnement de travail

Ouvez les couches peaks, eau, sol, IRIS et troncons_routes, fournies en bas de cet article, dans un projet QGIS vierge, ainsi qu'un fond de carte OpenStreetMap.

Vérifiez la bonne supersposition des couches (les projections), puis appliquez une symbologie évidente.

Enregistrez votre projet puis décochez les couches eau, sol, IRIS et troncons_routes.

Ouvrez également le Layouts Manager. De même pour les besoins de la démonstration, faîtes-en sorte que vos écrans vous permettent de visualiser en même temps votre projet QGIS, console et éditeur Python ouverts, le Layouts Manager, et encore un peu de place pour un futur layout.

Zoomez sur la couche peak. Nous allons exécuter les bribes de code suivantes de façon itérative, en les décryptant une par une. Au terme de ce chapitre, vous aurez un code complet vous permettant de générer une carte simple. Sauvegardez petit-à-petit l'intégralité du code dans un fichier texte.

Layout

Nous générons un layout vide, que vous allez voir apparaître dans le Layouts Manager.

project = QgsProject.instance()
manager = project.layoutManager()
layoutName = "PrintLayout1"
 
# Vérification de la non-existence d'un layout de même nom
layouts_list = manager.printLayouts()
for layout in layouts_list:
    if layout.name() == layoutName:
        manager.removeLayout(layout)
 
# Génération d'un layout vide
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName(layoutName)
 
manager.addLayout(layout)

Encart cartographique

Nous tirons maintenant une carte vide dans la layout.

# Charger une carte vide
map = QgsLayoutItemMap(layout)
map.setRect(20, 20, 20, 20)
 
# Mettre un canvas basique
rectangle = QgsRectangle(1355502, -46398, 1734534, 137094)
map.setExtent(rectangle)
layout.addLayoutItem(map)

Canvas courant

Nous plaçons le canvas courant dans la carte.

# Mettre finalement le canvas courant
canvas = iface.mapCanvas()
map.setExtent(canvas.extent())
 
layout.addLayoutItem(map)
 
# Redimensionner la carte
map.attemptMove(QgsLayoutPoint(5, 27, QgsUnitTypes.LayoutMillimeters))
map.attemptResize(QgsLayoutSize(220, 178, QgsUnitTypes.LayoutMillimeters))
 
map.setFrameEnabled(True)

Légende dynamique

Ci-après le code d'affichage de la légende complète (une légende fixe, statique). Cette portion du code est commentée (encadré dans une balise """, et donc inactive).

En effet juste en dessous, nous allons préférer ne retenir que les couches apparaissant dans la carte (cochées dans le panneau des couches).

"""
# Mettre une légende complète
legend = QgsLayoutItemLegend(layout)
legend.setTitle("Légende")
layout.addLayoutItem(legend)
legend.attemptMove(QgsLayoutPoint(246, 5, QgsUnitTypes.LayoutMillimeters))
"""
 
# Légende personnalisée
tree_layers = project.layerTreeRoot().children()
checked_layers = [layer.name() for layer in tree_layers if layer.isVisible()]
print(f"Je vais ajouter uniquement {checked_layers} à la légende." )
 
layers_to_remove = [layer for layer in project.mapLayers().values() if layer.name() not in checked_layers]
print(f"Je vais retirer {layers_to_remove} de la légende." )
 
legend = QgsLayoutItemLegend(layout)
legend.setTitle("Légende")
 
layout.addLayoutItem(legend)
 
legend.attemptMove(QgsLayoutPoint(240, 24, QgsUnitTypes.LayoutMillimeters))
 
# Cette ligne permettra de ne pas sortir les couches inutilisées de votre panneau des calques QGIS, mais uniquement de la légende
legend.setAutoUpdateModel(False)
 
m = legend.model()
g = m.rootGroup()
for l in layers_to_remove:
    g.removeLayer(l)
 
legend.adjustBoxSize()

Jouez avec les print et la console python avec votre layout ouvert pour bien comprendre ce qui se passe.

Algorithmique

Pour la génération de notre légende dynamique ci-dessus, nous avons commencer par lister 3 variables.

Une contenant toutes les couches (tree_layers)
Une autre contenant uniquement les couches affichées (checked_layers)
Enfin une dernière ne contenant que les couches non-cochées (layers_to_remove), grâce à une 1^ère soustraction des couches cochées

En algorithmique :

layers_to_remove = tree_layers - checked_layers

Ensuite nous affichons la légende complète puis en retirons layers_to_remove dans une 2^nd soustraction.

Toujours en algorithmique :

my_legend = legend - layers_to_remove

Note personnelle : on pourrait aller plus loin dans cette légende dynamique en masquant de la légende les couches dont aucune entité n'apparaît dans le canvas courant (la couche adresse notamment).

Titre

# Titre
title = QgsLayoutItemLabel(layout)
title.setText("Ma Jolie Carte")
title.setFont(QFont("Verdana", 28))
title.adjustSizeToText()
 
layout.addLayoutItem(title)
 
title.attemptMove(QgsLayoutPoint(5, 4, QgsUnitTypes.LayoutMillimeters))

Tout comme avec la légende, vous connaissez maintenant la méthode attemptMove qui vous permet de positioner les éléments dans votre layout.

Les chiffres en argument sont dans l'ordre : la distance à partir du bord gauche du layout, puis la distance à partir du bord haut.

Sous-titre

# Sous-titre
subtitle = QgsLayoutItemLabel(layout)
subtitle.setText("Est-elle belle ?")
subtitle.setFont(QFont("Verdana", 17))
subtitle.adjustSizeToText()
 
layout.addLayoutItem(subtitle)
 
subtitle.attemptMove(QgsLayoutPoint(5, 19, QgsUnitTypes.LayoutMillimeters))

Échelle

# Échelle
scalebar = QgsLayoutItemScaleBar(layout)
scalebar.setStyle('Single Box')
scalebar.setUnits(QgsUnitTypes.DistanceKilometers)
scalebar.setNumberOfSegments(2)
scalebar.setNumberOfSegmentsLeft(0)
scalebar.setUnitsPerSegment(5)
scalebar.setLinkedMap(map)
scalebar.setUnitLabel('km')
scalebar.setFont(QFont('Verdana', 20))
scalebar.update()
 
layout.addLayoutItem(scalebar)
 
scalebar.attemptMove(QgsLayoutPoint(10, 185, QgsUnitTypes.LayoutMillimeters))

Logo

# Logo
Logo = QgsLayoutItemPicture(layout)
Logo.setPicturePath("https://master-geomatique.org/templates/theme3005/images/logo-ucp-cyu.png")
 
layout.addLayoutItem(Logo)
 
Logo.attemptResize(QgsLayoutSize(40, 15, QgsUnitTypes.LayoutMillimeters))
Logo.attemptMove(QgsLayoutPoint(250,4, QgsUnitTypes.LayoutMillimeters))

Vous utilisez ici la méthode attemptResize pour redimensionner l'élément. Les chiffres en argument sont dans l'ordre la largeur, puis la hauteur.

Flèche nord

Maintenant que vous savez mettre une image, ajoutez donc une flèche vers le nord correctement positionnée.

Texte

# Texte
TextCustom = QgsLayoutItemLabel(layout)
TextCustom.setText("Le massif des Écrins est un grand massif montagneux des Alpes françaises situé dans les Hautes-Alpes et en Isère.\
\n\nIl abrite d'importants glaciers, tant en nombre qu'en taille et possède deux sommets de plus de 4 000 mètres.\
\n\nIl était autrefois également nommé massif du Pelvoux.\
\n\nL'Oisans (bassin de la Romanche) au nord-ouest, le Champsaur (haut-bassin du Drac) au sud-ouest, et le Briançonnais (bassin de la Guisane) au nord-est recouvrent une partie du massif.")
TextCustom.setFont(QFont("Verdana", 11))
 
layout.addLayoutItem(TextCustom)
 
TextCustom.attemptResize(QgsLayoutSize(60, 120, QgsUnitTypes.LayoutMillimeters))
TextCustom.attemptMove(QgsLayoutPoint(230, 80, QgsUnitTypes.LayoutMillimeters))

Source

Maintenant que vous savez mettre un texte (un simple label finalement, comme les titre et sous-titre), ajoutez la source.

Note personnelle : on pourrait aller plus loin en allant chercher les sources directement dans les méta-données de nos couches.

Export

Vous connaissez déjà la manœuvre pour un export. Pensez à adapter le chemin utilisé ci-dessous.

# Export PDF
manager = QgsProject.instance().layoutManager()
layout = manager.layoutByName("PrintLayout1")
exporter = QgsLayoutExporter(layout)
exporter.exportToPdf("C:/Users/georg/Downloads/Ma Belle Carte.pdf",QgsLayoutExporter.PdfExportSettings())

En ouvrant le PDF dans votre navigateur, en vous déplaçant sur la carte, en cochant-décochant certaines couches puis en re-lançant l'intégralité du code ci-dessus, vous verrez votre PDF exporté être mis à jour en fonction.

Note personnelle : on pourrait aller plus loin en ajoutant une vraie flèche nord, dynamique, suivant l'orientation du canvas.

Et maintenant ?

Maintenant vous savez :

Faire une boucle
Faire une sélection
Faire une sélection itérative (sélectionner dans une boucle)
Exporter une carte

Et bien que se passerait-il si vous mettiez la sélection et l'export dans une boucle ? La production de cartes en série, pardi🤘.

À vous de jouer !

Import de couches

Nous allons maintenant travailler à partir d'un nouveau projet QGIS vierge, et importer nos couches, et tout le bastring nécessaire, directement en Python ! Yep, on est comme ça ici.

(Attention au &)

...
# Main path
monCheminDeBase = r'C:\\Users\\Georges\\Downloads\\'
 
project = QgsProject.instance()
 
# For OSM
urlWithParams = "type=xyz&amp;url=http://tile.openstreetmap.org/{z}/{x}/{y}.png"
osm = QgsRasterLayer(urlWithParams, "OpenStreetMap", "wms")
 
# Path of layers
peaks = QgsVectorLayer(monCheminDeBase + 'peaks_selection/peaks_selection.shp', 'Sommets', 'ogr')
glaciers = QgsVectorLayer(monCheminDeBase + 'glaciers/glaciers.shp', 'Glaciers', 'ogr')
eau = QgsVectorLayer(monCheminDeBase + 'eau/eau.shp', 'Fleuves', 'ogr')
troncons_routes = QgsVectorLayer(monCheminDeBase + 'troncons_routes/troncons_routes.shp', 'Grandes routes', 'ogr')
iris = QgsVectorLayer(monCheminDeBase + 'iris/iris.shp', 'IRIS', 'ogr')
 
# Open layers
project.addMapLayer(peaks)
project.addMapLayer(glaciers)
project.addMapLayer(eau)
project.addMapLayer(troncons_routes)
project.addMapLayer(iris)
 
project.addMapLayer(osm)
osm.setOpacity(0.75)
 
iface.mapCanvas().refresh()
 
...

Attention : puisque nous ouvrons OSM, qui utilise une projection Pseudo-Mercator (3857), alors en fonction des réglages de notre QGIS (projection à la volée) et des projections des autres couches utilisées (ici 2154, RGF 93, pour nos glaciers et fleuves, et 4326, WGS 84, pour les sommets), et de leur ordre d'ouverture, et des éventuelles instructions de rafraîchissement du canvas, alors on pourrait observer des bugs. EN particulier lors de la 1^ère ouverture de QGIS, à la 1^ère exécution du code et sur la 1^ère carte.

Nous allons donc définir et forcer la projection et l'ellipsoïde de notre projet :

project = QgsProject.instance()
...
crs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
project.setCrs(crs)
project.setEllipsoid('EPSG:4326')

Symbologie

Bon allez, retour aux choses sérieuses🤓!

Symbologie simple sur des points

Pour appliquer une symbologie simple sur nos points, les sommets enregistrés dans la variable my_layer :

my_layer = QgsProject.instance().mapLayersByName("Sommets")[0]
 
symbol_peak = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': 'green', 'outline_color': 'black', 'size': '4'})
mylayer.renderer().setSymbol(symbol_peak)
mylayer.triggerRepaint()

Symbologie simple sur des lignes

Cette fois on précise la couleur en HTML.

myroads = QgsProject.instance().mapLayersByName("Routes")[0]
myrivers = QgsProject.instance().mapLayersByName("Fleuves")[0]
 
# Symbologie simple de la couche route
symbol = QgsLineSymbol.createSimple({'line_style': 'dash', 'line_width': '0.5', 'color': 'black'})
myroads.renderer().setSymbol(symbol)
myroads.triggerRepaint()
 
# Symbologie simple de la couche eau
symbol = QgsLineSymbol.createSimple({'line_style': 'solid', 'line_width': '0.75', 'color': '#0088CC'})
myrivers.renderer().setSymbol(symbol)
myrivers.triggerRepaint()

Symbologie simple sur des polygones

Cette fois on précise la couleur en RGB (50,165,77,75), afin d'ajouter un peu de transparence (ici 75, l'indice allant de 0 à 255), et on rafraîchit le panneau des couches afin que notre légende y soit aussi prise en compte :

symbol = QgsFillSymbol.createSimple({'line_style': 'solid', 'line_width': '0.2', 'color': '50,165,77,75'})
myground.renderer().setSymbol(symbol)
myground.triggerRepaint()
iface.layerTreeView().refreshLayerSymbology(myground.id())

Il est aussi possible d'examiner les propriétés de symbologie d'un layer :

print(my_layer.renderer().symbol().symbolLayers()[0].properties())

Symbologie catégorisée

Pour une symbologie catégorisée, il va nous falloir définir plusieurs symbologies, puis les attribuer en fonction de valeurs attributaires.

Ici nous mettons tous nos sommets en triangle bleu (▲, puis nous allons chercher un sommet unique en se basant sur son nom (Aiguille de la Gandolière, dans le champ NAME). Celui-ci passera en triangle bleu foncé et de plus grande taille (▲).

symbol_peak = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': '#0088CC', 'outline_color': 'black', 'size': '4'})
symbol_peak_selected = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': 'blue', 'outline_color': 'black', 'size': '7'})
 
color_peak = QgsRendererCategory(None, symbol_peak, 'Sommets', True)
color_peak_selected = QgsRendererCategory('Aiguille de la Gandolière', symbol_peak_selected, 'Aiguille de la Gandolière', True)
 
renderer = QgsCategorizedSymbolRenderer('NAME', [color_peak,color_peak_selected])
 
mylayer.setRenderer(renderer)
mylayer.triggerRepaint()

Les 2 premières lignes définissent les symbologies que nous utiliserons, puis dans les lignes 4 et 5, la méthode QgsRendererCategory prend 4 arguments, dans l'ordre :

La valeur recherchée dans le champ pour l'attribution de la symbologie
La symbologie
Le label à afficher dans la légende pour cette catégorie
Un booléen, pour l'application ou non de cette symbologie

Puis c'est la ligne 7 qui vient exécuter les conditions, avec 2 arguments :

Le nom de champ où rechercher les valeurs
Un tableau contenant les conditions à appliquer

Symbologie dynamique

Nous n'avons plus qu'à utiliser une variable pour gérer la symbologie dans la boucle de génération de nos cartes.

La variable contenant le nom du sommet sélectionné existe déjà dans notre code, et a déjà été passé en texte (layoutName). C'est que nous utiliserons, à la fois dans l'argument value (pour filtrer le champ NAME sur cette valeur) et dans l'argument label (pour l'afficher dans la légende).

...
for feat in mylayer.getFeatures():
...
    peak_name = feat['NAME']
...
    layoutName = str(peak_name)
...
    # Symbologie catégorisée dynamique
    symbol_peak = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': '#0088CC', 'outline_color': 'black', 'size': '4'})
    symbol_peak_selected = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': 'blue', 'outline_color': 'black', 'size': '7'})
 
    color_peak = QgsRendererCategory(None,symbol_peak,"Autres sommets",True)
    color_peak_selected = QgsRendererCategory(layoutName,symbol_peak_selected,layoutName,True)
 
    renderer = QgsCategorizedSymbolRenderer("NAME", [color_peak,color_peak_selected])
 
    mylayer.setRenderer(renderer)
    mylayer.triggerRepaint()
 
    # Charger une carte vide
    map = QgsLayoutItemMap(layout)
    map.setRect(20, 20, 20, 20)
...

Logique de programmation

Pour fonctionner, la symbologie dynamique doit donc se trouver :

Après la création de la variable peak_name (sinon la symbologie dynamique ne saura pas quel sommet prendre, puisque peak_name devient ensuite layoutName, variable utilisée dans l'une des conditions de notre symbologie)
Après la création de la variable layoutName (sinon elle ne saura pas l'utiliser en tant que chaîne de texte)
Avant la mise à jour du canvas (ce point n'est pas réellement nécessaire, mais il correspond à une logique de programmation et de lecture/compréhension de notre code)
Avant la mise à jour de la légende (puisque la symbologie apparaît dans la legénde. Là encore il n'est pas certain que cela soit nécessaire, mais ce n'est que pure logique itérative, procédurale)

Votre code doit pouvoir être aisément lu et compris par vous-même ou quelqu'un d'autre. Un jour vous souhaiterez peut-être le ré-utiliser et l'enrichir. Respecter une logique de prommation et bien commenter votre code n'est pas anodin.

Note personnelle : gérer les sauts de ligne et la taille de la légende quand le nom du sommet selectionné est trop long (déborde de la carte).

OK, mais nous pouvons aller encore un peu plus loin dans les contenus de nos cartes.

Générer des cartes

Nous allons maintenant générer 341 cartes en PDF, chacune zoomant sur nos 341 sommets (mais pour avancer plus efficacement, utilisez plutôt l'extraction nommée peaks_selection, en pièce jointe de cet article, seulement une dizaine de sommets choisis, qui permettra d'exécuter nos scripts plus rapidement).

Ces cartographies automatiques sont donc en quelque sorte d'une reproduction de la fonction Atlas de QGIS, mais nous allons pouvoir aller beaucoup plus loin : appliquer des zooms différents par cartes, des symbologies conditionnelles différentes, des classifications, ajouter des contenus web, utiliser des APIs, lancer des géo-traitements, jointures, intersections... Bref, c'est pas trop mal.

Et pour les plus pressés🧐, un exemple de code complet est disponible ici.

Sélection itérative

Précédemment nous sélectionnions notre sommet, star de notre carte, sur la base de son identifiant. Mais maintenant que nous souhaitons itérer sur toute la couche peaks, nous ne pouvons plus nommément utiliser les identifiants (en fait si : nous pourrions les écrire un par un dans le code, dans un IN, mais ce n'est pas le but du jeu, une entité ajoutée dans les données ne serait pas prise en compte, et une entité supprimée ferait bugger le code).

Il nous faut donc boucler sur les entités sommitales.

Inspirez vous de la Sélection sans fonction plus haut dans ce court, pour créer une sélection d'un sommet sans connaître son identifiant.

Votre sélection des sommets se trouvera donc dans une boucle. Nous ne pourrons donc plus zoomer ou nous déplacer manuellement sur la couche pour générer les canvas de nos cartes.

Pour plus de clarté commencez par mettre votre expression dans une variable à part (5^ème ligne ci-dessous). Utilisez des print pour visualiser progressivement le contenu de vos variables.

mylayer = QgsProject.instance().mapLayersByName("peaks_selection")[0]
 
for feat in mylayer.getFeatures():
    id_peak = feat['OSM_ID']
    expr = u"OSM_ID = '{}'".format(id_peak)
    myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( expr ))
 
    mylayer.selectByIds( [ f.id() for f in myselect ] )
 
    iface.mapCanvas().zoomToSelected(mylayer)

Bien, vous êtes capable de lister les différents sommets à l'intérieur d'une boucle.

Profitons-en pour ajouter une variable contenant leurs noms (champ NAME) et affichons là dans un print.

Contrôlons également le niveau de zoom avec la méthode zoomScale().

mylayer = QgsProject.instance().mapLayersByName("peaks_selection")[0]
 
for feat in mylayer.getFeatures():
    id_peak = feat['OSM_ID']
    peak_name = feat['NAME']
 
    expr = u"OSM_ID = '{}'".format(id_peak)
    myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( expr ))
    print(f"Voici {peak_name}" )
 
    mylayer.selectByIds( [ f.id() for f in myselect ] )
 
    iface.mapCanvas().zoomToSelected(mylayer)
    iface.mapCanvas().zoomScale(23000)

Notez que le code ci-dessous se termine en sélectionnant la dernière entité listée, mais il a bel et bien listé, sélectionné puis zoomé une par une sur chaque entité. Même si cela ne s'est pas vu !

QGIS en a même gardé une trace dans son cache : si vous jouez avec le bouton Zoom Précédent du menu de QGIS, puis examinez son nom en comparant avec votre liste print, vous vous en apercevrez.

Utiliser les données dans la carte

Maintenant vous allez adaptez dans cette boucle le code complet du chapitre Générer une carte, afin de nommer votre layout, votre carte, lui mettre un titre, l'id OSM en guise de sous-titre, et un export correctement nommé également.

Vous allez devoir sélectionner des blocs de code et jouer avec le bouton Tab (indentation) et l'alliance Shift+Tab (retrait d'indentation) afin d'indenter correctement votre code dans la nouvelle boucle.

Testez votre code de façon itérative et procédurale :

D'abord la création des layouts, allez ensuite vérifier qu'ils existent bien.
Ensuite le zoom de chaque layout sur chaque sommet, vérifiez-en quelques uns.
...

Layouts, canvas et légende

Il vous faudra passer votre variable peak_name en chaîne de texte pour éviter des erreurs d'interprêtation de la variable layoutName. Utilisez la méthode str().

N'oubliez pas non plus de tout mettre dans la boucle ! Soit une indentation à ajouter au code ci-dessous.

...
layoutName = str(peak_name)
 
project = QgsProject.instance()
manager = project.layoutManager()
 
# Vérification de la non-existence d'un layout de même nom
layouts_list = manager.printLayouts()
for layout in layouts_list:
    if layout.name() == layoutName:
        manager.removeLayout(layout)
 
# Génération d'un layout vide
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName(layoutName)
manager.addLayout(layout)
 
# Charger une carte vide
map = QgsLayoutItemMap(layout)
map.setRect(20, 20, 20, 20)
 
# Mettre un canvas basique
rectangle = QgsRectangle(1355502, -46398, 1734534, 137094)
map.setExtent(rectangle)
layout.addLayoutItem(map)
 
# Mettre finalement le canvas courant
canvas = iface.mapCanvas()
map.setExtent(canvas.extent())
layout.addLayoutItem(map)
 
# Redimensionner la carte
map.attemptMove(QgsLayoutPoint(5, 27, QgsUnitTypes.LayoutMillimeters))
map.attemptResize(QgsLayoutSize(220, 178, QgsUnitTypes.LayoutMillimeters))
map.setFrameEnabled(True)
 
# Légende personnalisée
tree_layers = project.layerTreeRoot().children()
checked_layers = [layer.name() for layer in tree_layers if layer.isVisible()]
# print(f"Je vais ajouter uniquement {checked_layers} à la légende." )
 
layers_to_remove = [layer for layer in project.mapLayers().values() if layer.name() not in checked_layers]
# print(f"Je vais retirer {layers_to_remove} de la légende." )
 
legend = QgsLayoutItemLegend(layout)
legend.setTitle("Légende")
 
legend.setLinkedMap(map)
 
layout.addLayoutItem(legend)
legend.attemptMove(QgsLayoutPoint(240, 24, QgsUnitTypes.LayoutMillimeters))
 
# Cette ligne permettra de ne pas sortir les couches inutilisées de votre panneau des calques QGIS, mais uniquement de la légende
legend.setAutoUpdateModel(False)
 
m = legend.model()
g = m.rootGroup()
for l in layers_to_remove:
    g.removeLayer(l)
 
# Ajuster la légende
legend.adjustBoxSize()
 
# Forcer la mise à jour de la légende
legend.updateLegend()
iface.mapCanvas().refresh()
...

Titres et sous-titres

La variable peak_name contenant le nom des sommets, à utiliser dans le titre de nos cartes, à déjà été passée en texte, avec la méthode str(). Elle peut donc être utilisée telle quelle.

Mais ce n'est pas le cas de l'identifiant OSM, que nous voulons ajouter dans le sous-titre. Nous allons le faire mais en y laissant un peu de texte, grâce à une chaîne formatée.

...
# Titre
title = QgsLayoutItemLabel(layout)
title.setText(layoutName)
title.setFont(QFont("Verdana", 28))
title.adjustSizeToText()
layout.addLayoutItem(title)
title.attemptMove(QgsLayoutPoint(5, 4, QgsUnitTypes.LayoutMillimeters))
title.attemptResize(QgsLayoutSize(220, 20, QgsUnitTypes.LayoutMillimeters))
 
# Sous-titre
subtitle = QgsLayoutItemLabel(layout)
subtitle.setText("Identifiant OSM : %s" % (str(id_peak)))
subtitle.setFont(QFont("Verdana", 17))
subtitle.adjustSizeToText()
layout.addLayoutItem(subtitle)
subtitle.attemptMove(QgsLayoutPoint(5, 19, QgsUnitTypes.LayoutMillimeters))
subtitle.attemptResize(QgsLayoutSize(220, 20, QgsUnitTypes.LayoutMillimeters))
...

Échelle

L'échelle doit être légèrement adaptée, puisque nous avons modifié le niveau de zoom :

...
# Échelle
scalebar = QgsLayoutItemScaleBar(layout)
scalebar.setStyle('Single Box')
scalebar.setUnits(QgsUnitTypes.DistanceKilometers)
scalebar.setNumberOfSegments(2)
scalebar.setNumberOfSegmentsLeft(0)
scalebar.setUnitsPerSegment(1)
scalebar.setLinkedMap(map)
scalebar.setUnitLabel('km')
scalebar.setFont(QFont('Verdana', 20))
scalebar.update()
layout.addLayoutItem(scalebar)
scalebar.attemptMove(QgsLayoutPoint(10, 185, QgsUnitTypes.LayoutMillimeters))
...

Logo

Le logo ne change pas :

...
# Logo
Logo = QgsLayoutItemPicture(layout)
Logo.setPicturePath("https://master-geomatique.org/templates/theme3005/images/logo-ucp-cyu.png")
layout.addLayoutItem(Logo)
Logo.attemptResize(QgsLayoutSize(40, 15, QgsUnitTypes.LayoutMillimeters))
Logo.attemptMove(QgsLayoutPoint(250,4, QgsUnitTypes.LayoutMillimeters))
...

Texte

Le texte fera l'objet d'un chapite dédié, ne le mettez pas pour l'instant.

Le bug de la 1^ère carte

Dépendant du niveau de zoom mentionné dans votre code, mais aussi de la projection utilisée dans votre projet, à la 1ère ouverture de QGIS et à la 1^ère exécution du code, avec le code tel quel, la 1ère carte exportée, a minima, ne va pas prendre le niveau de zomm souhaité, oups...

En fonction de votre situation (puisque dépendant de beaucoup de choses, mdr), à la 2^nd exécution du code, ce bug va être corrigé.

Mouais, pas très pro cette histoire... Nous allons donc forcer la mise à jour du niveau de zoom, non pas seulement dans le canvas, mais dans l'objet carte (map), cela avec un map.setExtent et un map.setScale :

...
# Zoom on single peak in canvas
extent = sommet_unique.extent()
iface.mapCanvas().setExtent(extent)
iface.mapCanvas().zoomScale(20000)
iface.mapCanvas().refresh()
 
# Zoom on single peak in the map (to force zoom on the 1st map)
map.setExtent(iface.mapCanvas().extent())
map.setScale(100000)
 
# Layout
layout.addLayoutItem(map)
...

Export

Là encore, chaîne formatée pour correctement nommer nos carte grâce à la variable layoutName.

Attention : le répertoire d'accueil de vos cartes doit exister !

...
# Export PDF
manager = QgsProject.instance().layoutManager()
layout = manager.layoutByName(layoutName)
exporter = QgsLayoutExporter(layout)
exporter.exportToPdf("C:/Users/georg/Downloads/cartes/%s.pdf" % (layoutName),QgsLayoutExporter.PdfExportSettings())

Le bug des layouts

Bien, on a déjà de belles cartes. Mais si vous allez checker les mises en page dans le gestionnaire d'impressions, vous verrez que malgré le bon export des carte, chaque impression a pris la légende et la symbologie du dernier sommet, arf...

Ceci n'est pas un bug, mais tient au comportement des renderers (renderer()). En effet à chaque fois qu'on lance un renderer sur une couche, celui-ci s'applique naturellement à toute la couche, puis vient mettre à jour toutes les impressions utilisant cette couche, yep 🤣

Logiquement donc, à chaque itération de la boucle, la carte est bien exporté à partir de l'impression sur le sommet itéré, mais toutes les impressions sont aussi mises à jour sur le sommet itéré.

Cela est donc sans conséquence sur nos exports et nous resteront ainsi. Si cela était un vrai problème, en fonction du contexte, il est possible de gérer les symbologies autrement, avec des presets ou des themes de cartes.

Date du jour

Nous allons ajouter la date du jour dans les noms de nos fichiers.

Ouvrez la console de votre machine (CMD si vous êtes sous Windows) et entrez simplement:

python

La console passe en mode Python et affiche quelques informations de version. Maintenant ce code pour récupérer la date du jour :

import datetime
date = datetime.date.today()
str(date)

Le même code va fonctionner dans la console Python de QGIS. La 1^ère ligne est une importation d'un module natif, à inclure en haut de votre fichier.

Le date peut être stockée hors de la boucle, nous la stockerons dans une variable today.

import datetime
mylayer = QgsProject.instance().mapLayersByName("peaks_selection")[0]
 
date = datetime.date.today()
today = str(date)
 
for feat in mylayer.getFeatures():
    id_peak = feat['OSM_ID']
    peak_name = feat['NAME']
 
...

Avant de l'inclure dans le nom de nos exports dans la chaîne formatée.

...
exporter.exportToPdf\
("C:/Users/georg/Downloads/cartes/%s-%s.pdf"\
% (today, layoutName),QgsLayoutExporter.PdfExportSettings())

Bien, pour y voir plus clair, supprimez manuellement toutes vos cartes PDF existantes et lancez votre code.

Ça ne marche pas trop mal😎, mais nous pouvons encore aller plus loin, et sortir ainsi d'un simple atlas.

Intersection

Maintenant vous savez :

Exporter des cartes en série
Intersecter des couches géographiques

Et bien ajoutez donc dans le sous-titre de vos cartes la commune où se trouve le sommet sélectionné.

À vous de jouer !

Vous allez devoir utiliser la couche des IRIS (disponible en pièce jointe de cet article) et l'intersecter avec vos sommets :

...
import os
import shutil
...
monCheminDeBase = r'C:\\Users\\georg\\Downloads\\'
...
project = QgsProject.instance()
...
peaks = QgsVectorLayer(monCheminDeBase + 'peaks_selection/peaks_selection.shp', 'Sommets', 'ogr')
iris = QgsVectorLayer(monCheminDeBase + 'iris/iris.shp', 'IRIS', 'ogr')
...
QgsProject.instance().addMapLayer(peaks)
QgsProject.instance().addMapLayer(iris)
...
#Intersect peaks and IRIS
peaks_iris_path = _peaks_iris + r'\\peaks_iris.shp'
 
processing.run('qgis:intersection', { \
"INPUT": peaks, \
"OVERLAY": iris, \
"INPUT_FIELDS": ["OSM_ID", "NAME", "OTHER_TAGS"], \
"OVERLAY_FIELDS": ["CODE_IRIS", "NOM_COM"], \
"OVERLAY_FIELDS_PREFIX": "", \
"OUTPUT": peaks_iris_path})
 
# Remove layers
project.removeMapLayer(peaks)
project.removeMapLayer(iris)
 
# Open the new intersected layer
peaks_iris = QgsVectorLayer(peaks_iris_path, "Sommets", "ogr")
project.addMapLayer(peaks_iris)
 
# Register layer
mes_sommets = project.mapLayersByName("Sommets")[0]

Puis appeler le nouveau champ disponible (NOM_COM) et modifier la chaîne formatée qui remplit le sous-titre :

...
    commune = feat['NOM_COM']
    ...
    subtitle.setText("Identifiant OSM : {} - Commune : {}".format((str(id_peak)), (commune)))
...

Notez qu'au lieu de l'ancien %, nous utilisons ici une méthode de formatage plus moderne (.format())

Attention

Puisque nous créons un shape, et que nous l'ouvrons dans QGIS, si vous souhaitez exécuter votre code plusieurs fois, il faudra retirer votre couche de votre projet après votre boucle.

Idéalement il faudrait aussi créer/recréer (et donc supprimer) le shape avant la boucle, inspirez-vous du chapitre Processing: Intersection de cet article.

Installer un module tiers

Pour installer un module Python dans QGIS il vous faudra passer par le shell d'OSGeo4W. Dans cet article se trouve un chapitre nommé Install a third-party module qui vous donnera de plus amples indications.

Exemples de fichiers XML

Contenus web

Vous vous rappelez l'encart de texte au sujet du Massif des Écrins, accompagnant la carte que nous avons générée plus haut ? En dur dans le code, nous l'avons volontairement mis de côté au moment de générer nos cartes en série. En effet afficher toujours le même texte dans nos 341 cartes n'avait pas d'intérêt.

Et bien nous allons maintenant enrichir nos cartes avec des contenus dynamiques, en lien avec les données cartographiées, qui iront remplir cet encart de texte. Et nous irons chercher ces contenus sur le web.

Plusieurs méthodes existent, APIs et web-scraping notamment, souvent avec l'utilisation de fichier JSON ou XML.

Il nous faudra donc faire un lien avec nos données, et ce de façon dynamique bien sûr.

Dans ce chapitre, nous ferons nos tests directement depuis la console Windows, ouvrez donc la console Python dans l'interface de commande Windows (CMD + python).

OSM

Exemple de lien vers le fichier XML de la base de données OSM, suffixé par l'id OSM (que justement, nous possèdons pour tous nos sommets) :

https://www.openstreetmap.org/api/0.6/node/582066938

import requests
 
MyOsmId = '582066938'
r = requests.get("https://www.openstreetmap.org/api/0.6/node/"+MyOsmId)
 
print(r.text)

Notez l'importation nécessaire du module en haut de fichier.

Clé-valeur

Nos sommets n'ont pas de donnée attributaire sur leur altitude. Mais il semble que les XML d'OSM la connaissent ! En effet le fichier contient un tableau associatif avec un ensemble de clé-valeur dont ele (elevation en anglais).

Nous allons parser le fichier XML afin de récupérer l'altitude dans une jolie variable.

Parsing

Grâce à une condition if, nous allons récupérer l'altitude, à condition que le fichier XLM contienne le tag ele, sinon nous n'affichons rien.

import requests
import xml.etree.ElementTree
 
MyOsmId = '582066938'
r = requests.get("https://www.openstreetmap.org/api/0.6/node/"+MyOsmId)
root = xml.etree.ElementTree.fromstring(r.content)
 
for child in root.iter('tag'):
    if child.attrib['k'] == 'ele':
        print (child.attrib['v']+" mètres")

Nous sommes donc capable de récupérer l'altitude, bien, maintenant utilisons-là dans le titre de nos cartes générées :

import datetime
import requests
import xml.etree.ElementTree
 
mylayer = QgsProject.instance().mapLayersByName("peaks")[0]
...
for feat in mylayer.getFeatures():
...
    # Récupérer altitude sur API OSM
    r = requests.get("https://www.openstreetmap.org/api/0.6/node/"+id_peak)
    root = xml.etree.ElementTree.fromstring(r.content)
 
    for child in root.iter('tag'):
        if child.attrib['k'] == 'ele':
            altitude = ", "+child.attrib['v']+" mètres"
 
    # Titre
    title = QgsLayoutItemLabel(layout)
    title.setText(layoutName+altitude)
    title.setFont(QFont("Verdana", 28))
    title.adjustSizeToText()
    layout.addLayoutItem(title)
    title.attemptMove(QgsLayoutPoint(5, 4, QgsUnitTypes.LayoutMillimeters))
...

Précédé d'une virgule, notre variable altitude n'a plus qu'à aller se concaténer à la suite de notre variable layoutName dans le titre de nos cartes. Ainsi les enregistrements pour lesquels nous ne connaissons pas l'atitude n'afficheront pas de caractère inopportun à la suite du titre.

Note personnelle : on pourrait aller plus loin en ajoutant dans les sources des cartes le nom du contributeur OSM ayant saisi ou modifié l'entité sommitale cible.

Encodage

Il serait bon d'encoder notre caractère spécial è dans le mot mètres, car ici il est affiché en tant que pure chaîne, et dans certains cas spéciaux nous pourrions avoir de mauvaises suprises.

(Attention au &)

import html
...
            altitude = ", "+child.attrib['v']+html.unescape(" m&amp;egrave;tres")
...

Note personnelle : faire aussi l'encodage du fichier lui-même, et parler de l'encodage des fichiers Python.

Wikidata

Certains sommets possèdent un id Wikidata (dans le champ OTHER_TAGS). Exemple de lien vers donnée Wikidata, suffixé par un id Wikidata :

Pour tester, allez chercher le contenu d'une page page HTML d'osm.wikidata.link :

import requests
r = requests.get("https://osm.wikidata.link/Q726652")
print(r.text)

Un peu lent mais bien fourni. Les page wikidata.org se prêtent tout de même mieux au web-scraping (les contenus y sont mieux organisés, à l'intérieur de jolies balises).

Récupérer des données Wikidata

Les contenus y sont si bien organisé que nous allons pouvoir récupérer le lien vers la page Wikipédia du sommet, et même cibler celle en français.

Dans la console Python de Windows :

import requests
from requests_html import *
 
session = HTMLSession()
 
MyWikiId = 'Q726652'
r = session.get("https://www.wikidata.org/wiki/"+MyWikiId)
 
my_content = r.html.find('.wikibase-sitelinkview-link-frwiki a', first=True)
#print(my_content)
 
my_link = my_content.xpath('//a/@href')
print(my_link)

Nous nettoyons ensuite le lien pour le rendre utilisable par le module requests :

import requests
from requests_html import *
 
session = HTMLSession()
 
MyWikiId = 'Q726652'
r = session.get("https://www.wikidata.org/wiki/"+MyWikiId)
my_content = r.html.find('.wikibase-sitelinkview-link-frwiki a', first=True)
 
#print(my_content)
my_link = my_content.xpath('//a/@href')
 
#print(my_link)
 
my_link_string = str(my_link).replace('[', "").replace(']', "").replace("'", "")
print(my_link_string)
 
r2 = session.get(my_link_string)
 
print(r2.text)

Mais c'est finalement les pages Wikipédia dont les contenus ne sont plus suffisament bien organisés... Ce serait sans doute tout de même possible en multipliant les conditions, mais pas très efficace.

L'API Wikipedia

Essayons avec la bibliothèque wikipedia de Python, faîte expressément pour ça ! Dans la console Python de l'invite de commande Windows :

import wikipedia
 
wikipedia.set_lang("fr")
 
summary_fr = wikipedia.summary("Aiguille Dibona", sentences=10)
my_page = wikipedia.page("Aiguille Dibona")
 
print(summary_fr)
print (my_page.url)

Héhé, c'est bien plus simple hein ?

Vérification des résultats

OK, ajoutons également une comparaison in afin de vérifier la pertinence du résultat. En effet nous cherchons ici à partir de simples mots-clés dans le moteur de recherche de Wikipédia, et les résulats peuvent parfois être surprenants.

import wikipedia
 
wikipedia.set_lang("fr")
myWikiContent = wikipedia.summary("L'Ourson", sentences=3)
 
# Verifier la pertinence du resultat
if str('Massif des Écrins').lower() in str(myWikiContent).lower():
    pass
 
else:
    myWikiContent = 'Oooouuuuppssss !!!'
 
print(myWikiContent)

Ceci n'est qu'un court exemple, on peut faire mieux, comme dans cet exemple complet.

Condition

Puis nous mettons le texte dans une condition if else avant d'afficher le texte. Peut-être pas nécessaire mais évite d'encombrer la carte avec des blocs vides.

Le contenu du texte est maintenant une concaténation de nos nouvelles variables :

  ...
  #Texte
    if myWikiContent is None:
        None
    else:
        TextCustom = QgsLayoutItemLabel(layout)
        TextCustom.setText(myWikiContent+"\n\n"+myWikiLink)
        TextCustom.setFont(QFont("Verdana", 11))
        layout.addLayoutItem(TextCustom)
        TextCustom.attemptMove(QgsLayoutPoint(230, 100, QgsUnitTypes.LayoutMillimeters))
        TextCustom.attemptResize(QgsLayoutSize(60, 100, QgsUnitTypes.LayoutMillimeters))
    ...

Astuce
Un exemple complet d'usage de l'API Wikipedia ici !

L'API Wikidata

https://pypi.org/project/Wikidata/

from wikidata.client import Client
 
client = Client()
entity = client.get('Q1617977', load=True)
 
print(entity)
print(entity.description)
print(entity.data)
 
url_fr = entity.data['sitelinks']['frwiki']['url']
 
print(url_fr)
 
image_prop = client.get('P18')
image = entity[image_prop]
 
print(image.image_resolution)
print(image.image_url)
 
print(entity.__dict__)

Code complet

Code complet pour générer les cartes à partir d'un projet QGIS vide (modifiez les chemins, vérifié sur QGIS 3.14.16 (Pi), préférez l'éditeur pour l'exécution).

Autres exemples d'APIs

Ce chapitre est une parenthèse dans ce tutoriel, pour présenter quelques exemples d'exploitation d'APIs.

Un XML (http://mob.u-strasbg.fr/)

L'université de Strasbourg entretient une API fournissant des POIs, exemple :

http://mob.u-strasbg.fr/cgi-bin/odudsApi.py?method=search&id=allpois

La forme du XML généré est plus classique que ceux d'OSM :

<rss version="2.0">
  <poi>
    <id>1</id>
    <name>Central</name>
    <latitude>48.5803337</latitude>
    <longitude>7.7655637</longitude>
    <floor/>
    <type>campus</type>
    <description/>
    <url>http://fr.wikipedia.org/wiki/Campus_central_de_Strasbourg</url>
    <url_image/>
    <tags/>
    <capacity>0</capacity>
  </poi>
 
  <poi>
    <id>2</id>
...

Pour extraire les données de ce type d'XML :

(Attention au &)

import requests
import xml.etree.ElementTree
 
r = requests.get('http://mob.u-strasbg.fr/cgi-bin/odudsApi.py?method=search&amp;id=allpois')
root = xml.etree.ElementTree.fromstring(r.content)
for my_poi in root.findall('poi'):
    my_name = my_poi.find('name').text
    print(my_name)

Afficher les points à partir de leurs coordonnées

Vous pouvez vous entraîner avec le code simplifié de cet article, chapitre Temporary layer.

Le code suivant, dans un projet QGIS vide, affiche un fond de carte OSM, crée puis affiche une couche virtuelle de points, avec 2 champs, peuplée des points issus d'un XML. Toujours à partir de l'API de l'université de Strasbourg :

(Attention au &)

import requests
import xml.etree.ElementTree
from qgis.PyQt.QtCore import QVariant
 
project = QgsProject.instance()
 
project.removeAllMapLayers()
iface.mapCanvas().refresh()
 
crs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
project.setCrs(crs)
project.setEllipsoid('EPSG:4326')
 
urlWithParams = "type=xyz&amp;url=http://tile.openstreetmap.org/{z}/{x}/{y}.png"
osm = QgsRasterLayer(urlWithParams, "OpenStreetMap", "wms")
 
# point_vector
point_vector = QgsVectorLayer("Point", "my_points", "memory")
 
project.addMapLayer(osm)
project.addMapLayer(point_vector)
 
root = project.layerTreeRoot()
root.setHasCustomLayerOrder (True)
order = root.customLayerOrder()
order.insert(0, order.pop(order.index(point_vector)))
order.insert(1, order.pop(order.index(osm)))
root.setCustomLayerOrder(order)
 
pr = point_vector.dataProvider()
pr.addAttributes([QgsField("id", QVariant.Int), QgsField("name", QVariant.String)])
point_vector.updateFields()
 
r = requests.get('http://mob.u-strasbg.fr/cgi-bin/odudsApi.py?method=search&amp;id=allpois')
root = xml.etree.ElementTree.fromstring(r.content)
 
for my_poi in root.findall('poi'):
    my_id = int(my_poi.find('id').text)
    my_name = str(my_poi.find('name').text)
    my_longitude = float(my_poi.find('longitude').text)
    my_latitude = float(my_poi.find('latitude').text)
    #print(my_id, my_name, my_longitude, my_latitude)
 
    f = QgsFeature()
    f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(my_longitude, my_latitude)))
    f.setAttributes([my_id, my_name])
    pr.addFeature(f)
    point_vector.updateExtents()
    project.addMapLayer(point_vector)
 
iface.mapCanvas().refresh()

Aller chercher l'image sous condition

Dans la boucle, ajoutez :

    my_image = str(my_poi.find('url_image').text)
    if my_image == "None":
        my_image = ""
    elif my_image != "None":
        my_image = "https://mob.u-strasbg.fr/geoloc/url_image/" + my_image
        print(my_image)

Un JSON simple (data.culture.gouv.fr, monuments historiques)

Le site data.culture.gouv.fr propose plusieurs APIs, dont les données sont fournis en JSON.

Exemple avec la liste des immeubles parisiens protégés au titre des monuments historiques :

https://data.culture.gouv.fr/api/records/1.0/search/?dataset=liste-des-immeubles-proteges-au-titre-des-monuments-historiques&q=&facet=reg&facet=dpt_lettre&refine.dpt_lettre=Paris

Quand vous êtes confronté à ce type de fichiers, dont la visualisation en ligne n'est pas toujours aisée, il existe ce type d'outil en ligne qui vous permettra d'en comprendre la structure ou de tester votre fichier :

http://jsonviewer.stack.hu
https://geojson.io
http://geojson.io/#data=data:text/x-url,[URL]

Sur jsonviewer utilisez le bouton Load JSON data et copiez-y l'URL (supprimez le « s » du https de votre URL si besoin). Ou alors vous pouvez coller l'intégralité du contenu du fichier dans l'onglet Text et cliquer sur Format.

Sur geojson.io copiez vos données GEOJSON dans l'encart prévu, ou ajoutez l'adresse de votre fichier GEOJSON comme dans l'exemple ci-dessus (à la place de [URL]).

Pour parser ce type de JSON, et en extraire les données :

(Attention au &)

import json
import requests
 
response = requests.get("https://data.culture.gouv.fr/api/records/1.0/search/?dataset=liste-des-immeubles-proteges-au-titre-des-monuments-historiques&amp;q=&amp;facet=reg&amp;facet=dpt_lettre&amp;refine.dpt_lettre=Paris")
my_file = json.loads(response.text)
 
my_json = my_file["records"]
 
for my_block in my_json:
    my_fields = my_block["fields"]
    print(str(my_fields["adresse_forme_editoriale"]).replace("'", ''))

Les lignes 4 et 5 en font un JSON lisible pour Python.

Les données qui nous intéressent sont dans le bloc records, nous l'isolons donc avec la ligne 7.

La ligne 9 boucle sur le bloc records, et pour chaque sous-blocs fields affiche les données du tag adresse_forme_editoriale.

Les autres APIs de data.culture.gouv.fr

Si vous allez sur l'onglet Les données, et que vous cliquez sur presque n'importe quelle API, vous verrez que leur organisation est à peu près similaire.

Exemple avec l'API Répertoire des Musées de France : base Muséofile, son onglet API permet de modifier les paramètres de l'API. L'URL par défaut (en bas du formulaire) est la suivante :

https://data.culture.gouv.fr/api/explore/v2.1/catalog/datasets/musees-de-france-base-museofile/records?limit=20

Les affichages sont bridés à un maximum de 100 enregistrements.
On contrôle le nombre d'enregistrements affichés avec le paramètre limit.
On a un paramètre total_count pour connaître le nombre d'enregistrements.
D'autres champs sont filtrables avec le paramètre where.

On peut donc construire des URLs et itérer à travers les pages JSON des APIs pour récupérer les bases de données dynamiquement.

Autre exemples d'URLs :

https://data.culture.gouv.fr/api/explore/v2.1/catalog/datasets/musees-de-france-base-museofile/records?limit=20&offset=20

https://data.culture.gouv.fr/api/explore/v2.1/catalog/datasets/musees-de-france-base-museofile/records?where=region%3D%22Nouvelle-Aquitaine%22&limit=50

Une autre solution est de lancer le téléchargement du CSV complet à chaque fois que nécessaire, puis l'utiliser, mais c'est moins souple et plus lourd.

Un GEOJSON

Pour ouvrir un fichier GEOJSON, ce ne sera pas compliqué (source : webgeodatavore) :

iface.addVectorLayer('https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_populated_places.geojson', 'Populated places', 'ogr')

Mais pour ouvrir des données GEOJSON brutes (construites à la volée par une API par exemple, ou simplement non-accessibles par URL), alors vous devrez d'abord les enregistrer dans un fichier, puis ouvir ce fichier. Vous pouvez vous inspirer du chapitre nommé Zones isochrones de cet article, qui utilise l'API openrouteservice.org pour créer des GeoJson de zones isochrones.

Standalone

Alors ça c'est toute une histoire... Essayez cette page ci-dessous (accès restreint) ou contactez-moi en privé !

https://hg-map.fr/astuces/85-pyqgis

Buffer

En guise d'entrainement, commencez par créer une zone tampon de 5 mètres autour de nos sommets.

import os
import shutil
 
monCheminDeBase = r'C:\\Users\\georg\\Downloads\\'
 
project = QgsProject.instance()
project.removeAllMapLayers()
project.clear()
iface.mapCanvas().refresh()
 
crs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
project.setCrs(crs)
project.setEllipsoid('EPSG:4326')
 
peaks = QgsVectorLayer(monCheminDeBase + 'peaks_selection/peaks_selection.shp', 'Sommets', 'ogr')
QgsProject.instance().addMapLayer(peaks)
 
# Directory for created layer
_peaks_buffer = monCheminDeBase + '_peaks_buffer'
if os.path.isdir(_peaks_buffer) == True:
    shutil.rmtree(_peaks_buffer)
if os.path.isdir(_peaks_buffer) == False:
    os.mkdir(_peaks_buffer)
 
# Buffer on peaks
peaks_buffer_path = _peaks_buffer + r'\\peaks_buffer.shp'
 
processing.run('native:buffer', { \
"INPUT": peaks, \
"DISTANCE": 5, \
"SEGMENTS": 10, \
"END_CAP_STYLE": 0, \
"DISSOLVE": False, \
"OUTPUT": peaks_buffer_path})
 
# Open the buffer layer
peaks_buffer = QgsVectorLayer(peaks_buffer_path, "Buffer area", "ogr")
project.addMapLayer(peaks_buffer)
 
# Register layer
my_buffer = project.mapLayersByName("Buffer area")[0]
 
# Buffer below
root = project.layerTreeRoot()
myBelowLayer = root.findLayer(peaks_buffer.id())
myClone = myBelowLayer.clone()
parent = myBelowLayer.parent()
parent.insertChildNode(-1, myClone)
parent.removeChildNode(myBelowLayer)

Mais le code ci-dessus fonctionne très mal n'est-ce pas ?

Essayez de tricher
Essayez de comprendre pourquoi ça ne marche pas
Inspectez les unités géométriques utilisées dans votre canvas
Contrôlez vos logs dans View / Panels / Log Messages

Utilisez plutôt le code suivant pour créer des buffer de 2 kms autour de nos 108 000 adresses des Hautes-Alpes.

import os
import shutil
 
monCheminDeBase = r'C:\\Users\\georg\\Downloads\\'
 
project = QgsProject.instance()
project.removeAllMapLayers()
project.clear()
iface.mapCanvas().refresh()
 
crs = QgsCoordinateReferenceSystem.fromEpsgId(2154)
project.setCrs(crs)
project.setEllipsoid('EPSG:2154')
 
address = QgsVectorLayer(monCheminDeBase + 'adresses/ADRESSE_05.shp', 'Adresses', 'ogr')
 
QgsProject.instance().addMapLayer(address)
 
# Directory for created layer
_address_buffer = monCheminDeBase + '_address_buffer'
if os.path.isdir(_address_buffer) == True:
    shutil.rmtree(_address_buffer)
if os.path.isdir(_address_buffer) == False:
    os.mkdir(_address_buffer)
 
# Buffer on address
address_buffer_path = _address_buffer + r'\\address_buffer.shp'
 
processing.run('native:buffer', { \
"INPUT": address, \
"DISTANCE": 2000, \
"SEGMENTS": 10, \
"END_CAP_STYLE": 0, \
"DISSOLVE": False, \
"OUTPUT": address_buffer_path})
 
# Open the buffer layer
address_buffer = QgsVectorLayer(address_buffer_path, "Buffer area", "ogr")
project.addMapLayer(address_buffer)
 
# Register layer
my_buffer = project.mapLayersByName("Buffer area")[0]
 
# Buffer below
root = project.layerTreeRoot()
myBelowLayer = root.findLayer(address_buffer.id())
myClone = myBelowLayer.clone()
parent = myBelowLayer.parent()
parent.insertChildNode(-1, myClone)
parent.removeChildNode(myBelowLayer)

Joindre un CSV

import os
import shutil
 
monCheminDeBase = r'C:\\Users\\georg\\Downloads\\'
 
project = QgsProject.instance()
project.removeAllMapLayers()
project.clear()
iface.mapCanvas().refresh()
 
crs = QgsCoordinateReferenceSystem.fromEpsgId(4326)
project.setCrs(crs)
project.setEllipsoid('EPSG:4326')
 
countries = QgsVectorLayer(monCheminDeBase + 'simple_countries/simple_countries.shp', 'Countries', 'ogr')
 
project.addMapLayer(countries)
myCountries = project.mapLayersByName('Countries')[0]
 
extent_countries = countries.extent()
iface.mapCanvas().setExtent(extent_countries)
iface.mapCanvas().refresh()
 
# IMPORT CSV
csv_path = 'file:///' + monCheminDeBase + 'World Stats.csv?delimiter=;'
my_csv = QgsVectorLayer(csv_path, 'Countries', 'delimitedtext')
project.addMapLayer(my_csv)
 
# JOIN 
shpField='COUNTRY_HB'
csvField='Country'
 
myJoin = QgsVectorLayerJoinInfo()
myJoin.setJoinFieldName(csvField)
myJoin.setTargetFieldName(shpField)
myJoin.setJoinLayerId(my_csv.id())
myJoin.setUsingMemoryCache(True)
myJoin.setJoinLayer(my_csv)
myJoin.setJoinFieldNamesSubset(['Population (2024)'])
myJoin.setPrefix('')
 
countries.addJoin(myJoin)
countries.dataProvider().forceReload()

Classification

Intervalles égaux

# Set layer name and desired params
ramp_name = 'Spectral'
# ramp_name = 'PuRd'
 
value_field = 'Population (2024)'
num_classes = 10
classification_method = QgsClassificationEqualInterval()
 
layer = QgsProject().instance().mapLayersByName('Countries')[0]
 
format = QgsRendererRangeLabelFormat()
format.setTrimTrailingZeroes(True)
 
default_style = QgsStyle().defaultStyle()
 
color_ramp = default_style.colorRamp(ramp_name)
 
renderer = QgsGraduatedSymbolRenderer()
renderer.setClassAttribute(value_field)
renderer.setClassificationMethod(classification_method)
renderer.setLabelFormat(format)
renderer.updateClasses(layer, num_classes)
 
renderer.updateColorRamp(color_ramp)
 
layer.setRenderer(renderer)
layer.triggerRepaint()

Écriture de points

Dans un shape

Le code ci-dessous est inspiré de l'excellent site opensourceoptions.com.

En partant du principe que vous avez ouvert un shape nommé my_points, en 4326, contenant un champ id et un champ text, le code suivant crée un point sur l'Université de Cergy-Pontoise :

(Attention à &)

layers = QgsProject.instance().mapLayersByName('my_points')[0]
layer = QgsVectorLayer(layers.dataProvider().dataSourceUri(), '', 'ogr')
caps = layer.dataProvider().capabilities()
 
if caps &amp; QgsVectorDataProvider.AddFeatures:
    feat = QgsFeature(layer.fields())
    feat.setAttributes([1, 'Université de Cergy-Pontoise'])
    feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(2.0616833, 49.0297879)))
    res, outFeats = layer.dataProvider().addFeatures([feat])

Dans un shape, éditer seulement un champ

with edit(mylayer):
    for feature in mylayer.getFeatures():
        feature['myfield'] = 'Nouvelle valeur'
        mylayer.updateFeature(feature)

Dans une couche virtuelle ou un GEOJSON

Le code ci-dessous est inspiré de l'excellent site opensourceoptions.com.

Nous créons une couche virtuelle de points, y ajoutons deux champs id et name, avant d'y écrire un 1^er point :

point_vector = QgsVectorLayer("Point", "my_points", "memory")
 
QgsProject.instance().addMapLayer(point_vector)
 
from qgis.PyQt.QtCore import QVariant
pr = point_vector.dataProvider()
pr.addAttributes([QgsField("id", QVariant.Int), QgsField("name", QVariant.String)])
point_vector.updateFields()
 
f = QgsFeature()
f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(49.0297879, 2.0616833)))
f.setAttributes([1,'Université de Cergy-Pontoise'])
pr.addFeature(f)
point_vector.updateExtents()
QgsProject.instance().addMapLayer(point_vector)

Fusion de couches

Un cas d'école : un collaborateur récupère un grand nombre de fichiers de données et souhaite les exploiter proprement. Il s'agit de plusieurs centaines de fichiers cadastraux allemands.

Téléchargez, dézippez et ouvrez la couche Koln gml.zip.

Le but est d'automatiser :

L'ouverture d'un grand nombre de fichiers GML
Leur conversion en shapes
Leur édition de structure
Enfin, leur fusion en un seul shape.

import os
import shutil
from os import listdir
from os.path import isfile, join
 
for g in globals():
    del g
 
monRepertoire = 'C:/Users/georg/Downloads/Koln/lines/'
 
# On liste les GML en triant d'éventuels fichiers non GML (GFS par exemple)
mesGML_ = [f for f in listdir(monRepertoire) if isfile(join(monRepertoire, f))]
mesGML = [val for val in mesGML_ if val.endswith('.gml')]
 
QgsProject.instance().removeAllMapLayers()
iface.mapCanvas().refresh()
 
NbreGmlValid = 0
NbreGmlNoValid = 0
 
# On va stocker les noms des GML correctement ouverts mais aussi leurs appels en layers
GmlOuvertsNoms = []
LayersGmlOuverts = []
 
# Boucle pour ouvrir tous les GML sur QGIS
for gml in mesGML :
    mon_gml = QgsVectorLayer(monRepertoire+gml, gml, 'ogr')
 
    # On précise la projection
    crs = mon_gml.crs()
    crs.createFromId(25832)
    mon_gml.setCrs(crs)
    QgsProject.instance().addMapLayer(mon_gml)
 
    # On compte et liste les GML correctement ouverts mais aussi leurs appels en layers
    if mon_gml.isValid():
        print(gml, 'est valide ! OKLM.')
        NbreGmlValid += 1
        GmlOuvertsNoms.append(gml)
        LayersGmlOuverts.append(mon_gml)
 
    # On compte les GML en erreur
    if not mon_gml.isValid():
        print(gml, 'est invalide ! Nous le supprimons.')
        QgsProject.instance().removeMapLayer(mon_gml.id())
        NbreGmlNoValid += 1
 
print('Nombre de GML correctement ouverts :',NbreGmlValid)
print('Nombre de GML invalides et supprimés :',NbreGmlNoValid)
 
# EXPORT EN SHAPE
repertoireShapes = 'shapes'
pathShape = os.path.join(monRepertoire, repertoireShapes)
 
# On va stocker les noms des shapes correctement ouverts mais aussi leurs appels en layers
ShapesOuvertsNoms = []
LayersShapesOuverts = []
 
# Création d'un répertoire pour accueillir nos shapes
if not os.path.exists(pathShape):
    os.makedirs(pathShape)
else:
    shutil.rmtree(pathShape)
    os.makedirs(pathShape)
 
# Boucle pour exporter nos GML en shape
for s,n in zip(LayersGmlOuverts, GmlOuvertsNoms):
    nom = n.replace('.gml', '')
    writer = QgsVectorFileWriter.writeAsVectorFormat(s, pathShape+'/' + nom + '.shp', 'utf-8', driverName='ESRI Shapefile', onlySelected=False)
 
    mon_shape = QgsVectorLayer(pathShape+'/' + nom + '.shp', nom, 'ogr')
 
    ShapesOuvertsNoms.append(nom)
    LayersShapesOuverts.append(mon_shape)
 
# OUVRIR LES SHAPES
QgsProject.instance().removeAllMapLayers()
for x in LayersShapesOuverts:
    QgsProject.instance().addMapLayer(x)
 
# On va supprimer tous les champs de tous nos shapes sauf 2 champs
# Lister les champs d'un layer si besoin et tester avec my_vectorlayer = iface.activeLayer()
# D'abord on crée la liste de champs à supprimer
for my_shape in LayersShapesOuverts:
    prov = my_shape.dataProvider()
    field_names = [field.name() for field in prov.fields()]
    listeChampSuppression = []
 
    for a in field_names:
        listeChampSuppression.append(a)
 
    listeChampSuppression.remove('gml_id')
    listeChampSuppression.remove('measuredHe')
 
    # Puis on supprime la liste de champs à supprimer
    for c in listeChampSuppression:
        with edit(my_shape):
            my_field = my_shape.fields().indexFromName(c)
            my_shape.dataProvider().deleteAttributes([my_field])
        my_shape.updateFields()
 
# Fusion des shapes
parameters = {'LAYERS': ShapesOuvertsNoms,'CRS': 'EPSG:25832','OUTPUT': monRepertoire+'results.shp'}
processing.run("qgis:mergevectorlayers", parameters)
maFusion = QgsVectorLayer(monRepertoire+'results.shp', 'Ma jolie fusion', 'ogr')
QgsProject.instance().addMapLayer(maFusion)

Calculatrice de champs pour Python

Ci-dessous un exemple d'utilisation de cette calculatrice Python à la documentation hasardeuse.

Dans Global expression :

def myfunction(myparam):
    if myparam == 'Aiguille Dibona':
        return 2

Puis dans Formula:

Observez que le champ NAME est appelé en le mettant entre chevrons <>.

value = myfunction(<NAME>)

Liens ou pièces jointes
chemin_de_fer.zip	[ ]	0 Ko
data_BDTOPO_V3_Dep05_adresse.zip	[ ]	3889 Ko
data_IRIS_2019.zip	[ ]	45905 Ko
decathlon_france.zip	[308 magasins Décathlon français depuis OSM le 27 décembre 2020]	11 Ko
eau.zip	[ ]	0 Ko
glaciers.zip	[ ]	231 Ko
iso_iris.zip	[Des zones isochrones à 15 minutes autour de 308 POIs.]	12125 Ko
Koln gml.zip	[ ]	2818 Ko
peaks.zip	[ ]	14 Ko
peaks_selection.zip	[ ]	1 Ko
simple_countries.zip	[ ]	1880 Ko
sol.zip	[ ]	0 Ko
troncons_routes.zip	[ ]	0 Ko
World Stats	[ ]	27 Ko

Menu