Plot a point at the boundary of all polygons facing the nearest road line in PostGIS

Question

I want to add/plot a point geometry at the boundary of all home polygons.

Conditions:

1. The point should be facing the road geometry. As it will considered as the approximate entry point for that building.
2. The point should not be at the corner of the polygon. It should be at the edge but not at the corners. (Not Mandatory but will be good if not at the corner)

Right now I am able to implement this by first calculating the centroid of each polygon using ST_Centroid and then using ST_ClosestPoint to detect the closest point of road from the centroid. Then I draw a line starting from centroid to the closest road point using ST_MakeLine and then detect the intersection point between the created line and building polygon using ST_Intersection and then detect the 1st point of intersection to plot it as the required Entry Point.

Below is the graphical representation of my implemented logic where blue line is the road, and green dot is the plotted Entry Point of which I managed to avoid the corners using this logic.

Below is the result I get, which is perfectly ok in looks of what I needed

So, yes I am getting the output I needed BUT the execution time for this process is very high and I can have around 10000 polygon home structures at once. So it will be really awesome if someone can help me implement this same logic or better in a very time efficient way. Like plotting 10000 entry points within few seconds.

Can you provide the solution with some sample queries which can be directly executed?

Note: Right now I am doing this 1 polygon at a time. So I am looping through each polygons and executing this query 10000 times for 10000 polygons. Cause when I try to execute it together it just takes forever.

Answer 1

While looping is unavoidable, using an index driven (spatial) (K)NN search will be most significant in boosting performance.

PostgreSQL implements the concept of looping over a table with the LATERAL statement (a JOIN extension), while PostGIS comes with a specialized 2D-distance operator <-> (among others; informally known as the (K)NN operators) that utilizes the GIST index of the lookup table when used in the ORDER BY expression.

I will use the same methodology (intersection with shortest line to centroid), but this does add quite some computational overhead. Needless to say, this approach may produce undesired results, which you seem to be aware of.

Try

SELECT  ST_GeometryN(ST_Intersection(p.geom, ln.geom), 1) AS geom
FROM    <polygons> AS p
CROSS JOIN LATERAL (
  SELECT ST_ShortestLine(l.geom, ST_PointOnSurface(p.geom)) AS geom
  FROM   <lines> AS l
  ORDER BY
         p.geom <-> l.geom
  LIMIT  1
) AS ln
;

Note that

• a GIST index on <lines>.geom is mandatory
• this may fail (with a NULL geometry), or return a LineString, in the rare edge case where the shortest line partially follows a Polygon boundary

Answer 2

I've faced a very similar problem, and the solution that best worked to me was the following one:

1. First, I split the building polygons into segments.
2. Second, I find the centroid for each segment (in order to find a point that is within the boundary of the building polygon but not in a corner, which would correspond to the start or to the end of segments).
3. Third, I find the distance to the closest roads for each centroid within a distance (in my case, they where 20 meters).
4. Since the centroids belonging to the same building keep the same id, I make a query to find the segment centroid in each building that has the minimum distance to a road. This will be the point I was looking for.

The PostGIS query that performs this calculation is as follows:

    create table entrance as  
  --subquery to calculate segment centroids (steps 1 and 2)
        WITH segments_centroid AS (
        SELECT id, st_centroid(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_Boundary(geometry)) AS pt FROM building) as dumps ),
  --subquery to calculate closest roads within 20 meters to the segment centroids (step 3)
       streets_within_20m as (select segments_centroid.id,geom, 
       st_distance(segments_centroid.geom, streets.geometry) as dist from
       segments_centroid,streets  where  st_dwithin (segments_centroid.geom, 
       streets.geometry,20)),
  --subquery to calculate the closest point, among the segment centroids of each building, to a road (step 4)
       closest_point as
       (select id, min(dist) as mindist from streets_within_20m group by id)
  --final query to obtain the points I was looking for
       SELECT streets_within_20m.id, streets_within_20m.geom 
       FROM streets_within_20m, closest_point 
       where streets_within_20m.id=closest_point.id and 
       streets_within_20m.dist=closest_point.mindist;

I used this query to calculate the point for almost 28000 buildings and it took 20 seconds. This is an image of the result:

Hope it helps