6

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:

sample scenario

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:

desired scenario

Can anyone suggest how to split circle into 12 sections?

  • 1
    A few years ago I've posted a question: gis.stackexchange.com/questions/81560/…. This could be interesting for you. You have to edit the function to get what you want. – Stefan May 23 '17 at 14:07
  • Thanks @Stefan. I just had a look at your code. I will try to modify it and see if it works! – khajlk May 23 '17 at 14:17
  • 2
    Pedantic note - those sectors aren't quadrants, because they do not divide the circle into four. – Toby Speight May 23 '17 at 17:14
9

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)

  • 1
    just brainstorming for the rounder circle: maybe use st_boundary to get the closure as linestring in conjunction with st_segmentize by the appropriate length (linelength divided by 12) to get 12 equal linestrings. then create a polygon with the lines vertices and the center point? – ThingumaBob May 23 '17 at 13:36
  • @tilt: Thanks a lot for your feedback. Your code worked but as you said I am not getting a rounded circle. I would still prefer to have a rounded circle with 12 quadrants but in any case your code is much helpful. Thanks a lot. Any suggestions for rounded circle? I want to give it a try before I mark your answer as a solution. – khajlk May 23 '17 at 14:15
  • 1
    I've added a version that does what you want. It's worth comparing to @Stefans version because his might be faster and more versatile. – tilt May 23 '17 at 14:56
  • Thanks @tilt. It worked and yielded what was required. Nevertheless, I would like to test and try Stefan's solution as well. – khajlk May 23 '17 at 15:09
  • @tilt: I just noticed that problem with your code is that it doesn't generate sections in increasing order from true north. For example, id (or gid) of 1st segment (angle 30 degree) should 1, 2nd segment's gid as 2 and so on.. in clockwise direction. Why is it so? – khajlk May 31 '17 at 19:16
4

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';
  • Thanks @Stefan. I am trying to test this code but so far unsuccessful. – khajlk May 23 '17 at 15:57
  • Which errors do you get? – Stefan May 23 '17 at 17:47
  • Well, I tested your function. The code runs but doesn't yield what was expected. There were some errors relating to return mismatch type. In CREATE TYPE you used ID as 'smallint' while in function it was 'integer' type. I modified CREATE TYPE with ID as 'integer' and it worked! However, slicing of wedges is not as expected. Maybe I need to play with the function parameters a bit more.In any case, thanks a lot! – khajlk May 27 '17 at 10:28
3

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.

  • Woahh!! your solution rocks and is wayyy faster, I must admit. Thanks a lot. – khajlk Aug 15 '17 at 20:00
2

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:

enter image description here

  • @xunik: Thanks a lot. Your code is useful and got my Upvote! but I would prefer to have a Postgis based solution. – khajlk May 27 '17 at 10:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.