Le portail francophone de la géomatique

white-shadow90 · Wed 18 October 2017 19:10

Bonjour,

Je cherche le moyen de tracer une ligne qui coupe un angle en 2 au niveau d'un sommet d'une polyligne. L'idée s'apparente à la création d'une bissectrice, mais je rencontre un problème pour les angles compris entre 179 et 180° donc j'ai envisagé une autre méthode qui me permet de récupérer les géométries des 2 segments qui composent l'angle sous forme de texte, ainsi que la géométrie du point de jonction de ces 2 segments et l'angle formé par les 2 segments en ce point.

Code:

CREATE TABLE public.p_sommets AS
SELECT
row_number() OVER() as id_sommet,
t1.id_angle,
t1.id_ilot,
t1.vertex,
t1.sp as angle_geom,
t1.ang,
t3.x,
t3.y,
t1.segment1,
t1.segment2,
t3.idu
FROM(
select 
        row_number() OVER () as id_angle,
        id_ilot,
       point_order as vertex,
       sp,
       --
       case when point_order = 1
         then last_value(ST_Astext(ST_Makeline(sp,ep))) over (partition by id_ilot)
         else lag(ST_Astext(ST_Makeline(sp,ep)),1) over (partition by id_ilot order by point_order)  end as segment1,
        ST_Astext(ST_Makeline(sp,ep)) as segment2,
       --
       abs(abs(
       case when point_order = 1
         then last_value(degrees(ST_Azimuth(sp,ep))) over (partition by id_ilot)
         else lag(degrees(ST_Azimuth(sp,ep)),1) over (partition by id_ilot order by point_order)
       end - degrees(ST_Azimuth(sp,ep))) -180 ) as ang
from (-- 2.- extract the endpoints for every 2-point line segment for each linestring
      --     Group polygons from multipolygon
      select id_ilot,
             generate_series(1, ST_Npoints(geom)-1) as point_order,
             ST_Pointn(geom, generate_series(1, ST_Npoints(geom)-1)) as sp,
             ST_Pointn(geom, generate_series(2, ST_Npoints(geom)  )) as ep
      from ( -- 1.- Extract the individual linestrings and the Polygon number for later identification
             select id_ilot,
                    (ST_Dump(ST_Boundary(geom))).geom as geom
              from cadastre.s_ilots ) as pointlist ) as segments) t1
INNER JOIN (
SELECT * FROM (
WITH t AS -- Transfor polygons in sets of points
    (SELECT idu,
            st_dumppoints(geom) AS dump
     FROM cadastre_qgis.geo_parcelle),
f AS -- Get the geometry and the indexes from the sets of points 
    (SELECT t.idu,
           (t.dump).path[1] AS part,
           (t.dump).path[3] AS vertex,
           (t.dump).geom AS geom
     FROM t)
-- Get all points filtering the last point for each geometry part
SELECT row_number() OVER () AS gid, -- Creating a unique id
       f.idu,
       f.part,
       f.vertex,
       ST_X(f.geom) as x, -- Get point's X coordinate
       ST_Y(f.geom) as y, -- Get point's Y coordinate
       f.geom::geometry('POINT',3949) as geom -- make sure of the resulting geometry type
FROM f 
WHERE (f.idu, f.part, f.vertex) NOT IN
      (SELECT f.idu,
              f.part,
              max(f.vertex) AS max
       FROM f
       GROUP BY f.idu,
                f.part)) t2) t3
ON ST_EQUALS(t1.sp, t3.geom)
WHERE ST_EQUALS(t1.sp, t3.geom) AND t1.sp && t3.geom
ORDER BY t3.idu

A partir de là, je ne vois pas comment écrire la requête qui permette de tracer une droite qui coupe l'angle en 2 partie égales, sur une distance donnée.

Je pense qu'il existe une solution avec ST_Rotate en demandant à un segment (comment l'identifier ?) de tourner (dans le sens des aiguilles) de la moitié de la valeur de l'angle, à partir de cette angle.
Mais j'imagine qu'il existe une autre solution (ou fonction) qui permette de couper facilement un angle en 2 parties égales.
Pour la méthode de la bissectrice, j'avais suivi la démarche suivante :
- Création d'un buffer à partir des sommets
- Récupération du point d'intersection (p1 et p2) entre le buffers et les segments qui composent l'angle
- Nouveau buffer de même rayon que le précédent, tracé depuis p1 et p2
- Récupération des points d'intersection (p3 et p4) de ces deux buffers
- Traçage d'une ligne entre p3 et p4
Mais comme indiqué précédemment, pour les angles plats, je n'ai pas d'intersection, au mieux une tangente (qui me convient) mais pour un certain nombre de cas, je n'ai pas de géométrie issue de ce calcul d'intersection.

Aussi, je suis preneur de n'importe quelle méthode me permettant de récupérer cette droite qui coupe mon angle en 2 parties égales.
Je vous remercie pour votre collaboration.

Post en semi doublon sur le forum SIG : http://www.forumsig.org/showthread.php/ … issectrice

Dernière modification par white-shadow90 (Wed 18 October 2017 19:12)

tumasgiu · Thu 19 October 2017 17:22

Salut,

très laborieusement:

Code:

SELECT
row_number() OVER () id,
st_multi(
st_union(ARRAY[s1,s2,
CASE WHEN st_pointn(bisect, 1)  && st_pointn(bisect, 2)
THEN st_makeline(st_pointn(s1, 1),
                  st_rotate(st_pointn(st_segmentize(s1, 1), 2), radians(90), st_pointn(s1, 1)))
ELSE bisect
END])) g
FROM
    (SELECT
    st_makeline(
        ST_LineInterpolatePoint(
            st_makeline(
                St_pointn(ST_Segmentize(s1, 1), 2),
                St_pointn(ST_Segmentize(s2, 1), 2)),
            0.5),
        St_intersection(s1, s2)) bisect
    , s1, s2
    FROM (SELECT
            CASE WHEN st_intersects(st_pointn(s1, 1), s2) 
            THEN s1
            ELSE st_reverse(s1) END s1,
            CASE WHEN st_intersects(st_pointn(s2, 1), s1) 
            THEN s2
            ELSE st_reverse(s2) END s2
            FROM ( SELECT
                        st_makeline('POINT(0 0)'::geometry,
                                          st_rotate('POINT(0 10)'::geometry, radians(r))) s1
                        ,st_makeline(st_rotate('POINT(0 10)'::geometry, radians(180+r)),
                                           'POINT(0 0)'::geometry) s2
                        FROM 
                        generate_series(0,  360, 20) r
                        union
                        SELECT
                        st_makeline('POINT(20 0)'::geometry,
                                          'POINT(20 10)'::geometry) s1
                        ,st_makeline(st_rotate('POINT(20 10)'::geometry, radians(r), 'POINT(20 0)'::geometry),
                                                        'POINT(20 0)'::geometry) s2
                        FROM 
                        generate_series(0,  360, 10) r ) foo
                        ) sample
        ) clean

EDIT la derniere sous requête (sample) génère simplement des lignes de test.

Dernière modification par tumasgiu (Thu 19 October 2017 17:45)

white-shadow90 · Mon 23 October 2017 10:24

Bonjour,

Je vous remercie. Je regarde ça de près et vous tiens au courant.

white-shadow90 · Mon 23 October 2017 16:25

Bonjour,

Je ne suis pas certain d'avoir compris toute la requête. Voici le code modifié lorsque je l'exécute à partir des segments obtenus dans ma table p_sommets.

Code:

CREATE TABLE public.test_bissectrice AS SELECT
row_number() OVER () id,
st_multi(
st_union(ARRAY[s1,s2,
CASE WHEN st_pointn(bisect, 1)  && st_pointn(bisect, 2)
THEN st_makeline(st_pointn(s1, 1),
                  st_rotate(st_pointn(st_segmentize(s1, 1), 2), radians(90), st_pointn(s1, 1)))
ELSE bisect
END])) g
FROM
    (SELECT
    st_makeline(
        ST_LineInterpolatePoint(
            st_makeline(
                St_pointn(ST_Segmentize(s1, 1), 2),
                St_pointn(ST_Segmentize(s2, 1), 2)),
            0.5),
        St_intersection(s1, s2)) bisect
    , s1, s2
    FROM (SELECT
            CASE WHEN st_intersects(st_pointn(s1, 1), s2) 
            THEN s1
            ELSE st_reverse(s1) END s1,
            CASE WHEN st_intersects(st_pointn(s2, 1), s1) 
            THEN s2
            ELSE st_reverse(s2) END s2
            FROM (  SELECT
        ST_GeomFromText(segment1) as s1,
        ST_GeomFromText(segment2) as s2
                   FROM public.p_sommets) foo) foo ) clean

Pour l'instant, je rencontre une erreur m'indiquant "Input geometries must be points or lines" et cela même lorsque je remplace le ST_GeomFromText par ST_MakeLine. Pourtant, quand je fais un test en créant mes segments, le type de géométrie des segments est ST_LineString

tumasgiu · Tue 24 October 2017 13:31

La fonction qui fait planter la requête est probablement st_intersection:
Il semblerait que l'intersection des segments de certains enregistrements
de votre table sommet ne soit ni une ligne ni un point, ce qui empêche
le st_makeline englobant de constituer la bissectrice.

Vous pouvez poser des questions sur les aspects de la requête que
vous ne comprenez pas bien.

white-shadow90 · Wed 25 October 2017 10:56

Bonjour,

Je ne comprends pas pourquoi le radians a une valeur de 90. j'imaginais plutôt une valeur qui serait variable.

EDIT : J'ai regardé l'intersection et effectivement j'avais des "GeometryCollection"

J'ai donc modifié la requête comme suit :

Code:

CREATE TABLE public.test_bissectrice AS SELECT
row_number() OVER () id,
st_multi(
st_union(ARRAY[s1,s2,
CASE WHEN st_pointn(bisect, 1)  && st_pointn(bisect, 2)
THEN st_makeline(st_pointn(s1, 1),
                  st_rotate(st_pointn(st_segmentize(s1, 1), 2), radians(90), st_pointn(s1, 1)))
ELSE bisect
END])) g
FROM
    (SELECT
    st_makeline(
        ST_LineInterpolatePoint(
            st_makeline(
                St_pointn(ST_Segmentize(s1, 1), 2),
                St_pointn(ST_Segmentize(s2, 1), 2)),
            0.5),

CASE
WHEN ST_GeometryType(St_intersection(s1, s2)) = 'ST_GeometryCollection' THEN CASE
        WHEN ST_GeometryType(ST_CollectionExtract(St_intersection(s1, s2), 1)) = 'ST_MultiPoint' THEN geom(ST_DUMP(ST_CollectionExtract(St_intersection(s1, s2), 1)))
        WHEN ST_GeometryType(ST_CollectionExtract(St_intersection(s1, s2), 2)) = 'ST_MultiLineString'  THEN geom(ST_DumpPoints(St_intersection(s1, s2))) END
    ELSE St_intersection(s1, s2) END ) as bisect,
     s1, s2
    FROM (SELECT
            CASE WHEN st_intersects(st_pointn(s1, 1), s2) 
            THEN s1
            ELSE st_reverse(s1) END s1,
            CASE WHEN st_intersects(st_pointn(s2, 1), s1) 
            THEN s2
            ELSE st_reverse(s2) END s2
            FROM (  SELECT
        ST_GeomFromText(segment1) as s1,
        ST_GeomFromText(segment2) as s2
                   FROM public.p_sommets) foo) foo ) clean

Elle s'exécute correctement. Je regarde le résultat et vous tiens informé

Dernière modification par white-shadow90 (Wed 25 October 2017 11:50)

tumasgiu · Wed 25 October 2017 12:00

Le st_rotate n'intervient qu'en cas d'angle plat.
Comme la valeur d'un angle plat est égal à 180°,
sa bissectrice sépare l'angle en deux parts égales,
180/2 = 90,
on "crée" la bissectrice en donnant un quart de tour
à l'un de deux segment qui compose l'angle par rapport
à l'intersection des deux segments.
Ce qui varie ce n'est pas le nombre d'unité
à donner à la rotation, mais son centre.

Pour vous représenter la chose, imaginez comment vous utiliseriez
un rapporteur pour dessiner la bissectrice d'un angle plat.

white-shadow90 · Wed 25 October 2017 12:49

OK, je vois, je n'avais pas procédé de la même manière que vous lorsque j'avais essayé de créer une bissectrice de la même manière qu'avec un compas.
Mon calcul de bissectrice posait problème pour des valeurs d'angles qui n'étaient pas forcément égales à 180, ces angles pouvaient être de 179,9 voire peut-être inférieurs. Je n'avais pas regardé la valeur minimale de l'angle pour lequel je n'arrivais pas à récupérer 2 points d'intersectionspour mes buffers de 2e génération.

J'avais procédé comme suit à partir de ma table p_sommets obtenue par la requête indiquée au début de ce fil de discussion :

Code:

-- Création de la couche des point d'origine de la bissectrice

CREATE TABLE public.p_obis AS                
SELECT
row_number() OVER () AS id_origine,
a.id_sommet,
a.idu as idu_1,
b.idu as idu_2,
a.x,
a.y,
a.angle_geom,
a.ang,
a.id_angle,
a.id_ilot,
a.segment1,
a.segment2
FROM 
public.p_sommets as a
INNER JOIN
public.p_sommets as b
ON (a.x = b.x and a.y = b.y)
WHERE a.idu != b.idu;


-- Défnition des SRID de p_obis

ALTER TABLE public.p_obis
ALTER COLUMN angle_geom TYPE geometry(point, 3949)
USING ST_SetSRID(angle_geom,3949);
          
              
-- Création de la couche des segments

CREATE TABLE public.l_segments AS                
SELECT
id_origine,
id_sommet,
id_angle,
ST_GeomFromText(segment1) as geom,
ST_Length(ST_GeomFromText(segment1)) as longueur
FROM public.p_obis t1
UNION SELECT 
id_origine,
id_sommet,
id_angle,
ST_GeomFromText(segment2) as geom,
ST_Length(ST_GeomFromText(segment2)) as longueur
FROM public.p_obis t2;


-- Défnition des SRID de l_segments

ALTER TABLE public.l_segments
ALTER COLUMN geom TYPE geometry(LineString, 3949)
USING ST_SetSRID(geom,3949);


-- Création de la couche des points d'intersection entre le buffer du point d'origine et la couche des segments

CREATE TABLE public.p_buf1 AS
WITH d AS 
    (SELECT id_origine,
            MIN(longueur) as rayon
     FROM public.l_segments
     GROUP BY id_origine)
SELECT
row_number() OVER () as id_buf1,
t1.id_origine,
d.rayon,
ST_INTERSECTION(ST_BOUNDARY(ST_BUFFER(t1.angle_geom, (0.95*d.rayon))),t2.geom) geom
FROM public.p_obis t1, public.l_segments t2, d
WHERE ST_INTERSECTS(ST_BUFFER(t1.angle_geom, (0.95*d.rayon)), t2.geom) IS TRUE AND t1.angle_geom && t2.geom    AND t1.id_origine = t2.id_origine AND d.id_origine = t1.id_origine;


-- Création des points d'intersection des deuxièmes buffers

CREATE TABLE public.p_buf2 AS SELECT
row_number() OVER () as id_buf2,
t1.id_origine,
ST_INTERSECTION(ST_BOUNDARY(ST_BUFFER(t1.geom,1.25*t1.rayon)),ST_BOUNDARY(ST_BUFFER(t2.geom,1.25*t1.rayon))) geom,
ST_GeometryType(ST_INTERSECTION(ST_BOUNDARY(ST_BUFFER(t1.geom,1.25*t1.rayon)),ST_BOUNDARY(ST_BUFFER(t2.geom,1.25*t1.rayon)))) as type_geom
FROM
public.p_buf1 as t1
INNER JOIN
public.p_buf1 as t2
ON (t1.id_origine = t2.id_origine)
WHERE t1.id_buf1 != t2.id_buf1;


-- Création de la ligne qui relie les 2 points composants les buf2

CREATE TABLE public.l_bissectrice AS SELECT
id_origine,
ST_MakeLine(t1.geom) geom FROM
p_buf2 t1
GROUP BY id_origine;

Au début, je pensais qu'il fallait conserver la même valeur de rayon du buffer entre la table buf1 et buf2 mais c'est ce qui me donnait des intersections de buf2 qui n'étaient pas forcément des points. Néanmoins avec le code que je viens d'indiquer j'avais un autre soucis, car le traçage de ce que je pense être une bissectrice (la ligne qui relie les couples de points créés par buf2) ne passe pas par mon point d'origine (EDIT : je n'arrive pas à uploader la capture d'écran alors qu'elle fait 15ko). A quoi pensez-vous que ce décalage puisse être dû ?

Sinon, toujours concernant votre requête. j'obtiens bien des linestrings mais elles sont composées des segments et d'une bissectrice qui n'est pas toujours orientée du côté à partir duquel j'avais calculé la valeur de mon angle (extérieur des parcelles).[

Dernière modification par white-shadow90 (Wed 25 October 2017 13:59)

white-shadow90 · Wed 25 October 2017 13:48

J'ai enlevé la partie

Code:

st_multi(
st_union(ARRAY[s1,s2,

Pour n'avoir plus que la bissectrice. Reste à trouver comment l'orienter
[EDIT] Et dans certains cas, j'ai des erreurs, le linestring qui reste correspond à un segment de l'angle et non à la bissectrice.

Dernière modification par white-shadow90 (Wed 25 October 2017 14:28)

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 Wed 18 October 2017 19:10

Couper un angle en 2 (équivalent bissectrice)

Code:

#2 Thu 19 October 2017 17:22

Re: Couper un angle en 2 (équivalent bissectrice)

Code:

#3 Mon 23 October 2017 10:24

Re: Couper un angle en 2 (équivalent bissectrice)

#4 Mon 23 October 2017 16:25

Re: Couper un angle en 2 (équivalent bissectrice)

Code:

#5 Tue 24 October 2017 13:31

Re: Couper un angle en 2 (équivalent bissectrice)

#6 Wed 25 October 2017 10:56

Re: Couper un angle en 2 (équivalent bissectrice)

Code:

#7 Wed 25 October 2017 12:00

Re: Couper un angle en 2 (équivalent bissectrice)

#8 Wed 25 October 2017 12:49

Re: Couper un angle en 2 (équivalent bissectrice)

Code:

#9 Wed 25 October 2017 13:48

Re: Couper un angle en 2 (équivalent bissectrice)

Code:

Pied de page des forums