7

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?

  • 1
    I don't really have time to look at this properly at the moment, but ST_ShortestLine and ST_LineCrossingDirection could hold some promise. – MickyT Nov 29 '17 at 22:40
  • @MickyT: Thank you for the suggestion. Care to explain a bit? – khajlk Nov 30 '17 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. – MickyT Nov 30 '17 at 21:15
  • @MickyT: Tested and provided feedback to your suggestion. – khajlk Nov 30 '17 at 21:51
6

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

  • 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. – khajlk Nov 30 '17 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. – MickyT Dec 3 '17 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? – khajlk Dec 3 '17 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. – MickyT Dec 4 '17 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. – khajlk Dec 4 '17 at 7:59
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;

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.