Determine the Winding of a Polygon in BigQuery to calculate Interior Angles
In this article I will show how to quickly determine the winding of a polygon. This can be used f.e. to calculate all interior angles for OpenStreetMap (OSM) polygons which- at the point of writing- arrange nodes both clockwise and counterclockwise. Finally, I will use the winding and calculated angles to find 'notches': concave indentations of polygon. I encountered the concept of notches in a paper on polygon complexity “ Measuring the Complexity of Polygonal Objects”.
A notch is characterized by an interior angle larger than 180º. In BigQuery we can easily calculate the angle between two vertices with the ST_ANGLE() function, but depending on if the vertices are ordered in a counterclockwise or clockwise direction, the function will either calculate the interior angle, or 360º- the interior angle. If you were to calculate the angles for a single or even a handful of polygons, visual inspection is of course enough to determine if the angles were calculated as interior angles of the polygon, but we want to do this at scale. ’12 million polygons at once’ kind of scale.
I am starting with a table 'ways' which already extracted the OSM ways I am interested in based on a few key-value combinations. We extract both the latitude and longitude columns as well as creating a geopoint from them.
get_geom as (
SELECT
ways.id,
ARRAY_LENGTH(collapse_id.nodes) as number_points,
node_data.longitude as lon,
node_data.latitude as lat,
ST_GEOGPOINT(node_data.longitude,node_data.latitude) as geog_point,
pos
FROM
ways
, ways.nodes AS node_id WITH OFFSET pos
LEFT JOIN `bigquery-public-data.geo_openstreetmap.planet_nodes` node_data
ON node_id.id = node_data.id
),The nodes making up the polygon are assigned to a position pos we got from the WITH OFFSET when unnesting the array. But we have to rearrange the data a bit so that we can calculate all angles. Starting at point 0, we are calculating the angle a at point 1, made up by the vertices 0–1 and 1–2. Repeating this procedure around the clock for each angle will bring us back to point 0. In order to calculate the angle at 0, we need the vertex 0–1 again, so we need to append both point 0 and point 1.
We calculated the number of points with ARRAY_LENGTH() in the previous step and use the result to assign positions to the repeated points (pos is 0-indexed) appended to our initial set of points.
append_point as (
SELECT
id,
number_points,
geo_point,
lon,
lat,
pos
FROM get_geom
UNION ALL
SELECT
id,
number_points,
geo_point,
lon,
lat,
number_points as pos
FROM get_geom
WHERE pos = 0
UNION ALL
SELECT
id,
number_points,
geo_point,
lon,
lat,
number_points +1 as pos
FROM get_geom
WHERE pos = 1
),As a next step, we use navigation functions to get the next (lead1) and the point following the next(lead2). This will allow us to calculate the angles with ST_ANGLE(geopoint, lead1_geogpoint, lead2_geogpoint) later on. We also pull the latitude and longitude of the next row.
get_points as(
SELECT
id,
geopoint,
number_points,
lon,
lat,
pos,
LEAD(geopoint,1) OVER(PARTITION BY id ORDER BY pos ASC) as lead1_geog_point,
LEAD(geopoint,2) OVER(PARTITION BY id ORDER BY pos ASC) as lead2_geog_point,
LEAD(lon,1) OVER(PARTITION BY id ORDER BY pos ASC) as next_lon,
LEAD(lat,1) OVER(PARTITION BY id ORDER BY pos ASC) as
next_lat
FROM append_point
)To determine if the nodes in a particular OSM polygon are arranged clockwise or counterclockwise, we use the concepts greatly explained here. It is the most intuitive explanation I found to why the formula below works.(Remember that our points contain point 0 and 1 twice and each row already includes the lead latitude and longitude. With the WHERE statement we stop before we come around to point 0 again).
get_winding as (
SELECT
id,
SUM((next_lon - lon)*(next_lat + lat)) as sum_vertices
FROM get_points
WHERE pos < number_points is not Null
GROUP BY id
)Now we can calculate all angles within the polygon and join the result of the winding formula.
get_angle as (
SELECT
id,
ARRAY_AGG(ST_ANGLE(geog_point, lead1_geog_point, lead2_geog_point) IGNORE NULLS) as angles
FROM get_points
WHERE pos < number_points
GROUP BY id
),
combine_angles_winding as(
SELECT
id,
angles,
IF(sum_vertices <0 , 'COUNTERCLOCKWISE', 'CLOCKWISE') as winding
FROM get_angle
LEFT JOIN get_winding
USING(id)
)With the information if the points in a polygon are arranged clockwise or counterclockwise, we can count the number of concave notches. Bigquery returns angles in radians, so we compare an angle to the arccosine of -1 to find all interior angles of more than 180º since:
ACOS(-1) = π rad = 180º
In a counterclockwise arrangement, an angle of more than 180º is a concave indentation.
SELECT
id,
COUNT(IF(winding = 'COUNTERCLOCKWISE', angle > ACOS(-1),angle < ACOS(-1))) as notches
FROM combine_angles_winding, UNNEST(angles) angle
GROUP BY
idIf counting notches is none of your interests, but you just want to display all interior angles correctly, you can use the IF logic to transform the angles for clockwise ordered polygons with
interior angle = ACOS(-1)*2 -angle
The nice thing about spatial data is that we can check our results plotted on a map. The picture below shows the desired result, we were able to find the 3 notches correctly (the title picture shows 759 concave notches of a clockwise arranged polygon).
Final comment:
In the context of polygon complexity I would also like to point to Jordan Ellenberg’s book “Shape: The Hidden Geometry of Information, Biology, Strategy, Democracy, and Everything Else”. In his last chapter, Ellenberg discusses gerrymandering and the attempt to determine if an electoral district has been tampered with excessively, based on geometric features of the districts shape. While he ultimately dismisses metrics like the ratio of the perimeter to the area of a polygon in his context, the chapter is certainly an interesting real-life example where understanding a polygon's complexity (or in his case randomness ) is of great interest.
