Identifying points on left or right side of street using PostGIS?

Question

The scenario:

I have points (addresses) and multiline strings (streets) i.e., ap and street tables in my PostgreSQL 9.5 db. Table ap contains 3 rows (sample points) while street contains 85 rows. The sample scenario is shown below.

What I need:

For each point ap, I need to identify whether it is on "left side" or "right side" of the street.

The code:

On GIS SE, I found some relevant questions already. For example:

• Detecting if point is on left or right side of line in PostGIS? and
• Determining what side of line point is on using PostGIS?

The following code is adapted from the first link.

With cp as
(
Select
    a.gid ap_id,
    a.geom ap,
    ST_SetSRID(ST_LineMerge(b.geom), 3044) line,
    ST_SetSRID(ST_ClosestPoint(b.geom, a.geom), 3044) p
from
    ap a
Left Join street b on
ST_Dwithin(ST_SetSRID(b.geom, 3044), ST_SetSRID(a.geom, 3044), 25)
Order by a.gid, ST_Distance(b.geom, a.geom)
)
, h as
(
Select
    ap_id, ap,
    line, p,
        ST_MakeLine(p, ap) vec,
        ST_MakeLine(p,
        ST_LineInterpolatePoint(line,
            ST_LineLocatePoint(line, p) * 1.01)) seg
From cp
Order by ap_id
)
Select
    ap_id,
    degrees(ST_Azimuth(ST_StartPoint(vec), ST_EndPoint(vec))) 
        - degrees(ST_Azimuth(ST_StartPoint(seg), ST_EndPoint(seg))) diff
from h
order by ap_id

The output:

The above code (for the sample scenario) works and the output is:

ap_id    diff
1       -269.999999988824
2        89.9999999946886
3        269.999999980956

Error:

However, the code for another point dataset (87 rows), throws following error:

ERROR:  line_interpolate_point: 2nd arg isn't within [0,1]

Also, the above code doesn't seem to work for cases when street is east-west having bearing of 90 degrees (Reference, here). I need to perform this task of identifying points for a large dataset. Thus, I am looking for suggestions to get through this problem. Is there any alternative way to identify whether points are on left side or right side of the street?

Answer 1

This makes use of your current query, however I have put ST_ShorttestLine of the ST_ClosestPoint and have added ST_LineCrossingDirection to the output. I would expect it to return either 1 or -1 for left or right. It may however need the shortest line extended to ensure that it crosses. On very simple geometries this appears to work unless the point is ahead of the line (eg closest point is the last point in street)

With cp as
(
Select
    a.gid ap_id,
    a.geom ap,
    ST_SetSRID(ST_LineMerge(b.geom), 3044) line,
    ST_SetSRID(ST_ShortestLine(b.geom, a.geom), 3044) vec
from
    ap a
Left Join street b on
ST_Dwithin(ST_SetSRID(b.geom, 3044), ST_SetSRID(a.geom, 3044), 25)
Order by a.gid, ST_Distance(b.geom, a.geom)
)
Select *,
    ST_LineCrossingDirection(vec,line) side
from cp;

Taking into account that floating point error is likely to start causing some issues and to try and address leading points I've adjusted the query to extend lines by adding points. The extensions are based on this answer. This appears to affect the points that are on the axis, but they may not be important for you.

I have left my test data in the query.

WITH ap AS (
SELECT gid
    ,ST_SetSRID(ST_MakePoint(x,y),3044) geom
    ,description
FROM (VALUES
     (1, -1, 0, 'Left Before')
    ,(2, 0, 1, 'Left Middle')
    ,(3, 2, 3, 'Left After')
    ,(4, 0, -1, 'Right Before')
    ,(5, 1, 0, 'Right Middle')
    ,(6, 3, 2, 'Right After')
    ,(7, -1, -1, 'On Before')
    ,(8, 1, 1, 'On Middle')
    ,(9, 3, 3, 'On After')
    ) ap(gid, x, y, description)
),
street AS (
SELECT ST_SetSRID(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(2,2)),3044) geom
),
--With 
       cp as
(
Select
    a.gid ap_id,
    a.geom ap,
    a.description,
    ST_SetSRID( -- Line with an extension added to it
        ST_AddPoint(
            ST_MakeLine(
                a.geom,
                ST_ClosestPoint(b.geom,a.geom)
            ),
            ST_Translate(
                a.geom, 
                sin(ST_Azimuth(a.geom, ST_ClosestPoint(b.geom,a.geom))) * ST_Distance(a.geom, ST_ClosestPoint(b.geom,a.geom)) * 1.01,
                cos(ST_Azimuth(a.geom, ST_ClosestPoint(b.geom,a.geom))) * ST_Distance(a.geom, ST_ClosestPoint(b.geom,a.geom)) * 1.01
            )
        )
    , 3044) vec,
    ST_SetSRID(
        ST_AddPoint(
            ST_LineMerge(b.geom),
            ST_Translate(
                ST_PointN(ST_LineMerge(b.geom),-2),
                sin(ST_Azimuth(ST_PointN(ST_LineMerge(b.geom),-2), ST_PointN(ST_LineMerge(b.geom),-1))) * (ST_Distance(ST_PointN(ST_LineMerge(b.geom),-1), ST_PointN(ST_LineMerge(b.geom),-2)) *1.1),
                cos(ST_Azimuth(ST_PointN(ST_LineMerge(b.geom),-2), ST_PointN(ST_LineMerge(b.geom),-1))) * (ST_Distance(ST_PointN(ST_LineMerge(b.geom),-1), ST_PointN(ST_LineMerge(b.geom),-2)) *1.1)
            )
        )
    , 3044) line
from
    ap a
Left Join street b on
ST_Dwithin(ST_SetSRID(b.geom, 3044), ST_SetSRID(a.geom, 3044), 25)
Order by a.gid, ST_Distance(b.geom, a.geom)
)
Select *,
    ST_LineCrossingDirection(line,vec) side
from cp;

Giving the following result

Answer 2

To check if the point is on the left or right side in respect to the direction of the street. You just need to create a left or right line buffer with the same distance to the point. If the point intersects with the buffer then it should be on the left, else right. This is the example code. You can replace intersects with any other functions.

SELECT  
 CASE
   WHEN ST_Intersects (
      ST_Buffer(line_geom, 
         ST_Distance(line_geom, point_geom)*1.1, 'side=left endcap=flat'),
      point_geom) = TRUE
  THEN 'T'
  ELSE 'F'
  END AS isleft
FROM some_table;