How to split circles in 12 sections in PostGIS?
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)
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';
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: