Skip to content
Merged
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 95 additions & 0 deletions docs/iris/example_code/graphics/anomaly_log_colouring.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
"""
Colouring anomaly data with logarithmic scaling
===============================================

In this example, we need to plot anomaly data where the values have a
"logarithmic" significance -- i.e. we want to give approximately equal ranges
of colour between data values of, say, 1 and 10 as between 10 and 100.

As the data range also contains zero, that obviously does not suit a simple
logarithmic interpretation. However, values of less than a certain absolute
magnitude may be considered "not significant", so we put these into a separate
"zero band" which is plotted in white.

To do this, we create a custom value mapping function (normalization) using
the matplotlib Norm class `matplotlib.colours.SymLogNorm
<http://matplotlib.org/api/colors_api.html#matplotlib.colors.BoundaryNorm>`_.
We use this to make a cell-filled pseudocolour plot with a colorbar.

NOTE: By "pseudocolour", we mean that each data point is drawn as a "cell"
region on the plot, coloured according to its data value.
This is provided in Iris by the functions :meth:`iris.plot.pcolor` and
:meth:`iris.plot.pcolormesh`, which call the underlying matplotlib
functions of the same names (i.e. `matplotlib.pyplot.pcolor
<http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.pcolor>`_
and `matplotlib.pyplot.pcolormesh
<http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.pcolormesh>`_).
See also: http://en.wikipedia.org/wiki/False_color#Pseudocolor.

"""
import iris
import iris.plot as iplt
import matplotlib.pyplot as plt
import matplotlib.colors as mcols
import matplotlib.ticker as mticks
import numpy as np
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unused import



def main():
# Load a sample air temperatures sequence.
file_path = iris.sample_data_path('E1_north_america.nc')
temperatures = iris.load_cube(file_path)

# Create a sample anomaly field for one year, by subtracting a time mean.
i_year = 122
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not keen on encouraging indexing for this type of operation. It can make brittle code that fails silently or unexpectedly. How about using coord_categorisation here to add a year coord? e.g.

import iris.coord_categorisation

iris.coord_categorisation.add_year(temperatures, 'time')
time_mean = temperatures.collapsed('time', iris.analysis.MEAN)
anomaly = temperatures - time_mean

year = 1982
specific_year_anomoly = anomolies.extract(iris.Constraint(year=year))

It's a simple opportunity to demonstrate and encourage this approach. You can then use year and time_mean.coord('year').bounds[0] to get the values for the title.

time_mean = temperatures.collapsed('time', iris.analysis.MEAN)
anomaly = temperatures[i_year] - time_mean

# Construct a plot title string explaining which years are involved.
times = temperatures.coord('time')
cube_years = [time.year for time in times.units.num2date(times.points)]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not needed if you use categorisation.

title = 'Temperature anomaly [{}, log scale]'.format(anomaly.units)
title += '\n{} differences from {}-{} average.'.format(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure I like full stops in plot titles...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, I wasn't sure either. Perhaps the newline is enough.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I reckon so!

cube_years[i_year], cube_years[0], cube_years[-1])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd consider getting the limits from time_mean cube.


# Setup scaling levels for the anomaly data.
minimum_log_level = 0.1
maximum_scale_level = 3.0

# Use a standard colour map which varies blue-white-red.
# For suitable options, see the 'Diverging colormaps' section in:
# http://matplotlib.org/examples/color/colormaps_reference.html
anom_cmap = 'bwr'

# Create a 'logarithmic' data normalization.
anom_norm = mcols.SymLogNorm(linthresh=minimum_log_level,
linscale=0,
vmin=-maximum_scale_level,
vmax=maximum_scale_level)
# Setting "linthresh=minimum_log_level" makes its non-logarithmic
# data range equal to our 'zero band'.
# Setting "linscale=0" maps the whole zero band to the middle colour value
# (i.e. 0.5), which is the neutral point of a "diverging" style colormap.

# Make a pseudocolour plot using this colour scheme.
mesh = iplt.pcolormesh(anomaly, cmap=anom_cmap, norm=anom_norm)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd be tempted to plot this on a LambertConformal projection given it's North America.
figure_1

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sweeet !


# Add a colourbar, with suitable "log" ticks, and end-extensions to show
# the handling of any out-of-range values.
tick_levels = [-3, -1, -0.3, 0.1, 0.3, 1, 3]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't the centre point be at 0.0 rather than 0.1? The asymmetry looks odd.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks odd, because -0.1 and +0.1 actually appear in the same place, so we must mark one only, to avoid the strings colliding.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine, then set it to zero which is between them and also in the same place. It's also independent of the value you select for minimum_log_level.

tick_labels = [-3, -1, -0.3, r'$\pm$0.1', 0.3, 1, 3]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd construct the centre tick label from minimum_log_level to make this more robust to change.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I couldn't really see the point of doing that, unless you were to generate all the ticks automatically
As we have it, the numbers and the labels are entirely explicit, so why generalise this part of it ?

I could indeed form the contour levels (and tick values) by calculation from the min and max with matplotlib.ticker.LogLocator. But that turns out to be a bit more involved than you might think, as I discovered when I tried it. This is my neatest result so far ...

    # Add a colourbar, with suitable "log" ticks, and end-extensions to show
    # the handling of any out-of-range values.

    # Calculate decade-spaced tick values from the limit values.
    tick_levels = mticks.LogLocator(subs=[1, 3]).tick_values(
        minimum_log_level,
        maximum_scale_level)
    tick_levels = [tick for tick in tick_levels
                   if minimum_log_level <= tick <= maximum_scale_level]
    tick_levels = [-tick for tick in tick_levels[::-1]] + tick_levels

    # Replace the central pair of ticks at -min, +min with a single tick, to
    # stop them overlapping.
    i_minval_tick = len(tick_levels) / 2 - 1
    tick_levels.pop(i_minval_tick)

    # Create tick labels, adding a +/- character to the (single) minimum label.
    tick_format = "{:.1g}"
    tick_labels = [tick_format.format(tick) for tick in tick_levels]
    tick_labels[i_minval_tick] = r'$\pm$' + tick_labels[i_minval_tick]

    colorbar_tick_formatter = mticks.FixedFormatter(tick_labels)

(just in case you thought I couldn't do it 😉 !)

But all this is, I think, just far too distracting.
By contrast, using manual selection has the advantage that that the ticks and labels are literals, which means you can see - without having to explicitly explain it - how the contour levels and tick values need to be constructed slightly differently.

@esc24 am I missing something here ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely far too distracting. I was just suggesting the middle label centre_label = r'$\pm${}'.format(minimum_log_level) so that if I change minimum_log_level to 0.15 or 0.2 for example (extending the central linear region), the colorbar tick label would reflect the change and indicate the range covered by the white colour. As it stands if I do that, the centre tick label become misleading. I suggest keeping the others all simple numeric literals for simplicity.

colorbar_tick_formatter = mticks.FixedFormatter(tick_labels)
bar = plt.colorbar(mesh, orientation='horizontal', extend='both',
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unused variable bar.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you consider labelling the colorbar? Looks odd to see numbers without a label and units although I see you've put it in the title.

ticks=tick_levels,
format=colorbar_tick_formatter)

# Add coastlines and a title.
plt.gca().coastlines()
plt.title(title)

# Display the result.
plt.show()


if __name__ == '__main__':
main()
37 changes: 37 additions & 0 deletions docs/iris/example_tests/test_anomaly_log_colouring.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
# (C) British Crown Copyright 2014, Met Office
#
# This file is part of Iris.
#
# Iris is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Iris is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Iris. If not, see <http://www.gnu.org/licenses/>.


# Import Iris tests first so that some things can be initialised before
# importing anything else.
import iris.tests as tests

import extest_util

with extest_util.add_examples_to_path():
import anomaly_log_colouring


class TestAnomalyLogColouring(tests.GraphicsTest):
"""Test the anomaly colouring example code."""
def test_anomaly_log_colouring(self):
with extest_util.show_replaced_by_check_graphic(self):
anomaly_log_colouring.main()


if __name__ == '__main__':
tests.main()
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.