Le portail francophone de la géomatique

with tmp as (
  -- on calcule les normales de chaque coté de chaque tronçon (distance de 300m prise au pif en fonction des données de test)
  SELECT
    c.dbid,
    st_startpoint(geom) AS geom,
    st_lineinterpolatepoint(geom, 0.01, 300) as pt_left,
    st_lineinterpolatepoint(geom, 0.01, -300) as pt_right
  FROM testcours c
  UNION ALL
  SELECT
    c.dbid,
    st_endpoint(geom) AS geom,
    st_lineinterpolatepoint(geom, 1, 300) as pt_left,
    st_lineinterpolatepoint(geom, 1, -300) as pt_right
  FROM testcours c
), tmp1 as (
    -- on regroupe les points identiques a partir desquels on va calculer la normale aux lignes en ce point
    -- pour certains points, on a deux lignes (lignes connectées), pour d'autres, une seule (extrémité de ligne)
    SELECT
      array_agg(t.dbid) AS dbids,
      array_agg(t.pt_left) AS lefts,
      array_agg(t.pt_right) AS rights,
      t.geom
    FROM tmp t
    GROUP BY t.geom
), tmp2 as (
  -- on duplique la ligne si une seule ligne en ce point, 
  --histoire de pouvoir tout le temps calculer une normale moyenne entre deux lignes (même identiques)
  select t.geom,
    case when array_length(t.dbids, 1) = 1 then dbids||dbids[1] else dbids end as dbids,
    case when array_length(t.lefts, 1) = 1 then lefts||lefts[1] else lefts end as lefts,
    case when array_length(t.rights, 1) = 1 then rights||rights[1] else rights end as rights
  from tmp1 t
), tmp3 as (
  -- on calcule les perpendiculaires moyennes en chaque point, qui donnent les lignes jaunes
  select t.dbids, st_makeLine(
        st_lineinterpolatepoint(st_makeLine(lefts[1], lefts[2]), 0.5),
        st_lineinterpolatepoint(st_makeLine(rights[1], rights[2]), 0.5)) as geom
  UNION ALL
  -- on ajoute l'exterieur du pg pour avoir un reseau de lignes (lignes jaunes + exterieur du PG) dont on veut reconstruire les surfaces
    select null::bigint[] as dbids, st_boundary(geom) as geom
    from testpg
  from tmp2 t
) , tmp4 as (
  -- on fabrique l'union des lignes et de l'exterieur du pg, avec st_union, histoire d'avoir un réseau de lignes connectées aux intersections,
  -- condition indispensable pour que st_polygonize marche
  select (st_dump(st_union(t.geom))).geom
  from tmp3 t
), tmp5 as (
  -- on construit tous les polygones issus des croisements des lignes
  -- précédentes avec st_polygonize, et on dump la collection obtenue en polygones individuels
  select (st_dump(st_polygonize(t.geom))).geom as geom
  from tmp4 t
), tmp6 as (
  -- on conserve les polygones contenus dans le fond:
  -- --pb de precision ici: st_contains ne renvoie pas les bons résultats: petit buffer négatif pour contourner le pb de précision
  select t.geom
  from tmp5 t, testpg p
  where st_contains(p.geom, st_buffer(t.geom, -0.1)) 
), tmp7 as (
  -- on associe les pg de fond issus du découpage aux tronçons, en prenant le milieu des troncons
  select c.dbid, t.geom
  from testcours c join tmp6 t
      on st_contains(t.geom, st_lineinterpolatepoint(c.geom, 0.5))
) select * from tmp7;

-- découpage du pg en deux berges:
drop table berge;
create table berge as (
    with tmp as (
    select st_union(geom) as geom
    from testcours
  ), tmp1 as (
    SELECT st_dump(st_difference(st_exteriorring(p.geom), t.geom)) as d
    FROM testpg p, tmp t
  ) select (t.d).path[1], (t.d).geom
  from tmp1 t
);

-- puis exploitation des berges et des troncons:

with tmp as (
  -- on calcule les points les plus proches des start et end points sur chaque berge (st_snapToGrid pour eviter des pb de précision)
  select c.dbid,
    st_startpoint(c.geom) as start,
    st_snapToGrid(st_lineinterpolatepoint(b1.geom, st_linelocatepoint(b1.geom, st_startpoint(c.geom))), 0.01) AS sb1,
    st_snapToGrid(st_lineinterpolatepoint(b1.geom, st_linelocatepoint(b1.geom, st_endpoint(c.geom))), 0.01) AS eb1,
    st_endpoint(c.geom) as end,
    st_snapToGrid(st_lineinterpolatepoint(b2.geom, st_linelocatepoint(b2.geom, st_endpoint(c.geom))), 0.01) AS eb2,
    st_snapToGrid(st_lineinterpolatepoint(b2.geom, st_linelocatepoint(b2.geom, st_startpoint(c.geom))), 0.01) AS sb2,
    -- on garde aussi la position du point d'intersection sur chaque berge, elle servira a st_lineSubString()
    st_linelocatepoint(b1.geom, st_startpoint(c.geom)) AS lsb1,
    st_linelocatepoint(b1.geom, st_endpoint(c.geom)) AS leb1,
    st_linelocatepoint(b2.geom, st_endpoint(c.geom)) AS leb2,
    st_linelocatepoint(b2.geom, st_startpoint(c.geom)) AS lsb2,
    b1.geom as b1geom,
    b2.geom as b2geom

  from testcours c, berge b1, berge b2
  WHERE b1.path = 1 AND b2.path = 2

), tmp1 as (
  -- on construit 4 arcs de ligne pour chaque troncon: ces arcs sont constitués
  -- des points calculés précédemment et des partie de berge interctées par les points sb1/eb1 et eb2/sb2
  select t.dbid,
    -- arc 1: le segment dans le fond de vallee, au niveau du startpoint
    st_snaptogrid(st_makeLine(array[t.sb2, t.start, t.sb1]), 0.01) as sub1,
    -- arc2: partie de berge 1
    st_snaptogrid(st_lineSubstring(b1geom, least(t.lsb1, t.leb1), greatest(t.lsb1, t.leb1)), 0.01) as sub2,
    -- arc 3: le segment dans le fond de vallee, au niveau du endpoint
    st_snaptogrid(st_makeLine(array[t.eb1, t.end, t.eb2]), 0.01) as sub3,
    -- arc4: partie de berge 2
    st_snaptogrid(st_lineSubstring(b2geom, least(t.lsb2, t.leb2), greatest(t.lsb2, t.leb2)), 0.01) as sub4
  from tmp t

), tmp2 as (
    -- on remet le tableau de géométries en ligne, car st_lineMerge prend une multilinestring à fusionner.
    -- pour chaque troncon, on obtient 4 lignes qu'on pourra GROUP BY dbid plus tard
    SELECT
      t.dbid,
      unnest(ARRAY [sub1, sub2, sub3, sub4]) AS lines
    FROM tmp1 t

), tmp3 as (
    -- on tente de reconstruire les polygones en mergeant les lignes par id de troncon
    SELECT
      t.dbid,
      st_lineMerge(st_union(t.lines)) as lines
    FROM tmp2 t
    GROUP BY t.dbid

)
-- on filtre les éventuels pg en erreur: on n'arrive pas a construire un pg a partir des 4 lignes
-- TODO: réordonner les points d'intersection sur les berges pour découper correctement 
select dbid, st_makePolygon(lines)
from tmp3
where geometryType(lines) = 'LINESTRING';