How to split circles in 12 sections in PostGIS?

Question

I have a point layer in my PostgreSQL db and I would like to split a large number of circles into 12 sections each. Sample scenario is:

The circle corresponds to a 50 meter buffer of a point layer, and for each point, I need to split the buffer in 12 sections (pie-like wedges). The desired scenario is shown in following figure:

Can anyone suggest how to split circle into 12 sections?

Answer 1

Interesting question.

My suggestion is to use a simplified circle (only 12 segments) and then calculate a delaunay triangulation on that. Here's a working example:

CREATE TABLE twelves AS

WITH points AS ( 
    SELECT 1::int as id, ST_MakePoint(100,100) as geom
)
,dump AS (
    SELECT id, (ST_DumpPoints(ST_Buffer(geom,50,3))).geom as geom 
    FROM points --insert your point table here
    UNION ALL 
    SELECT id, geom FROM points --same here
)

SELECT id, (ST_Dump(ST_DelaunayTriangles(ST_Collect(geom),0, 0))).geom geom
FROM dump
GROUP BY id;

If you want to have a nicer (rounder) circle, you will have a lot more work to do but there may be interesting solutions for that as well.

[un-asked-for advice:] By the way, if this is for a visualisation, there might be more convenient client-sides way to do this.

EDIT:

Here's a more complex version that adds rounded triangles, at the expense of some speed.

WITH points AS ( 
    SELECT 1::int as id, ST_MakePoint(100,100) as geom
)
,circle AS ( --first make a nice circle with a lot of segments
    SELECT id, ST_Buffer(geom, 50,25) as geom FROM points
)
,dump AS ( --make the pie segments, but make them bigger than the nice circle
    SELECT id, (ST_DumpPoints(ST_Buffer(geom,60,3))).geom as geom 
    FROM points --insert your point table here
    UNION ALL 
    SELECT id, geom FROM points --same here
)
,triangles AS (
    SELECT id, (ST_Dump(ST_DelaunayTriangles(ST_Collect(geom),0, 0))).geom geom
    FROM dump
    GROUP BY id
)
--now get the intersection between the nice circle and the segments
SELECT a.id, ST_Intersection(a.geom, b.geom) geom
FROM circle as a, triangles as b
WHERE a.id = b.id;

EDIT2:

Here's an even more complex example that gives the order of the triangles with a number

WITH points AS ( 
    SELECT 1::int as id, ST_MakePoint(100,100) as geom
)
,circle AS ( --first make a nice circle with a lot of segments
    SELECT id, ST_Buffer(geom, 50,25) as geom FROM points
)
,segments AS ( --create segments from a smaller circle so we can find out later wich triangle belongs to which segment
     SELECT id, (pt).path path, ST_MakeLine(lag((pt).geom, 1, NULL) OVER (PARTITION BY id ORDER BY id, (pt).path), (pt).geom) AS geom
      FROM (SELECT id, ST_DumpPoints(ST_Buffer(geom,40,3)) AS pt FROM points) as dumps
)
,dump AS ( --make the pie segments, but make them bigger than the nice circle
    SELECT id, (ST_DumpPoints(ST_Buffer(geom,60,3))).geom as geom
    FROM points --insert your point table here
    UNION ALL 
    SELECT id, geom FROM points --same here
)
,triangles AS ( --triangles will have a random order
    SELECT id, (ST_Dump(ST_DelaunayTriangles(ST_Collect(geom),0, 0))).geom geom
    FROM dump
    GROUP BY id
)
--now get the intersection between the nice circle and the segments, and add the ordernr of the triangle based on the segments we got earlier on
SELECT a.id, c.path[2]-1 path, ST_Intersection(a.geom, b.geom) geom
FROM circle as a, triangles as b
LEFT JOIN segments c ON ST_Intersects(b.geom,ST_Centroid(c.geom))
WHERE a.id = b.id
ORDER BY a.id, path;

(Thanks to Paul Ramsey for the excellent example of getting segments from a linestring: http://blog.cleverelephant.ca/2015/02/breaking-linestring-into-segments.html)

Answer 2

At the moment I don't work at the project, which includes the function I've posted as a comment.

But I post here the untested function, which could do what you need.

Same as in my question:

CREATE TYPE quadrant AS (id smallint,geom geometry(polygon,31468))

Change the SRID to your project SRID.

Usage:

SELECT (quadrant(20,0.0,90)).*

The first parameter stands for the point_id, the second one is the start angle and the third one is the step in how many wedges your circle is cutted (90°=4, 30°=12).

Now the untested function:

CREATE OR REPLACE FUNCTION quadrant(id integer,start double precision, stop integer) RETURNS SETOF quadrant AS $$
WITH centroid AS
(SELECT
    ST_Buffer(geom, 50) AS buffer,
    geom AS vertex,
    point_id
FROM your_point_layer
    WHERE point_id=$1
),
newline AS
(SELECT 
ST_SetSRID(ST_Translate(
                ST_Rotate(
                    ST_MakeLine(
                          ST_MakePoint(0.0,2000.0), --check this with your buffer distance (50m buffer vs. 2000m span (60m could be enough))
                          ST_MakePoint(0.0,0.0)),
                radians(($2+s.a)*-1)),
          ST_X(vertex), ST_Y(vertex)),
ST_SRID(vertex)) AS geom
FROM centroid, generate_series(0,$3,$3) AS s(a)
),
span AS
(SELECT
    centroid.point_id,
    ST_LineMerge(ST_Union(newline.geom)) AS geom
FROM newline, centroid
    GROUP BY point_id),
multiobject AS
(SELECT
    span.point_id,
    ST_Split(centroid.buffer,span.geom) AS geom,
    generate_series(1,100) AS n --check this regarding how many wedges you want to have
FROM span, centroid
    WHERE centroid.point_id=$1),
objects AS
(SELECT
    n,
    ST_GeometryN(multiobject.geom,n) AS geom
FROM multiobject
    WHERE n <= ST_NumGeometries(multiobject.geom))

SELECT
    point_id AS id,
    objects.geom
FROM objects, multiobject
    WHERE multiobject.n <= ST_NumGeometries(multiobject.geom)
$$ LANGUAGE 'sql';

Answer 3

I produced these 12 sections split circles by using PyQGIS; where they could be introduced into PostGIS database by using psycopg2 python module at the same script (not included here). PyQGIS script is:

import psycopg2
import numpy as np

bufferLength = 600
polygonSides = 12

registry = QgsMapLayerRegistry.instance()

layer = registry.mapLayersByName('point')

feat_points = [ feat for feat in layer[0].getFeatures() ] 
points = [ feat.geometry().asPoint() for feat in layer[0].getFeatures() ]

epsg = layer[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'buffer',
                           'memory')

prov = mem_layer.dataProvider()

for i, point in enumerate(points):
    outFeat = QgsFeature()

    outFeat.setGeometry(QgsGeometry.fromPolygon([[ QgsPoint(point[0] + np.sin(angle)*bufferLength, point[1] + np.cos(angle)*bufferLength)
                        for angle in np.linspace(0, 2*np.pi, polygonSides, endpoint = False) ]]))

    outFeat.setAttributes([i])
    prov.addFeatures([outFeat])

feats_mem = [ feat for feat in mem_layer.getFeatures() ]

mem_layer2 = QgsVectorLayer(uri,
                           'sections',
                           'memory')

prov = mem_layer2.dataProvider()

k = 0

for feat in feats_mem:
    geom = feat.geometry().asPolygon()[0]
    n = len(geom)
    new_pol = []
    for i in range(n-1):
        new_pol.append([[ points[k], geom[i], geom[i+1]]])

    feat = QgsFeature()

    buffer_geom = feat_points[k].geometry().buffer(500, -1)

    for i, element in enumerate(new_pol):
        feat = QgsFeature()
        geom = QgsGeometry.fromPolygon(element)
        new_geom = geom.intersection(buffer_geom)
        feat.setGeometry(new_geom)
        feat.setAttributes([i])
        prov.addFeatures([feat])

    QgsMapLayerRegistry.instance().addMapLayer(mem_layer2)

    k += 1

After running the script at the Python Console of QGIS, with a 3 featured point layer, I got:

Answer 4

Okay, I know this is ancient, but I somehow got back here today and had a couple of minutes to follow my brainstorming (see my comment above);

This query yields 4k circles in around 1.3 secs on my machine, with segment count starting at 12 o'clock. minor backdraw is the amount of vertices for the circle (same as most answers, though: between 4 to 6 per segment. its based on the buffer outline and I don´t know how the number of vertices is decided upon internally):

SELECT ROW_NUMBER() OVER() AS id,
       c.circle_id,
       c.seg_id,
       ST_MakePolygon(
         St_AddPoint(
           St_Addpoint(c.geom, c.center),
           ST_StartPoint(c.geom)
         )
       ) AS geom
FROM (
  SELECT b.circle_id,
         CASE
           WHEN n + 4 > 12
             THEN (n + 4) - 12
             ELSE n + 4
         END AS seg_id,
         b.center,
         ST_LineSubstring(
           b.geom,
           ((ST_Length(b.geom) / 12) * n) / ST_Length(b.geom),
           ((ST_Length(b.geom) / 12) * (n + 1))/ ST_Length(b.geom)
         ) AS geom
  FROM (    
    SELECT ROW_NUMBER() OVER() AS circle_id,
           a.geom AS center,
           ST_ExteriorRing(
             ST_Buffer(a.geom, 50)
           ) AS geom
    FROM <POINTS> AS a
  ) AS b
  CROSS JOIN
    generate_series(0, 11) AS n
  ORDER BY circle_id, seg_id
) AS c

Yay, fun question!!! Also, my milkshake is definetly better than yours.