My first QGIS map: the aid from Cyrene

Polar (circular) diagram

Orientation of geographic features can be an important piece of information, for instance when studying ancient field systems or analysing geological formations. The problem comes when we try to summarise such data graphically . Standard histograms are not ideal because they cut through the continuous distribution - one cannot easily understand the relationship between the beginning and the end of the diagram (which should stick together). Polar diagram is what we need, but it cannot be made without some fiddling around.


There already exists a plugin called Line direction histogram that does the job quite fine. However, QGIS comes with great matplotlib python library for data visualisation (windows version at least) which can be tweaked in all imaginable ways (colours, offsets etc.). The following script is simply passing data to matplotlib, directions have to be pre-calculated beforehand.

##Visualisation=group
##input_layer=vector
##azimuth_field=field input_layer
##bin_count=number 32
##hollow_radius=number 0

"""
This script is used for generating circular graphs 
using matplotlib

More details on: www.qgiswhisperer.com/2017/11/polar-diagram.html
"""

import numpy as np
import matplotlib.pyplot as plt

input = processing.getObject(input_layer)
provider = input.dataProvider()

fieldIdx = provider.fields().indexFromName(azimuth_field)
    
features = input.getFeatures()

out=[]

for f in features: 
    v = f[fieldIdx]
    if 0<=v<=360: out.append(f[fieldIdx])  

# the following code is mostly adapted from : 
# http://stackoverflow.com/questions/22562364/circular-histogram-for-python
N = int(bin_count)
theta = np.linspace(0,  2 * np.pi, N, endpoint=False)
radii, labels = np.histogram (np.array(out),
                              bins=np.linspace(0, 360, N+1, endpoint=True)
width = 2*np.pi  / N

ax = plt.subplot(111, polar=True)
#put north to top , arrange clockwise, give geographic labels (instead of angles)
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
ax.set_xticklabels(['N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'])

bars = ax.bar(theta, radii, width=width,bottom=hollow_radius)# 

# Use custom colors and opacity
for r, bar in zip(radii, bars):
    bar.set_facecolor(plt.cm.jet(r / 10.))
    bar.set_alpha(0.8)

plt.show()

NB. QGIS may crash when saving the image in .png format - .jpeg works better (https://hub.qgis.org/issues/13982)

Comments