11 Apr
2020
11 Apr
'20
11:41 a.m.
Message posté par : Nicolas Ribot
----------------------------------------
Ca peut facilement se mettre dans une fonction:
-----------------
Code :
create or replace function st_dividePolygon(pg geometry, num_parts int default 10,
num_points int default 4000)
returns geometry[]
language sql
as $$
with pg_pts as (
SELECT (ST_Dump(ST_GeneratePoints(pg, num_points))).geom AS geom
), pg_pts_clustered AS (
SELECT geom, ST_ClusterKMeans(geom, 10) over () AS cluster
FROM pg_pts
), pg_centers AS (
SELECT cluster, ST_Centroid(ST_collect(geom)) AS geom
FROM pg_pts_clustered
GROUP BY cluster
), pg_voronoi AS (
SELECT (ST_Dump(ST_VoronoiPolygons(ST_collect(geom)))).geom AS geom
FROM pg_centers
) SELECT array_agg(ST_Intersection(pg, b.geom))
FROM pg_voronoi b;
$$;
-----------------
Pour calculer le nombre de morceaux, st_area(geom)/max_area
La méthode ne fabrique pas des morceaux qui ont tous la meme surface, mais ont des
surfaces approchantes (ca peut qd meme varier du simple au double sur des pg a la forme
complexe, quand meme ;) )
Il existe des méthodes exactes (cf lien dans le blog de Paul Ramsey).
En image: les dept de france découpés pour avoir, pour chacun, des morceaux de surface ~
55 km2:
(seuls les dept qui seraient découpés sont calculés)
-----------------
Code :
create table deptcut_ea as
with tmp as (
select gid, round(st_area(geom) / 550000000) as num_parts, geom
from dept2154
) select gid, num_parts, unnest(st_dividePolygon(geom)) as geom
from tmp
where num_parts > 1;
-----------------
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=331999#p331999
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