Построение диаграммы Вороного в ПостГИС

12

Я пытаюсь построить диаграмму Вороного из сетки точек, используя модифицированный код отсюда . Это запрос SQL после моих изменений:

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

Ниже - результат моего запроса. введите описание изображения здесь

Как видите, я получаю «почти» правильную диаграмму вороной, за исключением выделенных точек, которые не имеют уникального многоугольника. Ниже описано, что производит алгоритм QGIS и что я хотел бы получить из запроса. Любые предложения, где проблема с моим кодом?

введите описание изображения здесь

DamnBack
источник
Возможно, вы могли бы также сравнить результат функции SpatiaLite "VoronojDiagram" gaia-gis.it/gaia-sins/spatialite-sql-latest.html и взглянуть на исходный код в gaia-gis.it/fossil/libspatialite/ индекс .
user30184
Хороший вопрос, я смотрел на тот же вопрос, на который вы ссылаетесь, с целью его ускорения, но время уходит все меньше. Я не знал об этой проблеме с внешними точками.
Джон Пауэлл
5
Для чего стоит ST_Voronoi, появившийся в PostGIS 2.3, Дэн Бастон скоро проверит его код - trac.osgeo.org/postgis/ticket/2259 выглядит довольно много, просто нужно его потянуть. Было бы здорово иметь ребята тестируют
LR1234567
Можете ли вы опубликовать набор точек, которые вы используете? Я не возражал бы сделать небольшое тестирование на этом
MickyT
@MickyT Вот ссылка на мои данные. Data SRID
2180.

Ответы:

6

Хотя это устраняет непосредственную проблему с запросом данных, я не доволен им как решением для общего пользования, и я вернусь к этому и предыдущему ответу, когда смогу.

Проблема заключалась в том, что исходный запрос использовал выпуклую оболочку по краям вороного, чтобы определить внешнее ребро для многоугольника вороного. Это означало, что некоторые из ребер вороного не достигли наружу, когда они должны были быть. Я смотрел на использование функциональности вогнутого корпуса, но на самом деле это не сработало, как я хотел.
В качестве быстрого исправления я изменил выпуклую оболочку, которая будет построена вокруг замкнутого набора ребер вороного, плюс буфер вокруг исходных ребер. Края voronoi, которые не закрываются, вытянуты на большое расстояние, чтобы попытаться убедиться, что они пересекают внешнюю сторону и используются для построения многоугольников. Вы можете поиграть с параметрами буфера.

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
источник
Спасибо за объяснение и быстрое решение проблемы! Он работает с моими данными (немного медленнее из-за ST_Union(ST_Buffer(geom))), но я продолжу тестирование с другими наборами точек. Пока что я буду ждать, как вы сказали, более общего решения. :)
DamnBack
у вас есть изображение, которое вы можете опубликовать на вашем окончательном выводе?
Джерил Кук
10

Следуя предложению @ LR1234567, чтобы опробовать новую функциональность ST_Voronoi , разработанную @dbaston, оригинальный удивительный ответ @MickyT (как указано в вопросе OP) и использование исходных данных теперь можно упростить до:

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

что приводит к такому выводу, идентично вопросу ОП.

введите описание изображения здесь

Однако это страдает от той же проблемы, которая теперь исправлена ​​в ответе MickyT, о том, что точки на вогнутом корпусе не получают вмещающий многоугольник (как и следовало ожидать от алгоритма). Я исправил эту проблему с помощью запроса с помощью следующей серии шагов.

  1. Рассчитайте вогнутую оболочку из входных точек - точки на вогнутой оболочке - это те, которые имеют неограниченные многоугольники на выходной диаграмме Вороного.
  2. Найдите исходные точки на вогнутом корпусе (желтые точки на диаграмме 2 ниже).
  3. Буферизуйте вогнутую оболочку (буферное расстояние произвольно и, возможно, может быть найдено более оптимально из входных данных?).
  4. Найдите самые близкие точки в буфере вогнутой оболочки, ближайшей к точкам в шаге 2. Они показаны зеленым цветом на диаграмме ниже.
  5. Добавьте эти точки в исходный набор данных
  6. Рассчитайте диаграмму Вороного из этого комбинированного набора данных. Как видно на 3-й диаграмме, точки на корпусе теперь имеют вмещающие многоугольники.

Диаграмма 2, показывающая точки на вогнутом корпусе (желтый) и ближайшие точки для буферизации на корпусе (зеленый). Диаграмма 2.

Очевидно, что запрос можно было бы упростить / сжать, но я оставил именно эту форму в виде серии CTE, поскольку таким образом легче последовательно выполнять шаги. Этот запрос выполняется для исходного набора данных в миллисекундах (в среднем 11 мс на сервере разработки), тогда как ответ MickyT с использованием ST_Delauney выполняется за 4800 мс на том же сервере. DBaston утверждает, что еще один порядок увеличения скорости может быть получен при построении по сравнению с текущей магистралью GEOS, 3.6dev, благодаря усовершенствованиям в процедурах триангуляции.

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;

Диаграмма 3, показывающая все точки, теперь заключенные в многоугольник диаграмма 3

Примечание. В настоящее время ST_Voronoi включает сборку Postgis из исходного кода (версия 2.3 или транк) и связывание с GEOS 3.5 или выше.

Изменить: я только что посмотрел на Postgis 2.3, как он установлен на Amazon Web Services, и кажется, что имя функции теперь ST_VoronoiPolygons.

Без сомнения, этот запрос / алгоритм может быть улучшен. Предложения приветствуются.

Джон Пауэлл
источник
@dbaston. Хотите знать, есть ли у вас какие-либо комментарии по этому подходу?
Джон Пауэлл
1
хорошо, все точки же получить охватывающий полигон, это просто , что это непропорционально большое. Если и как они должны быть уменьшены, это довольно субъективно, и, не зная точно, что нужно для внешних полигонов, трудно понять, каков «лучший» способ. Ваш, кажется, хороший метод для меня. Я использовал менее сложный метод, который по духу похож на ваш, отбрасывая дополнительные точки вдоль буферизованной границы выпуклой оболочки с фиксированным интервалом, определяемым средней плотностью точек.
dbaston
@dbaston. Спасибо, просто убедившись, что я не пропустил что-то очевидное. Алгоритм сокращения внешних полигонов до чего-то большего в соответствии с размером внутренних (в моем случае областей почтовых индексов) - это то, о чем мне придется подумать еще немного.
Джон Пауэлл
@ Джон Барса Спасибо за еще одно отличное решение. Скорость расчетов более чем удовлетворительна при таком подходе. К сожалению, я хотел бы использовать этот алгоритм внутри моего плагина QGIS, и он должен работать с PostGIS 2.1+ из коробки. Но я обязательно буду использовать это решение после официального выпуска PostGIS 2.3. В любом случае спасибо за такие исчерпывающие ответы. :)
DamnBack
@DamnBack. Добро пожаловать. Мне это нужно было для работы, и ваш вопрос на самом деле мне очень помог, так как я понятия не имел о выходе ST_Voronoi, а старые решения намного медленнее (как вы заметили). Это было весело, понять это тоже :-)
Джон Пауэлл
3

Если у вас есть доступ к PostGIS 2.3, попробуйте новую функцию ST_Voronoi, недавно совершенную:

http://postgis.net/docs/manual-dev/ST_Voronoi.html

Для Windows есть скомпилированные сборки - http://postgis.net/windows_downloads/

LR1234567
источник
Спасибо за информацию, что есть готовая встроенная функция ST_Voronoi - я проверю это. К сожалению, мне нужно решение, которое работает на версиях PostGIS 2.1+, поэтому запрос @MickyT наиболее близок к моим потребностям на данный момент.
DamnBack
@ LR1234567. Требуется ли для этого какая-либо конкретная версия GEOS. Завтра у меня есть время для тестирования 2.3 и ST_Voronoi.
Джон Пауэлл
Требуется GEOS 3.5
LR1234567