18 Feb
2021
18 Feb
'21
12:09 p.m.
Message posté par : Nicolas Ribot
----------------------------------------
Alors voila une méthode qui marche mieux:
Le principe est de travailler avec les arcs des deux couches (arc au sens topologique:
linestrings connectées entre elles aux extrémités) et de reconstruire les polygones à
partir de ces arcs:
• on veut tous les arcs des provinces (car contours de référence)
• et les arcs des communes internes aux provinces (ne suivant pas le contour des
provinces)
Pour obtenir ces arcs, st_linemerge, st_union et st_boundary sont utilisées (les
boundaries d'un pg sont des linestrings):
-----------------
Code :
-- prepa tables ajout clé primaire (id plus facile
alter table com add column id serial primary key ;
alter table province add column id serial primary key ;
select * from spatial_ref_sys where srid = 25830;
-- ajout clés primaires
alter table dts_commune add column id serial primary key ;
alter table dts_province add column id serial primary key ;
-- on tente de travailler avec les segments (arcs) connectés des commune:
-- union et dump.
select st_numgeometries(geom) from dts_province;
select st_numgeometries(geom) from dts_commune;
-- ok pg simple, on utilise st_geometryN(1) pour extraire le pg.
drop table dts_com_arcs;
create table dts_com_arcs as
select (st_dump(st_lineMerge(st_union(st_boundary(st_geometryN(geom, 1)))))).geom
from dts_commune;
alter table dts_com_arcs add column id serial primary key ;
-- et des provinces
drop table dts_prov_arcs;
create table dts_prov_arcs as
select (st_dump(st_lineMerge(st_union(st_boundary(st_geometryN(geom, 1)))))).geom
from dts_province;
alter table dts_prov_arcs add column id serial primary key ;
-----------------
Ensuite on trouve une méthode pour garder les arcs des communes internes aux provinces:
par ex avec une distance moyenne des points de l'arc aux contours de la province: les
arcs internes auront une grande distance moyenne par rapport aux arcs suivant le contour
des provinces:
-----------------
Code :
-- on veut faire la reconstruction des arcs des provinces avec les arcs internes des
communes
-- il faut identifier les arcs internes des communes:
-- plusieurs méthodes, par ex distance moyenne des points de l'arc au arcs des
communes
create table dts_com_arcs_pt as
select id as idarc, (st_dumppoints(geom)).geom
from dts_com_arcs;
-- pk
alter table dts_com_arcs_pt add column id serial primary key ;
-- spatial index:
create index on dts_com_arcs_pt using gist(geom);
-- et sur les arcs province
create index on dts_prov_arcs using gist(geom);
vacuum analyse dts_prov_arcs;
vacuum analyse dts_com_arcs_pt;
-- methode empirique ici: on determine quelle dist moyenne permet de choisir les arcs
internes:
-- disons 62m :)
drop table dts_com_internal_arc;
create table dts_com_internal_arc as
with tmp as (
select p.idarc, avg(t.dist) as avg_dist
from dts_com_arcs_pt p
cross join lateral (
select a.geom <-> p.geom as dist
from dts_prov_arcs a
order by a.geom <-> p.geom
limit 1
) as t
group by p.idarc
) select tmp.*, a.geom
from tmp join dts_com_arcs a on tmp.idarc = a.id
where avg_dist > 62;
-----------------
On corrige un peu la topologie en forcant la connexion des arcs internes des communes avec
les arcs des provinces, pour etre sur que notre réseau soit fermé: les fn de réf linéaires
marchent bien pour ca: on ajoute alors des points aux arcs internes, ce qui les connecte
aux arcs des provinces.
Il faut penser à cette étape à garder tous les arcs internes des communes: ceux prolongés
et ceux qui ne le sont pas
Une fois qu'on a les arcs internes prolongés et les arcs des provinces, on passe par
st_polygonize pour reconstruire les polygones.
A cette étape, les erreurs topo des communes peuvent engendrer des petits pg
supplémentaires:
pas grave, on fait une jointure spatiale entre les communes initiales et les centroid des
communes corrigées, on retrouve alors nos communes: un st_union finale group by commune
permet de nettoyer les petits polygones en erreurs:
-----------------
Code :
-- on corrige les erreurs topo éventuelles des arcs en prolongeant les arcs internes sur
les contours des provinces,
-- sauf si l'arc interne est connecté à un autre arc interne:
-- on cherche l'arc province le plus proche de chaque start/end de nos arcs internes
-- on reprend notre distance de 100m: si dist <, on snape ce point sur l'arc:
-- ce snap est fait avec les fn de référencement linéaire, qui sont rapides et précises
-- idx code le startpoint/endpoint: 0 pour start, -1 pour end: ca permettra de faire
-- un addpoint avec cet idx.
-- une fois les arcs internes prolongés, on peut construire le réseau de lignes
-- arcs internes + arcs province et reconstruire les pg communes avec les arcs province
-- on peut faire en 1 seule requete, ou détailler par étape, suivant le dataset
-- attention a l'etape tmp1: on filtre certains arcs => il faut etre sur qu'on
travaille avec tous les arcs
-- internes avant la reconstruction des pg:
drop table dts_newcom;
create table dts_newcom as
with tmp as (
select idarc, 0 as idx, st_startpoint(geom) as geom
from dts_com_internal_arc
UNION ALL
select idarc, -1, st_endpoint(geom) as geom
from dts_com_internal_arc
), tmp1 as (
select t.idarc, t.idx, t.geom,
st_lineinterpolatepoint(b.geom, st_linelocatepoint(b.geom, t.geom)) as newpt
from tmp t
cross join lateral (
select a.id as idarcprov, t.geom <-> a.geom as dist, geom
from dts_prov_arcs a
order by t.geom <-> a.geom
limit 1
) as b
where dist < 100
), tmp2 as (
select st_addPoint(a.geom, t.newpt, t.idx) as geom
from tmp1 t
join dts_com_arcs a on t.idarc = a.id
UNION ALL
select geom
from dts_prov_arcs
-- on ajoute les éventuels arcs internes filtrés à l'étape précédente:
UNION ALL
select geom
from dts_com_internal_arc a
where not exists (
select null
from tmp1
where a.idarc = tmp1.idarc
)
), tmp3 as (
select st_union(geom) as geom
from tmp2
) select (st_dump(st_polygonize(geom))).geom
from tmp3;
-- plus de pg que de communes:
-- les erreurs topo entrainent la création de petits pg: facile a merger avec leur voisin
-- pour les faire disparaitre: par ex en refaisant le lien spatial avec les communes et en
faisant l'union des pg par commune:
-- on peut faire un lien spatial (st_pointOnSurface pour retrouver les communes
initiales:
-- table finale des communes:
select cod_mun, mun_cod_prov, nombre_com, c.id, st_union(n.geom) as geom
from dts_commune c join dts_newcom n on st_contains(c.geom, st_pointOnSurface(n.geom))
group by 1, 2, 3, 4;
-----------------
Et en image avec la piece jointe.
Nico
----------------------------------------
Ce message est accompagné de fichiers, pour les télécharger, suivre le lien ci-dessous.
----------------------------------------
Le message est situé https://georezo.net/forum/viewtopic.php?pid=340426#p340426
Pour y répondre : geobd(a)ml.georezo.net ou reply de votre messagerie
Pour vous désabonner connectez-vous sur le forum puis Profil / Abonnement
--
Association GeoRezo - le portail géomatique
https://georezo.net