Le portail francophone de la géomatique

Vous pourriez faire quelquechose du genre :

SELECT agricole.id, st_intersection(agricole.geom,region.geom), (region)
WHERE st_intersects(agricole.geom,region.geom)
UNION
SELECT agricole.id, st_intersection(agricole.geom,departement.geom), (departement)
WHERE st_intersects(agricole.geom,departement.geom)
UNION
.
.
.

De cette façon, vous obtenez toutes les intersections de toutes les parcelles agricoles avec les différentes couches
et les informations associées.

Bonjour,

Merci à tous les deux pour vos réponses, je vais regarder ça attentivement dès que j'aurai un moment (je ne pensais pas avoir de réponses aussi précises si rapidement !!).

Je note cependant (au moins dans la deuxième réponse très détaillé de Nicolas) que le résultat ne semble pas correspondre totalement à mon objectif, cela est surement dû à une mauvaise explication de ma part. Le résultat que je souhaite n'est pas à la parcelle mais aux intersections de parcelles avec mes différents polygones (administratifs ou libre). Je peux donc avoir plusieurs résultats avec le même id de parcelle mais des id de polygones administratifs ou libres différents (ou sans id de polygone si jamais un bout ou la totalité de parcelle n'est pas compris dans une zone (cf pièce jointe) - les limites noires délimitant chaque entités souhaités en fin de traitement !

Je reviens vers vous rapidement dans tous les cas après le test de vos méthodes ..
Bonne journée,

Lib'

Il y a plusieurs façon de faire pour calculer toutes les parties des parcelles qui sont en intersection avec un autre polygone:
intersection et difference entre les objets, ou une nouvelle fonction de postgis 2.0: st_split.

Dans l'exemple proposé, le processus est le suivant:

• La première sous requete se base sur l'hypothese que les parcelles sont entièrement contenues dans les départements (pas de parcelles frontalieres...): l'intersection entre parcelle et département donne toujours la surface complète de la parcelle.
Si ce n'est pas le cas, il faut aussi utiliser st_split, comme dans la deuxième partie de la requête.

• La deuxième sous requete "splite" les morceaux de parcelles (intersection avec les département) avec le contour des polygones libres (st_split ne peut splitter un polygon (PG) qu'avec une linestring (LS). Utilisation de st_boundary pour renvoyer le contour extérieur d'un PG sous forme de LS).
L'avantage de st_split est de renvoyer une collection contenant directement les deux bouts découpés par le PG, sans avoir a faire une intersection et une différence entre les objets.

Dans st_split, un trick vraiment pratique pour gérer le cas suivant:
on veut à la fin lister TOUS les morceaux de parcelles, même ceux qui ne sont splités par aucun PG libre. Or pg_split, comme (presque) toutes les fonctions SQL, renvoie null si un des deux parametres est null, ce qui sera le cas pour un morceau de parcelle n'intersectant pas un PG libre.
la construction:

 coalesce(st_boundary(p.geom), st_setSRID('LINESTRING EMPTY'::geometry, 2154))

permet de transformer une valeur nulle en une géométrie vide, ce qui n'est pas du pareil: st_split(geom, geomVide) renvoie geom, exactement ce qu'on veut.

Dans cette partie, on utilise st_dump pour exploser le geometryCollection issue de st_split en objets individuels, c'est a dire en partie de parcelles découpées. (attention ici, suivant la topologie des données, st_split peut renvoyer autre chose que des surfaces: ca peut etre des points ou des lignes. Si vous n'etes intéressés que par les objets surfaciques, il faut utiliser st_collectionExtract avant st_dump pour en sortir les objets surfaciques.)

Là, suivant l'information que vous voulez, il y a plusieurs possibilités.
Dans la requête finale, on considère que l'info que vous voulez, c'est de savoir si une parcelle possède une surface en intersection avec un pg, et pas juste un bord. Or avec la partie st_split, si une parcelle est découpée en deux morceaux par un pg libre, cet id de pg libre est associé aux deux morceaux splités.


• Dans la troisième sous requête, on calcule ainsi la surface de l'intersection entre les bouts de parcelles et les pg libres qui les ont découpés (arrondie, histoire de contourner les problèmes de précision lors de la découpe): pour les bouts de parcelles qui ont une surface d'intersection = 0 , on sait que ces objets se touchent sur une frontiere mais ne s'intersectent pas.

Là encore on utilise coalesce avec une géométrie vide pour avoir tjs une valeur a st_area: 0 (pas d'intersection) ou la surface d'intersection.

• Enfin, dans la dernière partie, les attributs sont juste concaténés dans une colonne texte pour pouvoir afficher les id des polygones graphiquement.

row_number over () ...

permet de générer un id unique pour chaque partie de parcelle directement dans la requête.
La construction case when ... permet de tester des valeurs en SQL (sorte de IF), ce qui est super pratique, par ex pour formater notre chaine de texte comme on le veut.

Enfin, l'exemple suivant marche pour un cas simple: peu d'objets, simples: peu de chances d'avoir des erreurs topologiques lors des splits, intersections.
Dans la vraie vie (et avec 750 000 parcelles), je doute que ca se passe aussi bien, mais il existe pleins de solutions pour nettoyer, valider...

Pour des questions de perf, je ne vous conseille pas de faire cela en une seule requête, mais de générer des tables à chaque étape, que vous indexerez.
Ca permettra en outre de contrôler les différentes étapes si une se passe mal.

with tmp as (
  select p.id as parc_id,
    d.gid as dept_id,
    st_intersection(p.geom, d.geom) as inter_dept_geom
  from parcelle p left join departement d on st_intersects(p.geom, d.geom)
), tmp1  as (
    SELECT
      t.parc_id,
      t.dept_id,
      p.id                                                AS libre_id,
      st_dump(st_split(t.inter_dept_geom, coalesce(st_boundary(p.geom), st_setSRID('LINESTRING EMPTY'::geometry, 2154)))) AS dmp
    FROM tmp t LEFT JOIN pg_libre p ON st_intersects(t.inter_dept_geom, p.geom)
), tmp2 as (
    SELECT
      parc_id,
      dept_id,
      libre_id,
      round(st_area(st_intersection((dmp).geom, coalesce(p.geom, 'POLYGON EMPTY' :: GEOMETRY)))) AS area_inter,
      (dmp).geom
    FROM tmp1 t LEFT JOIN pg_libre p ON st_intersects((dmp).geom, p.geom)
) select row_number() over () as id,
    t.parc_id, 'did: '||t.dept_id||
  case when t.area_inter = 0 then '' else ' lid: '||t.libre_id end as ids,
  t.geom
from tmp2 t;

Nicolas