2

Is there a simple way to calculate a field in QGIS indicating if the vertices of each polygon feature in a layer are in clockwise or anti-clockwise direction?

For simplicity, lets assume the layer has multiple polygon features, with a single geometry per feature (no multifeatures/multipolygons), and without donuts or other complications.

I have found related questions, but none of them give an answer to my question.

These ones are ESRI specific:

This one resolves the problem from a math point of view, but not applied to QGIS

The same algorithm is implemented in a shapely function - see links provided by Andreas Muller:

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
jberrio
  • 796
  • 7
  • 23
  • If you can implement a little script in Python, The Shapely-Module has a function called signed_area: http://toblerity.org/shapely/shapely.algorithms.html. Polygons the must be read into a Python-Data-Structure and the first ring must be extracted. But anyway: all outer rings of polygons in a shapefile have the same orientation: clockwise. – Andreas Müller Oct 09 '18 at 09:24
  • Thanks for your answer. You are right, the shapefile standard indicates polygons must always go clockwise. However QGIS accepts polygons in any direction when the provider is not a shapefile. https://issues.qgis.org/issues/6283 . I am interested in a simple way calculate them in a field for those layers. The algorithm you suggest with shapely is the same one explained in the last link of my question. Yes, I may end up programming it but I am interested in checking first if there is a simpler solution... – jberrio Oct 09 '18 at 12:38

3 Answers3

3

This solution uses the Advanced Python field calculator in QGIS 2.18 described here and the function signed_area (i copied it from shapely source code)

In the dialog we choose the Input layer first and a Result field name (somename) and *Field type' (Float). In the Global expression field we insert code of a python function called signed_area:

def signed_area(pr):
    """Return the signed area enclosed by a ring using the linear time
algorithm at http://www.cgafaq.info/wiki/Polygon_Area. A value >= 0
indicates a counter-clockwise oriented ring."""
    xs, ys = map(list, zip(*pr))
    xs.append(xs[1])
    ys.append(ys[1])
    return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(pr)))/2.0

And in the field called Formula we write:

value = signed_area($geom.asPolygon()[0])

Because asPolygon returns a list like object, we make sure to access the outer ring by using index 0 on that list. The result can be saved to a temporary layer and shows up in QGIS after calculation. An disadvantage may be, the a totally new layer is created, instead of just a new field is created.

Andy
  • 1,365
  • 9
  • 16
Andreas Müller
  • 2,622
  • 13
  • 20
  • That is the kind of answer I was looking for, thanks. I made a small modification to the code so it delivers a text field with "clockwise" and "anti-clockwise" labels directly; you can see it in my post below. – jberrio Oct 09 '18 at 23:42
  • @Andreas Müller shouldn't it be xs.append(sx[0]) instead of 1? Because I suppose they're trying to close the ring by adding the first point to the list. – Daniel Marques Dec 23 '19 at 12:44
  • No, it's because the formula needs always an i+1 and at the last repetition of the list comprehension the second point is needed. Search for shoelace formula to get a visual explanation. – Andreas Müller Dec 23 '19 at 12:52
3

For anybody looking for the same answer, below is how I implemented the code based on Andreas Muller's answer. It delivers a text field with a "clockwise" or "anti-clockwise" label for each polygon.

  • In QGIS, go the the Processing Toolbox > Advanced Python field calculator
  • Select your input layer
  • Field type = String
  • Length = 30 (or choose whatever fits you)
  • Global expression:

(paste this)

def signed_area(pr2):
     """Return the signed area enclosed by a ring using the linear time
     algorithm at http://www.cgafaq.info/wiki/Polygon_Area. A value >= 0
     indicates a counter-clockwise oriented ring."""
     xs, ys = map(list, zip(*pr2))
     xs.append(xs[1])
     ys.append(ys[1])
     return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(pr2)))/2.0

def rotation_dir(pr):
     signedarea = signed_area(pr)
     if signedarea > 0:
         return "anti-clockwise"
     elif signedarea < 0:
         return "clockwise"
     else:
         return "UNKNOWN"
  • value: rotation_dir($geom.asPolygon()[0])

I verified this works well in QGIS for a vector layer (memory provider, not a shapefile) with polygons rotating in multiple directions (clockwise and anti-clockwise).

jberrio
  • 796
  • 7
  • 23
  • There is another interpretation of the return value. If signed_area returns 0, then the polygon ring is not closed. – Andreas Müller Oct 10 '18 at 07:06
  • The label in that case will be UNKNOWN. If I ever get that result it will lead me to have a look at the polygon and see what's causing the problem. – jberrio Oct 11 '18 at 07:56
3

If all you need is to determine whether a polygon has its points set CW or CCW, shapely objects have a property for that:

import shapely
polygon = shapely.geometry.Polygon([(0, 0), (1, 1), (1, 0)]) 
polygon.exterior.is_ccw # False
polygon = shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1)])
polygon.exterior.is_ccw # True
Ricardo
  • 141
  • 3