Datashader

[jbednar]

Speculative branch for exploring Datashader, HoloViews, and GeoViews. Fine to merge to the datashader branch, but not to master.

So far, includes changes to three files:

• conda_environment.yml
  • Added datashader, hvplot, and geoviews
  • Added the pyviz channel, as that's always the most up to date and is faster when used with conda's strict channel priority due to it being a small channel, but conda-forge alone probably works fine too.
  • Creating the environment is very slow with or without strict channel priority.
• NCL_leg_1.py
  • Very simple example of a plot with a legend, using HoloViews+Matplotlib
  • Basic plot is only one line of code (hv.Curve(vz, label="V") * hv.Curve(uz, label="U"))
  • Customized to look like the original takes four lines of code
• NCL_leg_1.ipynb
  • Jupyter Notebook version of NCL_leg_1.py
  • Shows off Bokeh-based interactive plots (with hover, zoom, etc.)
  • Launch with jupyter notebook NCL_leg_1.ipynb, then Run All
• issue272_case2i_ds.py
  • Extended with a context manager to do timings
  • Sent output to a PNG to allow accurate timing comparisons
  • Accepted a command-line argument factor for selecting the resolution; higher values make it quicker to run
  • Compares timings for Matplotlib-based and SciPy-based Delaunay triangulation (704 vs 529 seconds at 10km res)
  • Compares timings for Matplotlib-based and Datashader-based trimesh plotting (113 vs 3 seconds at 10km res)

The full timings for issue272_case2i_ds.py on my MacBook Pro are:

$ python issue272_case2i_ds.py 200
MPL Triangulation takes: 0.0015921592712402344 seconds
SciPy Triangulation takes: 0.002576112747192383 seconds
Datashader-based rendering 1 takes: 1.8307688236236572 seconds
Datashader-based rendering 2 takes: 0.3103179931640625 seconds
Matplotlib-based rendering takes: 0.17391610145568848 seconds

$ python issue272_case2i_ds.py 10
MPL Triangulation takes: 0.7670328617095947 seconds
SciPy Triangulation takes: 0.8238921165466309 seconds
Datashader-based rendering 1 takes: 1.839228868484497 seconds
Datashader-based rendering 2 takes: 0.40596699714660645 seconds
Matplotlib-based rendering takes: 1.5457789897918701 seconds

$ python issue272_case2i_ds.py 5
MPL Triangulation takes: 3.4231629371643066 seconds
SciPy Triangulation takes: 3.4271531105041504 seconds
Datashader-based rendering 1 takes: 1.807750940322876 seconds
Datashader-based rendering 2 takes: 0.47564196586608887 seconds
Matplotlib-based rendering takes: 4.945454120635986 seconds

$ python issue272_case2i_ds.py 1
MPL Triangulation takes: 704.985121011734 seconds
SciPy Triangulation takes: 529.0948584079742 seconds
Datashader-based rendering 1 takes: 9.263015031814575 seconds
Datashader-based rendering 2 takes: 3.3065507411956787 seconds
Matplotlib-based rendering takes: 113.213791847229 seconds

Note that in issue272_case2i_ds.py the np.arange() calls don't actually cover lon_max or lat_max, which you can see clearly in the factor=200 run (big whitespace on top and right of plot). Probably np.linspace() is more appropriate here, but even then you have to worry what happens at the latitude edge; presumably those values should be copied from one edge to the other.

Takeaways

• SciPy's Delaunay triangulation is a bit faster for large arrays, but maybe not that dramatically
• Datashader's first rendering is more than 10X faster than Matplotlib for large arrays. The first rendering includes JIT compilation of the code, which we will eventually be able to cache, but right now is done again each time the library is first used.
• Datashader is another 3X faster for subsequent renderings, which probably (untested) is what applies to interactive usage
• Triangulation dominates the overall time, but in practice I think most people with gridded data will create the grid once and render data on that grid many times, so the rendering step seems more important to speed up.

[jbednar]

I haven't gotten a chance to convert issue272_case1.py, and probably won't get a chance between now and Thanskgiving as I have a very busy week ahead. That's a pretty plot, though!

From a quick skim of the notes about it:

• Overlaying filled contours, streamlines, and vectors over the same map
  • Contours: hv.Contours(filled=True),
  • Streamlines: Supported as Contours, but not currently with arrows
  • Vectors: hv.VectorField(color=...)
  • Map: gf.ocean * gf.land * gf.coastline
  • Overlaying: Should be no issue for hv/gv once you have plots p1, p2, etc. of each of those items separately; just do (p1 * p2 * p3).
• Using inset_axes() to create additional axes for color bars
  • I haven't seen multiple colorbars on a single HoloViews plot, so I'd have to ask @jlstevens or @philippjfr if that's currently possible
• Creating custom label formats for colorbars
  • I think that's doable, but I haven't done it myself
• Creating a quiverkey
  • That looks like an hv.Arrow and a hv.Label above an hv.Rectangle
  • Not sure if those items can be placed in screen (not data) coordinates, though that only matters for a zoomable Bokeh plot
• Assigning a colormap to contour and quiver plots
  • Supported by both hv.Contour and hv.VectorField
• Add arrows to streamlines
  • Not currently supported as far as I know, but could be added.
• Using zorder to specify the order in which elements will be drawn
  • Order is determined by the order plots are specified, i.e. (p1 * p2 * p3) (back to front).

[erogluorhan]

@jbednar thank you very much for this demonstration. I am also cc'ing @clyne to continue our discussions here. I'd like to share a few thoughts and findings from my first pass:

1. NCL_leg_1.py example and its Jupyter notebook are just right to show the higher level implementation capabilities of Holoviews and Bokeh's interactive plotting convenience.

2. issue272_case1.py would be great to see above capabilities in a prettier plotting case including contours, streamlines, and vectors over the same map. There is no rush though for this; moreover, the GeoCAT team would work on implementing this case with Holoviews capabilities.

3. I have tried a number of resolutions (2000 km, 10 km, 3 km, or factor=200, 1, 0.3 respectively) for issue272_case2i_ds.py. Even though the Matplotlib-based and SciPy-based Delaunay triangulations have comparable performances over any resolution, the rendering performance of Datashader is very impressive (as you already pointed out 113 vs 3 seconds at 10km resolution, and even better for finer res)

4. Despite the great rendering performance of Datashader for most of the resolutions, I was unable to get an output successfully from Datashader for roughly 3 km resolution (factor=0.3) in my local Mac after running for several hours, which was also the case in matplotlib-based implementation. It would be interesting to see outputs for such or finer resolutions or to have further analysis of this (such as memory profiling etc.) to see the margins.

We'd have further analysis and discussions.

[clyne]

Just wanted to note that the 3km resolution is significant because that is what is used by many in the global numerical weather prediction research community. There are also groups running at ~1km resolution.

[jbednar]

Does a 3km or 1km resolution grid fit in memory on your machine? If not, you'd need to use a Dask-backed xarray, either out of core or distributed across multiple machines. Many of the people in Pangeo could help you get set up with an example of doing that, if that's the issue here...

[clyne]

@jbednar, can datashader.rasterize() be parallelized?

[jbednar]

Yes, if you provide a Dask array with your data; it will use whatever Dask workers are available to work on that array.

[clyne]

Thanks. Are there any examples you can point to? My attempts to daskize the TriMesh thus far result in a runaway process, but I could be doing something stupid. Ignoring chunking for now:

tris_da = da.from_array(tris)
nodes = hv.Points(points, vdims='z')
trimesh = hv.TriMesh((tris_da, nodes)).opts(filled=True, edge_color='z')

[jbednar]

Oh, if it's for trimesh, what's supported is a Dask DataFrame, not a Dask Array (not sure if that's what you've got in the code snippet). The matrix of what type of data object is supported for which type of plot is shown in a table at https://datashader.org/user_guide/Performance.html#Data-objects . As you can see, there are a lot of combinations that are supported, most of which are difficult to build and test on CI where multi-node and GPU jobs are painful and expensive. We have not done a good job of showing a tested example of each supported combination; instead for multi-node and GPU examples there are just some random notebooks linked from PRs, etc. It's definitely important to get good multi-node and GPU examples together, but it will take a lot of work that we're not funded for, so it's a hole in what we're currently offering. If just switching to a Dask DataFrame isn't enough to get you going with this example, let me know, and I can see if we can make the specific example use Dask. And if anyone knows of someone willing to volunteer or fund such work as a group, I'd be happy to supervise and guide it!

[clyne]

I'm afraid I'm flailing on this. I've not been able to figure out how to feed hv.TriMesh a Dask DataFrame that it likes. I'm hoping to discuss datashader+dask at a talk I'm giving at RMACC next week if I can get this sorted out. Any further guidance would be greatly appreciated. Thanks!

[jbednar]

I've asked one of my colleagues to adapt a trimesh example to use a Dask DataFrame for you to have as a guide. Hopefully he can do that soon!

[clyne]

Awesome. Thanks so much!

[jbednar]

Ok, here's the result: holoviz/holoviews#4927
Looks like there is some very unintuitive behavior involved, which we can fix, but the code in that issue should be runnable. Please chime in on that issue if you have trouble getting that example to work, and on this issue if you have trouble adapting that example to the one in this issue.

[clyne]

Thanks much, Jim. Greatly appreciated. I'll see if we can get it adapted ASAP and report back.

[clyne]

I was able to adapt the code here holoviz/holoviews#4927 to a data set with ~80M triangles.
However, I've tried multiple combination of cores, processes, and partitions, but have not been able to see any kind of speedup, or indication that more than one core is busy on 18-core skylake processor or 4-core Intel i7.

[jbednar]

Right; with Dask there's a first step of "get it working with Dask", and then a separate step of "once it works with Dask, configure a scheduler and set of workers to make good use of the available hardware". To get started with this second step, try https://distributed.dask.org/en/latest/quickstart.html, or if you've already set up the distributed scheduler and configured workers appropriately, please post the working (but not sped up) code that you have so far (and corresponding dataset) and we can see if we can spot something.

[jbednar]

That said, I'm surprised that the default configuration isn't automatically using all the available cores on a single-CPU machine. We'll look at that tomorrow using our examples and see if we can spot anything amiss.

[clyne]

I've attached the notebook here.

The data can be download from:

https://drive.google.com/drive/folders/17uxAJCm81qPA1fuwukzHaz4AnWDK7-OX?usp=sharing

Thanks!!

MPAS_Datashader_Test.ipynb.gz

[philippjfr]

Okay I had a look at it and discovered that the dask codepath in Datashader was indeed not being exercised by HoloViews. I've fixed that in this PR holoviz/holoviews#4935 and will release a new version of HoloViews with this fix (version 1.14.4) some time today. I've also taken the liberty to clean up the notebook:

MPAS_Datashader_Trimesh.ipynb.zip

Here's a gif of the rasterization in action:

[clyne]

Outstanding! Thank you much, @philippjfr and @jbednar. With the code restructuring even the serial version is performing much better. I look forward to getting your fix incorporated and will report back later. Is there a way to instrument the rasterizer to time the rendering updates (or at least the first rendering)?

[jbednar]

It would be easy to time Datashader when it's used as a standalone library rendering to an xarray object, but I'm not sure if there is any easy way to time Datashader as it's used inside HoloViews. If it's ok to include the entire process of rendering and plotting you can do something like:

I'd run that a few times to see if it varies, e.g. if something needed to be imported the first time.

[clyne]

Works great. Thanks!