ENH: Adding a contour fast-path to speed up contouring

[greglucas]

This adds an option to change the contouring method to use the coordinates in projection-space rather than data-space. It is much faster to project the points beforehand, rather than projecting the contours themselves.

Before adding tests, I'd like to get a feel for whether people even think this is a good idea to have in Cartopy or not? This update basically removes the "transform" keyword from matplotlib.axes.Axes.contourf after already transforming the points to projection-space.

Rationale

#1291 shows how slow filled contours can be in Cartopy because it is projecting the contour regions, not the points. This ends up slowing down the transformation pipeline because of the non-affine transforms of potentially many contour regions.

Implications

A new keyword argument to contourf, but no changes to the defaults.
Closes #1291

Timing example

Similar to standard MPL now, ~100x faster than the standard Cartopy for this example.

MPL: 0.4196739196777344
Cartopy-fast: 0.44753384590148926
Cartopy-standard: 35.90342998504639

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

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

# MPL
fig1 = plt.figure(figsize=(12, 9))
ax = fig1.add_subplot(111)
t0 = time.time()
ax.contourf(lons, lats, data)
fig1.savefig('contourf-mpl.png')
print("MPL:", time.time() - t0)

# Cartopy fast path
fig2 = plt.figure(figsize=(12, 9))
ax = fig2.add_subplot(111, projection=ccrs.Robinson())
t0 = time.time()
ax.contourf(lons, lats, data, transform=ccrs.PlateCarree(), fast=True)
ax.coastlines()
fig2.savefig('contourf-cartopy-fast.png')
print("Cartopy-fast:", time.time() - t0)

# Cartopy standard path
fig3 = plt.figure(figsize=(12, 9))
ax = fig3.add_subplot(111, projection=ccrs.Robinson())
t0 = time.time()
ax.contourf(lons, lats, data, transform=ccrs.PlateCarree())
ax.coastlines()
fig3.savefig('contourf-cartopy.png')
print("Cartopy-standard:", time.time() - t0)

[dopplershift]

I'm certainly interested just based on the performance. I'm trying to work through the mental lift here to decide if the new approach is somehow wrong...

[greglucas]

I'm not sure this is "wrong", so much as a different choice of when to calculate contour lines. I don't think this should be the default, but it may be OK as a good approximation to add for users who want a faster option.

The most obvious difference is the lack of boundary wrap in the coordinates. Here is the example using the fast-path. This will also make the lines straight in projected space rather than curved to follow the true projection, but I think that those are trade-offs users can decide on if they want this option.

[efiring]

It sounds like this option makes the behavior match Basemap, correct?

[Gwi7d31]

I would agree that the kwargs fast option should be defaulted to False as everyone's current scripting would then "break" and produce a potentially unwanted result.

[greglucas]

I think this is the behavior of Basemap, but I am not 100% sure and I don't have a local install on hand right now to test it out. I am curious if Basemap cuts off the edges of the example or if there is potentially a better way of accounting for the wrapped coordinates in projection-space within Basemap.

[greglucas]

I was able to verify that the images look very similar with Basemap. i.e. Basemap does not handle the wrapped coordinates either.
Following on from the previous example:

from mpl_toolkits.basemap import Basemap
fig4 = plt.figure(figsize=(12, 9))
m = Basemap(projection='robin', resolution='c', lon_0=0)
t0 = time.time()
m.contourf(lons, lats, data, latlon=True)
m.drawcoastlines()
fig4.savefig('contourf-basemap.png')
print("Basemap:", time.time() - t0)

Here is the same image with Basemap:

[greglucas]

I had another thought, maybe rather than fast as the keyword we should use something more descriptive that states what this is actually doing and could potentially even be used in other functions too if it is needed. Something to evoke the idea that we are doing the calculations in projection space rather than source space (and then projecting the result) would be valuable I think.

draw_in_projection_space, transform_early, transform_fast, ... Other ideas/suggestions?

[efiring]

I agree that the kwarg name could be improved. Among your suggestions, only "transform_early" works for me. (Drawing is always in projection space; and speed is a "why" not a "what is done".) Other possibilities:

transform_xy transform_grid transform_input

[efiring]

Given that cartopy is relatively immature, I think the choice of default behavior should be considered carefully, and not be assumed to match current behavior.

When the transform is done matters most when the grid is coarse. In that case, is it even clear that the contour vertex calculation is done better in the original coordinate system, in general? And, for the expected long-term user base of cartopy, how common are cases where the path transform approach is superior?

[greglucas]

I think you bring up some great points. The biggest difference between the approaches is with wrapped coordinates. If you have a model that produces output on a longitude grid from 160 - 200, the current behavior produces a patch that goes all the way out to the (-180, 180) bounds with a big gap in the middle.

This new method sorts the transformed points (required by contour) making the bounds go from (-179 to 180) in the contouring algorithm which then fills the internal gap where we really didn't have data originally.

I do think making the calculations in data-space is the better default, but I also think that there are good reasons for wanting to transform the points first and skip the slow non-affine projection of patches.

---

keyword names

The more I think about this, the more I wonder if we should make this an option that could be used by all of the plotting methods, similar to how we have the @_add_transform decorator. So, we should probably make this a generic name that could be used elsewhere too.

I like your transform_xy suggestion. Maybe rather than transform in the name we could use project_xy?

[mikecharles]

Any progress on this?

[greglucas]

I've just updated this PR for other folks that have been requesting it.

• I changed the name to transform_first, which I think gets across the idea that we are transforming the points right away, before doing anything else. Other suggestions on naming are still welcome!

• I pushed the logic up into a decorator and added it to contour as well, so now contour and contourf both get the new keyword argument option.

[QuLogic]

Should we name it something with stage in the name like the new interpolation_stage argument in Matplotlib?

[greglucas]

transform_stage could work. The only negative I have with that is that it would require a string answer to that instead of boolean True/False. I guess the string could be ("projection-space", vs. "data-space") transform_stage="data", but that may be more confusing to explain/use? I still slightly lean towards "transform_first=True" due to the boolean nature, but I'm not opposed to the "stage" suggestion if others have strong preferences.

[QuLogic]

Hmm, I just re-read the docstring, and I wonder if you were interested in writing an example that shows what it means?