Finding closest Intersection between line strings (PostGIS)

Question

I want to find the distance between one point (excavation sites) and a line (rivers) for the cardinal directions. What I thought would be the easiest way is to draw lines between the excavation site and a given maximum distance every 45° and find the intersection between the line and the river which isn't working properly yet (and maybe isn't the best idea to solve that?).

I have a table excavation_sites and a table dist_ex_sites_river (in which the lines and the distances should be stored). The dist_ex_sites_river table also has a foreign key for the excavation sites. I managed to create the lines already and write them in the db with the according step_nr (see figure - the red circles are the intersection points I want to find for measuring the distance) and the excavation site. Now I'd like to get the intersections to measure the distance which always returns None, maybe you guys can help me? This is the code I'm using right now:

SELECT ST_Intersection(dist_ex_sites_river.line_geometry,env_white_water.geom) FROM public.env_white_water, public.dist_ex_sites_river WHERE dist_ex_sites_river.step_nr_id_pk = 1 AND env_river.river_id_pk = 616;

I'm working with python and the psycopg2 package. I iterate through every id (in both tables) and want to append the distances to a list and choosing the minimum value in the end to insert that into the db in the end. Probably this is absolutely stupid and there is a much better way to do that. Thanks a lot for your help.

Answer 1

Given a point (excavation site), you can use an SQL-language function to project lines out in the cardinal directions, to a specified radius:

CREATE OR REPLACE FUNCTION cardinal_lines (g geometry, r float) 
RETURNS TABLE (dir varchar(2), geom geometry) AS
$$
    SELECT 'NE'::varchar(2), ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g) + sqrt(r/2), ST_Y(g) + sqrt(r/2)), ST_SRID(g)))
    UNION ALL
    SELECT 'SE'::varchar(2), ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g) + sqrt(r/2), ST_Y(g) - sqrt(r/2)), ST_SRID(g)))
    UNION ALL
    SELECT 'SW'::varchar(2), ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g) - sqrt(r/2), ST_Y(g) - sqrt(r/2)), ST_SRID(g)))
    UNION ALL
    SELECT 'NW'::varchar(2), ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g) - sqrt(r/2), ST_Y(g) + sqrt(r/2)), ST_SRID(g)))
    UNION ALL
    SELECT 'N'::varchar(2),  ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g)    , ST_Y(g) + r), ST_SRID(g))) 
    UNION ALL
    SELECT 'E'::varchar(2),  ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g) + r,  ST_Y(g)   ), ST_SRID(g)))
    UNION ALL
    SELECT 'S'::varchar(2),  ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g)    , ST_Y(g) - r), ST_SRID(g)))
    UNION ALL
    SELECT 'W'::varchar(2),  ST_MakeLine(g, ST_SetSRID(ST_MakePoint(ST_X(g) - r, ST_Y(g)    ), ST_SRID(g))) 
$$ language 'sql';

Using this function, you can write a query to find, for each excavation site, the closet river line in each of the cardinal directions:

SELECT DISTINCT ON (id, dir)
  id, 
  dir, 
  roads.linearid,
  ST_Distance_Sphere(
    ST_StartPoint(cardinal_lines.geom),
    ST_Intersection(cardinal_lines.geom, roads.geom)) 
FROM
  (SELECT id, (cardinal_lines(geom, 1)).* FROM ex_sites) cardinal_lines
  LEFT JOIN roads ON ST_Intersects(cardinal_lines.geom, roads.geom)
ORDER BY id, dir, st_distance_sphere ASC

(Note that I used a table of roads here - just replace this with your rivers table.)

The function above could be modified to return polygon slices instead of lines, so that "northeast" means 22.5 to 37.5 degrees, for example.