`Canvas().raster` issues with out-of-core computation

This is a follow up from [this on Discourse](https://discourse.holoviz.org/t/canvas-raster-failing-when-reading-with-dask/736) where @philippjfr suggested I open an issue.

At the time, I got the example I posted working but I'm not running against a wall with a larger dataset (177GB) and experiencing the same issues, so I thought I'd follow up.

Exploring issues, this should work after #553 and indeed what I'm trying to do is exactly the use case in #560: to generate a visualisation of the entire dataset that is coarse enough that fits in memory.

#### ALL software version info

- `datashader`: 0.11.1
- `xarray`: 0.16.1
- Python 3.8.6

#### Description of expected behavior and the observed behavior

I expect to be able to generate an aggregated `DataArray` object (`agg`) from a larger `xarray.DataArray` that does not fit in memory.

After several attempts with different chunk sizes, the result is always similar: the `compute()` load runs for a while, seems to complete all processing tasks but at the latest stages, I assume when it's merging pieces, it fails on a `KilledWorker` error.

#### Complete, minimal, self-contained example code that reproduces the issue

It is hard to provide a self-contained working example as I suspect the issue here is something that kicks in with a large dataset (indeed, the same code works like a charm in smaller ones).

The dataset I'm connecting to is:

`! rio info $mosaic_url | python -m json.tool`

Click to expand

```json

{
"blockxsize": 128,
"blockysize": 128,
"bounds": [
-222823.73719089525,
-213574.25107683009,
996789.2497132053,
1612237.380579703
],
"colorinterp": [
"gray",
"undefined",
"undefined",
"undefined"
],
"count": 4,
"crs": "EPSG:27700",
"descriptions": [
null,
null,
null,
null
],
"driver": "VRT",
"dtype": "uint16",
"height": 182437,
"indexes": [
1,
2,
3,
4
],
"lnglat": [
-2.2113098420427835,
56.186432587438965
],
"mask_flags": [
[
"nodata"
],
[
"nodata"
],
[
"nodata"
],
[
"nodata"
]
],
"nodata": 0.0,
"res": [
10.007902079383749,
10.007902079383749
],
"shape": [
182437,
121865
],
"tiled": true,
"transform": [
10.007902079383749,
0.0,
-222823.73719089525,
0.0,
-10.007902079383749,
1612237.380579703,
0.0,
0.0,
1.0
],
"units": [
null,
null,
null,
null
],
"width": 121865
}

```

The code I use to set up the attempt is:

```python
import xarray
import datashader as ds
from dask_kubernetes import KubeCluster
from dask.distributed import Client
import dask.array as da

# Set up cluster
cluster = KubeCluster.from_yaml('../../worker-spec.yml')
# Provision with up to X pods
cluster.scale(50)
# Connect Dask to the cluster
client = Client(cluster)

mosaic_url = f"http://{SERVER_IP}:8000/ghs_composite_s2/GHS-composite-S2.vrt"
r = xarray.open_rasterio(mosaic_url,
chunks={"x": 1280, "y": 1280}
)
r
```

The aggregation I'm trying is:

```python
ndvi = (r.sel(band=4) - r.sel(band=1)) / (r.sel(band=4) + r.sel(band=1))

cvs = ds.Canvas(plot_width=1212, plot_height=1824)
agg = cvs.raster(ndvi)
agg.compute()
```

#### Stack traceback and/or browser JavaScript console output

After running for a while and, apparently, completing most of the processing tasks (as monitored in the Dask dashboard), I receive a `KilledWorker` error

This may well be entirely of my fault through something I'm doing that is not a correct use of `xarray`/`datashader` but, since I posted on Discourse and was prompted to report here, I thought I'd give this a try. I'd really like to get this working, `datashader` is such a neat library and enabling us to do this kind of operations makes it feel like superpowers!