Construindo o diagrama de Voronoi no PostGIS

12

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. insira a descrição da imagem aqui

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?

insira a descrição da imagem aqui

DamnBack
fonte
Talvez você também possa comparar o resultado da função SpatiaLite " VoronojDiagram " gaia-gis.it/gaia-sins/spatialite-sql-latest.html e dar uma olhada no código fonte em gaia-gis.it/fossil/libspatialite/ índice .
user30184
Boa pergunta, estive analisando a mesma pergunta que você faz referência, com o objetivo de acelerá-la, mas continua ficando sem tempo. Eu não estava ciente desse problema com os pontos externos.
John John Powell
5
Pelo que vale a pena ter ST_Voronoi no PostGIS 2.3, Dan Baston logo verificará o código - trac.osgeo.org/postgis/ticket/2259 parece muito bem feito, só precisa ser puxado. Seria ótimo ter teste de pessoal
LR1234567
Você pode postar o conjunto de pontos que está usando? Eu me importaria de fazer um pouco de testes sobre isso mesmo
MickyT
@MickyT Aqui está um link para meus dados. Dados SRID é um 2180.
DamnBack

Respostas:

6

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.

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
fonte
Obrigado pela explicação e solução rápida para o problema! Funciona com meus dados (um pouco mais devidos ST_Union(ST_Buffer(geom))), mas continuarei testando com outros conjuntos de pontos. Enquanto isso, vou esperar como você disse - solução mais geral. :)
DamnBack 2/15/15
você tem uma imagem que pode postar em sua saída final?
precisa saber é o seguinte
10

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:

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

o que resulta nessa saída, idêntica à pergunta do OP.

insira a descrição da imagem aqui

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.

  1. Calcular o casco côncavo dos pontos de entrada - os pontos no casco côncavo são aqueles que possuem polígonos sem limites no diagrama de Voronoi de saída.
  2. Encontre os pontos originais no casco côncavo (pontos amarelos no diagrama 2 abaixo).
  3. Faça o buffer do casco côncavo (a distância do buffer é arbitrária e talvez possa ser encontrada de maneira mais otimizada a partir dos dados de entrada?).
  4. Encontre os pontos mais próximos no buffer do casco côncavo mais próximo dos pontos na etapa 2. Eles são mostrados em verde no diagrama abaixo.
  5. Adicione esses pontos ao conjunto de dados original
  6. Calcule o diagrama Voronoi deste conjunto de dados combinado. Como pode ser visto no terceiro diagrama, os pontos no casco agora têm polígonos anexos.

Diagrama 2 mostrando os pontos no casco côncavo (amarelo) e os pontos mais próximos para amortecer no casco (verde). Diagrama 2.

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.

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;

Diagrama 3 mostrando todos os pontos agora colocados em um polígono diagrama 3

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.

John Powell
fonte
@dbaston. Quer saber se você tem algum comentário sobre essa abordagem?
John Powell
1
bem, todos os pontos que conseguir um polígono colocando, é apenas que é desproporcionalmente grande. Se e como elas devem ser reduzidas é bastante subjetivo e, sem saber exatamente o que é desejado para os polígonos externos, é difícil saber qual é o "melhor" caminho. O seu parece um bom método para mim. Eu usei um método menos sofisticado, com espírito semelhante ao seu, eliminando pontos extras ao longo de um limite de casco convexo em um espaçamento fixo determinado pela densidade média de pontos.
dbaston
@dbaston. Obrigado, apenas certificando-me de que não estava perdendo algo óbvio. Um algoritmo para reduzir os polígonos externos para algo mais alinhado com o tamanho dos internos (no meu caso, áreas de código postal) é algo sobre o qual terei que pensar um pouco mais.
John Powell
@ John Barça Obrigado por outra ótima solução. A velocidade dos cálculos é mais do que satisfatória com essa abordagem. Infelizmente, eu gostaria de usar esse algoritmo dentro do meu plugin QGIS, e ele tem que funcionar com o PostGIS 2.1+ imediatamente. Mas certamente vou usar esta solução após o lançamento oficial do PostGIS 2.3. De qualquer forma, obrigado por respostas tão abrangentes. :)
DamnBack
@DamnBack. Você é muito bem-vindo. Eu precisava disso para trabalhar e sua pergunta realmente me ajudou muito, pois eu não tinha idéia do lançamento do ST_Voronoi, e as soluções mais antigas são muito mais lentas (como você notou). Foi divertido descobrir isso também :-)
John Powell
3

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/

LR1234567
fonte
Obrigado pela informação de que existe a função ST_Voronoi incorporada - eu vou dar uma olhada. Infelizmente, preciso de uma solução que funcione nas versões PostGIS 2.1+, portanto a consulta @MickyT é a mais próxima das minhas necessidades no momento.
precisa saber é o seguinte
@ LR1234567. Isso requer alguma versão específica do GEOS. Amanhã tenho tempo para testar o 2.3 e o ST_Voronoi.
John Powell
Requer GEOS 3.5
LR1234567