Identification des polygones «longs et étroits» avec PostGIS

10

J'ai un ensemble de polygones représentant de grandes zones, disent les quartiers de la ville. Je veux identifier les grandes zones qui se chevauchent entre elles.

Mais il y a un problème: parfois ces polygones se chevauchent le long de leur périmètre (car ils ont été dessinés avec peu de précision). Cela va générer des chevauchements longs et étroits dont je ne me soucie pas.

Mais d'autres fois, il y aura de grands chevauchements de polygones robustes, ce qui signifie de grandes zones où le polygone d'un quartier en chevauche un autre. Je souhaite sélectionner uniquement ceux-ci.

Voir l'image ci-dessous de juste les chevauchements. Imaginez que je voulais sélectionner uniquement le polygone bleu dans le coin inférieur gauche.

chevauchement

Je pourrais regarder des zones, mais parfois les zones étroites sont si longues qu'elles finissent par avoir des zones aussi grandes que le polygone bleu. J'ai essayé de faire un rapport surface / périmètre, mais cela a également donné des résultats mitigés.

J'ai même essayé de l'utiliser ST_MinimumClearance, mais parfois les grandes zones auront une partie étroite attachée, ou deux sommets très proches.

Des idées d'autres approches?


En fin de compte, ce qui a le mieux fonctionné pour moi, c'était d'utiliser un tampon négatif, comme suggéré par @Cyril et @FGreg ci-dessous.

J'ai utilisé quelque chose comme:

ST_Area(ST_Buffer(geom, -10)) as neg_buffer_area

Dans mon cas, les unités étaient en mètres, donc 10 m de tampon négatif.

Pour les polygones étroits, cette zone a renvoyé zéro (également, la géométrie serait vide). Ensuite, j'ai utilisé cette colonne pour filtrer les polygones étroits.

bplmp
la source
4
Certes, le rapport surface / périmètre pourrait être utilisé pour cela.
Vince
Il est difficile de dire où se trouvent les polygones distincts de l'image, mais faire quelque chose comme ça gis.stackexchange.com/a/265233/64838 pourrait fonctionner? Calculez la zone de délimitation pivotée minimale, puis jetez celles de faible largeur ou hauteur.
FGreg
Vous pouvez également essayer d'utiliser un tampon négatif comme décrit ici: Comment puis-je identifier des polygones vraiment fins dans mon fichier de forme?
FGreg

Réponses:

5

J'essaierais de créer un tampon négatif, s'il mange des polygones minces, alors c'est bien, s'il ne mange pas le polygone, alors c'est le mien ... :-)

exécutez ce script, après avoir défini 2/3 de la largeur des polygones linéaires ...

create table name_table as
SELECT ST_Buffer(
(ST_Dump(
(ST_Union(
ST_Buffer(
(geom),-0.0001))))).geom,
0.0001)) as geom from source_table

OS: -) ...

Cyril Mikhalchenko
la source
à la fin, votre suggestion est celle qui a le mieux fonctionné pour moi. J'ai fini par utiliser quelque chose comme ST_Area(ST_Buffer(geom, -10)), le -10 étant -10 mètres dans mon cas. Si quelque chose renvoyait 0 de cette expression, je pourrais le filtrer.
bplmp
9

Au lieu de zone / périmètre, il est préférable d'utiliser la zone divisée par le carré du périmètre (ou son inverse).

Ceci est également appelé "indice de forme". Le carré du périmètre divisé par l'aire a une valeur minimale de 4 * Pi () (dans le cas d'un disque, qui est la géométrie 2D la plus compacte), il peut donc être normalisé par 4 * Pi () pour une interprétation (des valeurs normalisées proches de 1 signifient alors que vous avez des objets très compacts et que les carrés ont des valeurs d'environ 1,27).

EDIT: Un seuil sur la zone serait utile pour éliminer les très petits artefacts, qui pourraient être compacts. L'indice de forme montrerait alors un meilleur contraste. EDIT: en plus de cette réponse, l'utilisation de ST_Snap pourrait vous aider à résoudre le problème avant qu'il ne survienne.

radouxju
la source
Merci! Mais je ne sais pas comment ST_Snap pourrait aider dans ce cas ... Si j'ai bien compris, vous proposez quelque chose comme (o.overlap_perimeter^2 / o.overlap_area) / (4 * Pi()) as overlap_ratio? Cela donne de moins bons résultats que la surface / périmètre.
bplmp
Utilisant maintenant o.overlap_perimeter / (4 * sqrt(o.overlap_area)) as overlap_ratioselon cet article, mais des résultats encore pires (bien qu'il soit difficile de quantifier ce que j'entends par pire) isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/I-7/135/… , page 183.
bplmp
2
Merci pour cela, je n'avais jamais entendu parler de "l'indice de forme". J'avais toujours pensé que l'utilisation d'un rectangle de délimitation minimale était la meilleure façon de répondre à ce genre de question. J'ai trouvé ceci, repository.asu.edu/attachments/111230/content/… , ce qui est intéressant.
John Powell
@JohnPowell papier intersting, merci. Je vois que ce que je connais comme un indice de forme est appelé indice de circularité dans le papier. Mon problème avec les rectangles de délimitation minimum est qu'il ne fonctionne pas avec des objets très concaves (par exemple en forme de U)
radouxju
@bplmp ST_Snap vous aiderait à accrocher les sommets des polygones "presque" adjacents afin qu'ils ne se chevauchent plus. Il n'y a pas d'échelle sur vos figures, mais votre artefact ressemble à des lignes, donc je suppose que vous pouvez utiliser une valeur de tolérance qui suffit pour éviter les artefacts mais n'affecte pas les grands polygones.
radouxju
5

Une option consisterait à utiliser le rapport entre l'aire du polygone et la ligne la plus longue pouvant être tracée à l'aide de ses extrémités. Identification de longs polygones étroits.

select * from polygons where ST_Length(ST_LongestLine(geom, geom)) < ST_Area(geom) * 4

Cela fonctionne assez bien pour les polygones de ruban. Vous pouvez ajuster le ratio (ce avec quoi vous multipliez la surface) en fonction de vos besoins et de votre projection.

HeikkiVesanto
la source
1

Il semble que cela puisse correspondre à votre cas d'utilisation: éliminer les polygones sélectionnés

Combine les polygones sélectionnés de la couche d'entrée avec certains polygones adjacents en effaçant leur frontière commune. Le polygone adjacent peut être soit celui qui a la plus grande ou la plus petite surface, soit celui qui partage la plus grande frontière commune avec le polygone à éliminer.

L'élimination est normalement utilisée pour se débarrasser des polygones de ruban, c'est-à-dire de minuscules polygones résultant de processus d'intersection de polygones où les limites des entrées sont similaires mais pas identiques.

Il semble que vous souhaitiez essayer l'option "Plus grande limite commune".

FGreg
la source
Je me rends compte maintenant que vous demandiez des solutions postgis et non qgis. Mes excuses, je ne pense pas que postgis ait une fonction équivalente mais je laisse cela à la postérité.
FGreg
0

Cela me semble être un cas d'utilisation parfait pour l' extension de topologie PostGIS . Le paramètre de tolérance de la topologie déterminera dans quelle mesure vous autorisez les sommets à s'accrocher à d'autres polygones existants, à gérer la faible précision des données source et à les nettoyer.

En bref, la stratégie est la suivante:

1. Activez l'extension de topologie

CREATE EXTENSION postgis_topology;

2. Créez une nouvelle topologie vide

SELECT topology.CreateTopology('neighborhoods_topo', 4326, 1e-7);

Le troisième paramètre est la tolérance, dans les unités du CRS; choisissez-le judicieusement. Idéalement, vous voulez un CRS où l'unité est en mètres. Si l'unité CRS n'est pas en mètres, comme avec WGS 84 aka 4326, utilisez ST_Transformpour reprojeter vos polygones.

3. Ajoutez une colonne TopoGeometry à la table des polygones

SELECT topology.AddTopoGeometryColumn('neighborhoods_topo', 'public', 'neighborhoods', 'topogeom', 'POLYGON');

Cela renvoie un nouveau layer_id. Enregistrez-le, il sera nécessaire plus tard. Il sera superposé 1si vous partez de zéro et incrémenté à chaque nouvel appel.

4. Ajoutez tous les polygones dans la topologie

UPDATE public.neighborhoods
SET topogeom = topology.toTopoGeom(geom, 'neighborhoods_topo', 1, 1e-7);

Cela peut prendre plusieurs heures pour un grand ensemble de données, soyez patient. 1est le layer_id retourné plus tôt.

5. Trouvez des visages apparaissant dans plusieurs quartiers

Trouvez tous les visages de la topologie présents dans 2 topogéométries ou plus. Je vais laisser la requête comme un exercice. Le plus simple est probablement avec la GetTopoGeomElementsfonction, puis regroupez par identifiant de visage et regardez ceux avec un nombre de 2 ou plus. Alternativement, vous pouvez créer une nouvelle table avec la géométrie nettoyée à partir de la colonne topogeom, simplement la convertir en géométrie standard topogeom::geometryet répéter ce que vous avez déjà, mais maintenant avec un jeu de données propre sans que le ruban ne se chevauche.

François
la source