J'essaie de construire un diagramme voronoi à partir d'une grille de points en utilisant du code modifié d' ici . Il s'agit d'une requête SQL après mes modifications:
DROP TABLE IF EXISTS example.voronoi;
WITH
-- Sample set of points to work with
Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM example."MeshPoints2d"),
-- Build edges and circumscribe points to generate a centroid
Edges AS (
SELECT id,
UNNEST(ARRAY['e1','e2','e3']) EdgeName,
UNNEST(ARRAY[
ST_MakeLine(p1,p2) ,
ST_MakeLine(p2,p3) ,
ST_MakeLine(p3,p1)]) Edge,
ST_Centroid(ST_ConvexHull(ST_Union(-- Done this way due to issues I had with LineToCurve
ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p1,p2),ST_MakeLine(p2,p3)))),'LINE','CIRCULAR'),15),
ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p2,p3),ST_MakeLine(p3,p1)))),'LINE','CIRCULAR'),15)
))) ct
FROM (
-- Decompose to points
SELECT id,
ST_PointN(g,1) p1,
ST_PointN(g,2) p2,
ST_PointN(g,3) p3
FROM (
SELECT (gd).Path id, ST_ExteriorRing((gd).geom) g -- ID andmake triangle a linestring
FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a -- Get Delaunay Triangles
)b
) c
)
SELECT ST_SetSRID((ST_Dump(ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, ST_ExteriorRing(ST_ConvexHull(v)))))))).geom, 2180)
INTO example.voronoi
FROM (
SELECT -- Create voronoi edges and reduce to a multilinestring
ST_LineMerge(ST_Union(ST_MakeLine(
x.ct,
CASE
WHEN y.id IS NULL THEN
CASE WHEN ST_Within(
x.ct,
(SELECT ST_ConvexHull(geom) FROM sample)) THEN -- Don't draw lines back towards the original set
-- Project line out twice the distance from convex hull
ST_MakePoint(ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * 2),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 2))
END
ELSE
y.ct
END
))) v
FROM Edges x
LEFT OUTER JOIN -- Self Join based on edges
Edges y ON x.id <> y.id AND ST_Equals(x.edge,y.edge)
) z
Ci-dessous - résultat de ma requête.
Comme vous pouvez le voir, j'obtiens un diagramme de voronoi «presque» correct, sauf les points en surbrillance qui n'ont pas de polygone unique. Voici ce que l'algorithme QGIS produit et ce que j'aimerais obtenir de la requête. Des suggestions où est le problème avec mon code?
postgis
sql
voronoi-thiessen
DamnBack
la source
la source
Réponses:
Bien que cela corrige le problème immédiat de la requête pour les données en question, je ne suis pas satisfait de cela comme solution pour une utilisation générale et je reviendrai sur cela et la réponse précédente lorsque je le pourrai.
Le problème était que la requête d'origine utilisait une coque convexe sur les bords de voronoi pour déterminer le bord extérieur du polygone de voronoi. Cela signifiait que certains des bords des voronoi n'atteignaient pas l'extérieur alors qu'ils auraient dû l'être. J'ai envisagé d'utiliser la fonctionnalité de coque concave, mais cela n'a pas vraiment fonctionné comme je le souhaitais.
Comme solution rapide, j'ai changé la coque convexe pour être construite autour de l'ensemble fermé des bords de voronoi plus un tampon autour des bords d'origine. Les bords de voronoi qui ne se ferment pas sont étendus sur une grande distance pour essayer de s'assurer qu'ils croisent l'extérieur et sont utilisés pour construire des polygones. Vous voudrez peut-être jouer avec les paramètres du tampon.
la source
ST_Union(ST_Buffer(geom))), mais je continuerai de tester avec d'autres ensembles de points. En attendant, j'attends, comme vous l'avez dit, une solution plus générale. :)Suite à une suggestion de @ LR1234567 d'essayer la nouvelle fonctionnalité ST_Voronoi qui a été développée par @dbaston, la réponse étonnante originale de @MickyT (comme indiqué dans la question OP) et l'utilisation des données originales peut maintenant être simplifiée pour:
ce qui donne cette sortie, identique à la question du PO.
Cependant, cela souffre du même problème, maintenant corrigé dans la réponse de MickyT, que les points sur la coque concave n'obtiennent pas de polygone englobant (comme on pourrait s'y attendre de l'algorithme). J'ai résolu ce problème avec une requête avec la série d'étapes suivante.
Diagramme 2 montrant les points sur la coque concave (jaune) et les points les plus proches du tampon sur la coque (vert).
La requête pourrait évidemment être simplifiée / compressée, mais j'ai laissé cette forme comme une série de CTE, car il est plus facile de suivre les étapes séquentiellement de cette façon. Cette requête s'exécute sur l'ensemble de données d'origine en millisecondes (11 ms en moyenne sur un serveur de développement) tandis que la réponse de MickyT utilisant ST_Delauney s'exécute en 4800 ms sur le même serveur. DBaston prétend qu'un autre ordre d'amplitude de vitesse peut être obtenu en construisant contre le tronc actuel de GEOS, 3.6dev, en raison des améliorations des routines de triangulation.
Diagramme 3 montrant tous les points maintenant inclus dans un polygone
Remarque: Actuellement, ST_Voronoi implique la construction de Postgis à partir de la source (version 2.3 ou tronc) et la liaison avec GEOS 3.5 ou supérieur.
Edit: je viens de regarder Postgis 2.3 tel qu'il est installé sur Amazon Web Services, et il semble que le nom de la fonction soit maintenant ST_VoronoiPolygons.
Nul doute que cette requête / cet algorithme pourrait être amélioré. Suggestions bienvenues.
la source
Si vous avez accès à PostGIS 2.3, essayez la nouvelle fonction ST_Voronoi, récemment validée:
http://postgis.net/docs/manual-dev/ST_Voronoi.html
Il existe des versions précompilées pour Windows - http://postgis.net/windows_downloads/
la source