Le portail francophone de la géomatique

Lucie D. · Mon 20 February 2017 10:17

Bonjour à tous,

Je travail sous postgresql 9.4, postgis 2.1.8, PG Admin III-1.20 et je souhaite découper une couche de polygones représentant le bâti d'un département avec une couche de polygones représentant l'occupation du sol.

Je voudrais affecter le type d'occupation du sol à chaque batiment là où il se situe.

- Première solution : requête de jointure du centroïde du bâtiment avec la géométrie du polygone d'occupation du sol...pas très pro ni précis car quid des bâtiments à cheval sur plusieurs polygones d'occupation du sol?

- Deuxième solution : je découpe mon bâtiment avec les objets de la couche d'occupation du sol, je calcule la superficie de chaque nouveaux petits bouts de bâtiments obtenus et je prends le type d'occupation du sol là où la superficie du petit bout de bâtiment est la plus importante par identifiant de bâtiment :

Code:

WITH tempa AS
(SELECT id_bati,max(surf_decoup) as s_max
FROM bati_decoup
GROUP BY id_bati)

SELECT id_bati,id_ocs
FROM bati_decoup,tempa
WHERE bati_decoup.id_bati=tempa.id_bati AND bati_decoup.surf_decoup=tempa.s_max

Cette requête suppose qu'un bâtiment n'est jamais à cheval exactement à 50%/50% sur deux polygones pour ne pas avoir de doublons d'identifiants id_bati.

Bon, le seul problème, c'est que je ne bloque pas sur cette requête que j'ai déjà effectué pour un autre travail mais sur la requête de découpage...

Voici ma requête :

Code:

DROP TABLE IF EXISTS bati_decoup;
CREATE TABLE bati_decoup
AS
SELECT bati.id as id_bati,ocs.id as id_ocs,st_intersection(bati.geom,ocs.geom) as geom_inter
FROM bati_decoup as bati,ocs
WHERE st_intersects(bati.geom,ocs.geom);

Le problème avec cette requête, c'est qu'après 4 jours, elle ne s'est toujours pas exécutée (je précise que ma couche d'OCS est uniquement sur l'agglo de Lyon tandis que mes bâtiments sont sur l'intégralité du département 69).

J'ai essayé la requête suivante:

Code:

DROP TABLE IF EXISTS bati_decoup;
CREATE TABLE bati_decoup
AS
SELECT bati.id as id_bati,ocs.id as id_ocs,st_intersection(bati.geom,ocs.geom) as geom_inter
FROM bati_decoup as bati
INNER JOIN ocs
ON st_intersects(bati.geom,ocs.geom) AND Not ST_IsEmpty(ST_Intersection(bati.geom,ocs.geom));

mais elle tourne elle aussi "dans le vide".

Bon, beaucoup me disent qu'il faut que je soit patiente mais le problème c'est que cette requête doit être reproductible sur l'ensemble du territoire nationale et je ne peux pas livrer une requête qui met une semaine à s'exécuter (surtout que ce découpage est loin d'être la seule requête de mon fichier .sql!Il s'agit en fait d'une étude complète sur le bâti).

Je précise que mes couches ont bien tous les index géo et attributaires nécessaires...

Avez-vous une idée pour optimiser ma requête?

Merci!!!

Lucie D.

Nicolas Ribot · Mon 20 February 2017 10:51

Bonjour,

Il faut regarder le plan de la requete pour savoir ce qu'elle fait et notamment si les index sont utilisés:
explain select...

Autre chose très importante pour ces cas d'intersection:
Il faut que les pg a intersecter soient petits par rapport à la zone totale, sinon l'index n'est pas très sélectif et un gros PG occ sol va se retrouver dans pleins d'intersections avec des petits batiments => long.

Avant, une solution consistait à créer une grille et à découper la couche avec de gros pg par rapport à cette grille. Maintenant, st_subdivide(geom) fait ca tout seul et TRES vite !

Si les PG occupation du sol sont gros, il faut créer une autre table occsol_cut en appliquant subdivide, puis index et utiliser cette couche.
Le gain de perf peut etre énorme.

Nicolas

Lucie D. · Mon 20 February 2017 11:23

Merci pour l'info!
Je vais tester comme ça, mais c'est vrai qu'il me semblait que les PG de cette couche d'occupation du sol étaient relativement petits (il s'agit de urban areas de l'agence européenne de l'environnement).

Mais en fait quand je regarde en zone périurbaine, les PG deviennent plus gros!

Je reviens sur le forum pour vous tenir au courant!

Merci!

Lucie D.

ChristopheV · Mon 20 February 2017 11:52

Bonjour,
Je vais peut-être dire une bêtise mais j'avais remarqué pour un sujet similaire qu'il est plus rapide de calculer l'intersection de deux multi-polygones.

L'utilisation de st_clusterwithin vous offrirait deux avantages :

Création d'un multi-polygone et possibilité en fonction de la distance de savoir si vous êtes dans du bâti dense (urbain) ou non (peri-urbain).
Comme on peut aussi récupérer le nombre de bâtiments dans une grappe on peut aussi gérer le bâti isolé.

L'idée est la suivante :

Je crée des grappes bâties (st_clusterwithin) avec une distance de 10 m, ie zone urbaine, je récupère la BBox de la grappe, dans une CTE je récupère l'occupation du sol de cette BBox et je calcule l'intersection des deux couches.

Lucie D. · Mon 20 February 2017 13:27

Bonjour ChristopheV,

Merci pour ce message qui tombe à pic car je suis justement en train de tester votre approche de la tache urbaine (https://georezo.net/forum/viewtopic.php?id=102129) qui est très intéressante! Je suis en phase test et du coup j'ai effectivement eu l'occasion d'utiliser ce fameux st_clusterwithin() pour la première fois! Bon le seul hic c'est que je travaille de plus en plus avec des fonctions qui nécessitent postgis 2.2 et je travaille sur un serveur distant avec un postgis 2.1.8...Donc j'ai installé un serveur postgresql sur mon propre PC avec une version 2.3.2 de postgis et comme ça je pourrai utiliser st_clusterwithin ET st_subdivide(geom) que je ne connaissais pas non plus mais qui ne fonctionne également qu'à partir de la 2.2...

Je vous tiens au courant!

Lucie D.

Dernière modification par Lucie D. (Mon 20 February 2017 14:10)

Lucie D. · Mon 20 February 2017 14:48

Nicolas,

Je suis en train de lancer la requête de découp du bati avec une couche d'occupation du sol sous_découpée avec st_subdivide(geom) [qui, au passage, s'exécute effectivement très rapidement].

Toutefois, comme vous me l'avez proposé, j'ai jeté un coup d’œil au plan de ma requête et en fait postgresql prend bien en compte mon index géographique sur ma couche bati, mais pas sur ma couche d'occupation du sol alors qu'elle a bien un index sur sa géométrie...?!

Est-ce normal? Je me suis demandé si ce n'était pas du au inner join. J'ai essayé de le remplacer avec un where mais le résultat est le même, seul l'index de ma couche bati est pris en compte... En revanche, quand je passe la souris sur la couche d'occupation du sol dans le plan de la requête, ça me dit : "seq scan" (voir pièce jointe).
Je ne sais pas ce que ça signifie (je ne maitrise pas les plans de requête ^^")?

Encore merci pour votre aide!

Bien cordialement,

Lucie D.

Dernière modification par Lucie D. (Mon 20 February 2017 14:48)

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

Lucie D. · Mon 20 February 2017 16:04

Et bien ma foi, le découpage grâce au st_subdivide(geom) de ma couche d'occupation du sol m'a effectivement permis de passer de plus de trois jours de requête à quelques minutes pour le découpage du bâti!

Je suis vraiment impressionnée! o_O

Merci beaucoup, du coup je dois légèrement adapter les requêtes qui suivent mais c'est parfait comme ça!

A bientôt pour de nouvelles aventures SQL!

Bien cordialement,

Lucie D.

Nicolas Ribot · Mon 20 February 2017 17:14

Ah cool, mieux en effet

Lucie D. a écrit:

Je suis vraiment impressionnée! o_O

Moi aussi: c'est vraiment un bel outil ce Postgis !

Lucie D. · Tue 30 May 2017 18:07

Bonjour,

Je reviens sur ce post car je me pose une question un peu bête.
J'avais lancé un st_subdivide(geom) sans préciser le nombre de sommets maxi...
Savez-vous quel est le nombre de sommets par défaut de la fonction quand on ne précise rien?

Merci par avance pour votre aide!

Bien cordialement,

Lucie D.

Nicolas Ribot · Tue 30 May 2017 18:36

http://postgis.net/docs/ST_Subdivide.html

(je vous invite à lire cette doc si vous utilisez postgis !)

NIcolas

Lucie D. · Wed 31 May 2017 09:20

Bonjour,

Merci.
Je m'attendais à cette réponse et c'est de ma faute, j'aurais du préciser que, bien évidemment j'ai lu et relu cette page postgis, mais soit je suis aveugle, soit mon niveau en anglais est trop bas, mais je ne trouve pas la réponse exacte à ma question.
Ok, le nombre minimum de sommets que l'on peut indiquer dans cette fonction est 8 sinon ça plante, mais ça ne me dit pas si c'est bien la valeur par défaut utilisée par st_subdivide quand on ne lui précise rien?

J'ai peut-être pas tout compris...?? =P

Merci beaucoup!

Lucie

ChristopheV · Wed 31 May 2017 09:53

Bonjour,

integer max_vertices=256

Cela signifie que la valeur par défaut de max_vertices est 256.

Lucie D. · Wed 31 May 2017 10:01

Bonjour,

Merci, c'est super, je ne l'avais pas du tout compris comme ça mais c'est super, merci beaucoup!

Bonne journée!

Lucie D.

ppluvinet · Wed 31 May 2017 10:19

Bonjour,
ravi de voir que ST_Subdivide a pu permettre une exécution plus rapide. C'est effectivement un réflexe qu'il faut avoir pour ce type de traitement !

Toutefois, je me suis interrogé sur ce croisement :
1. As-t-on besoin de faire un st_intersection : cela suppose qu'on peut avoir un même bâtiment dans plusieurs occupation différentes ? ca me semble conceptuellement peu approprié ?
2. Ne faut-il pas retirer les routes ? Ce sont souvent des polygones tentaculaires qui ne permettent pas d'utiliser les index correctement et un bâtiment ne peut pas appartenir à une route (à priori).

Bien sûr, je ne connais pas le sujet exact de travail de Lucie.

Bonne journée,

Dernière modification par ppluvinet (Wed 31 May 2017 10:22)

Lucie D. · Thu 01 June 2017 10:52

Bonjour Pascal,

Effectivement, il aurait été possible de retirer certains items de l'occupation du sol quand on est sûr que les bâtiments ne peuvent pas en faire partie...Et ça aurait été effectivement plus juste car je me retrouve malgré tout, avec ma méthode, avec quelques bâtiments dans l'eau (pas beaucoup heureusement)=> problèmes de précision géométrique? problème de méthode? Mais ça n'impacte pas vraiment mon travail.

En fait, je cherchais un moyen de récupérer le type d'occupation du sol dont la superficie de recouvrement par le bâti est la plus importante pour pouvoir affirmer que le bâtiment se trouve prioritairement sur telle ou telle occupation du sol (en cas de chevauchement sur plusieurs type d'OCS). J'ai donc fait des sommes de superficie de bouts de bâtiments par item OCS en prenant bien en compte le fait que ma couche d'OCS a été sous découpée par st_subdivide pour obtenir des objets avec moins de sommets (et donc doublons d'identifiants OCS).

Le résultat que j'ai obtenu avec ça est tout à fait cohérent, mais pour ce que je souhaitais faire, le simple croisement avec le centroide du bâtiment ne change quasiment pas les résultats à l'échelle d'un département et aurait suffit (j'ai fait le test).

Bonne journée!

Bien cordialement,

Lucie D.

ppluvinet · Thu 01 June 2017 12:31

Oui, j'aurais fait un simple update avec une condition d’intersection avec le centroïde de chaque bâtiment.
Attention tout de même avec l'urban areas de l'agence européenne de l'environnement. Elle est à prendre avec des pincettes sur certains postes de la nomenclature. Je me souviens d'avoir vu une carte de l'agglomération de Lyon qui était un peu surprenante ! Des zones commerciales un peu surreprésentées de souvenir...?
Bon courage,

ChristopheV · Thu 01 June 2017 14:28

Bonjour,
Nous ne sommes plus dans postgis mais il y a une chose qui me choque c'est cette notion

affirmer que le bâtiment se trouve prioritairement sur telle ou telle occupation du sol

Pour moi un bâtiment EST une occupation du sol. Il me semble que l'IGN considère que tout bâtiment est a minima une surface de 200 m² pour son ocs à grande échelle.

Avoir un bâtiment qui est situé sur plusieurs zone ocs c'est que ce bâtiment a été construit après l'OCS ou que l'OCS a été réalisée avec des données obsolètes.

Et nous nous retrouvons dans le cadre habituel, la problématique de la mise à jour des données à grande échelle ... bon il y a un organisme qui sait faire depuis 200 ans mais c'est pas à lui que l'on demande et c'est pas lui que l'on écoute, bon il faut dire que sa tutelle n'a pas vraiment envie.

Lucie D. · Thu 01 June 2017 16:04

Pour moi un bâtiment EST une occupation du sol.

J'en conviens totalement mais ça dépend quelle niveau de précision on souhaite. Urban areas a une nomenclature qui mélange occupation et usage des sols et avoir un bâtiment en zone agricole et un bâtiment dans une zone d'activités commerciales est pour moi différent.

Pour tout vous avouer :

1°/ je manipule des données que, il est vrai, je ne maitrise pas vraiment malgré mes recherches mais manque de temps et on fait tous un jour les braves exécutants ;-)

2°/ on ne sait pas si ces données vont être utilisées

3°/ Nous sommes en plein test alors merci pour vos conseils qui pourront peut-être faire pencher la balance...qui sait?

Lucie D.

ppluvinet · Thu 01 June 2017 16:05

Il y a différente définition/école et de règle appliquée à la production de cartographie d'occupation du sol. Le bureau d'études dans lequel je suis, travaille sur toutes ces notions. Selon les nomenclatures, un bâtiment peut être considéré à lui seul comme une occupation du sol, parfois ce sont les bâtiments et les espaces autour associés . Dans le second cas, on peut distinguer aisément par exemple les centres-ville et les zones pavillonnaires. D'ailleurs dans les zones pavillonnaires on a rarement des bâtiments de plus de 200m2, non ?
Par contre, un bâtiment est rarement à cheval sur 2 occupations du sol (quoique on a déjà vu des terrains de tennis sur des centre-commerciaux! )

Bon, on s'éloigne du sujet...

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Sites amis	Agenda Biblio Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

Annonce

#1 Mon 20 February 2017 10:17

Postgresql :Découpage couche batiments avec couche d'occupation du sol

Code:

Code:

Code:

#2 Mon 20 February 2017 10:51

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#3 Mon 20 February 2017 11:23

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#4 Mon 20 February 2017 11:52

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#5 Mon 20 February 2017 13:27

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#6 Mon 20 February 2017 14:48

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#7 Mon 20 February 2017 16:04

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#8 Mon 20 February 2017 17:14

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#9 Tue 30 May 2017 18:07

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#10 Tue 30 May 2017 18:36

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#11 Wed 31 May 2017 09:20

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#12 Wed 31 May 2017 09:53

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#13 Wed 31 May 2017 10:01

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#14 Wed 31 May 2017 10:19

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#15 Thu 01 June 2017 10:52

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#16 Thu 01 June 2017 12:31

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#17 Thu 01 June 2017 14:28

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#18 Thu 01 June 2017 16:04

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

#19 Thu 01 June 2017 16:05

Re: Postgresql :Découpage couche batiments avec couche d'occupation du sol

Pied de page des forums