Nous utilisons des cookies pour vous garantir la meilleure expérience sur notre site. Si vous continuez à utiliser ce dernier, nous considèrerons que vous acceptez l'utilisation des cookies. J'ai compris ! ou En savoir plus !.
banniere

Le portail francophone de la géomatique


Toujours pas inscrit ? Mot de passe oublié ?
Nom d'utilisateur    Mot de passe              Toujours pas inscrit ?   Mot de passe oublié ?

Annonce

Suite à un problème technique intervenu entre le 22 et le 23 mars, nous avons du procéder dans la soirée du 25 mars, à la restauration de la base de données du 24 mars (matinée).

En clair, nous avons perdu vos contributions et inscriptions du dimanche 24 et du lundi 25 mars.
Nous vous prions de nous excuser.

#1 Tue 03 November 2015 14:16

Librius
Participant actif
Lieu: Créteil
Date d'inscription: 9 Nov 2012
Messages: 67

[PostGis] Connaître les polygones qui intersectent mon polygone ..

Bonjour à tous,

Dans mon domaine : l'agriculture, j'ai un certain nombre de polygones classés comme suivants (1 table par ligne) :
  - France : 1 entité ;
  - Régions : 22 entités (pour le moment) ;
  - Départements : 96 ;
  - Polygones libres : 780 (suivant des courbes de relief mais pas toujours)
  - Parcelles agricoles : je n'ai plus le compte en tête mais je dois être dans l'ordre de grandeur de 750 000

Je dispose d'un certain nombre de parcelles agricoles et j'aimerais pour chacune d'entre elle récupérer l'information sur la ou les couches incluant partiellement ou totalement ma parcelle. J'ai mis un fichier joint pour présenter un exemple: je donne ma parcelle et je reçois l'information et les géométries (voir les découpages selon les zones).

Dans l'idée, je comptais faire une itération sur toutes les zones administratives + zones "polygones libres" mais cela me semble fastidieux.
1) Intersection "départements x parcelle" (je récupère à ce moment les régions associées) ; une parcelle pouvant être sur deux départements, il est possible qu'il y ait plusieurs sous-parcelles à ce premier niveau
2) Intersection "Polygones libres x résultat de l'intersection 1" mais je peux avoir des sous-parcelles qui se créent à chaque fois résultat d'intersection et il faut à chaque fois reprendre le résultat précédent.

Est-ce que cela vous semble cohérent ou connaissez une meilleure méthode pour arriver à cette finalité ? Il y aurait-il pas une astuce pour sélectionner déjà les polygones qui a priori matcheraient avec les contours de la parcelle (j'avais lu quelque chose sur les BBOX, je vais essayer de creuser de ce côté là ..)

En vous remerciant d'avance de toute l'aide que vous pourrez m'apporter .. smile
Bien cordialement,

Lib'


Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

Hors ligne

 

#2 Tue 03 November 2015 17:03

tumasgiu
Membre
Lieu: Ajaccio
Date d'inscription: 5 Jul 2010
Messages: 1129

Re: [PostGis] Connaître les polygones qui intersectent mon polygone ..

Vous pourriez faire quelquechose du genre :

Code:

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.

Hors ligne

 

#3 Wed 04 November 2015 00:02

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1534

Re: [PostGis] Connaître les polygones qui intersectent mon polygone ..

Bonsoir,

Le gros avantage de Postgis est de pouvoir faire de l'agregation spatiale facilement:
Ce n'est pas la peine de garder les couches france et région:
Faites les intersections sur les départements.
Pour avoir des info (attributs, géométries) au niveau de la region, faites un group by code_region.

Ensuite, le SQL va faire l'itération tout seul (il va se taper la partie fastidieuse).

Dans votre cas, vous voulez cumuler les informations d'intersection entre les parcelles et deux couches différentes:

Soit vous faites l'union des requetes faisant les intersections deux a deux, comme proposé par tumasgiu (mais ensuite, il faut refaire une requete d'agrégation pour ramener les infos à la parcelle), soit vous passez par une sous requete pour croiser les deux premieres couches, puis une requete pour croiser le résultat avec la derniere couche.

Pour cela, des index spatiaux sur vos 3 tables sont indispensables, sans quoi les requetes vont prendre des heures

Ensuite, pour agreger les résultats des intersections, vous pouvez utiliser des tableaux.
Ici, deux tableaux sont construits: dept_ids avec les identifiants des departements qui intersectent les parcelles et libre_ids contenant les identifiants des polygones libres qui intersectent les parcelles.

Code:

-- la premiere partie intersecte les parcelles avec les departements
with tmp as (
  select p.id as parc_id,
    array_agg(distinct d.gid order by d.gid) as dept_ids,
    p.geom
  from parcelle p left join departement d on st_intersects(p.geom, d.geom)
  group by p.id, p.geom

) select t.parc_id, -- la deuxieme partie croise les parcelles avec les polygones libres, en gardant les identifiants trouvés dans la partie précédente
    t.dept_ids,
    array_agg(distinct p.id order by p.id) as libre_ids
from tmp t left join pg_libre p on st_intersects(t.geom, p.geom)
  group by t.parc_id, t.dept_ids;

parc_id    dept_ids        libre_ids
1          {59,90}       {3}
2          {59}              {2}
3          {18,42,46}    {5,6}
4          {18,59}       {7}
5          {46,90}       {5}

(cf. image)

Dans la construction des tableaux (array_agg(...)) vous pouvez garder les identifiants identiques en enlevant le DISTINCT, si par exemple vous voulez savoir combien de fois un meme polygone intersecte une parcelle).

Vous pouvez aussi grouper tous les identifiants des couches dans un meme champ, suivant ce que vous voulez faire.

Nicolas


Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

Hors ligne

 

#4 Wed 04 November 2015 11:58

Librius
Participant actif
Lieu: Créteil
Date d'inscription: 9 Nov 2012
Messages: 67

Re: [PostGis] Connaître les polygones qui intersectent mon polygone ..

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 .. smile
Bonne journée,

Lib'


Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

Hors ligne

 

#5 Wed 04 November 2015 12:10

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1534

Re: [PostGis] Connaître les polygones qui intersectent mon polygone ..

Bonjour,

Oui, j'ai lu un peu rapidement votre message (c'est ca de répondre le soir...)

Vous voulez un peu plus que l'intersection en fait: il faut garder chaque bout de parcelle qui intersecte un objet: ca correspond a l'intersection parcelle/obj et a la difference parcelle/obj, ou encore au split de la parcelle par la linestring formant le polygone obj.

Je détaille les étapes dans un prochain message: une requete mettant en oeuvre quelques tricks postgis et postgreSQL bien utiles.

Nicolas

Hors ligne

 

#6 Wed 04 November 2015 13:06

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1534

Re: [PostGis] Connaître les polygones qui intersectent mon polygone ..

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:

Code:

 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.

Code:

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.

Code:

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


Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

Hors ligne

 

Pied de page des forums

Powered by FluxBB