Splitting polygon in equal parts based on polygon area?

Question

So, I came across this relevant question, but I think my problem is a bit different.

I have two tables i.e. polygon poly and points pts in my Postgres db. I am on Windows 10 (x64) based machine and using PostgreSQL 9.5.12/PostGIS 2.3. Some of the polygons contain points inside them while, some do not. The sample scenario is given as follows:

For each poly, the following query computes the area and counts number of pts inside polygon.

Select 
    a.gid poly_id, count(b.geom) pt_cnt,
    st_area(a.geom)::int poly_area
FROM
    poly a
LEFT JOIN
    pts b 
ON st_contains(a.geom, b.geom) 
GROUP BY a.gid
Order by a.gid;

The area of the polygons (from left to right) is 1079, 744, 340 meter-square respectively. Now, I would like to split these polygons' geometry based on the area i.e. if area is above 100 meter-square, then each polygon needs to be divided/split into approximate equal parts. Can somebody suggest me a workaround to fix the issue?

UPDATE

Here, is another example of the polygons in my data. I am aware that there can be much more complex examples (e.g. irregular polygons). In that case, I think exceptions can be added.

while this is reasonably straight forward for the (very strict) cases of your sample scenario's geometries (i.e. being simple rectangles), it actually becomes exponentially more complex if those polygons are irregular and the result set has to be even only approximately equal...how exactly do your real case polygons look like? — geozelot, Jun 13 '18 at 13:59
@ThingumaBob: Yes, you are right. In the real world scenario, polygons can be of irregular shapes. In my case, of course there are irregular polygons (I am going to add one more example scenario). Nevertheless, I would like to have, at least, a very rough approximation (not necessarily equal) in terms of split/sub-divided polygons. — khajlk, Jun 13 '18 at 14:37
I used two approaches for different cases with good enough results for me. in your case, see if you can use ST_ApproximateMedialAxis (you need postgis_sfcgal extension); you'd need to extend the resulting line using ST_ClosestPoint from it's start/end points to the polygons boundary and split it with that line as blade. might work, sry I can't get into detail right now. — geozelot, Jun 13 '18 at 17:46

score 2 · Accepted Answer · answered Jan 08 '19 at 10:13

@ThingumaBob I've tried your approach with different results. In this case I've found a rectangle not being split up equally, since the distance from the endpoint of the medial axis is equal to all corners. So I've adapted the approach based on "How to extend a straight line in postgis?"

Result closest point method:

Result new extended line method:

Query closest point:

 with medialaxis_multi as(
    select
        st_union(ST_ApproximateMedialAxis(geom)) as geom
    from
        polygons
),
medialaxis_single as(
    select
        (st_dump(geom)).geom as geom
    from
        medialaxis_multi
),
closestpoints as(
    select 
        st_closestpoint(st_boundary(polygons.geom),st_startpoint(medialaxis_single.geom)) as startpoint,
        st_closestpoint(st_boundary(polygons.geom),st_endpoint(medialaxis_single.geom)) as endpoint
    from 
        polygons,
        medialaxis_single
)
    select 
        st_makeline(startpoint, endpoint) 
    from 
        closestpoints
;

Working Query based on extended medialaxis: See also my Gist

-- blades part adapted from "https://gis.stackexchange.com/questions/104439/how-to-extend-a-straight-line-in-postgis"
with medialaxis_multi as(
    select
        polygons.id,
        st_union(ST_ApproximateMedialAxis(geom)) as geom
    from
        polygons
    group by 
        polygons.id
),
medialaxis_single as(
    select
        id,
        (st_dump(geom)).geom as geom
    from
        medialaxis_multi
),
blades as (
    select 
        id, 
        st_makeline(st_translate(a, sin(az1) * len, cos(az1) * len),st_translate(b,sin(az2) * len, cos(az2) * len)) as geom
    from (
            select 
                id,
                a,
                b, 
                st_azimuth(a,b) as az1, st_azimuth(b, a) as az2, st_distance(a,b) + 1 as len
            from (
                select 
                    id,
                    st_startpoint(geom) as a, 
                    st_endpoint(geom) as b
                from 
                    medialaxis_single
            )   as sub
        )   as sub2
)
select 
    blades.id,
    (st_dump(st_split(polygons.geom,blades.geom))).geom
from    
    polygons,
    blades
where
    polygons.id = blades.id
  ;

Splitting polygon in equal parts based on polygon area?

1 Answers1