Lessen the overlapping in Polygons created on a road in PostgreSQL

Question

I am working in an application of PostgreSQL. I have a list of points and want to create square of 500 sq. feet around these points. These points are in Geometry (EPSG:4326) format.Every point is center of square and final size of square should not be changed.

I have created a script in PostgreSQL to do this. But the problem is that on few places squares are more overlapped. I want to lessen the overlapping without changing there size and neither point should be outside of any square.Following data can be used for testing:

-84.26929 39.18387,-84.26929 39.18389,-84.26921 39.18399,-84.26613 39.18531,-84.26585 39.18542,-84.24252 39.19028,-84.24153 39.19059,-84.2408 39.19082,-84.23963 39.19123,-84.23803 39.19266,-84.23773 39.19336

Any idea?

so you have created an envelope of 500 sft around each point and if two points are within the respective threshold, their rectangles overlap. you want those rectangles evenly distributed along the line instead, covering every point but without any rectangles overlapping each other? does it matter if a point is on the edge of a polygon (so all is good as long as all points are somehow covered by a rectangle?) — geozelot, Aug 11 '17 at 11:14
Rectangle can overlap but they should not overlap more than 50 sqft. yes, it is true that all point should be covered by rectangle. If a few points are on edge that can be consider. But in every case the size of rectangle should be 500 sft. — M.Sharma, Aug 11 '17 at 12:09
Can you tell a bit about the exact use-case you are dealing with? I can't really understand the situation. You might also add a drawing that explains the situation. — tilt, Aug 11 '17 at 13:43
Yes, definitely I can show you a drawing. But just tell me how to add that in my post. — M.Sharma, Aug 11 '17 at 13:58
I think you would be better to start with a grid of 500sqft that covers the data. From there you can drop empty grid cells and maybe even slide grids lines to optimize the structure. — MickyT, Aug 14 '17 at 01:34
I could probably come up with a recursive query to do this, but it would be horribly slow. Basic steps would be create a square for each point. Keep the one with the most points. Create a square for every point outside. Keep the one with least overlap area that is >= 0. Repeat. The point count vs overlap are may want to considered as well. — MickyT, Aug 14 '17 at 02:09
The main problem I feel is that how to rotate a polygon to road side in my example. If I could do that , I believe in that case polygons will overlap less. — M.Sharma, Aug 14 '17 at 09:55
I think MickyT has a good point with starting with a grid and throwing out grid cells that are not overlapping any points. It will not be the best fit but it would fulfill your question. Here is a post that can bring you on track: https://gis.stackexchange.com/questions/16374/how-to-create-a-regular-polygon-grid-in-postgis — tilt, Aug 14 '17 at 13:14
@tilt I have seen the link. However, my problem is different. The problem in link is generating series to generate cell from polygon but I have points. I have successfully generated cells. The only thing left is to get azimuth and rotate cells along road. I have stucked here. — M.Sharma, Aug 14 '17 at 16:21
Why do you need to rotate the cells? Your original question just says the cells should not overlap, which MickyT's suggestion of using a grid will achieve. If you start with small grid cells, you can also use ST_Union to consolidate them into larger rectangles. — amball, Aug 16 '17 at 06:45
OK, I suggest editing your post to make this clear. Your diagram doesn't suggest that you want to match the bearing of the road. — amball, Aug 16 '17 at 07:52

MickyT · Accepted Answer · 2017-08-17T20:36:24.887

This is not a complete answer at the moment, but rather an indication of what could be done. I originally worked through this in SQL Server (for visualization and then ported it to PostGIS. I have put bothe versions in.

Disclaimers

First, you will need to make lines from the points that you have collected.
At the moment I haven't included that, but have concentrated on getting minimal overlapping rotated (as per comments) rectangles.

The following is based on a single line, but can be changed to handle multiple lines with out to many problems.

I have done it based on units only and the rectangle I am generating is 200 units as it suited the dummy line I made. Feel free to add some some data to your question.

Due to the way SQL Server STDifference works I have reversed the order that I segment the line in.

Explanation

This is a recursive query for a single line currently. It segments the the line from the end of the line to the beginning into sections of 200 units. This can be look at as a very specific point reduction. The last segment is not currently extended, but could quite easily be.

Each of these segments is then turned into a square of 200 x 200 units (40,000 sq units).

PostGIS

/*
 Altered to handle input from a geography datatype.
 This means chopping and changing between types to get around restrictions
*/
WITH RECURSIVE SEGMENT_LINE(x1,y1,Ring,X2,Y2,myLine) AS (
    SELECT 
        /*Initial coorindate for the line*/
        ST_X(ST_PointN(myLine::Geometry,1)) X1, ST_Y(ST_PointN(myLine::Geometry,1)) Y1,
        /*The ring for determining where to split line.  ** NOT REQUIRED ** */
        ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry)::Geography Ring,
        /*Case statements to deal with end of line*/
        CASE WHEN ST_Intersects(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry) THEN
            ST_X(ST_Intersection(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry))
        ELSE
            ST_X(ST_PointN(myLine::Geometry, -1))
        END X2,
        CASE WHEN ST_Intersects(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry) THEN
            ST_Y(ST_Intersection(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry))
        ELSE
            ST_Y(ST_PointN(myLine::Geometry, -1))
        END Y2,
        /*Split and remove the first portion of the line. */
        CASE WHEN ST_Intersects(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine) THEN
            ST_GeometryN(ST_Split(
                /*Nodes the line for the split at point.  Deals with floating point errors*/
                ST_LineMerge(ST_Split(myLine::Geometry,ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry))),
                ST_Intersection(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry)
            ),2)::Geography
        ELSE
            NULL
        END myLine
    FROM (VALUES (
        ST_GeographyFromText('SRID=4326;LINESTRING(-84.26929 39.18387,-84.26929 39.18389,-84.26921 39.18399,-84.26613 39.18531,-84.26585 39.18542,-84.24252 39.19028,
                             -84.24153 39.19059,-84.2408 39.19082,-84.23963 39.19123,-84.23803 39.19266,-84.23773 39.19336)'))
         ) AS Sample(myLine)
    UNION ALL
    SELECT 
        /*Initial coorindate for the line*/
        ST_X(ST_PointN(myLine::Geometry,1)) X1, ST_Y(ST_PointN(myLine::Geometry,1)) Y1,
        /*The ring for determining where to split line.  ** NOT REQUIRED ***/
        ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry)::Geography Ring,
        /*Case statements to deal with end of line*/
        CASE WHEN ST_Intersects(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry) THEN
            ST_X(ST_Intersection(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry))
        ELSE
            ST_X(ST_PointN(myLine::Geometry, -1))
        END X2,
        CASE WHEN ST_Intersects(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry) THEN
            ST_Y(ST_Intersection(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry))
        ELSE
            ST_Y(ST_PointN(myLine::Geometry, -1))
        END Y2,
        /*Split and remove the first portion of the line. */
        CASE WHEN ST_Intersects(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine) THEN
            ST_GeometryN(ST_Split(
                /*Nodes the line for the split at point.  Deals with floating point errors*/
                ST_LineMerge(ST_Split(myLine::Geometry,ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry))),
                ST_Intersection(ST_ExteriorRing(ST_Buffer(ST_PointN(myLine::Geometry,1)::Geography,500 * 0.3048)::Geometry),myLine::Geometry)
            ),2)::Geography
        ELSE
            NULL
        END myLine
    FROM SEGMENT_LINE
    WHERE myLine is not null
    )
SELECT  X1, Y1, X2, Y2, Ring,
    ST_MakeLine(ST_MakePoint( X1 , Y1 ),ST_MakePoint( X2 , Y2 )) seg,
    ST_MakePolygon(
        ST_MakeLine(ARRAY[
        ST_MakePoint( X1 - (Y2 - Y1) / 2 , Y1 + (X2 - X1) / 2 ),
        ST_MakePoint( X1 + (Y2 - Y1) / 2 , Y1 - (X2 - X1) / 2 ),
        ST_MakePoint( X2 + (Y2 - Y1) / 2 , Y2 - (X2 - X1) / 2 ),
        ST_MakePoint( X2 - (Y2 - Y1) / 2 , Y2 + (X2 - X1) / 2 ),
        ST_MakePoint( X1 - (Y2 - Y1) / 2 , Y1 + (X2 - X1) / 2 )]
        )
    ) Rectangle
FROM SEGMENT_LINE;

SQL Server 2012+

WITH /*RECURSIVE*/ SEGMENT_LINE AS (
    SELECT 
        Line.STPointN(Line.STNumPoints()).STX X1, Line.STPointN(Line.STNumPoints()).STY Y1,
        Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing() Ring,
        CASE WHEN Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersects(Line) = 1 THEN
            Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersection(Line).STPointN(1).STX
        ELSE
            Line.STPointN(1).STX
        END X2,
        CASE WHEN Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersects(Line) = 1 THEN
            Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersection(Line).STPointN(1).STY
        ELSE
            Line.STPointN(1).STY
        END Y2,
        CASE WHEN Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersects(Line) = 1 THEN
            Line.STDifference(Line.STPointN(Line.STNumPoints()).STBuffer(200))
        ELSE NULL END Line
    FROM (VALUES(Geometry::STGeomFromText('LINESTRING(0     0      ,100   500    ,500     600    ,750     300    ,900   700 )',0)))A(Line)
    UNION ALL
    SELECT 
        Line.STPointN(Line.STNumPoints()).STX X1, Line.STPointN(Line.STNumPoints()).STY Y1,
        Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing() Ring,
        CASE WHEN Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersects(Line) = 1 THEN
            Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersection(Line).STPointN(1).STX
        ELSE
            Line.STPointN(1).STX
        END X2,
        CASE WHEN Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersects(Line) = 1 THEN
            Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersection(Line).STPointN(1).STY
        ELSE
            Line.STPointN(1).STY
        END Y2,
        CASE WHEN Line.STPointN(Line.STNumPoints()).STBuffer(200).STExteriorRing().STIntersects(Line) = 1 THEN
            Line.STDifference(Line.STPointN(Line.STNumPoints()).STBuffer(200))
        ELSE NULL END Line
    FROM SEGMENT_LINE
    WHERE LINE IS NOT NULL
    )
SELECT  X1, Y1, X2, Y2, Ring,
    Geometry::STGeomFromText(CONCAT('LINESTRING(',X1,' ',Y1,', ',X2,' ',Y2,')'),0) Line,
    Geometry::STGeomFromText(CONCAT('POLYGON((',
    X1 - (Y2 - Y1) / 2, ' ' , Y1 + (X2 - X1) / 2, ', ',
    X1 + (Y2 - Y1) / 2, ' ' , Y1 - (X2 - X1) / 2, ', ',
    X2 + (Y2 - Y1) / 2, ' ' , Y2 - (X2 - X1) / 2, ', ',
    X2 - (Y2 - Y1) / 2, ' ' , Y2 + (X2 - X1) / 2, ', ',
    X1 - (Y2 - Y1) / 2, ' ' , Y1 + (X2 - X1) / 2,
    '))'),0) Rectangle
FROM SEGMENT_LINE;

For a line like this we get rectangles as followsIt's not perfect, but as stated above if this is the sort of thing you are looking for then it can be improved.

Yes, I was looking the same thing in PostGIS. I will try to convert your logic into PostGIS. However, this may be difficult for me to find all matching function , you have used in SQL SERVER. — M.Sharma, Aug 16 '17 at 06:45
I tried your code in SQL SERVER 2014 and get the following error :
Error converting data type varchar to float. Before this I change CONCAT function to plus because that is also not working in SQL SERVER 2014. — M.Sharma, Aug 16 '17 at 08:47
@M.Sharma CONCAT should work in versions of SQL Server above 2012. If you change it to +ing the string you will need to wrap CASTs around the coordinate expressions. Better still use FORMAT to ensure the number isn't munged. so something like 'POLYGON((' + FORMAT(X1+(Y2-Y1)/2,'#######.#######') + ... — MickyT, Aug 16 '17 at 19:52
@M.Sharma I'll look at converting this to PostGIS later today — MickyT, Aug 16 '17 at 19:53
T: I have given some sample data in my question. With that data your logic in PostGIS is not working. Your logic is not creating cells. It gives me a single entry after recursion. — M.Sharma, Aug 17 '17 at 11:12
@M.Sharma changed to handle Geography with Latitude /Longitude — MickyT, Aug 17 '17 at 20:37
I can not say why but your logic fails and cell size is not consistent, if I test with more data. I test this for 100 and 400 Points. Your algorithm is good. But did you check my logic to draw these cells.However, In my logic you will have to store cells in a table to do query to find intersected cells. — M.Sharma, Aug 18 '17 at 11:55
@M.Sharma That is because I used cartesian methods on a spherical coordinates to build the final square. You should find that the length of the segment created is correct. I would strongly suggest that you look at using a projection for your data. Or you could project your data do the gridding, then unproject it. — MickyT, Aug 18 '17 at 19:56
can you help me on this https://gis.stackexchange.com/questions/274424/alternative-of-st-translate-function-in-sql-server-2016 — M.Sharma, Mar 14 '18 at 07:43
@M.Sharma the answer given by DPSSpatial is the one I would have provided — MickyT, Mar 15 '18 at 17:38

Lessen the overlapping in Polygons created on a road in PostgreSQL

1 Answers1