Le portail francophone de la géomatique

Bonjour,

nous avons peut-être quelques éléments de réponse à vous apporter dans le cadre de la mise à jour de notre base de données occupation du sol.

Je tiens déjà à remercier Mathieu Bossaert, Responsable du système d'information du CEN LR, qui par son aide précieuse nous a permis d'obtenir les résultats attendus. Les éléments de réponse ci-dessous ont été rédigés par ma collègue Claire Séard, actuellement en formation.


L’objectif était à partir de deux couches d’occupation du sol de deux dates différentes de pouvoir représenter spatialement chaque changement d’occupation du sol en faisant les traitements sous Postgresql/PostGIS.

Les évolutions possibles entre les deux couches sont les suivantes :
- polygones de géométrie identique,
- fusion de polygones,
- découpage de polygones,
- remodelage complet des polygones.

Pour identifier spatialement chaque changement à la fois géométrique et attributaire, nous avons voulu découper les polygones des 2 couches l’une par rapport à l’autre pour obtenir autant de polygones qu’il y a de changements possibles.

Dans un premier temps, nous avons tenté de lier les deux couches d’occupation du sol avec des requêtes spatiales comme ST_intersects, ST_Overlaps en combinaison avec une modification des géométries avec les fonctions ST_Intersection et ST_Difference pour découper les polygones. Les résultats n’ont pas été concluants car la couche résultat possède différents types de géométrie : points, lignes, multi-lignes, polygones, multi-polygones et collection de géométrie, ce qui biaise ensuite les traitements sur cette couche. De plus, la couche résultat possède des zones vides sans polygones, notamment là où il y a des fusions de polygones.

Dans un second temps, pour éviter de découper les polygones entre eux, la réflexion s’est tournée sur la transformation des polygones en lignes. Cela permet d’obtenir la combinaison des limites des polygones de chacune des couches qui seront ensuite re-transformées en polygone.

La méthodologie, indiquée dans la requête ci-dessous, se décompose en plusieurs étapes de traitement :
- une première étape consiste donc à transformer les contours des polygones en polylignes avec la fonction ST_Boundary : 1ère sous-requête ci-dessous « polygones_ini » ;
- puis les lignes de l’association (UNION) de ces 2 couches sont unies avec la fonction ST_Union  pour éliminer les doublons : 2e sous-requête « limites » ;
- les polygones possibles sont ensuite crées à partir de ces lignes avec la fonction ST_Polygonize. Pour récupérer des polygones individuels, il faut combiner cette fonction avec ST_Dump : 3e  sous-requête « découpage ».

CREATE TABLE ocsol.evolution AS 
WITH polygones_ini as ( --- transformation des polygones en lignes
    SELECT st_boundary(the_geom) as geometrie
    FROM ocsol.ocsolannee1
    UNION
    SELECT st_boundary(the_geom) as geometrie
    FROM ocsol.ocsolannee2
)
, limites AS (--- union des lignes
    SELECT st_union(geometrie) as geometrie
    FROM polygones_ini 
)
, decoupage AS ( ---création des polygones
SELECT ((st_dump(ST_Polygonize(geometrie))).geom)::geometry(POLYGON,3943) AS geometrie 
---à rajouter pour verifier le type de géométrie : 
---,st_geometrytype((st_dump(ST_Polygonize(geometrie))).geom) as type_geom
FROM limites)

SELECT row_number() OVER() AS gid, geometrie
FROM decoupage;

Le résultat est une couche de polygones sans valeur attributaire avec un identifiant gid et une géométrie. Il est ensuite nécessaire d’attribuer à chacun de ces polygones les attributs des 2 couches d’occupation du sol pour évaluer les changements. Il faut rajouter des colonnes pour ces attributs et les mettre à jour en intersectant la couche obtenue avec les deux couches d’occupation du sol comparées. Étant donné la complexité des recoupements de géométries, une simple intersection des polygones entre eux ne suffit pas. La fonction ST_pointonsurface a donc été utilisée pour créer un point à l’intérieur de chaque polygone de notre couche résultat qui est ensuite intersecté avec chacune des 2 couches d’occupation du sol pour récupérer leurs attributs.

UPDATE ocsol.evolution SET (attribut1_annee1, attribut2_annee1) = (SELECT a. attribut1_annee1, a. attribut2_annee1 FROM ocsol.ocsolannee1 a WHERE ST_intersects(ST_pointonsurface(evolution.geometrie), a.the_geom));

UPDATE ocsol.evolution SET (attribut1_annee2, attribut2_annee2) = (SELECT a. attribut1_annee2, a. attribut2_annee2 FROM ocsol.ocsolannee2 a WHERE ST_intersects(ST_pointonsurface(evolution.geometrie), a.the_geom));

J'espère que ces éléments vous permettront d'apporter de l'eau à votre moulin,

Bonne journée.