créer une couche à partir de l’étendue de ses cartes à imprimer

Dans cet article, je vous parlerai d’un petit script QGIS python permettant de présenter sur votre carte l’étendue de vos cartes à imprimer. Ce sera aussi l’occasion de voir comment créer un traitement one-click pour la boîte à outils traitement.

J’habite une région magnifique, la Provence, et voilà un certain temps que je prépare mes randos du week end grâce à QGIS.

D’ailleurs, si vous ne le savez pas déjà, sachez qu’il est possible de bénéficier des cartes IGN dans QGIS sans clé avec cette adresse-là : 

https://wxs.ign.fr/pratique/geoportail/wmts?SERVICE=WMTS&REQUEST=GetCapabilities

(cf cet article sur le site de l’IGN)

J’ai une petite couche de données appelée tracés qui me permet de tracer mes futures randos, et de suivre celles que j’ai réalisées. J’ai aussi une couche de points qui me permet de marquer des secteurs intéressants.

Je pourrais utiliser geoportail.gouv.fr pour imprimer mes cartes, mais l’avantage que je vois avec QGIS est de vraiment pouvoir contrôler l’aspect de mes cartes, leur échelle d’impression. Je peux par exemple superposer une grille de 1 km de côté comme sur les cartes randos de l’IGN que l’on trouve chez le libraire. Cela s’avère très pratique une fois en balade.

Le souci, c’est que comme cela fait un certain temps que je crée mes cartes de randos sous QGIS, je me retrouve avec un nombre incalculable de composeurs d’impression (appelés composers sous QGIS 2 et print layouts sous QGIS 3)

Il s’avère que toutes les cartes de rando que j’ai exportées correspondent à des randos que j’ai réalisées. Autrement dit, si je pouvais représenter sur ma carte l’étendue de toutes ces cartes, cela me permettrait de me rappeler des secteurs où j’ai marché dans la région. 

Il y a d’autres cas où cela peut être utile, dans un cadre plus professionnel. Imaginez par exemple un agent de votre organisation qui aurait créé des cartes à tout va pour préparer une mission terrain, mais pour lequel il se révèlerait assez dur quelques mois plus tard de vous dire sur une carte les zones sur lesquelles il aurait travaillé.

En fait, il suffirait de créer une couche de polygones rectangulaires récupérés depuis l’étendue de mes cartes pour voir ces zones.

Le code python

S’il n’existe pas de fonction toute prête sous QGIS pour reporter dans une couche les polygones rectangulaires correspondant aux étendues des cartes à imprimer, nous avons quand même accès à  toutes les fonctions de l’API QGIS.

Si vous voulez tester votre code python, ce que vous pouvez tout simplement faire, c’ est d’ouvrir la console python et d’y créer un petit script (en allant dans le menu Extension > Console Python).

Voici le code que vous pourrez mettre dans la console ou dans votre script de console, pour récupérer dans une couche mémoire l’étendue de vos cartes (valable pour QGIS3).

def getCompositionExtents(crsCode=2154) :
    
    # layer
    vl = QgsVectorLayer("Polygon?crs=EPSG:%s"%(crsCode), "Extent Polygons", "memory")

    # defining CRS
    crs = vl.crs()
    crs.createFromId(crsCode)
    vl.setCrs(crs)
    
    # data provider
    pr = vl.dataProvider()

    # start editing
    vl.startEditing()
    
    # add attributes
    vl.addAttribute(QgsField("name", QVariant.String))


    # get compositions
    projectInstance= QgsProject.instance()
    projectLayoutManager = projectInstance.layoutManager()
    comps = projectLayoutManager.printLayouts()

    # iterate over compositions
    feats = list()
    for comp in comps :
        nom = comp.name()
        print(nom)
        for item in comp.items() :
            if (isinstance(item, QgsLayoutItemMap)) :
                print("toto")
                print(item)
                e = item.extent()
                feat = QgsFeature()
                feat.setGeometry(QgsGeometry.fromRect(e))
                feat.initAttributes(1)
                feat.setAttribute(0, QVariant(nom))
                feats.append(feat)

    # add features
    pr.addFeatures(feats)

    # commit changes
    vl.commitChanges()
    
    return(vl)
    
def addLayerToCanvas(vl) :
    # add layer to canvas
    QgsProject.instance().addMapLayer(vl, False)
    layerTree = iface.layerTreeCanvasBridge().rootGroup()
    layerTree.insertChildNode(-1, QgsLayerTreeLayer(vl))
    
    
vl = getCompositionExtents() # create layer from extents
addLayerToCanvas(vl) # add layer to canvas

Pour lancer le script, cliquez juste sur le petit triangle vert

Et voilà, cela devrait vous créer une couche avec l’étendue des cartes à imprimer de votre projet

Décodage

Voici les éléments clés du code :

projectInstance= QgsProject.instance()
projectLayoutManager = projectInstance.layoutManager()
comps = projectLayoutManager.printLayouts()

Ce code permet de récupérer les print layouts dans une variable appelée comps

feats = list()
for comp in comps :
    nom = comp.name()
    for item in comp.items() :
        if (isinstance(item, QgsLayoutItemMap)) :
            e = item.extent()
            feat = QgsFeature()
            feat.setGeometry(QgsGeometry.fromRect(e))
            feat.initAttributes(1)
            feat.setAttribute(0, QVariant(nom))
            feats.append(feat)

Ici, on itère à travers toutes les compositions. Pour chaque composition, on aura un objet, une feature.

feat.setGeometry(QgsGeometry.fromRect(e)) permet d’affecter une géométrie à la feature depuis le rectangle formé par l’étendue

feat.initAttributes(1)
feat.setAttribute(0, QVariant(nom))

permet d’attribuer la variable nom (qui comprend le nom de la composition) à la première colonne

def addLayerToCanvas(vl) :
    # add layer to canvas
    QgsProject.instance().addMapLayer(vl, False)
    layerTree = iface.layerTreeCanvasBridge().rootGroup()
    layerTree.insertChildNode(-1, QgsLayerTreeLayer(vl))

Cette fonction permet juste d’ajouter la couche créée artificiellement au canevas, soit, à la carte

Création d’un traitement python pour la boîte à outils Traitement

Exécuter un script depuis la console python n’est pas la chose la plus élégante à faire. Il y a deux externalités possibles à ce bout de code, voire même trois.

Les voici dans un ordre de difficulté d’implémentation croissante :

  1. Créer une macro à l’ouverture du projet : automatiquement, une couche se crée à l’ouverture du projet
  2. Créer un traitement pour la boîte à outils Traitement
  3. Créer un plugin

Dans mon cas, j’ai préféré créer un traitement QGIS. Un plugin serait un peu trop ambitieux pour ce cas d’usage, quoiqu’un bouton dans l’interface aurait pu être utile.

Pour créer le traitement, créez un nouveau script, par exemple depuis un modèle :

Vous trouverez sur ce lien le code complet pour créer un traitement pour la boîte à outils Traitement

Cela créera une rubrique My Scripts dans la boîte à outils Traitement.

Cela apparaîtra de cette façon-là quand vous l’exécuterez :

Décodage

Dans le script, les fonctions name, displayName, group, groupId permettent de paramétrer la façon dont le script apparaîtra dans la boîte à outils Traitement.

(sink, dest_id) = self.parameterAsSink(
    parameters,
    self.OUTPUT,
    context,
    flds,
    3,
   crs
)

permet de paramétrer les méta-data de la couche. 3 est le type géométrique de la couche. fldset crs sont les champs et système de projection de la couche

feats  = self.fGetExtentFeaturesFromPrintLayouts()
for feat in feats :
    sink.addFeature(feat, QgsFeatureSink.FastInsert)
    
# Return the results of the algorithm
return {self.OUTPUT: dest_id}

calcule les features associées aux étendues (voir fonction créée précédemment) et les met dans le sink.

Note : pour les fonctions qui sont personnelles, je rajoute le petit préfixe f afin de bien les dissocier des fonctions QGIS du modèle.

Il y a vraiment beaucoup à faire et de quoi s’amuser quand on plonge dans l’API et QGIS offre plein de voies possibles afin d’implémenter des traitements originaux.

Et après ?

Nous pourrions imaginer rapatrier le nom des communes intersectées par ces rectangles, voir quelle surface cela représente en tout et révéler les zones de la région qui resteraient encore à explorer !

Petit rappel des gists :

 

 

Datagistips

datagistips, axé sur l’opensource SIG et l’opendata des données spatialisées, essaie de rendre la donnée amusante, via des cas d’usage inspirés du quotidien. Il informe aussi de nouveautés, sur certaines technologies, à appliquer dans le domaine de la géomatique. Se voulant original, il tente de relier au domaine SIG, des sphères communément extérieures à celle-ci : celle des statistiques, de la récupération de données depuis Internet, ou du design génératif appliqué à la cartographie.

Parcourir le Monde

Quel est l’itinéraire idéal pour parcourir le Monde ? Combien de temps de vol cela prendrait ? Quels sont les deux pays les plus distants l’un de l’autre ?
Dans ce post, j’ai fait mumuse avec la carte du globe, et appris des choses au sujet des distances orthodromiques et des voyageurs de commerce. Bon voyage !


Dernièrement, j’ai eu un échange avec @matamix au sujet du calcul de trajets permettant de couvrir un ensemble de points de façon optimale.

Nous avons échangé au sujet du code, puis une question m’est venue : dans quel ordre devrais-je visiter les différents pays, dans l’éventualité d’ un tour du monde (rêvons un peu). C’est la problématique du voyageur de commerce, qui doit démarcher un certain nombre d’adresses en un minimum de temps. Il s’avère que la problématique du voyageur de commerce peut s’appliquer à des domaines non spatiaux, comme ordonner des séquences génétiques, voire regrouper des individus (au sens statistique) selon leur proximité : intéressant !

“Attrape-moi si tu peux !”

Tout d’abord, il faut se poser la question des distances entre les différents pays. L’espace dans lequel on doit travailler n’est pas euclidien, plan, mais sphérique. Ainsi, pour rejoindre le Japon depuis les Etats-Unis, inutile de traverser l’Atlantique, bien sûr. Les distances que l’on doit calculer sont les distances de grands cercles ou orthodromiques. Il existe apparemment différentes méthodes, la plus courante étant celle utilisant la loi sphérique des Cosinus

Sous R, on peut calculer 4 sortes de distances des grands cercles grâce au paquet geosphere.

J’ai récupéré les contours des pays depuis le magnifique Natural Earth, déterminé le centre des pays, et calculé les distances entre chaque, de sorte à réaliser une matrice de distance entre les différents pays. Cette matrice peut être vue comme un réseau constitué de noeuds : les pays, connectés entre eux par des arêtes de poids égal à la distance qui les sépare.

On obtient une matrice des distances orthodromiques, dont voici ci-dessous une illustration, avec 177 lignes et 177 colonnes (autant de pays que dans la couche des pays), soit 177^2 “cases”. Plus la couleur est claire et plus la valeur de distance est grande. On voit que la diagonale est obscure. Les poids sur celle-ci sont nuls, puisque la distance de la France à la France est nulle. On voit que la matrice est symétrique puisque la distance, à vol d’oiseau, est la même de la France à l’Espagne que de l’Espagne à la France.

Entrez dans la matrice…

D’ailleurs, à partir de cette matrice, je me suis amusé à déterminer les trajets les plus grands entre deux pays sur le globe. En voici la liste des 10 premiers, chaque ligne représentant un trajet, dans l’ordre d’importance :

En allant du Paraguay à Taïwan, on parcourt quasiment la moitié de la Terre, la circonférence de cette dernière étant de 40 075 km

Et voici la carte correspondante pour les 5 plus grands trajets :

De même, voici la liste des pays les plus “isolés” :

Revenons à nos moutons…ou plutôt à nos voyageurs de commerce. Maintenant que notre matrice a été constituée, nous pouvons calculer le chemin optimal qui parcourt tous les noeuds une seule fois en cumulant un poids minimal (une distance minimale). Il existe plusieurs méthodes permettant de résoudre la problématique du voyageur de commerce. Dans un premier temps, j’ai tenté la plus commune, la méthode des plus proches voisins, mais elle me trouvait un itinéraire assez bizarre.

Le voyageur de commerce qui avait mal planifié son voyage

En fait, c’est normal : la méthode des plus proches voisins prend un point au hasard, détermine le point le plus proche, et répète l’analyse de proche en proche jusqu’à ce que tous les points aient été visités. Une variante prend comme point de départ chaque noeud, compare les tournées, et choisit la meilleure. Cet algorithme est assez basique.

C’est avec l’algorithme Concorde que j’ai obtenu le meilleur résultat. Voici ce que dit Wikipédia à propos de Concorde :
“a state of the art implementation”
“one of the best exact TSP solvers currently available.”
“is widely regarded as the fastest TSP solver, for large instances, currently in existence”.
Bref , il s’agit du voyageur de commerce le plus performant du monde, apparemment. Il doit être milliardaire à l’heure qu’il est !

 Voici donc l’itinéraire que j’ai pu trouver pour parcourir les pays du monde entier, prenant en compte les distances sphériques entre pays. On voit que j’ai aussi inclus l’Antarctique : soyons un peu aventureux !

Carte avec les fuseaux horaires

Voici un extrait des pays, dans l’ordre, au moins pour les premiers (ou les derniers, selon le sens où l’on va) :

Bon, Corée du Nord, 16e dans la liste, on passe ?

En distance, en tout, cela représente un peu plus du tiers de la distance Terre-Lune, cette dernière étant de 384 400 km, et entre trois et quatre fois le tour de la Terre

Je me suis demandé combien de temps de vol cela représenterait pour parcourir tout ça. Il existe un site appelé travelmath qui donne le temps de vol entre tous les pays.

(et 11 secondes)

Avec un petit chouilla de scraping – rien de compliqué – j’ai pu récupérer ces temps tout le long de mon itinéraire, soit pour 176 couples origine-destination. J’ai fait le calcul : il faudrait 12 jours ! pour effectuer cet itinéraire. En bref, si vous projetez une année sabbatique, ça vous fait un jour à passer dans l’avion par mois en moyenne, pas très reposant…et peu écologique !

Dataventures

Bon, je dois vous avouer qu’utiliser le Concorde (l’algorithme) n’a pas été forcément une mince affaire. Dans un premier temps j’ai testé l’outil Concorde en interface graphique et me suis rendu compte qu’il fournissait des bons résultats :

Pas mal !

Mais le souci, c’était que visiblement, ce dernier ne prenait en compte que les fichiers de coordonnées, et pas les matrices de poids !..

J’ai essayé de convertir ma matrice de distances orthodromiques en coordonnées avec un échelonnage multi-dimensionnel et une fonction de R appelée cmdscale, mais cette méthode est loin d’être optimale car j’avais en coordonnées le Cameroun plus proche de la France que la Belgique (!). Ca m’a un peu trotté cette histoire jusqu’à ce que je conclue, bêtement, que cela venait du fait que l’espace dans lequel raisonne cmdscale ne peut prendre en compte la nature sphérique des distances (c’est comme si on était dans une autre dimension avec ma matrice).

Comme un petit problème..

J’ai donc dû me rabattre sur les exécutables source. J’ai téléchargé Concorde depuis un site universitaire canadien

Le site, qui n’a peut-être jamais crashé, contrairement à son outil éponyme

A noter que la licence ne prévoit qu’une utilisation académique, et non commerciale, malgré le nom de l’algorithme 😉 L’outil est exécutable sous Windows grâce à Cygwin, un émulateur Linux. Du coup, j’ai installé Cygwin en 64 bits. Hors, ça ne fonctionnait pas. Je me suis aperçu après moult recherches que Concorde.exe fonctionnait avec Cygwin 32 bits, et non 64 bits. Bon, finalement, lorsque concorde.exe a affiché l’aide dans l’invite de commandes, j’étais soulagé !

**Alleluïa**

Une fois tout bien configuré, j’ai pu utiliser Concorde dans R après avoir donné le chemin via concorde_path (youpi !) :

Cool, ça marche. Le rouge, c’est juste un avertissement
Allez, pour finir, voici des cartes que l’on trouve sur le web, qui donnent le chemin le plus court pour parcourir la planète (pas si éloignées que ça de la mienne) :
Aller d’une capitale à l’autre selon wolfram

La réponse était aussi sur reddit, sur le canal “Map Porn”

Crédits icônographie
Globe issu de Wikipédia
The Crown, by Oliviu Stolan from NounProject
The Wing by Alice Noir from NounProject

Geospatial fun