Sunday, 24 April 2016

Detecting if point is on left or right side of line in PostGIS?


I have a linestring table and a point table in postgis.


I know the closest line to any given point. What I need to know is on which "side" of that line is the point. I guess I have to do that by creating a perpendicular line from given point to the line (closest point on the line) and then compare the coordinates, but I don't know exactly how to do that, and if it's the proper way, since line changes it's direction.


I've made a picture to illustrate my task.


enter image description here


The line itself is black, its direction is shown with green arrows. I need to add a "side" column to the point table, so that red points should have value "right" and blue points should have value "left".


Can someone give an SQL code example of calculating a "side" value of a point?



Answer




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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...