From efeda78018a8b84b73e6febece6d116a30a405d7 Mon Sep 17 00:00:00 2001
From: Iury Simoes-Sousa <iuryt@pm.me>
Date: Fri, 13 Sep 2024 08:04:33 -0400
Subject: [PATCH] Interactive downsampling large vector fields (example
 gallery) (#1372)

Co-authored-by: maximlt <mliquet@anaconda.com>
---
 doc/reference/xarray/vectorfield.ipynb | 104 +++++++++++++++++++++++++
 1 file changed, 104 insertions(+)

diff --git a/doc/reference/xarray/vectorfield.ipynb b/doc/reference/xarray/vectorfield.ipynb
index 6b0c09de7..ec5761aa2 100644
--- a/doc/reference/xarray/vectorfield.ipynb
+++ b/doc/reference/xarray/vectorfield.ipynb
@@ -31,6 +31,7 @@
    "source": [
     "import numpy as np\n",
     "import xarray as xr\n",
+    "import holoviews as hv\n",
     "import cartopy.crs as ccrs"
    ]
   },
@@ -113,6 +114,109 @@
     "ds.hvplot.vectorfield(x='x', y='y', angle='angle', mag='mag',\n",
     "                      hover=False, geo=True, coastline=True)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Large Data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The visualization of vector fields from large datasets often presents a challenge.  Direct plotting methods can quickly consume excessive memory, leading to crashes or unresponsive applications. To address this issue, we introduce a dynamic downsampling technique that enables interactive exploration of vector fields without sacrificing performance. We first create a `sample_data` that contains 4,000,000 points."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xs, ys, U, V, crs = sample_data(shape=(2000, 2000))\n",
+    "\n",
+    "mag = np.sqrt(U**2 + V**2)\n",
+    "angle = (np.pi/2.) - np.arctan2(U/mag, V/mag)\n",
+    "\n",
+    "ds = xr.Dataset({'mag': xr.DataArray(mag, dims=('y', 'x'), coords={'y': ys, 'x': xs}),\n",
+    "                 'angle': xr.DataArray(angle, dims=('y', 'x'), coords={'y': ys, 'x': xs})}, \n",
+    "                attrs={'crs': crs})\n",
+    "ds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we just try to call `ds.hvplot.vectorfield` this is probably returning `MemoryError`. The alternative is to dynamically downsample the data based on the visible range. This helps manage memory consumption when dealing with large datasets, especially when plotting vector fields. We are going to use HoloViews to create a view (`hv.DynamicMap`) whose content is dynamically updated based on the data range displayed (tracked by the `hv.streams.RangeXY` stream), [find out more about these concepts](https://holoviews.org/user_guide/Custom_Interactivity.html)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def downsample_quiver(x_range=None, y_range=None, nmax=10):\n",
+    "    \"\"\"\n",
+    "    Creates a HoloViews vector field plot from a dataset, dynamically downsampling \n",
+    "    data based on the visible range to optimize memory usage.\n",
+    "\n",
+    "    Args:\n",
+    "        x_range (tuple, optional): Range of x values to include. Defaults to None (full range).\n",
+    "        y_range (tuple, optional): Range of y values to include. Defaults to None (full range).\n",
+    "        nmax (int, optional): Maximum number of points along each axis after coarsening. \n",
+    "                              Defaults to 10.\n",
+    "\n",
+    "    Returns:\n",
+    "        HoloViews DynamicMap: A dynamic vector field plot that updates based on the visible range.\n",
+    "    \"\"\"\n",
+    "\n",
+    "    if x_range is None or y_range is None:\n",
+    "        # No range provided, downsample the entire dataset for initial display\n",
+    "        xs, ys = ds.x.size, ds.y.size        # Get dataset dimensions\n",
+    "        ix, iy = xs // nmax, ys // nmax      # Calculate downsampling intervals\n",
+    "\n",
+    "        ix = max(1, ix)                     # Ensure interval is at least 1\n",
+    "        iy = max(1, iy)\n",
+    "\n",
+    "        sub = ds.coarsen(x=ix, y=iy, side=\"center\", boundary=\"trim\").mean()  # Downsample\n",
+    "    else:\n",
+    "        # Select data within the specified range\n",
+    "        sub = ds.sel(x=slice(*x_range), y=slice(*y_range))\n",
+    "\n",
+    "        # Downsample the selected data\n",
+    "        xs, ys = sub.x.size, sub.y.size\n",
+    "        ix, iy = xs // nmax, ys // nmax\n",
+    "        ix = max(1, ix)\n",
+    "        iy = max(1, iy)\n",
+    "        sub = sub.coarsen(x=ix, y=iy, side=\"center\", boundary=\"trim\").mean()\n",
+    "\n",
+    "    # Create the vector field plot\n",
+    "    quiver = sub.hvplot.vectorfield(\n",
+    "        x=\"x\",\n",
+    "        y=\"y\",\n",
+    "        mag=\"mag\",\n",
+    "        angle=\"angle\",\n",
+    "        hover=False,\n",
+    "    ).opts(magnitude=\"mag\")\n",
+    "\n",
+    "    return quiver"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create interactive plot components\n",
+    "range_xy = hv.streams.RangeXY()               # Stream to capture range changes\n",
+    "filtered = hv.DynamicMap(downsample_quiver, streams=[range_xy])  # Dynamic plot\n",
+    "filtered                                      # Display the plot"
+   ]
   }
  ],
  "metadata": {