Robin's Blog

Producing polar contour plots with matplotlib

In my field I often need to plot polar contour plots, and generally plotting tools don’t make this easy. In fact, I think I could rate every single graphing/plotting package in the world by the ease of producing a polar contour plot – and most would fail entirely! Still, I have managed to find a fairly nice way of doing this using my plotting package of choice: matplotlib.

I must warn you first – a Google search for matplotlib polar contour¬†or a similar search term will produce a lot of completely out-dated answers. The most commonly found answers are those such as this StackOverflow question and this forum post. In fact, the first question was asked by me last year – and got an answer which is now completely out of date. Basically, all of these answers tell you that you can’t do a polar contour plot directly in matplotlib, and you must convert your points from polar co-ordinates to cartesian co-ordinates first. This isn’t difficult, but is a pain to do, and of course you then end up with cartesian axes which doesn’t look great. The great news is that you can now do polar contour plots directly with matplotlib!

So, how do you do them? Simple really, you just create some polar axes and plot a contour plot on them:

fig, ax = subplots(subplot_kw=dict(projection='polar'))
cax = ax.contourf(theta, r, values, nlevels)

This produces a filled contour plot, as it uses the contourf function, using the contour function would give simple contour lines. The first three parameters which must be given to this function are all two-dimensional arrays containing: the radii, the angles (theta) and the actual values to contour. The final parameter is the number of contour levels to plot – you tend to want lower numbers for line contours and higher numbers for filled contour plots (to get a smooth look).

I never quite understood these two-dimensional arrays, and why they were needed. I normally had my data in the form of three lists that were basically columns of a table, where each row of the table defined a point and value. For example:

Radius Theta Value
10 0 0.7
10 90 0.45
10 180 0.9
10 270 0.23
20 0 0.5
20 90 0.13
20 180 0.52
20 270 0.98

Each of these rows define a point – for example, the first row defines a point with a radius of 10, an angle of 0 degrees and a value of 0.7. I could never understand why the contour function didn’t just take these three lists and plot me a contour plot. In fact, I’ve written a function that will do just that, which I will describe below, but first let me explain how those values are converted to two-dimensional arrays.

First of all, lets think of the dimensions: we obviously have dimensions here in our data, the radius and the angle. In fact, we could re-shape our values array so that it is two-dimensional fairly easily. We can see from the table above that we are doing all the azimuth angles for a radius of 10 degrees, then the same azimuth angles for a radius of 20 degrees, etc. Thus, rather than our values being stored in a one-dimensional list, we could put them in a two-dimensional table where the columns are the azimuth angles, and the rows are the radii:

0 90 180 270
10 0.7 0.45 0.9 0.23
20 0.5 0.13 0.52 0.98

This is exactly the sort of two dimensional array that we need to give to the contourf function. That’s not too hard to understand – but why on earth do the radii and angle arrays have to be two-dimensional too. Well, basically we just need two arrays like the one above, but with the relevant radii and angles in the cells, rather than the values. So, for the angles, we’d have:

0 90 180 270
10 0 90 180 270
20 0 90 180 270

And for the radii we’d have:

0 90 180 270
10 10 10 10 10
20 0 20 20 20

Then, when we take all three arrays together, each cell will define the three bits of information we need. So, the top left cell gives us an angle of 0, a radius of 10 and a value of 0.7. Luckily, you don’t have to make these arrays by hand – a handy NumPy function called meshgrid will do it for you:

>>> radii = np.arange(0, 60, 10)
>>> print radii
[ 0 10 20 30 40 50]
>>> angles = np.arange(0, 360, 90)
>>> print angles
[  0  90 180 270]
>>> np.meshgrid(angles, radii)
(array([[  0,  90, 180, 270],
       [  0,  90, 180, 270],
       [  0,  90, 180, 270],
       [  0,  90, 180, 270],
       [  0,  90, 180, 270],
       [  0,  90, 180, 270]]),
array([[ 0,  0,  0,  0],
       [10, 10, 10, 10],
       [20, 20, 20, 20],
       [30, 30, 30, 30],
       [40, 40, 40, 40],
       [50, 50, 50, 50]]))

One thing to remember is that the plotting function requires the angle (theta) in radians, not degrees, so if your data is in degrees (as it often is) then you’ll need to convert it to radians using the NumPy radians function.

After doing all of this you can get your data into the contour plotting function correctly, and you can get some polar axes for it to be plotted on. However, if you do this, you’ll find that your axes look something like this:

Example of polar axes with 0 at 3pm

You can see that zero degrees isn’t at the top, it’s at the ‘East’ or ’3 o’clock’ position, and the angles go round the wrong way! Apparently that’s how these things are often done in maths – but in my field particularly people want to have a polar plot like a compass, with zero at the top!

If you try and find how to do this, you’ll find a StackOverflow answer with a brilliant subclass of PolarAxes which does this for you. It’s brilliant that matplotlib allows you do this sort of customisation, but if you look below the accepted answer you’ll find a link to the matplotlib documentation for a function called set_theta_zero_location. This function very nicely takes a compass direction (“N” or “E” or “NE” etc) for where zero should be, and puts it there! Similarly, the function set_theta_direction sets the direction in which the angles will increase. All you need to do to use these is call them from the axes object:

ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)

The example above will set up the plot for a ‘normal’ compass-style plot with zero degrees at the north, and the angles increasing clockwise. If you find that these lines of code give an error you need to update your matplotlib version – these methods were only added in the latest version (v1.1.0).

So, now we’ve covered everything that I’ve gradually learnt about doing this, we can put it all together in a function. I use the function below whenever I want to plot a polar contour plot, and it works fine for me. It is documented through the docstring shown in the code below.

I can’t guarantee the code will work for you, but hopefully this post has been helpful and you’ll now be able to go away and create polar contour plots in matplotlib.

import numpy as np
from matplotlib.pyplot import *

def plot_polar_contour(values, azimuths, zeniths):
    """Plot a polar contour plot, with 0 degrees at the North.

    Arguments:

     * `values` -- A list (or other iterable - eg. a NumPy array) of the values to plot on the
     contour plot (the `z` values)
     * `azimuths` -- A list of azimuths (in degrees)
     * `zeniths` -- A list of zeniths (that is, radii)

    The shapes of these lists are important, and are designed for a particular
    use case (but should be more generally useful). The values list should be `len(azimuths) * len(zeniths)`
    long with data for the first azimuth for all the zeniths, then the second azimuth for all the zeniths etc.

    This is designed to work nicely with data that is produced using a loop as follows:

    values = []
    for azimuth in azimuths:
      for zenith in zeniths:
        # Do something and get a result
        values.append(result)

    After that code the azimuths, zeniths and values lists will be ready to be passed into this function.

    """
    theta = np.radians(azimuths)
    zeniths = np.array(zeniths)

    values = np.array(values)
    values = values.reshape(len(azimuths), len(zeniths))

    r, theta = np.meshgrid(zeniths, np.radians(azimuths))
    fig, ax = subplots(subplot_kw=dict(projection='polar'))
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)
    autumn()
    cax = ax.contourf(theta, r, values, 30)
    autumn()
    cb = fig.colorbar(cax)
    cb.set_label("Pixel reflectance")

    return fig, ax, cax

Categorised as: Programming, Python, Remote Sensing


5 Comments

  1. Thanks, this is a good write up. I’ve been struggling with plotting spectral wave data in a polar plot in matplotlib.

    I found your post after seeing your question on stackoverflow, might be worth linking them together.

  2. Laurie Wilkinson says:

    Thanks a lot for posting this. I can’t use it yet because my work computer won’t allow me to download matplotlib v1.1.0 but I will hopefully soon.
    To James Morrison: I’m working with wave spectral data too. If you haven’t heard of it, the WAFO toolbox could be good for you – http://www.maths.lth.se/matstat/wafo/. It has a python version too.

  3. John L says:

    I tried this out but it’s still quite a bit slower than what my previous method was: for each point along a ray, add an x,y,intensity to a scatter plot. This is still roughly 4x faster than a polar contour as above, but the output ends up looking pretty similar.

  4. Nils says:

    Thanks! This was super useful.

  5. Dusan says:

    Thank you very much for sharing the useful information! I have one question. Is there any chance to bring the polar contour plot to 3D (to make cylindrical plot)? I work with computational chemistry and performed a potential energy surface scan of a small 3 atomic molecule, where I was stretching one bond at different angles. The data which I produced would form sort of Mexican hat. I can plot them 3D in Cartesian coordinates, but the plot is not as easily readable as if it was plotted in cylindrical coordinates. Thanks!

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>