Estou tentando construir o diagrama voronoi da grade de pontos usando o código modificado daqui . Esta é a consulta SQL após minhas modificações:
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
Abaixo - resultado da minha consulta.
Como você pode ver, eu recebo o diagrama voronoi "quase" correto, exceto os pontos destacados que não têm polígono exclusivo. Abaixo está o que o algoritmo QGIS produz e o que eu gostaria de obter da consulta. Alguma sugestão onde está o problema com o meu código?
postgis
sql
voronoi-thiessen
DamnBack
fonte
fonte
Respostas:
Embora isso resolva o problema imediato com a consulta dos dados em questão, não estou satisfeito com isso como uma solução para uso geral e revisitarei esta e a resposta anterior quando puder.
O problema era que a consulta original estava usando um casco convexo nas bordas de voronoi para determinar a borda externa do polígono de voronoi. Isso significava que algumas das arestas voronoi não estavam atingindo o exterior quando deveriam estar. Eu olhei para usar a funcionalidade do casco côncavo, mas realmente não funcionou como eu queria.
Como solução rápida, mudei o casco convexo para ser construído em torno do conjunto fechado de arestas voronoi mais um amortecedor em torno das arestas originais. As arestas voronoi que não fecham são estendidas por uma grande distância para tentar garantir que elas cruzem o exterior e sejam usadas para construir polígonos. Você pode querer brincar com os parâmetros do buffer.
fonte
ST_Union(ST_Buffer(geom))
), mas continuarei testando com outros conjuntos de pontos. Enquanto isso, vou esperar como você disse - solução mais geral. :)Após uma sugestão de @ LR1234567 para experimentar a nova funcionalidade ST_Voronoi que foi desenvolvida por @dbaston, a incrível resposta original do @MickyT (conforme declarado na pergunta do OP) e o uso dos dados originais agora podem ser simplificados para:
o que resulta nessa saída, idêntica à pergunta do OP.
No entanto, isso sofre com o mesmo problema, agora corrigido na resposta de MickyT, de que os pontos no casco côncavo não recebem um polígono anexo (como seria de esperar do algoritmo). Corrigi esse problema com uma consulta com as seguintes séries de etapas.
Diagrama 2 mostrando os pontos no casco côncavo (amarelo) e os pontos mais próximos para amortecer no casco (verde).
Obviamente, a consulta poderia ser simplificada / compactada, mas deixei esse formato como uma série de CTEs, pois é mais fácil seguir as etapas sequencialmente dessa maneira. Essa consulta é executada no conjunto de dados original em milissegundos (média de 11ms em um servidor dev), enquanto a resposta de MickyT usando ST_Delauney é executada em 4800ms no mesmo servidor. O DBaston afirma que outra melhoria de velocidade de ordem de magnitude pode ser obtida com a construção contra o tronco atual do GEOS, 3.6dev, devido a melhorias nas rotinas de triangulação.
Diagrama 3 mostrando todos os pontos agora colocados em um polígono
Nota: Atualmente, o ST_Voronoi envolve a construção do Postgis a partir da fonte (versão 2.3 ou tronco) e a vinculação ao GEOS 3.5 ou superior.
Edit: Acabei de examinar o Postgis 2.3, pois ele está instalado no Amazon Web Services, e parece que o nome da função agora é ST_VoronoiPolygons.
Sem dúvida, essa consulta / algoritmo poderia ser melhorada. Sugestões são bem-vindas.
fonte
Se você tiver acesso ao PostGIS 2.3, tente a nova função ST_Voronoi, recentemente confirmada:
http://postgis.net/docs/manual-dev/ST_Voronoi.html
Existem compilações pré-compiladas para Windows - http://postgis.net/windows_downloads/
fonte