Plotting on map projection much slower on v0.6.1 than 0.6.0

[fmaussion]

The following code snippet produces an average of ERA-Interim temperature on a map:

import matplotlib.pyplot as plt 
import xray
import cartopy.crs as ccrs 
import time

netcdf = xray.open_dataset('ERA-Int-Monthly-2mTemp.nc')
t2_avg = netcdf.t2m.mean(dim='time')

start_time = time.time()
ax = plt.axes(projection=ccrs.Robinson())
if xray.__version__ == '0.6.0':
    t2_avg.plot(ax=ax, origin='upper', aspect='equal', transform=ccrs.PlateCarree()) 
else:
    t2_avg.plot(ax=ax, transform=ccrs.PlateCarree()) 
ax.coastlines()
plt.savefig('t_xray.png')
print("xray V{}: {:.2f} s".format(xray.__version__, time.time() - start_time))

I've been very careful to check that my environments are exact same (mpl 1.4.3, numpy 1.10.1, cartopy 0.13.0).

See the output for V0.6.0 and 0.6.1 (output from the latest master is similar to 0.6.1):

0.6.0:

python test_xray.py
/home/mowglie/.bin/conda/envs/climate/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):
xray V0.6.0: 3.21 s

0.6.1:

 python test_xray.py 
/home/mowglie/.bin/conda/envs/upclim/lib/python3.4/site-packages/numpy/lib/shape_base.py:431: FutureWarning: in the future np.array_split will retain the shape of arrays with a zero size, instead of replacing them by `array([])`, which always has a shape of (0,).
  FutureWarning)
/home/mowglie/.bin/conda/envs/upclim/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):
xray V0.6.1: 28.52 s

The first warning seems related to recent numpy. Note that a second warning appeared with xray V0.6.1.

It's interesting to mention that the bottleneck clearly is in the rendering (plt.savefig('t_xray.png')). Removing this line will make xray V0.6.1 faster than xray V0.6.0.

[shoyer]

We changed the default plot type from "imshow" in xray v0.6 to "pcolormesh" in xray v0.6.1. Does setting that manually make a difference? e.g., t2_avg.plot.imshow?

[fmaussion]

Thanks for the quick answer. changing to imshow indeed reduced the execution time back to normal. I wouldn't expect a factor 10 between imshow and pcolormesh, though...

[fmaussion]

If I remove the projection stuff imshow and pcolormesh are equally fast. Both functions must handle the projection differently, somehow?

[shoyer]

Yes, the difference in speed is definitely due to cartopy, which handles the projections. @pelson may be able to clarify whether imshow is always slower than pcolormesh when using cartopy, or if that is specific to particular projections.

We changed the default from imshow to pcolormesh in #608. I would not be opposed to changing the default back, given that we already product essentially equivalent plots with both methods. We did add a performance note to our (currently broken) docs: http://xray.readthedocs.org/en/latest/plotting.html#two-dimensions

@jhamman and/or @clarkfitzg may have some opinions here.

[fmaussion]

I guess it's OK to leave it as is, as long as the performance loss is documented. I like the automatic plot() function because it hides the boring mpl internals to my students and (until now) always did what I expected it to do.

30 seconds to generate a plot in a Notebook is definitely a fun killer. Maybe I should still upgrade to 0.6.1 and explain to my students that they should use imshow() instead of plot().

[pelson]

There is definitely scope for being smarter with cartopy's pcolormesh. There isn't an issue for it yet, but would be happy if you opened up a performance related issue in cartopy. pcolormesh will always be slower than imshow, but in most cases, not an order of magnitude slower!

[fmaussion]

Hi again,

I've made a self-contained example below. I have no clue about how cartopy works (blind xray user, sorry :-() and was not able to remove the dependency to xray without getting different plots for imshow() and pcolormesh()... If one of the xray gurus could help me to remove the xray part of the code I could open an issue in cartopy. I wonder if the problem comes from the fact that ERA-Interim lons are spanning 0-360?

import matplotlib.pyplot as plt 
import xray
import numpy as np
import cartopy.crs as ccrs 
import time

nlats, nlons = (241, 480)
lats = np.linspace(90, -90, nlats)
lons = np.linspace(0, 360-0.75, nlons)
l1, l2 = np.meshgrid(lons, lats)
data = xray.DataArray(l1 + l2, [('latitude',  lats), ('longitude', lons)])

start_time = time.time()
fig = plt.figure()
ax = plt.axes(projection=ccrs.Robinson())
data.plot.imshow(ax=ax, transform=ccrs.PlateCarree()) 
ax.coastlines()
plt.savefig('imshow.png')
print("imshow: {:.2f} s".format(time.time() - start_time))

start_time = time.time()
fig = plt.figure()
ax = plt.axes(projection=ccrs.Robinson())
data.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree()) 
ax.coastlines()
plt.savefig('pcolormesh.png')
print("pcolormesh: {:.2f} s".format(time.time() - start_time))

imshow: 3.09 s
pcolormesh: 27.53 s

[pelson]

Thanks @fmaussion - I've raised it in SciTools/cartopy#700. Thanks for putting together the self contained example - it's fine to have xray as a dependency on that. FWIW you would do something like plt.pcolormesh(lons, lats, l1 + l2, transform=ccrs.PlateCarree()) if you didn't have a DataArray.

[fmaussion]

Thanks @pelson ! So now this is when it becomes funny: I've been able to make three similar plots using xray's imshow, pcolormesh and cartopy's pcolormesh. Xray's pcolormesh takes twice the time as cartopy's:

import matplotlib.pyplot as plt 
from matplotlib.colors import Normalize
import xray
import numpy as np
import cartopy.crs as ccrs 
import time

nlats, nlons = (241, 480)
lats = np.linspace(90, -90, nlats)
lons = np.linspace(0, 360-0.75, nlons)
l1, l2 = np.meshgrid(lons, lats)
data = xray.DataArray(l1 + l2, [('latitude',  lats), ('longitude', lons)])

cmap = plt.get_cmap('viridis')
norm = Normalize(vmin=0, vmax=data.max().values)

start_time = time.time()
fig = plt.figure()
ax = plt.axes(projection=ccrs.Robinson())
data.plot.imshow(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap=cmap, vmin=0) 
ax.coastlines()
plt.savefig('imshow_xray.png')
print("imshow xray: {:.2f} s".format(time.time() - start_time))

start_time = time.time()
fig = plt.figure()
ax = plt.axes(projection=ccrs.Robinson())
data.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap=cmap, vmin=0) 
ax.coastlines()
plt.savefig('pcolormesh_xray.png')
print("pcolormesh xray: {:.2f} s".format(time.time() - start_time))

start_time = time.time()
fig = plt.figure()
ax = plt.axes(projection=ccrs.Robinson())
ax.pcolormesh(lons, lats, l1 + l2, transform=ccrs.PlateCarree(), cmap=cmap, norm=norm)
ax.coastlines()
plt.savefig('pcolormesh_cartopy.png')
print("pcolormesh cartopy: {:.2f} s".format(time.time() - start_time))

imshow xray: 3.06 s
pcolormesh xray: 27.50 s
pcolormesh cartopy: 12.14 s

[clarkfitzg]

This is surprising! But good to know.

[fmaussion]

Was someone able to reproduce this or is it just me?

[shoyer]

Yes, I'm seeing the same thing. Very weird -- I'll see if I can profile it.

On Wed, Nov 18, 2015 at 8:28 AM, Fabien Maussion notifications@github.com
wrote:

Was someone able to reproduce this or is it just me?

Reply to this email directly or view it on GitHub:
#657 (comment)

[jhamman]

I can't verify right now but it may have something to do with using masked arrays under the hood. There are no nan's in your example but xray still is converting the array to a masked_array before plotting. I bet plotting with pcolormesh is slower with masked arrays than with numpy arrays.

[fmaussion]

Changing the input to np.ma.asarray(l1 + l2) did not change anything as far as I can see...