Une fonction PL/pgsql pour l’analyse exploratoire de données spatiales avec des polygones de Voronoï

L’analyse exploratoire des données spatiales est une étape incontournable dans les processus d’analyse spatiale.

Un des outils pour cette exploration est la construction des polygones de Voronoï. (Voir Analyse exploratoire des données pour la géostatistique:les diagrammes de Voronoï)

Si l’extension Geostatistical Analyst d’ArcGis permet une grande variété d’affichages d’une couche de polygones de Voronoï, avec QGis on est limités au simple affichage des polygones et des valeurs de l’attribut analysé. Ce qui manque c’est les différents affichages correspondant à des statistiques locales.

Dans un précédent article (Calculer des moyennes locales sur des polygones de Voronoï avec Qgis), qui commence à dater, nous avions donné les moyens d’afficher la vue « Moyenne » des polygones: Pour chaque polygone on recherche les polygones qui ont un côté ou un vertex commun avec lui, puis on calcule la moyenne d’un attribut pour cet ensemble (la moyenne du polygone central plus les valeurs des polygones adjacents). C’est cette moyenne qui est affectée au polygone au lieu de la valeur originale de l’attribut..

C’est vrai que dans les traitements disponibles dans QGis, ou bien utilisant des scripts R, on peut calculer ces statistiques. mais le processus et long et fastidieux.

Pour les utilisateurs de Postgresql/Postgis, il faut noter que Postgis ne possède pas du tout ce type de fonctions.

La solution présentée ici implique la disponibilité d’une base Postgresql/postgis et de QGis:

  • On utilise QGis pour calculer les polygones de Voronoï en créant une table Postgresql en sortie.
  • On exécute la fonction dans la base de données pour créer les champs statistiques dans la couche de VoronoÏ
  • On utilise QGis pour visualiser les différents champs statistiques

Les statistiques locales des polygones de Voronoï

Dans l’analyse exploratoire il y a trois types de statistiques utilisées avec les polygones de Voronoï:

Des fonctions de lissage qui amenuisent les fortes
variations entre polygones voisins, de manière à mieux voir les
phénomènes globaux sur la zone d’étude.

En calculant la moyenne de chaque polygone avec ses voisins, la variation entre voisins devient moins brutale. On obtient alors une carte plus lissée des données. Ceci est utile quand il y a une trop forte variation entre les points voisins ce qui a pour résultat sur la carte globale de masquer ou de rendre plus difficile de voir les phénomènes globaux.

-Des fonctions pour mettre en évidence la variabilité locale. Si les outils de lissage (moyenne, mode, médiane) sont des outils qui s’intéressent à ce que l’on peut appeler la tendance centrale d’une distribution, les outils comme l’entropie, les écarts-type et les écarts inter-quantile s’intéressent à la dispersion des distributions.

L’interprétation de ces types de carte dépendent de la connaissance des données, car la variabilité sera toujours exprimée en cinq classes, avec des bornes différentes selon les données. Par contre, la carte d’entropie ne se présente pas de la même manière. Elle a toujours 5 classes mais les bornes des classes ne dépendent pas des données traitées. Elles sont fixes, la valeur d’entropie prenant une valeur comprise entre 0 et 2.3219 .

  • Des fonctions de recherche de valeurs anormales, telle que le cluster. Il est important d’identifier les valeurs aberrantes pour deux raisons: elles peuvent être des anomalies réelles du phénomène, ou la valeur peut avoir été mesurée ou enregistrée incorrectement.

Les valeurs d’entropie fournissent une mesure de dissimilarité entre polygones voisins. Dans la nature, vous vous attendez à ce que les choses plus rapprochées soient plus  semblables que des choses plus éloignées. Par conséquent, les valeurs aberrantes locales peuvent être identifiées par des zones d’entropie élevée.

La méthode cluster identifie les polygones qui sont dissemblables à leurs voisins environnants. Vous vous attendez à ce que la valeur enregistrée dans un polygone particulier soit similaire à au moins un de ses voisins. Par conséquent, cet outil peut être utilisé pour identifier d’éventuels valeurs aberrantes locales. Le cluster prend l’ensemble des points et classe toutes les valeurs en cinq classes. Pour chaque polygone on affiche sa classe, si et seulement si, au moins un polygone voisin est de même classe. Si tous les polygones voisins sont de classe différente, on affiche le polygone en gris.

Création des polygones de Voronoï avec QGis 3.8

Dans les traitements de QGis 3.8 vous avez plusieurs algorithmes permettant de créer les polygones de Voronoï à partir d’une couche de points. Pour nos besoins, nous allons utiliser le traitement Polygones de Voronoï de la Géométrie Vectorielle de QGis.

Cet outil permet d’enregistrer la couche en sortie directement dans une base Postgresql/Postgis.

Le résultat du traitement s’affiche dans QGis:

Avec PgAdmin vous pourrez voir la définition de la table créée:

Vous aurez besoin des noms de champs correspondant :

  • à l’identifiant
  • à la géométrie
  • au paramètre numérique que vous voulez analyser

Attention aux noms des colonnes importées à partir de la table attributaire de la couche d’origine. Remarquez que les noms sont compris entre guillemets. Ces guillemets font partie du nom. Vous devrez donc les ajouter dans l’appel de la fonction:

Dans cet exemple, si vous rentrez comme nom de l’attribut ‘DEBIT_MOY’ vous aurez une erreur disant que la relations DEBIT_MOY n’existe pas.

Il faut rentrer ‘  » DEBIT_MOY  » ‘ pour que ça fonctionne.

La fonction PL/pgsql

La fonction ci-dessous prend une couche de polygones de Voronoï et calcule trois statistiques locales.

Les quatre paramètres dans l’appel de la fonction sont:

  • le nom de la table contenant les polygones de Voronoï
  • le nom du champ identifiant de cette table
  • le nom de l’attribut numérique à analyser
  • le nom du champ géoémtrie de cette table

 En retour, la fonction crée trois attributs:

  • moyenne, qui contient la valeur de la moyenne locale (somme des valeurs des voisins et du polygone central, divisée par le nombre de voisins +1)
  • entrop, qui contient la valeur de l’entropie calculée avec la formule (voir l’article avec la description)

  • cluster, qui contient la valeur 99 si aucun des voisins appartient à la même classe que le polygone central ou la valeur de classe du polygone central si cette valeur est présente dans au moins un polygone voisin.

Code de la fonction
–Cette fonction s’applique à une table contenant déjà des polygones de voronoi et au moins un attribut numérique –dont on souhaite analyser sa distribution spatiale –Le résultat sont trois colonnes permettant de représenter — voronoi moyenne -> colonne « moyenne » — voronoi cluster -> colonne « cluster » — voronoi entropie -> colonne « entrop » –tablex est le nom de la table en entrée –gid est le nom de la colonne identifiant de la table –attr est le nom de la colonne (obligatoirement numérique) à traiter –‘||geom||’ est le nom du cham géométrie de la table –* CREATE OR REPLACE FUNCTION public.af_voronoi(tablex character varying, gid character varying,attr character varying,geom character varying) –* RETURNS text LANGUAGE ‘plpgsql’ — AS $BODY$ — DECLARE rec record; rec2 record; nb integer; L1 numeric; L2 numeric; L3 numeric; L4 numeric; L5 numeric; LL integer; Pi1 integer :=0; Pi2 integer:=0; Pi3 integer:=0; Pi4 integer:=0; Pi5 integer:=0; BEGIN –Création des colonnes de résultat et temporaires EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS smartcl’; EXECUTE ‘alter table ‘||tablex||’ add column smartcl integer’; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS entrop’; EXECUTE ‘alter table ‘||tablex||’ add column entrop float’; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS Nv’; EXECUTE ‘alter table ‘||tablex||’ add column Nv integer’; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P1′; EXECUTE ‘alter table ‘||tablex||’ add column P1 integer’ ; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P2′; EXECUTE ‘alter table ‘||tablex||’ add column P2 integer’ ; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P3′; EXECUTE ‘alter table ‘||tablex||’ add column P3 integer ‘; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P4′; EXECUTE ‘alter table ‘||tablex||’ add column P4 integer’ ; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P5′; EXECUTE ‘alter table ‘||tablex||’ add column P5 integer’ ; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS moyenne’; EXECUTE ‘alter table ‘||tablex||’ add column moyenne float ‘; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS cluster’; EXECUTE ‘alter table ‘||tablex||’ add column cluster integer’ ; EXECUTE ‘select count(*) as nb from ‘||tablex|| » INTO rec; — –calcul du voronoi « moyenne » EXECUTE ‘update ‘||tablex||’ vor set moyenne= (SELECT (sum(a.’||attr||’)::real+b.’||attr||’::real)/(count(a.’||attr||’)::real+1) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ group by b.’||gid||’)’; — –classement de l’attribut en cinq classes de même nombre d’échantillons — calcul des limites des classes EXECUTE ‘select max(‘||attr||’) as L1 from (select ‘||attr||’ from ‘||tablex||’ order by ‘||attr||’ asc limit ((2597/5)*1)::int)as a’ INTO rec; L1:=rec.L1; — EXECUTE ‘select max(‘||attr||’) as L2 from (select ‘||attr||’ from ‘||tablex||’ order by ‘||attr||’ asc limit ((2597/5)*2)::int)as a’ INTO rec; L2:=rec.L2; — EXECUTE ‘select max(‘||attr||’)as L3 from (select ‘||attr||’ from ‘||tablex||’ order by ‘||attr||’ asc limit ((2597/5)*3)::int)as a’ INTO rec; L3:=rec.L3; — EXECUTE ‘select max(‘||attr||’)as L4 from (select ‘||attr||’ from ‘||tablex||’ order by ‘||attr||’ asc limit ((2597/5)*4)::int)as a’ INTO rec; L4:=rec.L4; — EXECUTE ‘select max(‘||attr||’)as L5 from (select ‘||attr||’ from ‘||tablex||’ order by ‘||attr||’ asc limit ((2597/5)*5)::int)as a’ INTO rec; L5:=rec.L5; — –affectation d’une classe à chaque enregistrement de la table FOR rec IN EXECUTE ‘select ‘||gid||’ as gid,’||attr||’ as attr from ‘||tablex|| » LOOP — CASE WHEN rec.attr <=L1 THEN EXECUTE ‘update ‘||tablex||’ set smartcl=1 where ‘||gid||’=’||rec.gid; WHEN rec.attr <=L2 and rec.attr > L1 THEN EXECUTE ‘update ‘||tablex||’ set smartcl=2 where ‘||gid||’=’||rec.gid; WHEN rec.attr <=L3 and rec.attr > L2 THEN EXECUTE ‘update ‘||tablex||’ set smartcl=3 where ‘||gid||’=’||rec.gid; WHEN rec.attr <=L4 and rec.attr > L3 THEN EXECUTE ‘update ‘||tablex||’ set smartcl=4 where ‘||gid||’=’||rec.gid; ELSE EXECUTE ‘update ‘||tablex||’ set smartcl=5 where ‘||gid||’=’||rec.gid; END CASE ; END LOOP; –* –calcul du voronoi « entropie » –* EXECUTE ‘update ‘||tablex||’ vor set Nv= (SELECT (count(a.smartcl)) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ group by b.’||gid||’)’; EXECUTE ‘update ‘||tablex||’ vor set P1= (SELECT (count(a.smartcl)) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ and a.smartcl=1 group by b.’||gid||’)’; EXECUTE ‘update ‘||tablex||’ vor set P2= (SELECT (count(a.smartcl)) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ and a.smartcl=2 group by b.’||gid||’)’; EXECUTE ‘update ‘||tablex||’ vor set P3= (SELECT (count(a.smartcl)) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ and a.smartcl=3 group by b.’||gid||’)’; EXECUTE ‘update ‘||tablex||’ vor set P4= (SELECT (count(a.smartcl)) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ and a.smartcl=4 group by b.’||gid||’)’; EXECUTE ‘update ‘||tablex||’ vor set P5= (SELECT (count(a.smartcl)) FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ and a.smartcl=5 group by b.’||gid||’)’; EXECUTE ‘update ‘||tablex||’ set entrop=0′; EXECUTE ‘update ‘||tablex||’ set entrop= entrop +((P1::real/Nv::real)*(ln(P1::real/Nv::real))) where P1 IS NOT NULL’; EXECUTE ‘update ‘||tablex||’ set entrop= entrop+((P2::real/Nv::real)*(ln(P2::real/Nv::real))) where P2 IS NOT NULL’; EXECUTE ‘update ‘||tablex||’ set entrop= entrop+((P3::real/Nv::real)*(ln(P3::real/Nv::real))) where P3 IS NOT NULL’; EXECUTE ‘update ‘||tablex||’ set entrop= entrop+((P4::real/Nv::real)*(ln(P4::real/Nv::real))) where P4 IS NOT NULL’; EXECUTE ‘update ‘||tablex||’ set entrop= entrop+((P5::real/Nv::real)*(ln(P5::real/Nv::real))) where P5 IS NOT NULL’; EXECUTE ‘update ‘||tablex||’ set entrop=entrop*(-1) where entrop <>0′; — –calcul du voronoi « cluster » EXECUTE’update ‘||tablex||’ vor set cluster= (SELECT a.smartcl FROM ‘||tablex||’ as a JOIN ‘||tablex||’ as b ON st_touches(a.’||geom||’,b.’||geom||’) where b.’||gid||’ =vor.’||gid||’ and a.smartcl=b.smartcl group by b.’||gid||’,a.smartcl )’; — EXECUTE ‘update ‘||tablex||’ set cluster=99 where cluster IS NULL’; — –effacement des colonnes temporaires EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS Nv’; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P1′; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P2′; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P3′; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P4′; EXECUTE ‘alter table ‘||tablex||’ DROP COLUMN IF EXISTS P5′; RETURN ‘OK’; END; $BODY$;

Vous pouvez télécharger directement le fichier texte du code SQL en cliquant ici.

Représentation « standard »

Pour avoir une représentation de même type que celle de Geostatistical Analyst d’ArcGis, vous pouvez télécharger deux styles, un pour cluster et un autre pour l’entropie, en cliquant ici.

Vous aurez pour la variable cluster une représentation de type:

Et pour l’entropie, de type:

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *