Finding opposite side of a polygon

Question

I have a polygon and known 2 points. One is the center of the polygon and another is outside the polygon. I'm trying to find the opposite point of the polygon where the line string of the 2 known points intersects the polygon.

Below is the SQL I'm using to find that point by taking the line string and extending it out past the polygon then getting the last point of intersection.

st_astext(st_endpoint(ST_Intersection(pg.poly_geog::geometry,
    (
    select
        st_makeline(a,
        st_translate(b,
        sin(az2) * len,
        cos(az2) * len))
    from
        (
        select
            a,b,
            ST_Azimuth(b,a) as az2,
            ST_Distance(a,b) + 0.00001 as len
        from
            (select
                st_makepoint(c1.longitude, c1.latitude)::geography::geometry as a,
                ST_Centroid(st_astext(pg.poly_geog))::geography::geometry as b )a )b))))

This seemed to work until I started running more tests and found that it only works if the linestring between the points is in a certain heading. If the linestring has a heading of right to left it works fine and I can get the first intersection point and the above calculation of the opposite intersection point. If the line string is heading in the oposite direction for another polygon then I'm getting the the point returned from the same side as the intitial intersection.

If I switch around st_endpoint to st_start point then it just changes around the direction that the function works for.

Ive also tried ST_GeometryN as well but that only works in a certain direction also.

How can I get this to function regardless of the heading on the linestring or how can I determine whether I need to use startpoint or endpoint

does it work if you swap the geometries in st_intersection? Otherwise compute the intersection, then extract both start and end points then keep the one that is the further away from the source — JGH, Oct 10 '20 at 11:24
Thank you! Of course its something straight forward like that. I spent the best part 8 hours looking at it before I posted the question — briansexton, Oct 10 '20 at 12:05
check out this gisse post. It has a function for extending lines that might be useful — jbalk, Oct 14 '20 at 01:57

Cyril Mikhalchenko · Answer 1 · 2020-11-30T06:50:08.840

I'm not sure if I understood the auto-translation of your question text correctly, but judging by the SQL code and the picture,

run the script as a CTE

(if necessary, you can convert it into a subquery form, who's more comfortable for you)

WITH
    tbla AS (SELECT ST_MakeLine(a.geom, ST_Centroid(b.geom)) geom FROM pnt a, polygon b),
    tblb AS (SELECT ST_Azimuth(ST_StartPoint(geom), ST_EndPoint(geom)) as azimuth, ST_Distance(ST_StartPoint(geom), ST_EndPoint(geom)) + 0.00001 as length FROM tbla),
    tblc AS (SELECT ST_MakeLine(a.geom, ST_Translate(a.geom, sin(azimuth)*length, cos(azimuth)*length)) geom FROM tbla a, tblb b),
    tbld AS (SELECT ST_Intersection(a.geom, ST_ExteriorRing(b.geom)) geom FROM tblc a JOIN polygon b ON ST_Intersects(a.geom, b.geom))
             SELECT (ST_Dump(geom)).geom geom FROM tbld;

Check your result,

The script is called: ST_PolygonPunchedThroughArrow...

If you are looking for only 1 opposite point, just continue the geoinstrument logic a bit further, for example, like this:

WITH
tbla AS (SELECT ST_MakeLine(a.geom, ST_Centroid(b.geom)) geom FROM pnt a, polygon b),
tblb AS (SELECT ST_Azimuth(ST_StartPoint(geom), ST_EndPoint(geom)) as azimuth, ST_Distance(ST_StartPoint(geom), ST_EndPoint(geom)) + 0.00001 as length FROM tbla),
tblc AS (SELECT ST_MakeLine(a.geom, ST_Translate(a.geom, sin(azimuth)*length, cos(azimuth)*length)) geom FROM tbla a, tblb b),
tbld AS (SELECT ST_Intersection(a.geom, ST_ExteriorRing(b.geom)) geom FROM tblc a JOIN polygon b ON ST_Intersects(a.geom, b.geom)),
tble AS (SELECT (ST_Dump(geom)).geom geom FROM tbld)
         SELECT (a.geom) geom, ST_Distance(a.geom, b.geom) dist FROM tble a, tbla b ORDER BY dist DESC LIMIT 1;

Translated with www.DeepL.com/Translator (free version)

Finding opposite side of a polygon

1 Answers1