Index de l'article

Qgis3PythonInteropérabilité Python/QGIS - Principaux concepts de Python, indentation, itération, variable, méthode, boucle, fonction, condition, tuple, paramètre, argument, chaîne formatée, importation - Fonctions géométriques - Écriture de fichiers - Fonctions SIG - Génération de cartes - Standalone - Milieu alpin

Sources et liens divers :

Tutoriel écrit sur QGIS 3.14 'Py' et Python 3.8, mais j'ai pu constater qu'il fonctionne de façon assez similaire sur des versions précédentes, ainsi que sur QGIS 3.16.

 


Découverte de la console Python dans QGIS

Ce 1er 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). Avec Python nul besoin de déclarer ses variables (contrairement à PHP par exemple, où c'est conseillé). 

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. Vous pouvez aussi interroger les quelques 500 méthodes disponibles pour cet objet :

dir(layer)

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.

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 :

  1. layer = iface.activeLayer()
  2. for f in layer.getFeatures():
  3. 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 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 :

Inconsistent use of tabs and spaces in indentation.

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).
  • Note personnelle : ici prévoir exposé sur les langages interprêtés/compilés et Python

É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érer l'emplacement du dossier 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, dans ce cas le seul nom du fichier suffira) :

  1. with open('C:/Users/Georges/Downloads/SommetsEcrins.txt', 'w') as file:
  2. for f in layer.getFeatures():
  3. geom = f.geometry()
  4. line = '{};{};{:.4f};{:.4f}\n'.format(f['name'], f['osm_id'], geom.asPoint().y(), geom.asPoint().x())
  5. 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, un peu particulière. Dans cette 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 1er 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. Ensuite seulement nous y indiquons ce que les zones vides doivent contenir, en respectant le même ordre (le 1er champ mentionné sera contenu dans la 1ère zone vide, 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.

  • Note personnelle : ici prévoir exposé sur les chaînes formatées en Python

É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.

  1. out_file = open('C:/Users/Georges/Downloads/SommetsEcrins.csv', 'w', encoding='cp1252')
  2. out_file.write("name" + "," + "osm_id" + "," + "y" + "," + "x" + "\n")
  3.  
  4. for f in layer.getFeatures():
  5. geom = f.geometry()
  6. line = '{},{},{:.4f},{:.4f}\n'.format(f['name'], f['osm_id'], geom.asPoint().y(), geom.asPoint().x())
  7. out_file.write(line)
  8.  
  9. 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. Si vous procédez ainsi, observez d'ailleurs d'évidents gains de performance par rapport à l'écriture du fichier texte. En effet il semble qu'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 2nd 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).

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

La 2nd 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/Georges/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/Georges/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)

Export de PDF

Pour exporter un PDF il faut d'abord avoir créer 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("../MaJolieCarte1.pdf",QgsLayoutExporter.PdfExportSettings())

Sélection et fonction

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 1er plan et la couche linéaire en 2nd plan.

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

  1. mylayer = QgsProject.instance().mapLayersByName("ADRESSE_05")[0]
  2. myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"CODE_POST" = \'05000\'' ) )
  3. mylayer.selectByIds( [ f.id() for f in myselect ] )
  4. 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 2nd 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 1er 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 !

  • Note personnelle : ici prévoir exposé sur les différents usages possibles de Python.

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

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)

Sélection intelligente

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).

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.

Mais certaines de ces zones se superposent, et les IRIS, tout ou partie, existent autant de fois qu'il appartiennent à une zone. Après avoir mis une transparence sur les zones isochrones, 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.

actionPour 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.

  1. layer_selected = QgsProject.instance().mapLayersByName("decathlon_france")[0]
  2. layer_to_select = QgsProject.instance().mapLayersByName("iso_iris")[0]
  3.  
  4. def SelectionAuto():
  5. selected_features = layer_selected.selectedFeatures()
  6. for i in selected_features:
  7. if layer_selected.selectedFeatureCount() == 1:
  8. id_mag = i['id']
  9. myselect = layer_to_select.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"id_mag" = \'%s\'' % id_mag ) )
  10. layer_to_select.selectByIds( [ f.id() for f in myselect ] )
  11. if layer_selected.selectedFeatureCount() > 1:
  12. layer_to_select.removeSelection()
  13.  
  14. 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).

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, eausolIRIS 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.

  1. project = QgsProject.instance()
  2. manager = project.layoutManager()
  3. layoutName = "PrintLayout1"
  4.  
  5. # Vérification de la non-existence d'un layout de même nom
  6. layouts_list = manager.printLayouts()
  7. for layout in layouts_list:
  8. if layout.name() == layoutName:
  9. manager.removeLayout(layout)
  10.  
  11. # Génération d'un layout vide
  12. layout = QgsPrintLayout(project)
  13. layout.initializeDefaults()
  14. layout.setName(layoutName)
  15.  
  16. manager.addLayout(layout)

Encart cartographique

Ouvrez le layout vide. Nous tirons maintenant une carte vide dans la layout.

  1. # Charger une carte vide
  2. map = QgsLayoutItemMap(layout)
  3. map.setRect(20, 20, 20, 20)
  4.  
  5. # Mettre un canvas basique
  6. rectangle = QgsRectangle(1355502, -46398, 1734534, 137094)
  7. map.setExtent(rectangle)
  8. layout.addLayoutItem(map)

Canvas courant

Nous plaçons le canvas courant dans la carte.

  1. # Mettre finalement le canvas courant
  2. canvas = iface.mapCanvas()
  3. map.setExtent(canvas.extent())
  4.  
  5. layout.addLayoutItem(map)
  6.  
  7. # Redimensionner la carte
  8. map.attemptMove(QgsLayoutPoint(5, 27, QgsUnitTypes.LayoutMillimeters))
  9. map.attemptResize(QgsLayoutSize(220, 178, QgsUnitTypes.LayoutMillimeters))

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).

  1. """
  2. # Mettre une légende complète
  3. legend = QgsLayoutItemLegend(layout)
  4. legend.setTitle("Légende")
  5. layout.addLayoutItem(legend)
  6. legend.attemptMove(QgsLayoutPoint(246, 5, QgsUnitTypes.LayoutMillimeters))
  7. """
  8.  
  9. # Légende personnalisée
  10. tree_layers = project.layerTreeRoot().children()
  11. checked_layers = [layer.name() for layer in tree_layers if layer.isVisible()]
  12. print(f"Je vais ajouter uniquement {checked_layers} à la légende." )
  13.  
  14. layers_to_remove = [layer for layer in project.mapLayers().values() if layer.name() not in checked_layers]
  15. print(f"Je vais retirer {layers_to_remove} de la légende." )
  16.  
  17. legend = QgsLayoutItemLegend(layout)
  18. legend.setTitle("Légende")
  19.  
  20. layout.addLayoutItem(legend)
  21.  
  22. legend.attemptMove(QgsLayoutPoint(240, 24, QgsUnitTypes.LayoutMillimeters))
  23.  
  24. # Cette ligne permettra de ne pas sortir les couches inutilisées de votre panneau des calques QGIS, mais uniquement de la légende
  25. legend.setAutoUpdateModel(False)
  26.  
  27. m = legend.model()
  28. g = m.rootGroup()
  29. for l in layers_to_remove:
  30. g.removeLayer(l)
  31.  
  32. 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 2nd 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

  1. # Titre
  2. title = QgsLayoutItemLabel(layout)
  3. title.setText("Ma Jolie Carte")
  4. title.setFont(QFont("Verdana", 28))
  5. title.adjustSizeToText()
  6.  
  7. layout.addLayoutItem(title)
  8.  
  9. 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

  1. # Sous-titre
  2. subtitle = QgsLayoutItemLabel(layout)
  3. subtitle.setText("Est-elle belle ?")
  4. subtitle.setFont(QFont("Verdana", 17))
  5. subtitle.adjustSizeToText()
  6.  
  7. layout.addLayoutItem(subtitle)
  8.  
  9. subtitle.attemptMove(QgsLayoutPoint(5, 19, QgsUnitTypes.LayoutMillimeters))

Échelle

  1. # Échelle
  2. scalebar = QgsLayoutItemScaleBar(layout)
  3. scalebar.setStyle('Single Box')
  4. scalebar.setUnits(QgsUnitTypes.DistanceKilometers)
  5. scalebar.setNumberOfSegments(2)
  6. scalebar.setNumberOfSegmentsLeft(0)
  7. scalebar.setUnitsPerSegment(5)
  8. scalebar.setLinkedMap(map)
  9. scalebar.setUnitLabel('km')
  10. scalebar.setFont(QFont('Verdana', 20))
  11. scalebar.update()
  12.  
  13. layout.addLayoutItem(scalebar)
  14.  
  15. scalebar.attemptMove(QgsLayoutPoint(10, 185, QgsUnitTypes.LayoutMillimeters))

Logo

  1. # Logo
  2. Logo = QgsLayoutItemPicture(layout)
  3. Logo.setPicturePath("https://master-geomatique.org/templates/theme3005/images/logo-ucp-cyu.png")
  4.  
  5. layout.addLayoutItem(Logo)
  6.  
  7. Logo.attemptResize(QgsLayoutSize(40, 15, QgsUnitTypes.LayoutMillimeters))
  8. 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.

Petit exercice

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

Vous pouvez faire autrement si vous le souhaitez.

Texte

  1. # Texte
  2. TextCustom = QgsLayoutItemLabel(layout)
  3. TextCustom.setText("Le massif des Écrins est un grand massif montagneux des Alpes françaises situé dans les Hautes-Alpes et en Isère.\
  4. \n\nIl abrite d'importants glaciers, tant en nombre qu'en taille et possède deux sommets de plus de 4 000 mètres.\
  5. \n\nIl était autrefois également nommé massif du Pelvoux.\
  6. \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.")
  7. TextCustom.setFont(QFont("Verdana", 11))
  8.  
  9. layout.addLayoutItem(TextCustom)
  10.  
  11. TextCustom.attemptResize(QgsLayoutSize(60, 120, QgsUnitTypes.LayoutMillimeters))
  12. TextCustom.attemptMove(QgsLayoutPoint(230, 80, QgsUnitTypes.LayoutMillimeters))

Petit exercice

Maintenant que vous savez mettre un texte (un simple label en réalité, 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.

  1. # Export PDF
  2. manager = QgsProject.instance().layoutManager()
  3. layout = manager.layoutByName("PrintLayout1")
  4. exporter = QgsLayoutExporter(layout)
  5. exporter.exportToPdf("C:/Users/Georges/Downloads/qgis/scripts/Ma Jolie 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.

Exercice

Maintenant vous savez :

  • Faire une boucle
  • Faire une sélection
  • 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 !


Générer des cartes

Nous allons maintenant générer 341 cartes en PDF, chacune zoomant sur nos 341 sommets.

Il s'agit 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, ajouter des contenus web, lancer des géo-traitements...

* Pour les plus pressés (🧐) un exemple de code complet est disponible ici.

Sélection en aveugle

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")[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().

  1. for feat in mylayer.getFeatures():
  2. id_peak = feat['OSM_ID']
  3. peak_name = feat['NAME']
  4. expr = u"OSM_ID = '{}'".format(id_peak)
  5. myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( expr ))
  6. print(f"Voici {peak_name}" )
  7.  
  8. mylayer.selectByIds( [ f.id() for f in myselect ] )
  9.  
  10. iface.mapCanvas().zoomToSelected(mylayer)
  11. 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 boutons QGIS dans le code Python

Petite parenthèse : notez qu'il existe une façon rigolote de zoomer sur les entités sélectionnées, en utilisant la liste des boutons QGIS :

  1. for feat in mylayer.getFeatures():
  2. id_peak = feat['OSM_ID']
  3. peak_name = feat['NAME']
  4. expr = u"OSM_ID = '{}'".format(id_peak)
  5. myselect = mylayer.getFeatures( QgsFeatureRequest().setFilterExpression ( expr ))
  6.  
  7. mylayer.selectByIds( [ f.id() for f in myselect ] )
  8.  
  9. eMenu = iface.viewMenu()
  10. eMenu.actions()[12].trigger()
  11.  
  12. iface.mapCanvas().zoomScale(23000)

Ci-dessus notre variable eMenu contient le menu QGIS, et la méthode actions()[12] exécute le douzième bouton (dans mon cas c'est le bouton Zoom sur la sélection).

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

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().

...
layoutName = str(peak_name)
...

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.

  1. ...
  2. # Titre
  3. title = QgsLayoutItemLabel(layout)
  4. title.setText(layoutName)
  5. title.setFont(QFont("Verdana", 28))
  6. title.adjustSizeToText()
  7. layout.addLayoutItem(title)
  8. title.attemptMove(QgsLayoutPoint(5, 4, QgsUnitTypes.LayoutMillimeters))
  9.  
  10. # Sous-titre
  11. subtitle = QgsLayoutItemLabel(layout)
  12. subtitle.setText("Identifiant OSM : %s" % (str(id_peak)))
  13. subtitle.setFont(QFont("Verdana", 17))
  14. subtitle.adjustSizeToText()
  15. layout.addLayoutItem(subtitle)
  16. subtitle.attemptMove(QgsLayoutPoint(5, 19, QgsUnitTypes.LayoutMillimeters))
  17. ...

Échelle

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

...
scalebar.setUnitsPerSegment(1)
...

Logo, légende et texte

Le logo ne change pas. La legende non plus pour l'instant (plus tard nous appliquerons une symbologie dynamique qui impactera la légende).

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

Export

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

  1. ...
  2. # Export PDF
  3. manager = QgsProject.instance().layoutManager()
  4. layout = manager.layoutByName(layoutName)
  5. exporter = QgsLayoutExporter(layout)
  6. exporter.exportToPdf("C:/Users/Georges/Downloads/qgis/scripts/exports/%s.pdf" % (layoutName),QgsLayoutExporter.PdfExportSettings())

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 :

  1. import datetime
  2. date = datetime.date.today()
  3. 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 dans la boucle, nous la stockerons dans une variable today.

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

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

...
exporter.exportToPdf\
("C:/Users/Georges/Downloads/qgis/scripts/exports/%s-%s.pdf"\
% (today, layoutName),QgsLayoutExporter.PdfExportSettings())

Bien, mais nous pouvons encore aller plus loin.


Symbologie dynamique

Symbologie unique

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

  1. symbol_peak = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': 'green', 'outline_color': 'black', 'size': '4'})
  2. mylayer.renderer().setSymbol(symbol_peak)
  3. mylayer.triggerRepaint()

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 ().

  1. symbol_peak = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': '#0088CC', 'outline_color': 'black', 'size': '4'})
  2. symbol_peak_selected = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': 'blue', 'outline_color': 'black', 'size': '7'})
  3.  
  4. color_peak = QgsRendererCategory(None,symbol_peak,"Sommet",True)
  5. color_peak_selected = QgsRendererCategory("Aiguille de la Gandolière",symbol_peak_selected,"Aiguille de la Gandolière",True)
  6.  
  7. renderer = QgsCategorizedSymbolRenderer("NAME", [color_peak,color_peak_selected])
  8.  
  9. mylayer.setRenderer(renderer)
  10. 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).

  1. ...
  2. for feat in mylayer.getFeatures():
  3. ...
  4. peak_name = feat['NAME']
  5. ...
  6. layoutName = str(peak_name)
  7. ...
  8. # Symbologie catégorisée dynamique
  9. symbol_peak = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': '#0088CC', 'outline_color': 'black', 'size': '4'})
  10. symbol_peak_selected = QgsMarkerSymbol.createSimple({'name': 'Triangle', 'color': 'blue', 'outline_color': 'black', 'size': '7'})
  11.  
  12. color_peak = QgsRendererCategory(None,symbol_peak,"Autres sommets",True)
  13. color_peak_selected = QgsRendererCategory(layoutName,symbol_peak_selected,layoutName,True)
  14.  
  15. renderer = QgsCategorizedSymbolRenderer("NAME", [color_peak,color_peak_selected])
  16.  
  17. mylayer.setRenderer(renderer)
  18. mylayer.triggerRepaint()
  19.  
  20. # Charger une carte vide
  21. map = QgsLayoutItemMap(layout)
  22. map.setRect(20, 20, 20, 20)
  23. ...

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)

 

  • Note personnelle : ici prévoir exposé sur les grandes logiques de programmation existantes.

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.

 


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.

  • Note personnelle : ici prévoir exposé sur les APIs

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).

Pip

Petite parenthèse : Pip est un installateur d'extensions pour Python, vous allez sans doute devoir vous y intéresser, et lui-même l'installer pour la suite de ce tutoriel. Je vous laisse faire.

Module requests

Si besoin d'installer requests (dans le shell Windows, pas Python) :

pip install requests

XML 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 sommet) :

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 :

  1. import datetime
  2. import requests
  3. import xml.etree.ElementTree
  4.  
  5. mylayer = QgsProject.instance().mapLayersByName("peaks")[0]
  6. ...
  7. for feat in mylayer.getFeatures():
  8. ...
  9. # Récupérer altitude sur API OSM
  10. r = requests.get("https://www.openstreetmap.org/api/0.6/node/"+id_peak)
  11. root = xml.etree.ElementTree.fromstring(r.content)
  12.  
  13. for child in root.iter('tag'):
  14. if child.attrib['k'] == 'ele':
  15. altitude = ", "+child.attrib['v']+" mètres"
  16.  
  17. # Titre
  18. title = QgsLayoutItemLabel(layout)
  19. title.setText(layoutName+altitude)
  20. title.setFont(QFont("Verdana", 28))
  21. title.adjustSizeToText()
  22. layout.addLayoutItem(title)
  23. title.attemptMove(QgsLayoutPoint(5, 4, QgsUnitTypes.LayoutMillimeters))
  24. ...

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.

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 Wikidatasuffixé par un id Wikidata :

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

  1. import requests
  2. r = requests.get("https://osm.wikidata.link/Q726652")
  3. 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).

Module requests-html

Si besoin, installez requests-html (dans le shell Windows, pas Python) :

pip install requests-html
Pour une installation à l'université de Cergy (Python 2) :
py -m pip install requests-html

Récupérer des données Wikipédia

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 :

  1. import requests
  2. from requests_html import *
  3.  
  4. session = HTMLSession()
  5.  
  6. MyWikiId = 'Q726652'
  7. r = session.get("https://www.wikidata.org/wiki/"+MyWikiId)
  8.  
  9. my_content = r.html.find('.wikibase-sitelinkview-link-frwiki a', first=True)
  10. #print(my_content)
  11.  
  12. my_link = my_content.xpath('//a/@href')
  13. print(my_link)

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

  1. import requests
  2. from requests_html import *
  3.  
  4. session = HTMLSession()
  5.  
  6. MyWikiId = 'Q726652'
  7. r = session.get("https://www.wikidata.org/wiki/"+MyWikiId)
  8. my_content = r.html.find('.wikibase-sitelinkview-link-frwiki a', first=True)
  9.  
  10. #print(my_content)
  11. my_link = my_content.xpath('//a/@href')
  12.  
  13. #print(my_link)
  14.  
  15. my_link_string = str(my_link).replace('[', "").replace(']', "").replace("'", "")
  16. print(my_link_string)
  17.  
  18. r2 = session.get(my_link_string)
  19.  
  20. 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.

La bibliothèque (library) wikipedia de Python

Essayons avec la bibliothèque wikipedia de Python, faîte expressément pour ça !

Si besoin de l'installer (dans le shell Windows, et non Python) :

pip install wikipedia
Pour une installation à l'université de Cergy (Python 2) :
py -m pip install wikipedia

Et maintenant dans la console Python de l'invite de commande Windows :

  1. import wikipedia
  2.  
  3. wikipedia.set_lang("fr")
  4. w = wikipedia.summary("Aiguille Dibona", sentences=3)
  5.  
  6. my_page = wikipedia.page("Aiguille Dibona")
  7.  
  8. #print (my_page.url)
  9. print (w)

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

Installation de modules Python dans QGIS

Mais nous avons un problème. En effet la bibliothèque wikipedia de Python n'est pas forcément accessible depuis QGIS, même si vous l'avez déjà installé sur Windows.

Testez la seule importation du module wikipedia depuis la console Python de QGIS pour en avoir le cœur net.

import wikipedia

Cela tient à ce que la version de Python utilisée par QGIS n'est pas exactement la même, puisque sous Windows QGIS embarque sa propre installation de Python.

Pour installer un module tierce, il vous faudra passer par le shell d'OSGeo4W. Suivez ce tutoriel pour installer un module sur la version de Python utilisée par votre QGIS sous Windows.

Dans le shell d'OSGeo4W :

py3_env
python -m pip install wikipedia

Pour une installation à l'université de Cergy (droits administrateur)
python -m pip install wikipedia --user

Une fois la bibliotheque wikipedia installée, les codes précédents exécuté depuis QGIS fonctionneront également.

Bien, nous n'avons plus qu'à l'adapter à notre code de génération de cartes.

Ignorer les exceptions

Nous allons ajouter un mot-clé à notre recherche, pour éviter les contenus sans rapport à certains noms de sommet (recherche textuelle oblige), mais surtout nous mettons la recherche dans un try, afin d'ignorer les erreurs quand la recherche ne trouve rien.

Nous prendrons les 1ères phrases de la page (sentences) ainsi que le lien URL lui-même.

Nous 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.

  1. import datetime
  2. import requests
  3. import xml.etree.ElementTree
  4. import wikipedia
  5.  
  6. mylayer = QgsProject.instance().mapLayersByName("peaks")[0]
  7. ...
  8. for feat in mylayer.getFeatures():
  9. ...
  10. peak_name = feat['NAME']
  11. ...
  12. layoutName = str(peak_name)
  13. ...
  14. # Récupérer le résumé Wikipédia
  15. try:
  16. wikipedia.set_lang("fr")
  17. myWikiContent = wikipedia.summary("massif des Écrins "+layoutName, sentences=2)
  18.  
  19. my_page = wikipedia.page(layoutName)
  20. myWikiLink = my_page.url
  21.  
  22. # Verifier la pertinence du resultat
  23. if layoutName.lower() in str(myWikiContent).lower():
  24. None
  25. else:
  26. myWikiContent = ""
  27. myWikiLink = ""
  28.  
  29. except Exception:
  30. pass
  31. ...
  • Note personnelle : la recherche et la comparaison peuvent être améliorées, elle semble ne pas marcher dans certains cas (sommet de l'Ourson par exemple...).

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 :

  1. ...
  2. #Texte
  3. if myWikiContent is None:
  4. None
  5. else:
  6. TextCustom = QgsLayoutItemLabel(layout)
  7. TextCustom.setText(myWikiContent+"\n\n"+myWikiLink)
  8. TextCustom.setFont(QFont("Verdana", 11))
  9. layout.addLayoutItem(TextCustom)
  10. TextCustom.attemptMove(QgsLayoutPoint(230, 100, QgsUnitTypes.LayoutMillimeters))
  11. TextCustom.attemptResize(QgsLayoutSize(60, 100, QgsUnitTypes.LayoutMillimeters))
  12. ...

Code complet pour générer les cartes à partir d'un projet QGIS vide (modifiez les chemins, fonctionne dans l'éditeur Python exclusivement)

  • Note personnelle : ajouter une image d'illustration dans les cartes depuis les données Wikidata.

 


Autres exemples d'APIs

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

Un XML classique

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

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 :

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)

Un JSON

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 :

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 :

Sur jsonviewer utilisez le bouton Load JSON data et copiez-y l'URL (supprimez le « s » du https de votre URL si besoin). 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 :

  1. import json
  2. import requests
  3.  
  4. 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")
  5. my_file = json.loads(response.text)
  6.  
  7. my_json = my_file["records"]
  8.  
  9. for my_block in my_json:
  10. my_fields = my_block["fields"]
  11. print(my_fields["wadrs"])

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 wadrs.

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.

Exemple avec l'excellente API de zones isochrones d'openrouteservice (votre obtiendrez une clé gratuite après une rapide inscription) :

import requests
 
body = {"locations":[[8.681495,49.41461],[8.686507,49.41943]],"range":[300,200]}
headers = {
    'Accept': 'application/json, application/geo+json, application/gpx+xml, img/png; charset=utf-8',
    'Authorization': 'VOTRE CLE API',
    'Content-Type': 'application/geo+json; charset=utf-8'
}
call = requests.post('https://api.openrouteservice.org/v2/isochrones/driving-car', json=body, headers=headers)
#print(call.status_code, call.reason)
#print(call.text)
 
with open('C:/Users/Georges/Downloads/myfile.geoson', 'w', encoding='utf-8') as outfile:
    outfile.write(call.text)
 
iface.addVectorLayer('C:/Users/Georges/Downloads/myfile.geoson', 'Isochrones', 'ogr')

Afficher des points à partir de leurs coordonnées issues d'une API

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 :

  1. import requests
  2. import xml.etree.ElementTree
  3.  
  4. # VIDER PROJET
  5. QgsProject.instance().removeAllMapLayers()
  6. iface.mapCanvas().refresh()
  7.  
  8. # PARAMETRES OSM
  9. urlWithParams = "type=xyz&amp;url=http://tile.openstreetmap.org/{z}/{x}/{y}.png"
  10. osm = QgsRasterLayer(urlWithParams, "OpenStreetMap", "wms")
  11.  
  12. # CREER UNE COUCHE VIRTUELLE DE POINTS
  13. point_vector = QgsVectorLayer("Point", "my_points", "memory")
  14.  
  15. # AFFICHER LES COUCHES
  16. QgsProject.instance().addMapLayer(osm)
  17. QgsProject.instance().addMapLayer(point_vector)
  18.  
  19. # AJOUTER CHAMPS DANS point_vector
  20. from qgis.PyQt.QtCore import QVariant
  21. pr = point_vector.dataProvider()
  22. pr.addAttributes([QgsField("id", QVariant.Int), QgsField("name", QVariant.String)])
  23. point_vector.updateFields()
  24.  
  25. # PARSER LE FICHIER
  26. r = requests.get('http://mob.u-strasbg.fr/cgi-bin/odudsApi.py?method=search&amp;id=allpois')
  27. root = xml.etree.ElementTree.fromstring(r.content)
  28.  
  29. for my_poi in root.findall('poi'):
  30. my_id = int(my_poi.find('id').text)
  31. my_name = str(my_poi.find('name').text)
  32. my_longitude = float(my_poi.find('longitude').text)
  33. my_latitude = float(my_poi.find('latitude').text)
  34. #print(my_id, my_name, my_longitude, my_latitude)
  35.  
  36. # EDITER COUCHE VIRTUELLE
  37. f = QgsFeature()
  38. f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(my_longitude, my_latitude)))
  39. f.setAttributes([my_id, my_name])
  40. pr.addFeature(f)
  41. point_vector.updateExtents()
  42. QgsProject.instance().addMapLayer(point_vector)
  43.  
  44. # RAFRAICHIR
  45. iface.mapCanvas().refresh()
  46.  
  47. # CRS
  48. my_crs=QgsCoordinateReferenceSystem(3857)
  49. QgsProject.instance().setCrs(my_crs)

Exercice

Testez le lien Wikipédia que nous avons mis dans nos cartes. Ils ne fonctionnent pas toujours. en effet nos PDF interprêtent ces liens en tant que chaînes de texte, et les sauts de ligne viennent parfois casser le lien.

Vous avez déjà eu des cours d'HTML l'année dernière, et bien sûr vous vous passionnez pour les solutions web.

Et bien à la place d'afficher le lien complet, qui ne marche même pas, mettez-nous un joli nom de lien fonctionnel !

  • Note personnelle : mettre le code d'affichage propre du lien en correction.

Standalone

Ce chapitre est une parenthèse dans ce tutoriel pour découvrir les possibilités standalone d'utilisation de QGIS grâce à Python.

Depuis la console Python accessible depuis python-qgis.bat (sous Windows : dans OSGEO4W / bin), vous pouvez tester vos scripts en mode standalone.

Le but est de ne plus ouvrir QGIS pour réaliser nos traitements.

La console Python de votre machine Windows, ou même celle du shell d'OSGeo4W, sera insuffisante. Testez par exemple cette importation :

from qgis.core import *

Elle ne fonctionne pas, car la version Python de votre machine ne sait pas où trouver QGIS. Nous y reviendrons.

Ouvrir des shapes

Auparavant nous avons travailler sur des projets QGIS avec nos shapes préalablement ouverts. En mode standalone il vous faudra aller les chercher directement, et eventuellement gérer les projections ainsi que leur ordre d'affichage.

Vérifions déjà que nous sommes bien capables de le faire depuis QGIS. Observez l'URL d'appel d'OpenStreetMap, vous pouvez la retrouver dans ses propriétés depuis QGIS.

  1. QgsProject.instance().setCrs(QgsCoordinateReferenceSystem(4326))
  2.  
  3. urlWithParams = "type=xyz&amp;zmin=0&amp;zmax=19&amp;url=http://tile.openstreetmap.org/{z}/{x}/{y}.png"
  4.  
  5. rlayer = QgsRasterLayer(urlWithParams, "OpenStreetMap", "wms")
  6. troncons_routes = QgsVectorLayer("C:/Users/Georges/Downloads/qgis/troncons_routes/troncons_routes.shp", "Routes", "ogr")
  7. sommets = QgsVectorLayer("C:/Users/Georges/Downloads/qgis/peaks/peaks.shp", "Mes sommets", "ogr")
  8. addresses = QgsVectorLayer("C:/Users/Georges/Downloads/qgis/data_BDTOPO_V3_Dep05_adresse/data_BDTOPO_V3_Dep05_adresse/ADRESSE_05.shp", "Addresses", "ogr")
  9.  
  10. QgsProject.instance().addMapLayer(addresses)
  11. QgsProject.instance().addMapLayer(troncons_routes)
  12. QgsProject.instance().addMapLayer(sommets)
  13. QgsProject.instance().addMapLayer(rlayer)
  14.  
  15. # Passer les sommets en haut des calques
  16. root = QgsProject.instance().layerTreeRoot()
  17. root.setHasCustomLayerOrder (True)
  18. order = root.customLayerOrder()
  19. order.insert(0, order.pop( order.index( sommets ) ) )
  20. root.setCustomLayerOrder( order )

Mais pour travailler en standalone, il vous faudra faire quelques importations et prendre des précautions supplémentaires :

  1. import os
  2. import datetime, requests, html, wikipedia
  3. import xml.etree.ElementTree
  4. from qgis.core import *
  5. from qgis.gui import *
  6. from qgis.utils import *
  7. from qgis.PyQt.QtGui import *
  8.  
  9. # Initialiser QGIS
  10. qgs = QgsApplication([], False)
  11.  
  12. # Cette ligne semble inutile, chercher pourquoi
  13. qgs.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
  14. qgs.initQgis()
  15. # Ces lignes également
  16. #for alg in QgsApplication.processingRegistry().algorithms():
  17. # print(alg.id(), "->", alg.displayName())
  18. if len(QgsProviderRegistry.instance().providerList()) == 0:
  19. raise RuntimeError('No data providers available.')
  20. ...
  21. qgs.exitQgis()

Essayer d'adapter votre code de génération de carte en mode standalone. En commençant par la génération d'une seule carte ce sera plus facile à déboguer. Une chose particulièrement subtile sera de parvenir à agir sur l'interface QGIS alors même qu'il n'y a plus d'interface QGIS (les fonctions de zoom par exemple).

Code complet de génération d'une carte en standalone

Vous pouvez aller plus loin en vous intéressant à la modification des variables d'environnement de votre machine. Cela vous permettra par exemple d'exécuter vos scripts depuis la version Python exécutable directement par Windows, et donc par un simple double-clic sur vos fichiers Python. Nous y reviendrons également.

  • Note personnelle : ici prévoir exposé sur les variables d'environnement et les possibilités standalone et serveur de QGIS)
  • Note personnelle : approfondir, variables d'environnement sous Windows et Ubuntu, centrer les layouts, PyGIS.

Fonctions géométriques

Nous allons lister dans un fichier CSV les adresses des Hautes-Alpes se trouvant à 2kms des limites d'un glacier.

Au terme de cet exercice, il nous faut un code complet permettant de générer ce fichier à partir d'un projet QGIS vide, et d'y inclure le nombre d'adresses concernées.

Buffer

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

  1. vlayer1 = QgsVectorLayer("C:/Users/Georges/Downloads/qgis/peaks/peaks.shp", "Mes sommets", "ogr")
  2. QgsProject.instance().addMapLayer(vlayer1)
  3.  
  4. bufferPath = "C:/Users/Georges/Downloads/qgis/buffers/buffer1.shp"
  5.  
  6. processing.run('native:buffer', {"INPUT": vlayer1, "DISTANCE": 5, \
  7. "SEGMENTS": 10, "END_CAP_STYLE": 0, "DISSOLVE": False, "OUTPUT": bufferPath})
  8.  
  9. vbuffer1 = QgsVectorLayer(bufferPath, "Buffer 1", "ogr")
  10. QgsProject.instance().addMapLayer(vbuffer1)
  11.  
  12. # Passer Buffer en bas
  13. root = QgsProject.instance().layerTreeRoot()
  14. myAlayer = root.findLayer(vbuffer1.id())
  15. myClone = myAlayer.clone()
  16. parent = myAlayer.parent()
  17. parent.insertChildNode(-1, myClone)
  18. parent.removeChildNode(myAlayer)

Les 2 premières lignes appellent la couche peaks, la 4ème stocke le chemin de notre future zone tampon.

C'est les lignes 6 et 7 qui à elles-seules créent le buffer grâce à l'instruction native/buffer. D'autres paramètres sont disponibles.

Les lignes 9 et 10 appellent ensuite notre couche contenant la zone tampon dans le canvas. Le reste du code est là pour passer les tampons sous les sommets, afin d'y voir clair.

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.

  1. vlayer2 = QgsVectorLayer("C:/Users/Georges/Downloads/qgis/data_BDTOPO_V3_Dep05_adresse/data_BDTOPO_V3_Dep05_adresse/ADRESSE_05.shp", "Addresses", "ogr")
  2. QgsProject.instance().addMapLayer(vlayer2)
  3.  
  4. bufferPath = "C:/Users/Georges/Downloads/qgis/buffers/buffer1.shp"
  5.  
  6. processing.run('native:buffer', {"INPUT": vlayer2, "DISTANCE": 2000, \
  7. "SEGMENTS": 10, "END_CAP_STYLE": 0, "DISSOLVE": False, "OUTPUT": bufferPath})
  8.  
  9. vbuffer1 = QgsVectorLayer(bufferPath, "Buffer 1", "ogr")
  10. QgsProject.instance().addMapLayer(vbuffer1)
  11.  
  12. # Passer Buffer en bas
  13. root = QgsProject.instance().layerTreeRoot()
  14. myAlayer = root.findLayer(vbuffer1.id())
  15. myClone = myAlayer.clone()
  16. parent = myAlayer.parent()
  17. parent.insertChildNode(-1, myClone)
  18. parent.removeChildNode(myAlayer)

Glaciers Alpins

Nous allons compter combien d'adresses se trouvent à 2 kms de distance des limites d'un glacier.

Téléchargez les glaciers contenus dans l'emprise de nos adresses Hautes-Alpines en utilisant l'extension QuickOSM (filtre natural=glaciers). Une couche sur ce modèle est disponible dans les pièces jointes de cet article.

Enregistrez-en un shape dans le même système de projection que vos adresses/tampons.

Pour contrôler la projection d'un layer.

{layer}.crs()

Intersection

L'algorithme de calcul d'intersection se présente de façon assez similaire que celle du buffer. D'autres paramètres sont disponibles.

  1. processing.run('qgis:intersection', {\
  2. "INPUT": {LayerInput},\
  3. "OVERLAY": {LayerOverlay},\
  4. "INPUT_FIELDS": {[list]},\
  5. "OVERLAY_FIELDS": {[list]},\
  6. "OVERLAY_FIELDS_PREFIX": {string},\
  7. "OUTPUT": {File}})

Introduisez l'algorithme ci-dessus dans votre code, exécutez-le puis zoomez, ci-dessous, sur le hameau de la Bérarde, au cœur de la zone de glaciers des Écrins, pour contrôler vos résultats.

  1. # Zoom sur la Berarde
  2. myselect = layerAdresse.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"NOM_1" = \'LA BERARDE\'' ) )
  3. layerAdresse.selectByIds( [ f.id() for f in myselect ] )
  4. iface.mapCanvas().zoomToSelected(layerAdresse)
  5. iface.mapCanvas().zoomScale(100000)

Vos intersections doivent exister pour chaque adresse intersectées (présence de doublons opportuns).

Si votre vérification est bonne (même nombre d'entités intersectées dans un groupe de zones tampons proches que d'adresses au sein de ce groupe), alors il ne vous reste plus qu'à compter le nombre total d'intersections.

{Layer}.featureCount()

Export CSV

Vous savez déjà exporter un fichier CSV (voir chapitre Exporter des fichiers plus haut dans ce court), à vous de jouer !

N'oubliez pas d'inclure le nombre d'adresses concernées dans le nom du fichier.

Nettoyage

La couche de nos zone tampons et celle de nos intersections sont à présent inutiles. Nous avons fait le job, nous pouvons donc supprimez ces couches mais aussi fermer QGIS.

  1. from osgeo import ogr
  2. ...
  3. # Nettoyage
  4. QgsProject.instance().removeMapLayer(layerBuffer)
  5. QgsProject.instance().removeMapLayer(intersection)
  6.  
  7. driver = ogr.GetDriverByName("ESRI Shapefile")
  8.  
  9. driver.DeleteDataSource(bufferPath)
  10. driver.DeleteDataSource(interstectPath)
  11.  
  12. os._exit(0)

Il semble toute fois que la suppression totale des shapes soit impossible à réaliser à la suite du code complet. Sans doute une mémoire cache à vider dans QGIS ou Python.

Code complet de la génération des buffers, des intersections et du fichier CSV (modifiez les chemins, fonctionne dans l'éditeur Python exclusivement)

  • Note personnelle : aller plus loin en se passant de la génération des shapes des buffers (couche virtuelle à la place), des glaciers (API Overpass d'OSM) et des intersections (couche virtuelles). Afin d'obtenir le simple fichier CSV des adresses à 2 kms d'un glacier, et leur nombre.

Exercice

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 :

  1. ...
  2. iris = QgsVectorLayer("C:/Users/Georges/Downloads/qgis/data_IRIS_2019/data_IRIS_2019/iris.shp", "IRIS", "ogr")
  3. ...
  4. QgsProject.instance().addMapLayer(iris)
  5. ...
  6. # Intersecter les sommets et les IRIS
  7. peaks_iris_path = "C:/Users/Georges/Downloads/qgis/exports/peaks_iris.shp"
  8.  
  9. processing.run('qgis:intersection', {\
  10. "INPUT": peaks,\
  11. "OVERLAY": iris,\
  12. "INPUT_FIELDS": ["OSM_ID", "NAME", "OTHER_TAGS"],\
  13. "OVERLAY_FIELDS": ["CODE_IRIS", "NOM_COM"],\
  14. "OVERLAY_FIELDS_PREFIX": "",\
  15. "OUTPUT": peaks_iris_path})
  16.  
  17. peaks_iris = QgsVectorLayer(peaks_iris_path, "Sommets", "ogr")
  18. QgsProject.instance().addMapLayer(peaks_iris)
  19.  
  20. # Retirer ancienne couche des sommets et celle des IRIS
  21. QgsProject.instance().removeMapLayer(peaks)
  22. QgsProject.instance().removeMapLayer(iris)

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()

Code complet de l'export des cartes avec les noms de commune en sous-titre et les liens aux APIs (modifiez les chemins, testez ce code de préférence dans l'éditeur Python, pas la console).


É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 :

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(49.0297879, 2.0616833)))
    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 1er 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)

Exposés

Durant le cours je demanderai aux étudiants de courts exposés (entre 3 et 6 minutes), à présenter au début des séances suivantes.

  • Les langages interprêtés/compilés et Python
  • Les chaînes formatées en Python
  • Les différents usages possibles de Python (y compris hors SIG)
  • Les grandes logiques de programmation (tous langages)
  • Les grands processus de développement (le mode projet)
  • Les APIs
  • Python et ArcGIS

 

 

Liens ou pièces jointes
Accéder à cette adresse URL (https://hg-map.fr/extern/data/shapes/france/chemin_de_fer.zip)chemin_de_fer.zip[ ]0 Ko
Télécharger ce fichier (data_BDTOPO_V3_Dep05_adresse.zip)data_BDTOPO_V3_Dep05_adresse.zip[ ]3889 Ko
Télécharger ce fichier (data_IRIS_2019.zip)data_IRIS_2019.zip[ ]45905 Ko
Télécharger ce fichier (decathlon_france.zip)decathlon_france.zip[308 magasins Décathlon français depuis OSM le 27 décembre 2020]11 Ko
Accéder à cette adresse URL (https://hg-map.fr/extern/data/shapes/france/eau.zip)eau.zip[ ]0 Ko
Télécharger ce fichier (glaciers.zip)glaciers.zip[ ]231 Ko
Télécharger ce fichier (iso_iris.zip)iso_iris.zip[Des zones isochrones à 15 minutes autour de 308 POIs.]12125 Ko
Télécharger ce fichier (MasterGeom2_ProgPython_DevoirMaison.pdf)Devoir-maison[ ]422 Ko
Télécharger ce fichier (peaks.zip)peaks.zip[ ]14 Ko
Accéder à cette adresse URL (https://hg-map.fr/extern/data/shapes/france/troncons_routes.zip)troncons_routes.zip[ ]0 Ko