iris.plot does not always plot cloud base/top extent correctly

[mo-avdjturp]

🐛 Bug Report

iris.plot does not always plot cloud base/top extent correctly. The plotted region is sometimes less than it should be.

How To Reproduce

Steps to reproduce the behaviour:

1. cube=iris.load(cloud_base_file)
2. contour_levels = np.linspace(0.0, 8000.0, 9.)
3. cf = iplt.contourf(cube[0], contour_levels)
4. plt.show()

(Example cloud base file is in /scratch/fpdh/CB/20210819T0600Z/20210819T0000Z_T+06_CB_base_test1_as_height.pp;
example cloud top file is in /scratch/fpdh/CB/20210819T0600Z/20210819T0000Z_T+06_CB_top_test1_as_height.pp)

Expected behaviour

Cloud base field extent should be the same as cloud top field extent. (Independent testing and plotting with matplotlib.pyplot confirms that this should be the case).

Actual behaviour
Cloud base field extent is less than cloud top field extent.
When plotted with matplotlib.pyplot, the extents are the same.
(See example plots; plt=plotted with matplotlib.pyplot, iplt=plotted with iris.plot)

Screenshots

Environment

• OS & Version: [e.g., Ubuntu 20.04 LTS]
• Iris Version: [e.g., From the command line run python -c "import iris; print(iris.__version__)"]
  RHEL7
  python-3.4.10
  iris-2.4.0 (from scientific software stack - default-current environment)

Additional context

Click to expand this section...

Please add additional verbose information in this section e.g., code, output, tracebacks, screenshots etc


[trexfeathers]

Thanks so much for taking the time @mo-avdjturp!

[mo-avdjturp]

Issue can be avoided by using matplotlib.pyplot to plot the field, e.g.:
lon_grid, lat_grid = np.meshgrid(cube.coord('longitude').points, cube.coord('latitude').points)
cf = plt.contourf(lon_grid, lat_grid, cube.data, levels= contour_levels)

[wjbenfold]

From a bit of investigation, I think that this is being caused mostly (if not completely) by the antialiasing in iris.plot.contourf. I'm not too surprised that this issue is present in iris-2.4.0 given significant improvements were made in iris-3.1 (see #4127, #4150, #4086 and #4263).

However, I'm able to reproduce this in iris 3.2.dev0, suggesting there's an edge case that this data is hitting. Note that the matplotlib and "fixed" iris plots are not identical - there are ~1 pixel differences in various places.

matplotlib.pyplot.contourf

iris.plot.contourf

iris.plot.contourf (with undocumented kwarg antialiased=False)

@rcomer, as the person who last fixed some antialiasing, any advice before I look further as to what might be the culprit?

[rcomer]

Thanks @wjbenfold for passing me the script. I think part of the problem is that we have a dodgy pp-file. Looking at the data itself, there are a heckofalot of grid cells with value -1.e-30. Since this is (I think) cloud fraction, negative values are impossible, so presumably -1.e-30 is supposed to be a missing data indicator. However, if we look at the pp-header by

fields = iris.fileformats.pp.load(filename)
field = next(fields)
print(field)

we see that bmdi is set to -1e+30

I tried making a new cube with

new_cube = cube.copy(ma.masked_where(cube.data == -1.e-30, cube.data))

Plotting this with iplt.contourf gave me something more sensible:

Though it still doesn't look the same as the matplotlib version

I agree that the antialiasing seems to be the culprit, but I'm not currently seeing anything that relates to the previous fixes we put in.

[wjbenfold]

@mo-avdjturp, based on my (shallow) understanding of the initial issue, if anything @rcomer's iplt version now looks better than the plt one - should we keep digging into the differences or does this seem resolved to you?

[rcomer]

I'm also wondering whether a contour plot is the best choice for this data set. Taking the version plotted with matplotlib, if we zoom in on the Maritime Continent, we can't see much in the way of contours and it's mostly just clusters of the same colour. Other regions are similar. @mo-avdjturp have you considered using pcolormesh instead for these plots? I think contours work best when the you have a smoothly varying field.

[mo-avdjturp]

Thanks for your efforts and insights, @rcomer and @wjbenfold . Yes, the 1.0E-30 value is intended to be a missing data indicator (perhaps I've got this incorrect and it should be -1.0E+30; the file was generated by some other python code I'm writing). I hadn't thought of plotting using pcolormesh and agree that would be a better tool to use. I'll try looking at the masking the field too. I'm not likely to get time to look at this more for a week or so though.

[github-actions]

In order to maintain a backlog of relevant issues, we automatically label them as stale after 500 days of inactivity.

If this issue is still important to you, then please comment on this issue and the stale label will be removed.

Otherwise this issue will be automatically closed in 28 days time.

[github-actions]

This stale issue has been automatically closed due to a lack of community activity.

If you still care about this issue, then please either:

• Re-open this issue, if you have sufficient permissions, or
• Add a comment stating that this is still relevant and someone will re-open it on your behalf.