[GIS] Identifying points on left or right side of street using PostGIS

postgis

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.

enter image description here

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:

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?

Best Answer

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

Results Pane

Related Question