Géométries découpées et détection des intersections #31

JulienCorny · 2023-07-27T08:38:09Z

Contexte

Lorsqu’on créé une nouvelle ZH (depuis l’interface de l’application), si son tracé chevauche une ZH existante, le tracé de la nouvelle ZH est automatiquement découpé pour qu’il corresponde aux limites de la ZH existante. Les 2 polygones doivent donc au final avoir un ou plusieurs côté en commun (voir exemple ci-dessous).

Donc dans l’application, lorsqu’on retourne dans le formulaire de la nouvelle ZH (onglet Carte), si on ne touche pas au tracé, on devrait pouvoir passer à l’onglet suivant sans avoir le message : « La géométrie a été découpée car elle intersectait une autre zone humide. ».

Description du problème

Mais comme repéré par Colas Geier dans #30 :

" Lorsque la géométrie est découpée et enregistrée dans la base, elle continue d'intersecter les géométries qui l'entourent.
Le message d'enregistrement indiquera toujours La géométrie a été découpée car elle intersectait une autre zone humide."

En fait, ce message est affiché parce que la fonction PostGIS ST_Intersects détecte une intersection même lorsque la ZH est déjà découpée. Ce qui est intéressant (ou problématique ! ) c’est que les ST_Intersects ne détecte pas systématiquement une intersection car il doit probablement se produire une imprécision au moment du découpage des ZH (avec ST_Difference), ce qui a pour conséquences plusieurs cas de figures représentés ci-dessous et qui se présentent aléatoirement selon les découpages :

cas 1 : découpage attendu, une ligne commune (type LineString).
cas 2 : le découpage provoque un positionnement légèrement décalé (dans ce cas ST_Intersects=False)
cas 3 : le découpage provoque une zone d’intersection de type Polygon (ST_Intersects=True).

Pistes de résolution

Il faudrait donc trouver une solution pour réussir à détecter et différencier les cas où il s’agit de 2 ZH qui s’intersectent vraiment (avant découpage) de ceux où il s’agit d’une limite commune (déjà découpées).

A noter que :

des alternatives PostGIS à ST_Intersects comme ST_Touches pourraient fonctionner dans le cas 1 mais pas dans les cas 2 et 3.
dans le cas 3, PostgreSQL interprète les intersections comme un type LineString, alors que SQLAlchemy détecte bien un polygone à partir des mêmes coordonnées exactement. Il doit y avoir une différence d'implémentation de ST_Intersection dans PostgreSQL et SQLAlchemy ? D’ailleurs, toujours avec exactement les mêmes coordonnées, SQL Server interprète aussi le ST_Intersection comme un polygone. Voir ci-dessous un test fait sur PostgreSQL, SQLAlchemy et SQL Server. Donc il n'est pas possible non plus de détecter une ZH vraiment découpée ou non à partir du type d'intersection.

Une solution testée qui semble bien fonctionner est de calculer le pourcentage que représente la zone d’intersection par rapport à la zone totale de chaque ZH intersectée. On peut considérer que si le pourcentage est inférieur à 0,001 % alors il s’agit d’une ZH déjà découpée. Qu'en pensez-vous ?

Une autre piste serait de trouver en amont une alternative au découpage des ZH (autre que ST_Difference de PostGIS) pour améliorer la précision mais je ne suis pas sûr que ce soit possible.

Exemple de différence d’interprétation d’une intersection entre PostgreSQL, SQLAlchemy et SQL Server :

PostgreSQL :

SELECT st_geometrytype(ST_Intersection(
ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
,
'MULTIPOLYGON (((3.25007 45.546124, 3.247643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry
));

=> on obtient un type LineString

SQLAlchemy :

DB.session.query(func.st_geometrytype(func.ST_Intersection(polygon1, polygon2))).scalar()

=> polygon1 et 2 sont exactement les même WKT qu'au-dessus pour PostgreSQL. On obtient un type Polygon

SQL Server :

DECLARE @g1 geometry = 'POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))'

DECLARE @G2 geometry= 'POLYGON ((3.25007 45.546124, 3.247643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124))'

DECLARE @g geometry;
SET @g = @G2.STIntersection(@g1)

SELECT @g1
UNION ALL
SELECT @G2

SELECT @g.STGeometryType()

=> on obtient un type Polygon

A noter que pour les 3, que ce soit des polygones ou multipolygones, ou mélange des 2, on obtient exactement le même résultat.

camillemonchicourt · 2023-07-27T14:47:50Z

Ce qui m'interroge le plus, c'est d'intersecter les ZH automatiquement entre elles.

Ça me semble très discutable de modifier une ZH existante quand on en crée une nouvelle.

Mais si on veut garder ce fonctionnement d'intersection automatique, peut-être regarder autre chose que le ST_intersects, par exemple le ST_Dwithin ou autre du genre qui ne poserait par de soucis avec les frontières communes.

JulienCorny · 2023-07-27T14:53:33Z

Salut Camille,

C'est l'inverse, on ne modifie pas l'existante, on découpe la nouvelle quand le tracé recouvre une partie d'une ZH existante.

cen-cgeier · 2023-07-28T09:25:45Z

Salut à tous,

ST_intersects renverra toujours True dès lors qu'à minima un sommet entre deux géométries sont communs.
Lorsque j'avais rencontré ce cas de figure dans le cadre de traitements de données, j'avais opté pour la théorie suivante, mêler ST_Intersects avec ST_Buffer dans la condition de filtrage.

PostgreSQL :

SELECT ST_Difference(new_geom,old_geom) 
WHERE ST_Intersects(ST_Buffer(new_geom,-0.00000001), old_geom)

SQLAlchemy :

q = (
DB.session.query(func.ST_Difference(new_geom, old_geom))
.filter(
    func.ST_Intersects(func.ST_Buffer(new_geom,-0.00000001), old_geom))
)

Cette solution devrait permettre d'identifier les Cas 3 non souhaitable et discriminer les cas 1 et cas 2.
Je ne sais pas si c'est une bonne solution, mais jusqu'à présent ça a fonctionné.

camillemonchicourt · 2023-07-28T10:23:11Z

A priori on peut faire "ST_Intersects() AND NOT ST_Touches()".

Vu sur https://gis.stackexchange.com/questions/273993/st-relate-intersect-ignore-shared-boundary

cen-cgeier · 2023-07-28T12:07:21Z

Je ne suis pas certain que ST_Touches() soit mieux adapté. Cette fonction revoit True dans des cas plutôt particulier.
J'ai testé plusieurs cas en repartant de l'exemple de @JulienCorny :

t1 : une frontière commune entre la geom1 et la geom2, mais pas de sommet commun (exemple d'origine proposé par @JulienCorny)
t2 : Aucun sommet commun entre les deux géométries, mais 1 sommet de geom2 situé à l’intérieur de la geom1
t3 : Deux sommets communs reliés par une frontière commune entre la geom1 et la geom2
t4 : Un sommet commun entre geom1 et geom2 mais recouvrement de geom1 sur geom2

Résultats obtenus :

id	geom1	geom2	st_touches	st_intersects
t1	POLYGON ((3.243927 45.544396, ...	MULTIPOLYGON (((3.25007 45.546124, ...	False	True
t2	POLYGON ((3.243927 45.544396, ...	MULTIPOLYGON (((3.25007 45.546124, ...	False	True
t3	POLYGON ((3.243927 45.544396, ...	MULTIPOLYGON (((3.25007 45.546124, ...	True	True
t4	POLYGON ((3.243927 45.544396, ...	MULTIPOLYGON (((3.25007 45.546124, ...	False	True

Code testé :

with t1 as 
(SELECT 
	'POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))'::geometry geom1,
	'MULTIPOLYGON (((3.25007 45.546124, 3.247643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry geom2,
	ST_Touches(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.25007 45.546124, 3.247643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry
	),
	ST_Intersects(
		ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry,
	'MULTIPOLYGON (((3.25007 45.546124, 3.247643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry)
), t2 as (
	SELECT
	'POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))'::geometry geom1,
	'MULTIPOLYGON (((3.25007 45.546124, 3.245643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry geom2,
	ST_Touches(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.25007 45.546124, 3.245643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry
	),
	ST_Intersects(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.25007 45.546124, 3.245643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry
	)
),t3 as (
	SELECT
	'POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))'::geometry geom1,
	'MULTIPOLYGON (((3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.2469716663316515 45.54457224155132)))'::geometry geom2,
	ST_Touches(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.2469716663316515 45.54457224155132)))'::geometry
	),
	ST_Intersects(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.2469716663316515 45.54457224155132)))'::geometry
	)
), t4 as (
SELECT 
	'POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))'::geometry geom1,
	'MULTIPOLYGON (((3.25007 45.546124, 3.242082 45.545568, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry geom2,
	ST_Touches(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.25007 45.546124, 3.242082 45.545568, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry
	),
	ST_Intersects(
	ST_Multi('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))')::geometry
	,
	'MULTIPOLYGON (((3.25007 45.546124, 3.242082 45.545568, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry
	)
)
select 't1' id, t1.* from t1
	union all 
SELECT 't2' id, t2.* from t2
	union all 
SELECT 't3' id, t3.* from t3
	union all 
SELECT 't4' id, t4.* from t4
;

camillemonchicourt · 2023-08-01T05:37:29Z

OK merci pour ces tests et retours sur ST_Touches.

Avec ces éléments, la solution la plus pertinente me semble celle du buffer négatif avant de vérifier l'intersection.

JulienCorny · 2023-08-01T07:55:18Z

En fait les différences de type de geométrie entre PostgreSQL, SQLAlchemy et SQLServer pour les intersections sont dues à des différences d’arrondis pour le traitement des coordonnées géographiques à un moment du process :

pour obtenir le même résultat avec PostgreSQL, il faut passer le string du WKT dans la fonction ST_AsText qui a par défaut le paramètre "maxdecimaldigits = 15" : là on obtient un type 'ST_Polygon' pour l'intersection comme pour SQLAlchemy et SQLServer. Si on passe maxdecimaldigits = 16, là on a bien un 'ST_LineString' :

SELECT ST_GeometryType(ST_Intersection(
ST_AsText('POLYGON ((3.243927 45.544396, 3.242082 45.545568, 3.242575 45.547957, 3.246328 45.547567, 3.246693 45.546139, 3.247232 45.544988, 3.2471258888892076 45.54481822222273, 3.2469716663316515 45.54457224155132, 3.246157 45.545162, 3.243927 45.544396))'::geometry, 16)::geometry
,
ST_AsText('MULTIPOLYGON (((3.25007 45.546124, 3.247643 45.545643, 3.246569 45.54393, 3.248631 45.5416, 3.249855 45.541495, 3.251874 45.5439, 3.25007 45.546124)))'::geometry, 16)::geometry
));

Donc finalement, le découpage des géométries avec ST_Difference ne pose pas de problème et est précis. C'est l'approximation de l'intersection calculée avec SQLAlchemy qui posait problème (arrondi des float en Python ??). Etant donné que l'implémentation de ST_AsText dans SQLAlchemy n'a pas le paramètre 'maxdecimaldigits', passer directement par un format binaire avec des WKB à la place des WKT permet de résoudre ce problème. Et à ce moment là une combinaison de ST_Intersects = True et ST_Touches = False permet de ne garder que les géométries pas encore découpées.

Et on a le même souci de précision pour la fonction ST_Intersection : si on ne modifie pas dans SQLAlchemy la valeur par défaut du paramètre gridSize = -1 (float), on obtient pour certaines intersections complexes un type d’intersection ‘ST_GeometryCollection’ incorrect (car mélange de polygones et lignes) alors qu'on devrait avoir un type ‘ST_MultiLineString’.

This was referenced Jul 27, 2023

fix (backend): save geom when is intersected #30

Merged

Pr30 review #32

Closed

camillemonchicourt mentioned this issue Aug 2, 2023

ref_geo point & cor area linear PnX-SI/RefGeo#12

Merged

camillemonchicourt added the enhancement New feature or request label Feb 22, 2024

edelclaux mentioned this issue May 21, 2024

Pr30 review #67

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Géométries découpées et détection des intersections #31

Géométries découpées et détection des intersections #31

JulienCorny commented Jul 27, 2023

camillemonchicourt commented Jul 27, 2023

JulienCorny commented Jul 27, 2023

cen-cgeier commented Jul 28, 2023

camillemonchicourt commented Jul 28, 2023 •

edited

Loading

cen-cgeier commented Jul 28, 2023

camillemonchicourt commented Aug 1, 2023

JulienCorny commented Aug 1, 2023

Géométries découpées et détection des intersections #31

Géométries découpées et détection des intersections #31

Comments

JulienCorny commented Jul 27, 2023

camillemonchicourt commented Jul 27, 2023

JulienCorny commented Jul 27, 2023

cen-cgeier commented Jul 28, 2023

camillemonchicourt commented Jul 28, 2023 • edited Loading

cen-cgeier commented Jul 28, 2023

camillemonchicourt commented Aug 1, 2023

JulienCorny commented Aug 1, 2023

camillemonchicourt commented Jul 28, 2023 •

edited

Loading