Slow contour plots

[Luis García Carreras]

Hi,

I'm not sure if this a matplotlib or iris issue (or just unavoidable), but haven't found anything online addressing this issue. I'm doing global contour plots for a dataset of 1537x2048 points, and it's incredibly slow. I haven't timed it, but it's sometimes over 30 minutes for a single plot. I have been plotting contours from different runs, and then difference plots (I just subtract the cubes, and plot). For some strange reason the difference plots are much slower to plot.

Is this normal for a dataset this size? I'm not sure how to time code to see exactly where the bottleneck is, but it gets stuck at the simple iris.plot.contourf command. Interactive plotting is turned off. This is all run on a remote server (Jasmin), although I don't think that would make a difference.

Thanks!

[Andrew Dawson]

1537x2048 is a lot of data points. Try using pcolormesh instead of contour it should be much faster since it doesn't have to do any contouring calculations. If you want contours I'd suggest down-sampling the data first. It is unlikely anyone will be able to spot pixel level differences in the raw data unless the plot is going to be huge anyway!

[Luis García Carreras]

Thanks a lot for the tip! It's a lot of points, but I was surprised by how incredibly slow it was. The pcolormesh tip worked a treat though, once I figured out how to get it to work like contourf. Takes a minute or so now, not fast but acceptable.

Is there an inbuilt iris routine for downscaling data? Or is it just a question of indexing the cube at regular intervals? As you say, the high resolution is unnecessary in global plots really, unless you're planning on covering a living room wall with it, and it's not exciting enough for that.

Cheers!

[Andrew Dawson]

  Takes a minute or so now, not fast but acceptable.

 

I'm quite surprised it takes this long. I tried this short test to see how fast matplotlib alone can plot a data set of this size:

 

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-10*np.pi, 10*np.pi, 2048)
y = np.linspace(-10*np.pi, 10*np.pi, 1537)
X, Y = np.meshgrid(x, y)
Z = np.cos(X * Y) * 0.25 * np.sin(X + Y)
plt.contourf(x, y, Z)
plt.contour(x, y, Z)
plt.show()

 

Is that example also slow on your environment? It takes about 10 seconds on mine. Of course it may depend on the complexity of the data too, would you mind doing a quick test to see if it is just matplotlib being slow or if iris is making it worse? Could you make the contour plot using matplotlib functions and just the raw data from your cube, something along the lines of:

 

x = cube.coord('latitude').points
y = cube.coord('longitude').points
z = cube.data
plt.contourf(x, y, z)

 

substituting in the correct names for your data set. Is that faster or the same?

[Niall Robinson]

Further to that, its worth separating out the data loading. Iris defers loading data until it has to. It might be that the plots are talking a long time because this is actually the first time that iris has had to load the data. If you do a

>>> your_cube.data

Before you plot, you will the know the data is loaded.

Niall

[Andrew Dawson]

I was going to suggest that but it still shouldn't take that long, and it doesn't really explain the vast difference between pcolormesh and contourf. It could still be a factor though I suppose.

[Niall Robinson]

> the vast difference between pcolormesh and contourf.

Yeh - I suspect you are right in that that is just from matplotlib having to calculate contours from a lot of data.

[Luis García Carreras]

Hi,

thanks for all the feedback, and sorry about my slowness. Really need to set up email alerts.

Niall, the loading data issue did occur to me, but it's not a contributing factor here, as the data is already post-processed (actually time averages of daily data).

Andrew, I've done as you suggested, and tested plotting my data both with matplotlib and iris, with both contourf and pcolormesh. I used plt.close() between calls. I tried two different arrays (absolute values, and a difference plot), as I'd noticed the latter was taking significantly longer.

Matplotlib:

plt.contourf(x, y, z): 10-15s

plt.pcolormesh(x, y, z): 10s

Iris:

iplt.contourf(cube): 4 minutes - >10 minutes

iplt.pcolormesh(cube): 30-35s

So, iplt.pcolormesh is not as bad as I stated earlier (I hadn't timed it, and tidying the plot up may add some extra time), but still a significant overhead in iris compared to matplotlib - although it does plot it on a map, so probably fair enough.

For iplt.contourf on the other hand, the difference does seem very large. Also there's a big difference between the two datasets - the difference plot I gave up on when it clocked 10 minutes.

I did also try your simpler matplotlib test. It takes me 20s as opposed to 10s for you, but calling plt.clf() then replotting it cuts the time to 10s. I thought it was a bit strange half the time was just setting up the figure.. not sure if it makes a difference to the above results.

Cheers,

Luis

[Luis García Carreras]

One more thing. I just used iplt.contourf by mistake, and had to halt the code. The traceback was exactly the same as last time, so I'm guessing this is where the code gets stuck. Might give an indication of the source of the slowdown? Full traceback below in case it's useful (sorry it's a bit long!)

/home/users/lgarciacarreras/python_scripts/gl_general/irisutils.pyc in mapcontour(varplot, levels, title, cbar, cmap, fout, verbose, extend, cbtitle, min0, boxplot)

    208         pl = iplt.pcolormesh(varplot, cmap=colmap, norm=norm)

    209     else:

--> 210         pl = iplt.contourf(varplot, levels=levels, cmap=cmap, extend=extend)

    211 

    212     #title

/usr/lib/python2.7/site-packages/iris/plot.pyc in contourf(cube, *args, **kwargs)

    575     coords = kwargs.get('coords')

    576     kwargs.setdefault('antialiased', True)

--> 577     result = _draw_2d_from_points('contourf', None, cube, *args, **kwargs)

    578 

    579     # Matplotlib produces visible seams between anti-aliased polygons.

/usr/lib/python2.7/site-packages/iris/plot.pyc in _draw_2d_from_points(draw_method_name, arg_func, cube, *args, **kwargs)

    281         result = _map_common(draw_method_name, arg_func,

    282                              iris.coords.POINT_MODE, cube, plot_defn,

--> 283                              *args, **kwargs)

    284     else:

    285         # Obtain data array.

/usr/lib/python2.7/site-packages/iris/plot.pyc in _map_common(draw_method_name, arg_func, mode, cube, plot_defn, *args, **kwargs)

    534 

    535     # Draw the contour lines/filled contours.

--> 536     return draw_method(*new_args, **kwargs)

    537 

    538 

/usr/lib/python2.7/site-packages/cartopy/mpl/geoaxes.pyc in contourf(self, *args, **kwargs)

    968                         sub_trans.force_path_ccw = True

    969 

--> 970         result = matplotlib.axes.Axes.contourf(self, *args, **kwargs)

    971         self.autoscale_view()

    972         return result

/usr/lib/python2.7/site-packages/matplotlib/axes.pyc in contourf(self, *args, **kwargs)

   7724         if not self._hold: self.cla()

   7725         kwargs['filled'] = True

-> 7726         return mcontour.QuadContourSet(self, *args, **kwargs)

   7727     contourf.__doc__ = mcontour.QuadContourSet.contour_doc

   7728 

/usr/lib/python2.7/site-packages/matplotlib/contour.pyc in __init__(self, ax, *args, **kwargs)

   1292         are described in QuadContourSet.contour_doc.

   1293         """

-> 1294         ContourSet.__init__(self, ax, *args, **kwargs)

   1295 

   1296     def _process_args(self, *args, **kwargs):

/usr/lib/python2.7/site-packages/matplotlib/contour.pyc in __init__(self, ax, *args, **kwargs)

    828                                      transform=self.get_transform(),

    829                                      zorder=zorder)

--> 830                 self.ax.add_collection(col)

    831                 self.collections.append(col)

    832         else:

/usr/lib/python2.7/site-packages/matplotlib/axes.pyc in add_collection(self, collection, autolim)

   1489         if autolim:

   1490             if collection._paths and len(collection._paths):

-> 1491                 self.update_datalim(collection.get_datalim(self.transData))

   1492 

   1493         collection._remove_method = lambda h: self.collections.remove(h)

/usr/lib/python2.7/site-packages/matplotlib/collections.pyc in get_datalim(self, transData)

    172 

    173         if not transform.is_affine:

--> 174             paths = [transform.transform_path_non_affine(p) for p in paths]

    175             transform = transform.get_affine()

    176         if not transOffset.is_affine:

/usr/lib/python2.7/site-packages/matplotlib/transforms.pyc in transform_path_non_affine(self, path)

   2190             return path

   2191         elif not self._a.is_affine and self._b.is_affine:

-> 2192             return self._a.transform_path_non_affine(path)

   2193         else:

   2194             return self._b.transform_path_non_affine(

/usr/lib/python2.7/site-packages/cartopy/mpl/geoaxes.pyc in transform_path_non_affine(self, src_path)

    163         # in matplotlib, but is in Cartopy).

    164         geoms = cpatch.path_to_geos(src_path,

--> 165                                     getattr(self, 'force_path_ccw', False))

    166 

    167         for geom in geoms:

/usr/lib/python2.7/site-packages/cartopy/mpl/patch.pyc in path_to_geos(path, force_ccw)

    201                 isinstance(collection[-1][0], Polygon) and

    202                 isinstance(geom, Polygon) and

--> 203                 collection[-1][0].contains(geom.exterior)):

    204             collection[-1][1].append(geom.exterior)

    205         else:

/usr/lib/python2.7/site-packages/shapely/geometry/base.pyc in contains(self, other)

    384     def contains(self, other):

    385         """Returns True if the geometry contains the other, else False"""

--> 386         return bool(self.impl['contains'](self, other))

    387 

    388     def crosses(self, other):

/usr/lib/python2.7/site-packages/shapely/predicates.pyc in __call__(self, this, other, *args)

      9         self._validate(this)

     10         self._validate(other, stop_prepared=True)

---> 11         return self.fn(this._geom, other._geom, *args)

     12 

     13 class RelateOp(Delegating):

/usr/lib/python2.7/site-packages/shapely/geos.pyc in errcheck_predicate(result, func, argtuple)

    183     return retval

    184 

--> 185 def errcheck_predicate(result, func, argtuple):

    186     if result == 2:

    187         raise PredicateError("Failed to evaluate %s" % repr(func))

KeyboardInterrupt: 

On Wednesday, 23 April 2014 11:25:04 UTC+2, Luis García Carreras wrote: