Colorier une carte (avec postgis)

Afin de créer une carte “Police/Justice” montrant les ressorts des Tribunaux de Grande Instance, j’ai besoin de colorier ceux-ci pour les rendre plus visibles qu’avec juste leurs frontières.

Colorier une carte avec un nombre limité de couleurs en faisant en sorte que 2 polygones contigus soient toujours de couleur différente est un problème qui peut sembler simple, mais après avoir un peu cherché je n’ai pas trouvé d’algorithme “tout fait” pour faire cela.

La théorie veut que 4 couleurs suffisent pour colorier une carte (sauf cas particuliers) et y parvenir à coup sûr est toutefois assez complexe. On doit sûrement pouvoir s’en sortir avec pas mal de récursion.

Par contre si on ne se limite pas à 4 couleurs mais à 5 voire 6, ça devient relativement plus facile à implémenter.

Premier essai

Voilà la première méthode que j’ai appliqué en utilisant postgis

1. ajouter un champ “couleur” à ma table contenant mes polygones
2. attribuer la couleur “1” à un polygone arbitrairement (notre point de départ)
3. chercher un polygone (p1) de couleur nulle touchant un (ou plusieurs) polygone (p2) déjà coloré

• lister les couleurs déjà utilisées sur les polygones (p2) adjacents à p1
• attributer à p1 la première couleur non utilisée par ces polygones p2
• boucler sur 3 jusqu’à plus soif

Voici la requête utilisée en 3:

with u as (
select p1.id as u_id,
substring(regexp_replace(‘123456’,’[‘||string_agg(distinct(p2.couleur),’’)||’]*’,’’),1,1) as u_couleur
from poly p1
join poly p2 on (st_touches(p1.geom,p2.geom) and p2.couleur is not null)
where p1.couleur is null
group by p1.id
limit 1
)
update poly set couleur=u_couleur from u where id=u_id;

Explication de la ligne 3…

- string_agg sert à combiner les couleurs utilisées en une chaîne de texte

- ce texte est utilisé pour construire une expression régulière qui sert à supprimer les couleurs trouvées de la liste des couleurs possibles ‘123456’

Et voilà, en 80s tout était colorié pour mes 160 TGI !

Un défaut est que les premières couleurs sont attribuées en priorité et donc plus fréquentes. Pour ré-équilibrer, on peut réattribuer des couleurs 5 peu fréquentes aux polygones de couleur 1 très fréquents et ne touchant pas un polygone de couleur 5.

Ceci se fait par exemple avec la requête suivante:

with u as (
select p1.id as u_id
from poly p1
left join poly p2 on (st_touches(p1.geom, p2.geom) and p2.couleur=’5')
where
p1.couleur=’1' and p2.id is null
limit 15
)
update poly set couleur=’5' from u where id=u_id;

En ligne 7 on peut régler le nombre de polygones à basculer de la couleur 1 à 5

Cet algorithme simpliste m’a permis d’obtenir un coloriage des communes françaises avec 7 couleurs… mais on doit pouvoir faire mieux.

On vise les 4 couleurs ?

Sur l’animation visualisant le remplissage, on voit que celui-ci n’est pas fait en “tâche d’huile”.

L’intuition me dit qu’un remplissage en tâche d’huile permettrait sûrement de limiter l’usage des dernières couleurs voire à arriver sur des cas optimaux à 4 couleurs.

Le problème tourne donc au choix du prochain polygone à colorier.

Avec la première requête, on laissait postgis nous retourner le premier polygone trouvé ce qui fait que l’ordre n’est pas du tout contrôlé et franchement aléatoire.

Après plusieurs expérimentations, il semble qu’en choisissant le polygone ayant le plus grand nombre de voisins déjà colorés et avec un voisin comportant à la couleur la plus “élevée”, on arrive à faire un coloriage plus optimal.

Le coloriage de mes TGI se fait maintenant avec seulement 4 couleurs !

La nouvelle requête utilisée:

with u as (
select p1.id as u_id,
substring(regexp_replace(‘123456’,’[‘||string_agg(distinct(p2.couleur),’’)||’]*’,’’),1,1) as u_couleur
from poly p1
join poly p2 on (st_touches(p1.geom, p2.geom) and p2.couleur is not null)
where p1.couleur is null
group by p1.id
order by count(p2.*) DESC, max(p2.couleur) DESC
limit 1
)
update poly set couleur=u_couleur from u where id=u_id;

L’ajout du ORDER BY en ligne 8 est le seul changement et sert à sélectionner le polygone ayant le plus de voisins déjà colorés (count) et utilisant la couleur la plus élevée (max).

L’inconvénient c’est que les requêtes sont beaucoup plus longues, on passe à 4 minutes au lieu de 80s.

Optimisons le temps de calcul

L’utilisation de ST_Touches de postgis est coûteuse en calculs, l’idée est donc de faire ce calcul une seule fois, de stocker dans une table “voisins” les voisinages entre polygones et de n’utiliser que cette table par la suite.

On peut la créer ainsi:

CREATE TABLE voisins AS (SELECT p1.id as a, p2.id as b FROM poly p1 JOIN poly p2 ON (ST_Touches(P1.geom, p2.geom) and st_geometrytype(st_intersection(p1.geom, p2.geom))!=’ST_Point’)));

J’ai aussi ajouté un test pour ne pas considérer comme voisins des polygones qui ne se touchent que par un point.

Un index sur la première colonne n’est pas du luxe car les recherches se feront toujours sur elle.

CREATE INDEX voisins_a ON voisins (a);

La requête principale peut donc être ré-écrite comme suit:

with u as (
select p1.id as u_id,
substring(regexp_replace(‘123456’,’[‘||string_agg(distinct(p2.couleur),’’)||’]*’,’’),1,1) as u_couleur
from poly p1
join voisins v on (v.a=p1.id)
join poly p2 on (p2.id=v.b and p2.couleur is not null)
where p1.couleur is null
group by p1.id
order by count(p2.*) DESC, max(p2.couleur) DESC
limit 1
)
update poly set couleur=u_couleur from u where id=u_id;

Le coloriage des 165 TGI descend ainsi à 10s… soit 24 fois plus rapide :)

Colorier la carte de France des communes ?

Passer de quelques centaines de polygones à plus de 36000 crée forcément des cas où cet algorithme simpliste va diverger. Trouver le coloriage optimal avec 4 couleurs sur un graphe aussi important n’est pas gagné, on va au moins tenter de réduire le nombre de couleurs nécessaires: objectif 5.

• Premier constant la table des voisinages est assez volumineuse (plus de 200000 voisinages).
• Deuxième constat, l’effet tâche d’huile n’est pas optimal car l’algo se disperse assez facilement.

Pour palier à ces deux problèmes, plutôt que de vouloir colorier toute la France d’un coup, je le fais département par département. Ainsi, la table voisinage n’est remplie qu’avec les voisinages des communes d’un département, puis on fait tourner la boucle de coloriage, on passe au département limitrophe suivant et on colorie et ainsi de suite.

Le temps de coloriage d’un département est d’environ une minute, et au final 7 communes seulement se sont vues attribuées une sixième couleur.

Éliminer les tâches

Ces 7 tâches je vais les éliminer en remettant leur couleur à 0 puis en faisant de même pour les 3 niveaux de communes limitrophes… et ensuite je relance le coloriage de ces trous en croisant les doigts…

Bingo ! 5 couleurs suffisent !

Un petit peu de ré-équilibrage (aléatoire) entre les couleurs et voici le résultat après à peine plus d’une heure de calcul:

Les

Les scripts et requêtes utilisés pour colorer les communes de France sont sur https://github.com/cquest/coloriage

J’ai essayé de le rendre assez générique pour pouvoir l’adapter rapidement au coloriage de tout ensemble de polygones.

Le résultat pour les 35888 communes au 1/1/2016 est disponible en “ressource communautaire” (bas de page) sur https://www.data.gouv.fr/fr/datasets/decoupage-administratif-communal-francais-issu-d-openstreetmap/

Et la version au 1/1/2018 est elle aussi disponible sur https://www.data.gouv.fr/fr/datasets/decoupage-administratif-communal-francais-issu-d-openstreetmap/#resource-0693014d-0a4a-4690-930a-0a7d5013ba56