contourf method significantly slower than basemap

[crazyapril]

Description

I'm migrating my codes from Basemap to Cartopy. Generally Cartopy is always faster, however I find contourf method can take me seconds. Here's my example code.

import time

import cartopy
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage as sn
from mpl_toolkits.basemap import Basemap

#-------------------
#  Create data
#-------------------
x, y = np.mgrid[60:140:0.25, 0:45:0.25]
r = np.random.randint(0, 255, size=(320, 180), dtype='uint8')
r = sn.filters.gaussian_filter(r, (3, 3), mode='constant')

#-------------------
#  Cartopy
#-------------------
ax = plt.axes(projection=cartopy.crs.LambertConformal(
    central_longitude=100.0, central_latitude=35.0))
ax.set_extent((60,140,0,45))
for step in (20, 10, 5):
    tic = time.time()
    col = ax.contourf(x, y, r, transform=cartopy.crs.PlateCarree(),
        levels=np.arange(100,161,step))
    toc = time.time()
    points = sum(len(e.get_paths()) for e in col.collections)
    print('[Cartopy] Step: {} Points: {}  Time: {:.0f}ms'.format(
        step, points, (toc - tic) * 1000))

#-------------------
#  Basemap
#-------------------
_map = Basemap(projection='lcc', resolution=None,
    lon_0=100.0, lat_0=35.0, llcrnrlat=0,
    urcrnrlat=45, llcrnrlon=60, urcrnrlon=140)
for step in (20, 10, 5):
    tic = time.time()
    col = _map.contourf(x, y, r, levels=np.arange(100,161,step))
    toc = time.time()
    points = sum(len(e.get_paths()) for e in col.collections)
    print('[Basemap] Step: {} Points: {}  Time: {:.0f}ms'.format(
        step, points, (toc - tic) * 1000))

The result:

[Cartopy] Step: 20 Points: 123  Time: 729ms
[Cartopy] Step: 10 Points: 273  Time: 1433ms
[Cartopy] Step: 5 Points: 528  Time: 2887ms
[Basemap] Step: 20 Points: 123  Time: 118ms
[Basemap] Step: 10 Points: 273  Time: 151ms
[Basemap] Step: 5 Points: 528  Time: 179ms

which means Cartopy can be up to 16x slower on my old laptop. After some profiling work I find the culprit is https://github.com/SciTools/cartopy/blob/master/lib/cartopy/mpl/geoaxes.py#L1505, it seems this step is to calculate new datalim and autoscale the extent. This is how far I can get to since I don't have enough time to dig in.

• Why does this step take so much time?
• As a quick workaround, can we just check if the extent is already set (which is the most scenario imo) before autoscaling?

Environment

• OS: Windows 10
• cartopy: 0.17.0
• basemap: 1.2.0
• matplotlib: 3.0.2
• shapely: 1.6.4.post1

Thanks!

[efiring]

This is a critical problem. In trying to translate one of my oceanography class notebooks from basemap to cartopy, a 2-second basemap plot took over 2 minutes with cartopy. That's prohibitive. As @crazyapril notes, the problem seems to be coming from a calculation added in 2016: https://github.com/SciTools/cartopy/pull/730/commits/d74eeddb22cc98c7a13f0f69b3c8949f0bd00d3a @pelson, do you remember why this was needed, and why it is done only for filled contours, not for line contours? Is there a good test case that shows the need for this calculation?

[pelson]

Fundamentally cartopy's contour implementation is doing a lot more work than Basemap needs to do, as we are projecting the contour lines once they have been calculated. The commit referenced is ensuring that the data limits are respected by transforming the coordinates from source to target projection.

If this calculation is a prohibitive/expensive operation, then it is always possible to fall back to the equivalent Basemap behaviour - namely, you transform your data yourself, and pass the reprojected data to be contoured in the target project. It isn't ideal as the whole point of cartopy is that you shouldn't have to worry about this stuff (like datelines and poles), but it does highlight that cartopy is a superset of the capabilities available with Basemap.

To illustrate my meaning:

Current code:

ax = plt.axes(projection=target)
ax.contourf(src_x, src_y, data, transform=source)

Becomes:

ax = plt.axes(projection=target)
target_x, target_y, target_data = whatever_used_to_be_done_in_basemap(src_x, src_y, data)
ax.contourf(target_x, target_y, target_data, transform=target)

[efiring]

In the case of contouring, whatever_used_to_be_done_in_basemap is simply applying the projection transformation from lon, lat to the projected coordinates (meters?) via the __call__ method of the Basemap object, which invokes the corresponding proj4 function. The cartopy equivalent would be via the CRS.transform_points method, correct?

There is still the question of why the expensive calculation to check data limits is invoked only for contourf and not for contour.

And, I suspect the calculation is itself much more expensive than it might need to be, because it is having to navigate a tangle of python steps to do many transform calls instead of just transforming a very large number of points in one call into proj4.

I appreciate the convenience, generality, and elegance of cartopy's design. I'm hoping that performance in some of the worst cases can be greatly improved with some optimizations.

[Laura-Guillory]

I've experienced this issue as well, generating a lot of contourf maps in Basemap takes me 15 minutes but in Cartopy it takes 1 hour 4 minutes.

I tried working around it as suggested above with CRS.transform_points, but got ValueError: x and y arrays must have the same length. I think I'm missing something here, since surely you should be able to transform longitudes and latitudes that don't form a perfect square (same length)?

Anyway, hoping that whatever causes this slowdown can be fixed, or a workaround appears.

[dopplershift]

To do it manually with lat/lon arrays that each define 1 dimension of a 2D region, you need to do the work of assembling them into the same size arrays. You can accomplish this using np.meshgrid.

[crazyapril]

So one year later I restarted the migration based on the suggestion above and in case someone may need it, just replace ax.contourf(x, y, data, transform=source_proj) with:

ret = ax.projection.transform_points(source_proj, np.array(x),
    np.array(y), np.array(data)) # This method only accept ndarray!
xx = ret[..., 0]
yy = ret[..., 1]
data = ret[..., 2]
ax.contourf(xx, yy, data)

Still, I think in this case the slowdown it brings somewhat outweigh the convenience. It would be better if this behavior is explained in the documentation. Also, I find if latitudes cross the dateline, using only np.meshgrid as @dopplershift suggested would cause a blank output. A bug?

[dopplershift]

Unless your data actually represent a height above ground, you shouldn't be passing that to transform_points, so I would try the above as:

ret = ax.projection.transform_points(source_proj, np.array(x),
    np.array(y)) # This method only accept ndarray!
xx = ret[..., 0]
yy = ret[..., 1]
ax.contourf(xx, yy, data)

[karlwx]

I don't think the data transformation is the bottleneck. It may make it slower, but the problem here as noted above is the calculation of the data limits. After doing some testing, even without a data transformation cartopy is MUCH slower than matplotlib with contouring.

I use cartopy in several operational scripts and it has gotten to the point where the plotting speed (specifically when using contourf) is a major limiting factor; the plots just can't keep up with the data coming in (unless we had 100's of cpu's available). As the community continues to migrate to cartopy, I think others will run into the same problem. I think this is an issue which should get some more attention to find out where exactly the bottlenecks are and what can be done. I would think that all of the contour math calculations should be vectorized, if possible, as a start.

[greglucas]

I remember profiling all of this a while ago and the bottleneck is in the transform_non_affine function because we are projecting the contour lines from the source projection to target projection as mentioned above. There is a slow for loop to go through all of the geometries and project them into the new projection and currently there really isn't a great way to get around that. There is ongoing work to merge pygeos and Shapely capabilities so that you can do more vectorized calculations which I think would directly apply to this situation.

If you need it operationally fast, did you try the workaround mentioned above?

[karlwx]

Greg, I'm not sure I'm understanding that workaround correctly. In my test, the data are already in axes coordinates. Cartopy is slower even when the coordinates do not need to be transformed. To explain what I mean, here's an example - I am not passing any transform kwarg to contourf - but the contourf call is far slower with a cartopy axes instance than matplotlib:

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

lat = np.arange(-10,90,0.1)
lon = np.arange(0,100,0.1)
lons,lats = np.meshgrid(lon,lat)
data = np.random.rand(len(lat), len(lon))

fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) # cartopy axes
#ax = fig.add_subplot(111) # contouring on regular matplotlib axes is much faster

ax.contourf(lons, lats, data) # no data transform

Here's a screenshot of the profile for contouring on the cartopy axes:

What takes ~15 seconds to do w/ cartopy takes only ~5 seconds with matplotlib alone.

[greglucas]

Thanks @karlwx! I did not realize that was the case even without transforms. I also don't get why the same accurate bounding box calculations aren't computed for the contour case, but only for the filled version here.

I agree this seems like a pretty big speed problem, especially if someone is already setting the extent manually and doesn't care about the limits being updated accurately. I think the lowest hanging fruit would be to add a keyword argument to contourf to tell it to ignore the accurate limit calculations, leave it True by default and allow a user to override it if they want the fast path. https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/mpl/geoaxes.py#L1574-L1580

[karlwx]

After digging into this, it appears the issue runs deeper than we thought. While changing the above code improves speed of the contourf call itself, It turns out there are some serious slow-downs when the image is drawn for display. I added plt.savefig() to the code and here is the result: It looks like something to do with the matplotlib Axes draw wrapper, but I'm not entirely sure how it all fits together.

[Gwi7d31]

I have run into this issue as well trying to use contourf via MetPy. Right now, an RTMA for 2mtr temperatures across the CONUS takes ~5min to plot the contour on an X5650. contour is just fine and I don't see any speed degradations.

@param [String] passed_run_hour

@param [String] passed_final_out_directory

@param [String] passed_netcdf_file_path

def makeTheGraphics(self, passed_run_hour, passed_final_out_directory, passed_netcdf_file_path): data = self.getNetCDFData(passed_netcdf_file_path) print(list(data.variables)) self.logging.info(f'Establishing variables...')

Establish the variables

lat = data['latitude'] lon = data['longitude'] times = data['time'] temp2mtr = data['TMP_2maboveground'][0]

Convert temp to °F

temp2mtr.metpy.convert_units('degF') dp2mtr = data['DPT_2maboveground'][0] visibility = data['VIS_surface'][0]

THIS IS WHERE MULTIPROCESSING SHOULD START

self.logging.info(f'Beginning creation of plot...')

Establish the projection

crs = ccrs.Mercator()

Based upon 100dpi as default 1500x1000 Matplotlib

fig = plt.figure(figsize=(15., 10.), frameon=False)

Matplotlib Adds an Axes to the figure as part of a subplot arrangement

ax = fig.add_subplot(1, 1, 1, projection=crs)

Set the lat/long bounds of the image

region = "conus" bounds = [(-127, -65, 23, 50)] ax.set_extent(*bounds, crs=ccrs.PlateCarree())

Add details to the map like borders and counties

ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.75) ax.add_feature(cfeature.STATES.with_scale('50m'), linewidth=0.5) # edgecolor='#0099ff' ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=1.0)

ax.add_feature(USCOUNTIES.with_scale('5m'))

Apply Title

plt.suptitle('2mtr Temp Test', x=0.25, y=0.25, horizontalalignment='left', verticalalignment='top', fontsize = 16)

ax.set_title(f'Valid: {passed_run_hour}Z', loc='left') self.logging.info(f'Adding contours to plot...')

Contour Plot Of 2mtr Temp

cf = ax.contourf(lon, lat, temp2mtr, range(-40, 120, 5), cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree()) c32f = ax.contour(lon, lat, temp2mtr, [32], transform=ccrs.PlateCarree(), colors='r') cb = fig.colorbar(cf, orientation='horizontal', extend='max') self.logging.info(f'Saving file -- {passed_final_out_directory}{passed_run_hour}/{region}-2mTemp.jpg') fig.savefig(f'{passed_final_out_directory}{passed_run_hour}/{region}-2mTemp.jpg', pad_inches=0.0, format="jpg")

[karlwx]

I had some success improving the speed of hi-res plots like that by using xarray coarsen to resample the data to a lower resolution (did this with National Blend of Models), but it kind of defeats the purpose of the hi-res guidance in the first place. The bottom line is, the plotting time is still excessive.

[Gwi7d31]

I'll mess around with things a little more this evening. I pulled up some code from when I was playing with GOES16 data a year ago or so. In that code, I used pcolormesh to generate the image. I can't remember the difference between that and contourf right now.

[karlwx]

So the difference is pcolormesh does filled grid boxes, while contourf is filled contour lines. At that resolution on a national map, pcolormesh might look okay. It will also be MUCH faster. The disadvantage of pcolormesh is it produces a "blocky" look with low resolution data. contourf is always going to be more visually-appealing, but the trade-off is how slow it is.

[Gwi7d31]

I'll run a benchmark between the two and also drop the images to show comparisons in their quality. This could be an alternative at least when dealing with hi-res data being displayed in CONUS.

[greglucas]

Putting a few cross-referenced issues here with other libraries dealing with the non-affine transforms being slow.

MPL issue with slow non-affine contours: https://github.com/matplotlib/matplotlib/issues/11464

Astropy's fix: https://github.com/astropy/astropy/pull/7568

I think there will still be issues with the contours if one of them gets wrapped and needs to be cut along a boundary or something like that though, so the Astropy fix may not be totally suitable for our use-case.

[greglucas]

@karlwx and anyone else interested, I opened up PR #1690 with a faster contour option. Feel free to take a look at that and comment there if you have other thoughts/ideas on it, or if it doesn't work for your use-case.

[johnrobertlawson]

@karlwx and anyone else interested, I opened up PR #1690 with a faster contour option. Feel free to take a look at that and comment there if you have other thoughts/ideas on it, or if it doesn't work for your use-case.

Hi greglucas: was this ever implemented? I've hunted around for any speed optimisations recently and I don't see anything. I've just benchmarked a whole suite of plotting jobs in NCL, basemap, and cartopy, and NCL beats cartopy in over half the time, and I'm banging my head on a wall trying to rectify this. Any help - or a lay-of-the-land in case I missed something - would be greatly appeared, thanks.

[Gwi7d31]

Hi John,

I can tell you that the code changes do indeed work. I helped to verify the code when this PR came up because I had the exact same issue you are describing. With the issue still being open, the changes were not merged in. What I've personally done is make an Anaconda build with the latest MetPy build and then take the code changes from the PR and apply them to your build.

Or if you're not an Anaconda type person and you just need the cartopy build, then clone the PR'd branch here: https://github.com/greglucas/cartopy/tree/contour_fastpath

Best,

Garrett GwiazdaMeteorologist Senior Software Engineer @.** (802)399-4157wxyard.com* http://www.wxyard.com Queensbury, NY 12804

On Wed, Jun 16, 2021 at 11:38 AM John Lawson @.***> wrote:

@karlwx https://github.com/karlwx and anyone else interested, I opened up PR #1690 https://github.com/SciTools/cartopy/pull/1690 with a faster contour option. Feel free to take a look at that and comment there if you have other thoughts/ideas on it, or if it doesn't work for your use-case.

Hi greglucas: was this ever implemented? I've hunted around for any speed optimisations recently and I don't see anything. I've just benchmarked a whole suite of plotting jobs in NCL, basemap, and cartopy, and NCL beats cartopy in over half the time, and I'm banging my head on a wall trying to rectify this. Any help - or a lay-of-the-land in case I missed something - qouls be greatly appeared, thanks.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/SciTools/cartopy/issues/1291#issuecomment-862485178, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABK6TXB3IH37G2QWKRFELITTTDARFANCNFSM4HDJAQWA .

[johnrobertlawson]

Hi John, I can tell you that the code changes do indeed work. I helped to verify the code when this PR came up because I had the exact same issue you are describing. With the issue still being open, the changes were not merged in. What I've personally done is make an Anaconda build with the latest MetPy build and then take the code changes from the PR and apply them to your build. Or if you're not an Anaconda type person and you just need the cartopy build, then clone the PR'd branch here: https://github.com/greglucas/cartopy/tree/contour_fastpath Best, Garrett GwiazdaMeteorologist

Well, that is highly appreciated, thank you! I am a "conda person" so I'll crack on with this. Many thanks.

[greglucas]

Would you two mind commenting on the PR with your thoughts? If this is a feature that would help you it would be great to get feedback from users on it and try to get it in a release! :)

[Gwi7d31]

Just a nice extra verification, here is my RabbitMQ generating thousands of HRRR graphics an hour with the upload rate at 2 images a second on a rickety old server (2 - X5650 using 14 of the 24 threads) :

[image: image.png]

Garrett GwiazdaMeteorologist Senior Software Engineer @.** (802)399-4157wxyard.com* http://www.wxyard.com Queensbury, NY 12804

On Wed, Jun 16, 2021 at 11:53 AM John Lawson @.***> wrote:

Hi John, I can tell you that the code changes do indeed work. I helped to verify the code when this PR came up because I had the exact same issue you are describing. With the issue still being open, the changes were not merged in. What I've personally done is make an Anaconda build with the latest MetPy build and then take the code changes from the PR and apply them to your build. Or if you're not an Anaconda type person and you just need the cartopy build, then clone the PR'd branch here: https://github.com/greglucas/cartopy/tree/contour_fastpath Best, Garrett GwiazdaMeteorologist Senior Software Engineer @.** (802)399-4157wxyard.com http://wxyard.com* http://www.wxyard.com Queensbury, NY 12804

… <#m-7520066267390444570> On Wed, Jun 16, 2021 at 11:38 AM John Lawson @.***> wrote: @karlwx https://github.com/karlwx https://github.com/karlwx and anyone else interested, I opened up PR #1690 https://github.com/SciTools/cartopy/pull/1690 <#1690 https://github.com/SciTools/cartopy/pull/1690> with a faster contour option. Feel free to take a look at that and comment there if you have other thoughts/ideas on it, or if it doesn't work for your use-case. Hi greglucas: was this ever implemented? I've hunted around for any speed optimisations recently and I don't see anything. I've just benchmarked a whole suite of plotting jobs in NCL, basemap, and cartopy, and NCL beats cartopy in over half the time, and I'm banging my head on a wall trying to rectify this. Any help - or a lay-of-the-land in case I missed something - qouls be greatly appeared, thanks. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#1291 (comment) https://github.com/SciTools/cartopy/issues/1291#issuecomment-862485178>, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABK6TXB3IH37G2QWKRFELITTTDARFANCNFSM4HDJAQWA .

Well, that is highly appreciated, thank you! I am a "conda person" so I'll crack on with this. Many thanks.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/SciTools/cartopy/issues/1291#issuecomment-862496295, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABK6TXDBPW2T4ZCGPLOCZ3TTTDCHTANCNFSM4HDJAQWA .

[wxninja]

Hi all--

Is this PR still open and ongoing? I'll admit it's been an exciting find. We're currently in the middle of a (long overdue) migration from basemap, but are really struggling to engineer a solution to this problem. Contourf just takes too long, especially the dozens--if not hundreds--of images we create routinely through one programming cycle.

I'd love to test this, if it's still an option. I'm working in a Conda framework, which I readily admit is not my specialty, so I'm a little out of my element. I've attempted to PIP Greg's clone for the fast_transform (https://github.com/greglucas/cartopy.git) using 'git' inside the Conda framework, but I must have fouled up the process. When I attempt to employ the 'fast' or 'fast_transform' in code, it fails, telling me it can't recognize that kwarg:

ax.contourf(lons, lats, data, transform=ccrs.PlateCarree(), fast=True) UserWarning: The following kwargs were not used by contour: 'fast'

Anyway, I'm sure I'm doing something wrong...but if anybody can give me a kick in the pants so I can give this a try, I'd appreciate it.

As always, thanks!

[greglucas]

@wxninja, you'll have to be able to install from source to get it to work. But, you will want to checkout the specific branch contour_fastpath, not the master branch. If you go over to the PR https://github.com/SciTools/cartopy/pull/1690, you should be able to keep up-to-date with what is going on.

Here are the GitHub suggested commands you need to check that PR out locally and test it.

git checkout -b greglucas-contour_fastpath master
git pull https://github.com/greglucas/cartopy.git contour_fastpath

Feel free to chime in on the PR over there with comments/suggestions!

[wxninja]

@greglucas Thanks! I'm up and running now. Turns out I just kept getting the syntax wrong when trying to pull 1690. Duh.

Just wanted to chime in and say thank you for engineering this potential solution. I think it'll be a great help, and expand the 'usability' of Cartopy for those who need to make large amounts of graphics, or make graphics quickly.

[johnrobertlawson]

Greg is an absolute treasure to this group and goodness knows how many people he has saved in the more “operational” world. I mean, you’re all treasures really…

On a serious note, I must add yet another two cents that this speed issue is critical for so many users told to migrate from basemap to cartopy. I appreciate the lack of time and resources. I wonder if there’s any ideas or desire over writing a grant proposal to get some real labour devoted to optimisation of cartopy. I’ve only had experience leading proposals in academia, strictly on weather science and not things like this, but the people volunteering their free time on here deserve support! And more help! (I’m running at 100% myself so it’s frustrating I can’t offer more help more often).

John

On Sep 16, 2021, at 06:17, wxninja @.***> wrote:

 @greglucas Thanks! I'm up and running now. Turns out I just kept getting the syntax wrong when trying to pull 1690. Duh.

Just wanted to chime in and say thank you for engineering this potential solution. I think it'll be a great help, and expand the 'usability' of Cartopy for those who need to make large amounts of graphics, or make graphics quickly.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[greglucas]

Thanks for the kind words, @johnrobertlawson!

Note to everyone in this thread: The PR addressing this issue has just been merged and you will be able to use the keyword argument transform_first=True to get the faster behavior in the next release. Be warned that it does not handle wrapped coordinates as well, but on limited domains it should provide a significant speed improvement!

[Gwi7d31]

Greg,

I would also like to say thank you again for your fix. I have been using it since November and it has worked like a charm. It was so desperately needed! And with my best knowledge being in other software libraries, you really did come in and save the day...month...year!?!?!

Best,

Garrett GwiazdaMeteorologist Senior Software Engineer @.***

wxyard.com http://www.wxyard.com Queensbury, NY 12804

On Thu, Sep 16, 2021 at 4:52 PM Greg Lucas @.***> wrote:

Thanks for the kind words, @johnrobertlawson https://github.com/johnrobertlawson!

Note to everyone in this thread: The PR addressing this issue has just been merged and you will be able to use the keyword argument transform_first=True to get the faster behavior in the next release. Be warned that it does not handle wrapped coordinates as well, but on limited domains it should provide a significant speed improvement!

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/SciTools/cartopy/issues/1291#issuecomment-921236305, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABK6TXERTHCMJUEH27UCV53UCJKIDANCNFSM4HDJAQWA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[karlwx]

Thanks for the kind words, @johnrobertlawson!

Note to everyone in this thread: The PR addressing this issue has just been merged and you will be able to use the keyword argument transform_first=True to get the faster behavior in the next release. Be warned that it does not handle wrapped coordinates as well, but on limited domains it should provide a significant speed improvement!

Would it be possible to implement something similar with barbs? I am experiencing significant slowdowns with the barb transformation with high-resolution data. Appears I'm not the only one: https://stackoverflow.com/questions/71578820/use-transform-first-keyword-in-barbs