Construire un diagramme de Voronoi dans PostGIS

12

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. entrez la description de l'image ici

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?

entrez la description de l'image ici

DamnBack
la source
Peut-être pourriez-vous aussi comparer le résultat de la fonction SpatiaLite "VoronojDiagram" gaia-gis.it/gaia-sins/spatialite-sql-latest.html et jeter un œil au code source dans gaia-gis.it/fossil/libspatialite/ index .
user30184
Belle question, j'ai regardé la même question que vous référencez, en vue d'accélérer, mais continuez à manquer de temps. Je n'étais pas au courant de ce problème avec les points extérieurs.
John Powell
5
Pour ce que ça vaut, nous avons ST_Voronoi à venir dans PostGIS 2.3, Dan Baston va bientôt vérifier le code pour cela - trac.osgeo.org/postgis/ticket/2259 semble quasiment terminé, il suffit d'être tiré. Ce serait génial d'avoir tests de personnes
LR1234567
Pouvez-vous publier l'ensemble des points que vous utilisez? Cela me dérangerait de faire un peu de test sur ce sujet moi
MickyT
@MickyT Voici un lien vers mes données. Le SRID de données est un 2180.
DamnBack

Réponses:

6

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.

WITH 
    -- Sample set of points to work with
    Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM 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, (SELECT ST_ExteriorRing(ST_ConvexHull(ST_Union(ST_Union(ST_Buffer(edge,20),ct)))) FROM Edges))))))).geom, 2180) geom
INTO 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)) * 200),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 200))
            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;
MickyT
la source
Merci pour l'explication et la solution rapide pour le problème! Cela fonctionne avec mes données (un peu plus lentement en raison de 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. :)
DamnBack
avez-vous une image que vous pouvez publier sur votre sortie finale?
Jeryl Cook
10

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:

WITH voronoi (vor) AS 
     (SELECT ST_Dump(ST_Voronoi(ST_Collect(geom))) FROM meshpoints)
SELECT (vor).path, (vor).geom FROM voronoi;

ce qui donne cette sortie, identique à la question du PO.

entrez la description de l'image ici

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.

  1. Calculez la coque concave des points d'entrée - les points sur la coque concave sont ceux qui ont des polygones non bornés dans le diagramme de Voronoi en sortie.
  2. Trouvez les points d'origine sur la coque concave (points jaunes sur le schéma 2 ci-dessous).
  3. Tamponner la coque concave (la distance du tampon est arbitraire et pourrait peut-être être trouvée de manière plus optimale à partir des données d'entrée?).
  4. Trouvez les points les plus proches sur le tampon de coque concave les plus proches des points de l'étape 2. Ils sont affichés en vert dans le diagramme ci-dessous.
  5. Ajouter ces points à l'ensemble de données d'origine
  6. Calculez le diagramme de Voronoi de cet ensemble de données combiné. Comme on peut le voir sur le 3ème diagramme, les points sur la coque ont maintenant des polygones englobants.

Diagramme 2 montrant les points sur la coque concave (jaune) et les points les plus proches du tampon sur la coque (vert). Diagramme 2.

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.

WITH 
  conv_hull(geom) AS 
        (SELECT ST_Concavehull(ST_Union(geom), 1) FROM meshpoints), 
  edge_points(points) AS 
        (SELECT mp.geom FROM meshpoints mp, conv_hull ch 
        WHERE ST_Touches(ch.geom, mp.geom)), 
  buffered_points(geom) AS
        (SELECT ST_Buffer(geom, 100) as geom FROM conv_hull),
  closest_points(points) AS
        (SELECT 
              ST_Closestpoint(
                   ST_Exteriorring(bp.geom), ep.points) as points,
             ep.points as epoints 
         FROM buffered_points bp, edge_points ep),
  combined_points(points) AS
        (SELECT points FROM closest_points 
        UNION SELECT geom FROM meshpoints),
  voronoi (vor) AS 
       (SELECT 
            ST_Dump(
                  ST_Voronoi(
                    ST_Collect(points))) as geom 
        FROM combined_points)
 SELECT 
     (vor).path[1] as id, 
     (vor).geom 
 FROM voronoi;

Diagramme 3 montrant tous les points maintenant inclus dans un polygone diagramme 3

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.

John Powell
la source
@dbaston. Vous vous demandez si vous avez des commentaires sur cette approche?
John Powell
1
bien, tous les points n'obtenir un polygone qui entoure, il est juste que c'est disproportionnée. Si et comment les réduire, c'est assez subjectif, et sans savoir exactement ce que l'on souhaite pour les polygones externes, il est difficile de savoir quelle est la "meilleure" manière. La vôtre me semble une belle méthode. J'ai utilisé une méthode moins sophistiquée qui ressemble à la vôtre, en déposant des points supplémentaires le long d'une limite de coque convexe tamponnée à un espacement fixe déterminé par la densité de points moyenne.
dbaston
@dbaston. Merci de m'assurer de ne rien manquer d'évident. Un algorithme pour réduire les polygones externes à quelque chose de plus en ligne avec la taille des polygones internes (dans mon cas, les zones de code postal) est quelque chose auquel je devrai réfléchir davantage.
John Powell
@John Barça Merci pour une autre excellente solution. La vitesse des calculs est plus que satisfaisante avec cette approche. Malheureusement, je voudrais utiliser cet algorithme dans mon plugin QGIS, et il doit fonctionner avec PostGIS 2.1+ prêt à l'emploi. Mais j'utiliserai sûrement cette solution après la sortie officielle de PostGIS 2.3. Quoi qu'il en soit, merci pour ces réponses complètes. :)
DamnBack
@DamnBack. Vous êtes les bienvenus. J'en avais besoin pour le travail et votre question m'a beaucoup aidé, car je n'avais aucune idée de la sortie de ST_Voronoi, et les anciennes solutions sont beaucoup plus lentes (comme vous l'avez remarqué). C'était amusant de le découvrir aussi :-)
John Powell
3

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/

LR1234567
la source
Merci pour l'information qu'il y a une fonction intégrée ST_Voronoi à venir - je vais la vérifier. Malheureusement, j'ai besoin d'une solution qui fonctionne sur les versions PostGIS 2.1+, donc la requête @MickyT est la plus proche de mes besoins pour le moment.
DamnBack
@ LR1234567. Cela nécessite-t-il une version particulière de GEOS? J'ai le temps demain de tester 2.3 et ST_Voronoi.
John Powell
Nécessite GEOS 3.5
LR1234567