[GIS] Determining what side of line point is on using PostGIS

postgis

I want this in order to generate house numbers, odd or even depending on wether a feature is on one side of a street line or not . Thanks

In the end i think the most portable way is to use this( BTW line is actually 'linestring' not segment):

ST_LineCrossingDirection( linefromAtoB,streetline) <0 where A is the point , B is ST_ShortestPoint(thepoint,streetline); if <0 then linefromAtoB enters through the left of The streetline; if >0 the it enters from right side of streetline). You have to determine what street line is closest to each point first.

I made a mistake:
C is the symmetric of A in relation to B= ST_ShortestPoint(thepoint,streetline).
C's coordinates are XC=2*(XB-XA) -XA and for y YC=2*(YB-YA) -YA.
ST_LineCrossingDirection( linefromAtoC,streetline) <0 where A is the point , ; if <0 then linefromAtoC enters through the left of The streetline; if >0 the it enters from right side of streetline). You have to determine what street line is closest to each point first.

Best Answer

If you are doing this in code, start by getting the closest point on the line to the point:

http://paulbourke.net/geometry/pointlineplane/

Note you will need to test each segment of a LINESTRING, and make sure the point falls on the segment as the raw computation is for an infinite line.

Then get the vector from this point to the line (pt.x - intersection.x, pt.y - intersection.y)

This will be a perpendicular to the line, within floating point rounding errors.

A line with slope x, y (end.x - start.x, end.y - start.y) has two perpendiculars.

-y, x is to the left, y, -x is to the right. Compare to your vector. You will only need to compare the signs of the x and y components.

There are other ways to do this and if you need very high performance this may not be the fastest.

Related Solutions

PostGIS – Detecting Point Location Relative to Line: Left or Right Side

select (ST_Azimuth(h.vec) - ST_Azimuth(h.seg))
from (
    select 
        ST_MakeLine(cp.p, point.geom) vec,
        ST_MakeLine(cp.p, 
            ST_LineInterpolatePoint(
                line.geom, 
                ST_LineLocatePoint(line.geom, cp.p) * 1.01)
        ) seg
        from (
            select 
                ST_ClosestPoint(line.geom, point.geom)
        ) p as cp
    ) as h

So the idea is to calculate angle between closest line segment, and vector from closest point on the line to your point.

get a closest point on a line

select ST_ClosestPoint(line.geom, point.geom)

create the vector from closest point to your point

ST_MakeLine(cp.p, point.geom) vec

create a vector among your line

ST_MakeLine(
    --original point
    cp.p, 
    --find a point next to the closest point on line
    ST_LineInterpolatePoint(line.geom, 
         ST_LineLocatePoint(line.geom, cp.p) * 1.01)) seg

get the difference between directions

ST_Azimuth(h.vec) - ST_Azimuth(h.seg)

So right and left will be greater than zero and lower than zero.

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

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

Best Answer

Related Solutions

PostGIS – Detecting Point Location Relative to Line: Left or Right Side

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

Related Question