1

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:

|improve this question
  • If you can implement a little script in Python, The Shapely-Module has a function called signed_area: 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 9 '18 at 9: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. 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 9 '18 at 12:38
1

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.

|improve this answer
  • 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 9 '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
1

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).

|improve this answer
  • 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 7: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 7:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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