2

This question Calculate field of view / angle of a point to a shapefile was attempting to achieve a similar outcome to my need (calculating the angle of view between observer locations and a polygon). The image from the question is shown below for ease, with the angle of interest being alpha:

Angle between observer point and a polygon

My solution is to use a virtual layer to calculate the bearing between each viewpoint and each of the nodes on the polygon layer and then use aggregate expressions to give me the information I need.

I have the following code at present but testing it throws a syntax error (1 - near ".": syntax error) that I cannot work out (there may be other errors and I'm not convinced it's the most optimal code but I would like to use a virtual layer).

SELECT 
  vp_union.*,
  minimum(vp_union.bearing) AS bearing_min,
  median(vp_union.bearing) AS bearing_med,
  maximum(vp_union.bearing) AS bearing_max,
  range(vp_union.bearing) AS bearing_range,
  distance(
    vp_union.geometry,
    closest_point(
      st_union(redline_outer.geometry),
      vp_union.geometry
    ),
    vp_union.geometry
  ) AS distance
FROM
  (SELECT
    vp.*,
    degrees(
      azimuth(
        vp.geometry,
        redline_nodes.geometry
      )
    ) AS vp.bearing
  FROM
    "Viewpoint Layer 2D" AS vp
  UNION
    (SELECT
      nodes_to_points(
        st_union(redline.geometry)
      ) AS geometry
    FROM
      "Redline" AS redline
    ) AS redline_nodes
  ) AS vp_union,
  "Redline" AS redline_outer
GROUP BY
  vp_union.name
7
  • 1
    Try changing this ) AS vp.bearing into ) AS bearing
    – Taras
    Commented Nov 26, 2021 at 14:05
  • @Taras - good spot. Now there's an issue with an opening bracket somewhere... I'm wondering if it's the way I've set up the nested SELECT query Commented Nov 26, 2021 at 14:14
  • @BERA correct - one point layer ("Viewpoint Layer 2D") and one polygon layer ("Redline"). Commented Nov 26, 2021 at 14:15
  • @BERA every point for every polygon node. Commented Nov 26, 2021 at 14:22
  • 1
    @babel I had just realised there are floors in my logic (there's no way to determine which segment between the smallest angle and biggest angle is the correct, also there are also issues where the site lies either side of 0 degrees). So I had come out at using a convex hull but instead of testing all the points, create a convex hull that includes the viewpoint and then extracting just the two line sections that share a node with the viewpoint is all that's needed. Commented Nov 27, 2021 at 9:51

1 Answer 1

1

So I was well out on my original idea. In the end I realised that there is a possible case where a viewpoint could lie within the convex hull of the polygons, so I went back to testing the individual vertices. It's long, convoluted and a bit slow (I'm sure someone could optimise).

I'm aware this doesn't quite give you alpha (but it did what I needed it too) and it's relatively trivial to calculate that angle using the start and end points of the linestring for each viewpoint.

WITH RECURSIVE
  -- generate table with all the points from the polygons layer
  polygons_points AS(
  SELECT
    1 AS fid,
    1 AS n,
    point_n(polygons.geometry,1) AS geometry
  FROM
    "Redline" AS polygons
  UNION ALL
  SELECT
    ogc_fid AS fid,
    n + 1 AS n,
    point_n(polygons.geometry,n+1)
  FROM
    "Redline" AS polygons
  INNER JOIN
    polygons_points ON polygons.ogc_fid = polygons_points.fid
  WHERE
    n < ST_NPoints(polygons.geometry)
  ),
  
  -- generate table with spokes from polygons_points to the points layer
  spokes AS(
  SELECT
    points.VP AS pid,
    polygons_points.fid,
    polygons_points.n,
    MakeLine(points.geometry,polygons_points.geometry) AS geometry
  FROM
    "Viewpoint Locations 2D" AS points
  CROSS JOIN
    polygons_points
  ),
  -- select sightlines from spokes that do not cross the original polygons
  sightlines AS(
  SELECT spokes.*
  FROM
    spokes
  INNER JOIN
    "Redline" AS polygons
  ON
    polygons.ogc_fid = spokes.fid
  AND
    ST_Crosses(spokes.geometry,polygons.geometry) <> 1
  ),
  -- select and join the polygon segments that intersect with the sightlines
  polygons_segments AS(
  SELECT
    sightlines.pid,
    sightlines.fid,
    sightlines.n,
    ST_AsText(sightlines.geometry) AS point_txt,
    ST_Union(ST_GeometryN(ST_DissolveSegments(ST_ExteriorRing(polygons.geometry)),n)) AS geometry
  FROM
    sightlines
  INNER JOIN
    "Redline" AS polygons
  ON
    polygons.ogc_fid = sightlines.fid
  AND
    ST_Intersects(sightlines.geometry, polygons.geometry)
  GROUP BY
    sightlines.pid
  )

SELECT * FROM polygons_segments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.