Mentrek, etc: octobre 2015

Dans Postgres, avec le SQL spatial et des statistiques simples, on pose des faits et des questions politiques

Les inondations, les submersions de portions du littoral, l'exposition à des risques industriels sont des questions qui préoccupent, à juste titre, les citoyens. Elles traitent au fond du voisinage, de l'environnement proche. Là où les gens vivent dit beaucoup sur leurs vies : les réseaux de relations, les emplois, les services auxquels ils peuvent prétendre, les équipements publics accessibles, les déplacements contraints, les sources d'agrément ou de nuisance, les risques encourus...

Du côté des aménageurs, des géographes on peut voir ces phénomènes divers à travers la question de la répartition de la population. Parce que, dans les decennies passées, les effets de l'étalement urbain et l'artificialisation rapide des sols ont joué à plein, cette question profondément politique a pris une acuité particulière - et c'est là où la population rejoint le peuple. Qui, traitant des affaires de la cité, peut ignorer où il se trouve ?

Open data pour le débat public

Eurostat publie sur son site un jeu de données carroyées sur la population à l'échelle de l'Europe. C'est une grille de polygones (carrés) de 1 km de côté et un tableau de chiffres dont chaque ligne est identifiée par par une coordonnée de grille et un code de pays. Comme il existe des coordonnées identiques entre pays et que les éléments de la grille ne sont pas liés à un pays, il n'est pas possible de relier les deux de façon rigoureuse de prime abord. Mais on va faire ça avec PostgreSQL et des fonctions spatiales de SQL.

Après avoir créé dans Postgres une base et ajouté les extensions spatiales (voir le manuel de Postgis), on y importe le tableau de population et la grille d'éléments. Pour ce dernier, on emploie shp2pgsql avec une projection EPSG 3035, faite pour représenter l'ouest du continent.

Les jointures spatiales pour relier l'information

On importe ensuite le tableau des entités administratives des pays européens NUTS. L'idée est de faire une jointure spatiale entre la table des éléments de la grille et la table des entités NUTS et de récupérer le code nuts en guise de code pays. À ceci près que les deux ne sont pas toujours identiques. Dans le cas du Royaume-Uni, il existe un code pays pour l'Angleterre, l'Écosse, le Pays de Galles et l'Irlande du nord. Au préalable, on créé, dans la table NUTS, un attribut pays (cntr_code) reprenant le code nuts, sauf pour ces cas particuliers.

D'ordinaire, pour enrichir des données dans une table, on joint une autre table sur un ou plusieurs champs communs. Ici, il n'y a pas de champ commun mais c'est une relation spatiale, qui peut être d'inclusion, d'intersection, de distance (moyenne, minimale...), qui permet de relier entre elles les informations et d'ajouter du sens.

L'autre idée est de simplifier la géométrie de la grille avec la fonction st_centroid(). Un ensemble de points, centres de chaque carré, est équivalent en information. Il sera notamment plus simple de déterminer à quel pays rattacher un point qu'un carré qui pourrait être à cheval sur une frontière.

On créé donc la nouvelle table, qui comprend une géométrie de type point et dont chaque élément est identifié par une coordonnée de grille et un code de pays.

create table eupop_q as
select c.geom, c.code as grid_id, c.cntr_code, c.nuts_id
from
    (select a.code,b.cntr_code,b.nuts_id,a.geom
      from
        (select * from eustat.eeu
          where (niveau = 0 and nuts_id != 'UK') or (niveau = 1 and nuts_id like 'UK%') ) as b
      join eustat.eupop_grid as a
      on st_within(st_centroid(a.geom),b.geom)
    ) as c
with data
;

Une jointure spatiale sur une table de près de deux millions d'enregistrements prend un certain temps, surtout sur une machine plus jeune du tout. Pour l'exploitation, il existe quelques moyens pour accélérer les traitements. Mais, il faut l'avoir à l'esprit, une requête mal réfléchie (je suis sûr que c'est le cas ici) prendra un temps considérable sur un matériel ultra rapide, dès lors que les tables ont beaucoup de lignes.

La table obtenue par la requête ci-dessus n'est pas parfaite, loin s'en faut. On se rend compte qu'elle omet tous les points situés un peu en dehors des limites administratives, sur les côtes, côté mer. Le carroyage d'Eurostat n'est précis qu'à 700 m près, et parfois l'imprécision est à dessein, par souci de préserver les personnes (certains carrés contiennent un seul individu). On a trouvé plus de 20 000 points dans ce cas. Comment rattacher ces points au pays le plus proche et leur attribuer un code ?

Logique des ensembles

La méthode consiste à chercher, pour tout point situé en dehors des limites administratives, la distance à chaque pays et d'autre part de rechercher, pour chacun de ces point, la distance au plus proche pays. Les deux sous-requêtes ne diffèrent que par le jeu des groupements (group by). La jointure entre les deux tables résultantes donne le pays le plus proche, avec ses attributs. On réintègre le résultat dans le tableau de points (avec insert into).

select x.geom,x.code as grid_id,x.cntr_code,x.nuts_id
from (
select min(st_distance(st_closestpoint(a.geom,b.geom),b.geom)) as mini,b.id,b.geom,b.code,a.cntr_code,a.nuts_id
from (select * from eustat.eeu where (niveau = 0 and nuts_id != 'UK') or (niveau = 1 and nuts_id LIKE 'UK%')) as a,
(select id,st_centroid(geom) as geom,code from eustat.eupop_r) as b
group by b.id,b.geom,b.code,a.cntr_code,a.nuts_id
) as x
join (
select min(st_distance(st_closestpoint(c.geom,d.geom),d.geom)) as mini,d.id,d.code
from (select * from eustat.eeu where (niveau = 0 and nuts_id != 'UK') or (niveau = 1 and nuts_id LIKE 'UK%')) as c,
(select id,st_centroid(geom) as geom,code from eustat.eupop_r) as d
group by d.id,d.code
) as y
on x.id = y.id
where x.mini = y.mini
;

Ya peut-être une meilleure façon mais je l'ai pas trouvée. La table est complète. Manquent quelques paradis fiscaux et le Vatican, ce qui est négligeable.

Une jointure plus loin, avec le tableau de population, sur les attributs de grille et de pays, et on obtient la table qui sera exploitée dans l'outil de cartographie QGIS.

Après un petit travail sur la visualisation dans QGIS (cercles proportionnels à la population) le résultat est spectaculaire et éloquent. La présence ou l'absence de population révèle bien sûr les agglomérations mais aussi les fleuves, les routes, les lacs, les montagnes, les vallées de montagne.

Mais la contemplation a ses limites. Et le SQL spatial donne des résultats chiffrés, par exemple pour évaluer le voisinage d'une centrale nucléaire. Par exemple, combien de gens vivent à proximité de nos centrales, disons 20 km ?

Nous avons préparé un fichier de points avec quelques établissements industriels. La requête combine les fonctions spatiales et les statistiques, ici une simple somme.

select c.nom_centrale,sum(p.population) as pop_evacues
from eustat.centralesnuc as c join eustat.eupop_q as p
on st(intersects(st_buffer(c.geom,20000),p.geom)
;