Identifying points on left or right side of street using 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.

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:

``````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?

• I don't really have time to look at this properly at the moment, but `ST_ShortestLine` and `ST_LineCrossingDirection` could hold some promise. Commented Nov 29, 2017 at 22:40
• @MickyT: Thank you for the suggestion. Care to explain a bit? Commented Nov 30, 2017 at 7:31
• I have added an answer as an explanation, but please note I have not tested this theory and it may not work. Commented Nov 30, 2017 at 21:15
• @MickyT: Tested and provided feedback to your suggestion. Commented Nov 30, 2017 at 21:51

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_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_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

• Thank you! I just tested your suggested code. It returned 0, -1, 0 for all three points. Means, shortest lines of points 1 and 3 did not cross. So, probably, your idea "may need to extend the shortest line to ensure it crosses" is true. Commented Nov 30, 2017 at 21:48
• @JibranKhan I have tested and updated the query to use extensions on the lines to ensure crossings. Seems to work OK except points which sit directly on the extensions of the street line. Commented Dec 3, 2017 at 19:28
• I appreciate your efforts. Thank you. I just executed the whole same thing suggested by you and 'line' and 'side' columns are empty. I understand this is for test data. But, there should have been something in the output - at least - in the 'side' column. Any idea what's happening? Commented Dec 3, 2017 at 22:17
• @JibranKhan which versions are you working with? I have Postgis 2.3.3 running on PostGreSQL 9.6. The query above has a single null for vec and side when I run it. I have not encountered nulls on the line column. Commented Dec 4, 2017 at 0:18
• That's strange. I am using PostgreSQL 9.5 with PostGIS 2.2. I have Postgres 9.6 installed at other machine. Will check and get back. Commented Dec 4, 2017 at 7:59

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;
``````