Last active
December 10, 2024 14:14
-
-
Save jeaneric/fbc2a0a322cd4587e2f27d14d53cd227 to your computer and use it in GitHub Desktop.
PostGIS function to get the side of a point or geometry from a line
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
-- PostGIS function to get the side of a point or geometry from a line | |
-- ***** side function | |
-- ***** does not support m values. Only the digitizing "orientaiton" of the line. | |
-- -- inspired by: https://gis.stackexchange.com/a/156585 | |
CREATE OR REPLACE FUNCTION st2_side ( | |
geom_line geometry, | |
geom_other geometry | |
) | |
RETURNS text AS | |
$BODY$ | |
with shortestline as( | |
select ST_ShortestLine(geom_line, geom_other) as geom | |
), vec_start as ( | |
select ST_StartPoint((select shortestline.geom from shortestline)) as pt | |
), dist_along as( | |
select ST_LineLocatePoint(geom_line, vec_start.pt) as dist_along | |
from vec_start | |
), new_line_pt1 as( | |
SELECT | |
case when dist_along > 0.999 then | |
ST_LineInterpolatePoint(geom_line, 0.999) | |
when dist_along < 0.001 then | |
ST_StartPoint(geom_line) | |
else | |
vec_start.pt | |
end as pt | |
from dist_along, vec_start | |
), new_line_pt2 as( | |
SELECT | |
case when dist_along > 0.999 then | |
vec_start.pt | |
when dist_along < 0.001 then | |
ST_LineInterpolatePoint(geom_line, 0.001) | |
else | |
ST_LineInterpolatePoint(geom_line, dist_along + 0.001) | |
end as pt | |
from dist_along, vec_start | |
), new_line as( | |
select ST_MakeLine(new_line_pt1.pt, new_line_pt2.pt) as new_line | |
from new_line_pt1, new_line_pt2 | |
), h as( | |
select | |
shortestline.geom vec, | |
new_line seg | |
from shortestline, vec_start, new_line | |
), az as( | |
select (ST_Azimuth(ST_StartPoint(h.vec), ST_EndPoint(h.vec)) - ST_Azimuth(ST_StartPoint(h.seg), ST_EndPoint(h.seg))) as az | |
from h | |
) | |
select --az::text | |
case when az = 0 or az is null THEN | |
'middle' | |
when az > PI() or (az<0 and az > (-1*PI())) THEN | |
'left' | |
ELSE | |
'right' | |
end | |
from az | |
$BODY$ | |
LANGUAGE sql; | |
-- *** Testing: | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 2 2)'), ST_GeomFromEWKT('POINT(1 1)')); -- middle | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 2 2)'), ST_GeomFromEWKT('POINT(1 0)')); -- right | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 2 2)'), ST_GeomFromEWKT('POINT(0 1)')); -- left | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 2 2)'), ST_GeomFromEWKT('POINT(-1 0)')); -- left | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 2 2)'), ST_GeomFromEWKT('POINT(0 -1)')); -- right | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 2 2)'), ST_GeomFromEWKT('POINT(4 0)')); -- right | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 1 0)'), ST_GeomFromEWKT('POINT(0.5 1)')); -- left | |
SELECT st2_side( ST_GeomFromEWKT('LINESTRING(0 0, 1 0)'), ST_GeomFromEWKT('POINT(0.5 -1)')); -- right |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment