On the 10th of April, I created a map to see the evolution of COVID cases in France as spatialized sparklines.
I first saw it like a good use case on how to use some QGIS related functions (QGIS is an opensource geospatial tool)
Another experiment on rendering #COVID-19 evolution in France in #qgis using geometry generator Tonight, I'll put the code or write a tuto ▶️Use of relation_aggregate, some blur effects and a pinch of array "magic" 🧙♂️#datavizpic.twitter.com/Ugr7svXNx9
“Les Artisans cartographes” mentioned the map on twitter. Thereafter, many people started retweeting it. I didn’t expect it to make so much “noise”
For me, this representation was an efficient way to see evolutions on a map because you could observe the dynamics at a single glance.
Maybe this instantaneous look can help in emergency cases because you don’t need to press play and wait for every frame to complete in order to get “the big picture”.
For me, there was nothing much innovative or disruptive about this map. It had been a long time the idea was in my head and I had already experimented it around census data.
I thought someone had already created something similar. But Alberto Cairo, a world-reknown dataviz expert, seemed to have never seen something like this before, which surprised me.
Well before I put the QGIS code to create the map on github and explained it on this blog, some people on twitter started to create their own. Later, I saw some similar maps emerge.
That’s a gallery of these initiatives you’ll discover here.
You can see that the Washington Post map is more sophisticated, with additional use of thickness to illustrate the current state (actual number of cases) and color to support the evolution
Note the Washington Post map was showcased on flowingdata, which is a reference media relaying dataviz experiments (also tutorials)
This week-end, because of containment, I had lots of time to experiment around data and covid-19 cases, trying to figure out which graphical representation could draw the dynamics of the pandemic.
I saw a lot of graphics. The curves and their dynamics are crucial with the pandemic. Not only it is important to know what is the actual situation, but also how it was before, and how we can project ourselves to the future.
The above map is the map I made and shared on twitter :
Another experiment on rendering #COVID-19 evolution in France in #qgis using geometry generator Tonight, I'll put the code or write a tuto ▶️Use of relation_aggregate, some blur effects and a pinch of array "magic" 🧙♂️#datavizpic.twitter.com/Ugr7svXNx9
It had quite a lot of success. It’s a mix of a line plot and a map. It gives a sense, at a single glance, of the dynamics, and looks good when animated.
After some feedback on twitter, I decied to improve the map with labels :
I only used opensource tools for it : qgis and imagemagick when animated. What’s great with opendata and open tools is that everyone can potentially integrate the data and experiment with no cost. You just have to let the ideas come, select them and find the ways to transcribe in the constrained environment of a machine : here, my computer. Somehow, the screen is my canvas painting, my mouse my brush and my keyboard my palette.
The goal of this article is to explain the method behind. It will give you a nice tutorial covering 3 qgis techniques :
1-n relations
geometry generators
the use of arrays
Some next articles could get deeper into the subject with styling.
Here’s how the covid fr dataset looks like. It is in a long format with one row by day, unlike world John Hopskins worldwide dataset which is wide (with each column being a day)
Relation
We’ll add a relation between the departements dataset and the covid dataset. We use the department code INSEE_DEPT as a key connecting the parent table (departments layer) and the child table (covid dataset). Go to Project > Properties > Relation
Generative design
Next, go to the department layer and use a centroid symbol, but for this centroid symbol, you’ll use geometry generator to define the linestring.
As we want lines, choose line / polyline features :
And add this code in the expression (I’ll explain it just afterwards) :
To understand this code, it’s important to split it into key parts. Obviously, I didn’t write it in a go. It was a trial-and-error process. Let me help you.
First, this line uses the relation we configured to store all the cases for each feature (here, each department), in an array
Before constructing the linestring object, let’s consider this :
the height of the segment will be defined by the number of cases
and the spacing between nodes, simply, by the index of the node. One index = one node.
the line will start from the departement centroid
So, once we have this array of number of cases, we’ll iterate over it to get the consecutive nodes of the lines. We’ll just need the centroid coordinates, the indexes and the number of cases for each of them to build the lines
For the X coordinates of the nodes it’s this formula :
To get the index (@element), we use generate_series
generate_series(1, @atlasid),
-- X Coordinates
to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
-- Y coordinates
to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
)
To get the number of cases, we get the array element corresponding to the index :
@nb_cas[@element]
The iteration over indexes is made with array_for_each.
array_foreach(
generate_series(1, @atlasid),
-- X Coordinates
to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
-- Y coordinates
to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
)
Once we have the array of coordinates, we paste them with array_to_string with the comma as delimiter, so that they can be interpreted as WKT string for a linestring
array_to_string(
array_foreach(
generate_series(1, @atlasid),
-- X Coordinates
to_string(x(centroid($geometry)) + (@element/array_length(@nb_cas)*@width)) ||' '||
-- Y coordinates
to_string(y(centroid($geometry)) + @nb_cas[@element]*@height)
),
',')
Our final geometry, transformed from a WKT string, will be :
To make this the most programmatic as possible, I configured variables for the total width of the lines (@width), and the @height factor for each number of case. This way, it’s easy to change the appearance of the map
Also, I added a variable for the atlas feature id (@atlas_id) because I wanted to get the individual frames I could assemble as a GIF with ImageMagick.
For the animation, the trick consists in using an atlas and generating a series which will stop at the atlas current feature id, this one going from 0 to (n-1), n being the number of days.
So the first image will have one node (the centroid), the second one the centroid + the seconde node, the third image, the centroid + the two following nodes.
Note : why did I say 0 and (n-1) insteand of 1 and n ? It’s because the first element of an array, in qgis, is like in python : it has number 0.
Maybe later I’ll write about styling the data underneath, but now, you have the mechanism and logic behind !
There’s a github repo 😉
You struggled to follow all these steps ? Wait ! I put all the project on a github repo with all the data you’ll need and a totally reproducible map. You’ll just have to follow the steps on the github page (actually, one single click on the project after downloadgin qgis)
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à :
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
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
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 :
Créer une macro à l’ouverture du projet : automatiquement, une couche se crée à l’ouverture du projet
Créer un traitement pour la boîte à outils Traitement
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.
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 !