1

I am attempting to find the nearest neighbour to a geometry in a certain direction but am becoming stuck at the planning stage. I understand the basic tools that I need but am failing to find anything that can limit on a direction range.

I am finding the nearest neighbour and getting the coordinates of the two closest points as so:

SELECT
    tb1.degrees
    --taking ST_StartPoint and ST_EndPoint to get my closest points
    ,ST_ShortestLine(
        tb1.geom
        ,(
            SELECT geom
            FROM table2 tb2
            ORDER BY tb1.geom <-> tb2.geom ASC
            LIMIT 1
        )
    )
FROM table1 tb1;

But I do not want to find the nearest neighbour to the 'left' of this geometry, only right. When I say left, I actually mean (degrees-90) of the bearing for tb1.geom.

So how do I find the nearest neighbour but only in directions n°+180? Or, is this possible?

geozelot
  • 30,050
  • 4
  • 32
  • 56
Phish
  • 165
  • 6
  • I am airing on drawing a line from centroid of geom 1 at n°+90 for a distance that I know will intersect some and then finding the nearest neighbour from the returned list but it feels wildly unhygienic. – Phish Mar 22 '21 at 20:38
  • So why don't you add filter by X coordinate of centroids coords of other geometries greater than the X coordinate of centroid of your base geom ? – J. Monticolo Mar 22 '21 at 21:15
  • I may be being stupid but the left and right is only relevant to the direction/bearing of the geometry. – Phish Mar 22 '21 at 21:22
  • Your geometry is in geodesic system like WGS84 ? So yes, my solution will work (I think) only with projected coordinate system. – J. Monticolo Mar 22 '21 at 21:32
  • add the azimuth to your query https://gis.stackexchange.com/questions/77353/how-to-calculate-the-azimuth-of-a-line-of-a-multilinestring-in-postgis – Mapperz Mar 22 '21 at 22:08
  • @J.Monticolo No I am in Cartesian. I do use things like st_project but transforming to web mercator short of doing the trig later on (when I have learnt how) – Phish Mar 22 '21 at 22:12
  • @Mapperz I am not sure how this will aid the question but I will – Phish Mar 22 '21 at 22:13
  • I calculate the bearing using Azimuth in tb1 – Phish Mar 22 '21 at 22:14
  • @J.Monticolo The bearing may be 10 degrees, I need to know what is left (which in this case left will be 10-191 degrees). I may be going in the opposite direction so 190 degrees (which in this case left will be 190-10 degrees). – Phish Mar 22 '21 at 22:16
  • So, the only solution I saw is yours, to construct an oriented box to contains "right" geometries as a pre-filter before the nearest neighbour. Or draw a line, do a ST_BUFFER only on one side (see https://postgis.net/docs/ST_Buffer.html) and add a filter to contains geometries. – J. Monticolo Mar 23 '21 at 08:30
  • @J.Monticolo Thanks for the response. Buffer side wont work as I cannot let PG to decide what is left and right but the line may work. Centroid of polygon is going to potentially be too far "left" so I will draw the line from the centroid, interpolate when it intersects the exterior ring of the polygon and then draw a line from this point at the degree+90 bearing again. Then I can see what geometries this line intersects with and then get the closest. As I mentioned previously I feel like I need a shower, that the solution is not clean but perhaps as a first solution it will give something. – Phish Mar 23 '21 at 10:00
  • I am extremely lucky that the geometries I am trying to find are long travel network LineStrings so I KNOW that it will intersect with what I need. If this solution was required to find anything within the bearing range this solution would not work. Perhaps at a known likley distance you could draw the line and buffer the new line left and right so not buffering at the line points as you mentioned? I have not been in this world a long time but it is a lot of fun, sadly, for every bit of fun it seems to create 2 bits of stress; perhaps I am doing it wrong. – Phish Mar 23 '21 at 10:04

1 Answers1

0

If you skip the part where it is highly ambiguous to determine a positional geometric predicate on anything but an infinite line, that check itself is nothing but a 2D cross product. The issue here is that it has to be applied to geometric primitives, i.e. an actual two-vertex vector representation, of which a LineString may have plenty (each segment between two consecutive vertices). In effect, for each closest Point candidate you'd need to find the closest segment on the given LineString, then check for the positional predicate between their vector representations.

This could potentially be implemented in the core math of, say, ST_DWithin, and maybe even the index lookup of the (K)NN operator <->, since somewhere in their algorithms they have to include these per-segment operations already.

But it isn't, it might never be due to ambiguity and lack of usefulness, and in SQL this is bound to be somewhat dirty and inefficient...


...but not entirely useless - and rather surprisingly performant; tailor yourself some handy functions, and use a classic (K)NN query

  1. Create a set of functions to loop over a given LineString vertices, find the closest two-vertex segment to a given point, and check if it is left or right of it:
    -- Utility function; returns the cross product of three Points
    -- if cp < 0 the point is to the right, cp > 0 to the left, cp = 0 on the line
    CREATE OR REPLACE FUNCTION _ST_DeterminantSign (
      IN  pt GEOMETRY(POINT),
      IN  v1 GEOMETRY(POINT),
      IN  v2 GEOMETRY(POINT),
      OUT FLOAT
    ) LANGUAGE SQL AS
      $BODY$
        SELECT ( ST_X($3)-ST_X($2) ) * ( ST_Y($1)-ST_Y($2) ) - ( ST_Y($3)-ST_Y($2) ) * ( ST_X($1)-ST_X($2) );
      $BODY$
    ;
    

    -- loop over line segments of a given LineString, find its closest segment to a given point and check if it is to the right CREATE OR REPLACE FUNCTION ST_IsRightOf ( IN pt GEOMETRY(POINT), IN ln GEOMETRY(LINESTRING), OUT BOOLEAN ) LANGUAGE SQL AS $BODY$ SELECT _ST_DeterminantSign($1, ST_PointN($2, n), ST_PointN($2, n+1)) < 0 FROM GENERATE_SERIES(1, ST_NPoints($2)-1) AS n ORDER BY ST_MakeLine(ST_PointN($2, n), ST_PointN($2, n+1)) <-> $1 LIMIT 1; $BODY$ ;

    -- loop over line segments of a given LineString, find its closest segment to a given point and check if it is to the left CREATE OR REPLACE FUNCTION ST_IsLeftOf ( IN pt GEOMETRY(POINT), IN ln GEOMETRY(LINESTRING), OUT BOOLEAN ) LANGUAGE SQL AS $BODY$ SELECT _ST_DeterminantSign($1, ST_PointN($2, n), ST_PointN($2, n+1)) > 0 FROM GENERATE_SERIES(1, ST_NPoints($2)-1) AS n ORDER BY ST_MakeLine(ST_PointN($2, n), ST_PointN($2, n+1)) <-> $1 LIMIT 1; $BODY$ ;

    These functions here only work on LineString and Point geometries!
  2. Run
    SELECT t.id, t.geom
    FROM   <table1> AS t1
    CROSS JOIN LATERAL (
      SELECT t2.id, t2.geom
      FROM   <table2> AS t2
      WHERE  ST_IsRightOf(ST_EndPoint(ST_ShortestLine(t1.geom, t2.geom)), t1.geom)
      ORDER BY
             t1.geom <-> t2.geom
      LIMIT  1
    ) AS t
    ;
    
    to get the t2.id & t2.geom that is the closest to the right of each line in t1.

The core operation is fully index driven (ordering by distance via <->) and the filter applied only on the NN index subset returned by the <-> operation; the function execution time itself increases linearly with the amount of vertices in the given LineString.


Related:

geozelot
  • 30,050
  • 4
  • 32
  • 56
  • 1
    @Phish did you try this? While this has not attracted much positive resonance - or rather none at all you ungrateful lot ,) - I am extremely pleased with the performance of the directionally restricted KNN search using the functions - I am going to wrap this into a custom operator and put it into a GitHub Gist this weekend. – geozelot Apr 07 '21 at 08:51
  • Make a long thin triangle from your origin centred on a radial. Find all shapes that interect the triangle. Take the closest one to the origin. Modify the angle of the triangle as a tolerance measure. – wingnut May 09 '21 at 04:09