diff --git a/CITATION.cff b/CITATION.cff index 3cc5003..f8d05f3 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -11,6 +11,21 @@ authors: orcid: https://orcid.org/0000-0002-2666-8493 website: https://github.com/anissa111 affiliation: UCAR/NCAR + - family-names: Eroglu + given-names: Orhan + orcid: https://orcid.org/0000-0003-3099-8775 + website: https://github.com/erogluorhan + affiliation: UCAR/NCAR + - family-names: Chmielowiec + given-names: Philip + orcid: + website: https://github.com/philipc2 + affiliation: UCAR/NCAR + - family-names: Clyne + given-names: John + orcid: https://orcid.org/0000-0003-2788-9017 + website: https://github.com/clyne + affiliation: UCAR/NCAR - name: "Advanced Visualization Cookbook contributors" # use the 'name' field to acknowledge organizations website: "https://github.com/ProjectPythia/advanced-viz-cookbook/graphs/contributors" title: "Advanced Visualization Cookbook" diff --git a/README.md b/README.md index 035cf74..0077713 100644 --- a/README.md +++ b/README.md @@ -9,13 +9,13 @@ This Project Pythia Cookbook covers advanced visualization techniques building u ## Motivation -The possibilities of data visualization in Python are almost endless. Already using `matplotlib` the workhorse behind many visualization packages, the user has a lot of customization options available to them. `cartopy`, `metpy`, `seaborn`, `geocat-viz`, and `datashader` are all also great packages that can offer unique additions to your Python visualization toolbox. +The possibilities of data visualization in Python are almost endless. Already using `matplotlib` the workhorse behind many visualization packages, the user has a lot of customization options available to them. `cartopy`, `metpy`, `seaborn`, `geocat-viz`, and `datashader` are all also great packages that can offer unique additions to your Python visualization toolbox. -This Cookbook will house various visualization workflow examples that use different visualization packages, highlight the differences in functionality between the packages, any noteable syntax distinctions, and demonstrate combining tools to achieve a specific outcome. +This Cookbook will house various visualization workflow examples that use different visualization packages, highlight the differences in functionality between the packages, any noteable syntax distinctions, and demonstrate combining tools to achieve a specific outcome. ## Authors -[Julia Kent](@jukent), [Anissa Zacharias](@anissa111) +[Julia Kent](@jukent), [Anissa Zacharias](@anissa111), [Orhan Eroglu](@erogluorhan), [Philip Chmielowiec](@philipc2), [John Clyne](@clyne) ### Contributors @@ -43,6 +43,10 @@ In this section we will demonstrate how to visualize data that is on a structure Animated plots are great tools for science communication and outreach. We will demonstrate how to make your plots come to life. In this book, we use "animated plots" to refer to stable animations, such as the creation of gifs or videos. +### Interactivity + +Dynamically rendering, animating, panning & zooming over a plot can be great to increase data fidelity. We will showcase how to use Holoviz technologies with Bokeh backend to create interactive plots, utilizing an unstructured grid data in the Model for Prediction Across Scales (MPAS) format. + ## Running the Notebooks You can either run the notebook using [Binder](https://binder.projectpythia.org/) or on your local machine. diff --git a/_toc.yml b/_toc.yml index a3e1f89..113b725 100644 --- a/_toc.yml +++ b/_toc.yml @@ -6,9 +6,9 @@ parts: - file: notebooks/how-to-cite - caption: Basics of Geoscience Visualization chapters: - - file: notebooks/comparison - - file: notebooks/good-viz - - file: notebooks/plot-elements + - file: notebooks/comparison + - file: notebooks/good-viz + - file: notebooks/plot-elements - caption: Specialty Plots chapters: - file: notebooks/taylor-diagrams @@ -16,9 +16,9 @@ parts: - caption: Visualization of Structured Grids chapters: - file: notebooks/spaghetti - - caption: Interactive Visualization - chapters: - - file: notebooks/mpas-datshader - caption: Animation chapters: - file: notebooks/animation + - caption: Interactivity + chapters: + - file: notebooks/interactive-holoviz-mpas diff --git a/environment.yml b/environment.yml index a1e61ec..cdb9cb0 100644 --- a/environment.yml +++ b/environment.yml @@ -18,6 +18,7 @@ dependencies: - uxarray - datashader - geocat-datafiles + - geoviews - tropycal - pip - pip: diff --git a/notebooks/animation.ipynb b/notebooks/animation.ipynb index 65b20ac..11b4a37 100644 --- a/notebooks/animation.ipynb +++ b/notebooks/animation.ipynb @@ -469,7 +469,9 @@ "\n", "Creating animations in `matplotlib` might seem intimidating, but is easier when you know the options and purpose of each method. These visualizations can be a powerful tool to display and understand time-dependent data.\n", "\n", - "### What's next?\n" + "### What's next?\n", + "\n", + "In the final section of this cookbook, let’s look at [interactive plotting with Holoviz](interactive-holoviz-mpas) tools." ] }, { @@ -504,7 +506,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.11.6" } }, "nbformat": 4, diff --git a/notebooks/images/hv-dynamic-zoomed.png b/notebooks/images/hv-dynamic-zoomed.png new file mode 100644 index 0000000..d195187 Binary files /dev/null and b/notebooks/images/hv-dynamic-zoomed.png differ diff --git a/notebooks/images/hv-non-dynamic-zoomed.png b/notebooks/images/hv-non-dynamic-zoomed.png new file mode 100644 index 0000000..7218a42 Binary files /dev/null and b/notebooks/images/hv-non-dynamic-zoomed.png differ diff --git a/notebooks/interactive-holoviz-mpas.ipynb b/notebooks/interactive-holoviz-mpas.ipynb new file mode 100644 index 0000000..7459536 --- /dev/null +++ b/notebooks/interactive-holoviz-mpas.ipynb @@ -0,0 +1,737 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "# HoloViz Visualization" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "The ability to dynamically render, pan, zoom, animate and perform other dynamic operations on data can provide many benefits, such as providing greater data fidelity within the same plot. [HoloViz](https://holoviz.org/) provides high-level tools (such as Holoviews, Datashader, Geoviews, etc.) to visualize even the very large datasets efficiently.\n", + "\n", + "This notebook explores interactively plotting using an unstructured grid dataset in the [MPAS](https://mpas-dev.github.io/) format with Holoviews, Datashader, and Geoviews." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{attention}\n", + "If you are in this notebook to actually learn about and/or visualize unstructured grids, we highly recommend checking out the [UXarray Cookbook](https://projectpythia.org/unstructured-grid-viz-cookbook/README.html) that is a comprehensive showcase of workflows & techniques for visualizing Unstructured Grids using [UXarray](https://uxarray.readthedocs.io/). UXarray is a Python package that:\n", + "- Provides Xarray-styled functionality for working with unstructured grids built around the UGRID conventions\n", + "- Supports not only MPAS but also several other formats such as UGRID, SCRIP, and Exodus\n", + "- Enables significant data analysis and visualization functionality for unstructured grid research, which makes working with unstructured grids a breeze\n", + " - e.g. UXarray internally handles majority of the utility functions and data preparation steps mentioned throughout this notebook and provides user with convenient visualization and analysis functions\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook demonstrates:\n", + "1. Use of HoloViz tools for interactive plotting\n", + "2. Different interactivity schemes\n", + "3. Use of the MPAS format's connectivity information to render data on the native grid (hence avoiding costly Delaunay triangulation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The flow of the content is as follows:\n", + "\n", + "1. Package imports\n", + "2. MPAS preprocessing for visualization\n", + " - Utility functions\n", + " - Data loading\n", + " - Triangular mesh generation using MPAS's cell connectivity array from the primal mesh\n", + "3. Interactive Holoviz Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Prerequisites\n", + "\n", + "| Concepts | Importance | Notes |\n", + "| --- | --- | --- |\n", + "| [Holoviews](https://holoviews.org/) | Necessary | |\n", + "| [Geoviews](https://geoviews.org/) | Useful | Not necessary for plotting but useful for adding features |\n", + "| [Matplotlib](https://matplotlib.org/) | Useful | |\n", + "| [MPAS](https://mpas-dev.github.io/) | Useful | Not necessary for interactive plotting but useful for understanding the data used |\n", + "| [Xarray](https://xarray.pydata.org/) | Useful | |\n", + "\n", + "- **Time to learn**: 60 minutes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import math as math\n", + "\n", + "import cartopy.crs as ccrs\n", + "import dask.dataframe as dd\n", + "import geocat.datafiles as gdf # Only for reading-in datasets\n", + "import geoviews.feature as gf # Only for displaying coastlines\n", + "import holoviews as hv\n", + "import numpy as np\n", + "import pandas as pd\n", + "from holoviews import opts\n", + "from holoviews.operation.datashader import rasterize as hds_rasterize\n", + "from numba import jit\n", + "from xarray import open_mfdataset\n", + "\n", + "n_workers = 1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## MPAS Preprocessing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The MPAS data requires some preprocessing to get it ready for visualization such as implementation \n", + "of a few utility functions, loading the data, and creating triangular mesh out of the data to rasterize." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Utility functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def unzipMesh(x, tris, t):\n", + " \"\"\"Splits a global mesh along longitude.\n", + "\n", + " Examine the X coordinates of each triangle in 'tris'. Return an array of 'tris' where\n", + " only those triangles with legs whose length is less than 't' are returned.\n", + "\n", + " Parameters\n", + " ----------\n", + " x: numpy.ndarray\n", + " x-coordinates of each triangle in 'tris'\n", + " tris: numpy.ndarray\n", + " Triangle indices for each vertex in the MPAS file, in counter-clockwise order\n", + " t: float\n", + " Threshold value\n", + " \"\"\"\n", + " return tris[\n", + " (np.abs((x[tris[:, 0]]) - (x[tris[:, 1]])) < t)\n", + " & (np.abs((x[tris[:, 0]]) - (x[tris[:, 2]])) < t)\n", + " ]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def triArea(x, y, tris):\n", + " \"\"\"Computes the signed area of a triangle.\n", + "\n", + " Parameters\n", + " ----------\n", + " x: numpy.ndarray\n", + " x-coordinates of each triangle in 'tris'\n", + " y: numpy.ndarray\n", + " y-coordinates of each triangle in 'tris'\n", + " tris: numpy.ndarray\n", + " Triangle indices for each vertex in the MPAS file\n", + " \"\"\"\n", + " return ((x[tris[:, 1]] - x[tris[:, 0]]) * (y[tris[:, 2]] - y[tris[:, 0]])) - (\n", + " (x[tris[:, 2]] - x[tris[:, 0]]) * (y[tris[:, 1]] - y[tris[:, 0]])\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def orderCCW(x, y, tris):\n", + " \"\"\"Reorder triangles as necessary so they all have counter clockwise winding order.\n", + " CCW is what Datashader and MPL require.\n", + "\n", + " Parameters\n", + " ----------\n", + " x: numpy.ndarray\n", + " x-coordinates of each triangle in 'tris'\n", + " y: numpy.ndarray\n", + " y-coordinates of each triangle in 'tris'\n", + " tris: numpy.ndarray\n", + " Triangle indices for each vertex in the MPAS file\n", + " \"\"\"\n", + " tris[triArea(x, y, tris) < 0.0, :] = tris[triArea(x, y, tris) < 0.0, ::-1]\n", + " return tris" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def createHVTriMesh(x, y, triangle_indices, var, n_workers=1):\n", + " \"\"\"Create a Holoviews Triangle Mesh suitable for rendering with Datashader\n", + "\n", + " This function returns a Holoviews TriMesh that is created from a list of coordinates, 'x'\n", + " and 'y', an array of triangle indices that addresses the coordinates in 'x' and 'y', and\n", + " a data variable 'var'. The data variable's values will annotate the triangle vertices\n", + "\n", + " Parameters\n", + " ----------\n", + " x: numpy.ndarray\n", + " Projected x-coordinates of each triangle in 'tris'\n", + " y: numpy.ndarray\n", + " Projected y-coordinates of each triangle in 'tris'\n", + " triangle_indices: numpy.ndarray\n", + " Triangle indices for each vertex in the MPAS file, in counter-clockwise order\n", + " var: numpy.ndarray\n", + " Data variable from which the triangle vertex values are read.\n", + " n_workers: int\n", + " Number of workers, for Dask\n", + " \"\"\"\n", + " # Declare verts array\n", + " verts = np.column_stack([x, y, var])\n", + "\n", + " # Convert to pandas\n", + " verts_df = pd.DataFrame(verts, columns=[\"x\", \"y\", \"z\"])\n", + " tris_df = pd.DataFrame(triangle_indices, columns=[\"v0\", \"v1\", \"v2\"])\n", + "\n", + " # Convert to dask\n", + " verts_ddf = dd.from_pandas(verts_df, npartitions=n_workers)\n", + " tris_ddf = dd.from_pandas(tris_df, npartitions=n_workers)\n", + "\n", + " # Declare HoloViews element\n", + " tri_nodes = hv.Nodes(verts_ddf, [\"x\", \"y\", \"index\"], [\"z\"])\n", + " trimesh = hv.TriMesh((tris_ddf, tri_nodes))\n", + " return trimesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "@jit(nopython=True)\n", + "def triangulatePoly(verticesOnCell, nEdgesOnCell):\n", + " \"\"\"Triangulate MPAS dual mesh:\n", + "\n", + " Triangulate each polygon in a heterogenous mesh of n-gons by connecting\n", + " each internal polygon vertex to the first vertex. Uses the MPAS\n", + " auxilliary variables verticesOnCell, and nEdgesOnCell.\n", + "\n", + " The function is decorated with Numba's just-in-time compiler so that it is translated into\n", + " optimized machine code for better peformance\n", + "\n", + " Parameters\n", + " ----------\n", + " verticesOnCell: numpy.ndarray\n", + " Connectivity array that stores the vertex indices that surround a given cell\n", + " nEdgesOnCell: numpy.ndarray\n", + " Number of edges on a given cell.\n", + " \"\"\"\n", + "\n", + " # Calculate the number of triangles. nEdgesOnCell gives the number of vertices for each cell (polygon)\n", + " # The number of triangles per polygon is the number of vertices minus 2.\n", + " nTriangles = np.sum(nEdgesOnCell - 2)\n", + "\n", + " triangles = np.ones((nTriangles, 3), dtype=np.int64)\n", + " nCells = verticesOnCell.shape[0]\n", + " triIndex = 0\n", + " for j in range(nCells):\n", + " for i in range(nEdgesOnCell[j] - 2):\n", + " triangles[triIndex][0] = verticesOnCell[j][0]\n", + " triangles[triIndex][1] = verticesOnCell[j][i + 1]\n", + " triangles[triIndex][2] = verticesOnCell[j][i + 2]\n", + " triIndex += 1\n", + "\n", + " return triangles" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load data and coordinates" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{info}\n", + "The \"dyamond_1\" global datasets used in this example are courtesy of NCAR's Falko Judt and were produced as part of the \n", + "[DYAMOND](http://dx.doi.org/10.1186/s40645-019-0304-z) initiative. They are all from the same experiment but run at several \n", + "resolutions from 30km to 3.75km.\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Currently, the 30-km resolution dataset is used in this example and is available at [geocat-datafiles](https://github.com/NCAR/geocat-datafiles).\n", + "However, the other resolutions of these data are stored on NCAR's Glade data storage because of their size. Due to their size, \n", + "the higher resolution data sets are only distributed with two variables in them:\n", + "\n", + "+ relhum_200hPa: Relative humidity vertically interpolated to 200 hPa\n", + "+ vorticity_200hPa: Relative vorticity vertically interpolated to 200 hPa\n", + "\n", + "The \"relhum_200hPa\" variable is computed on the MPAS 'primal' mesh, while \"vorticity_200hPa\" is computed on the MPAS\n", + "'dual' mesh." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Load data\n", + "datafiles = (\n", + " gdf.get(\n", + " \"netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/diag.2016-08-20_00.00.00_subset.nc\"\n", + " ),\n", + " gdf.get(\"netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/x1.655362.grid_subset.nc\"),\n", + ")\n", + "\n", + "primalVarName = \"relhum_200hPa\"\n", + "dualVarName = \"vorticity_200hPa\"\n", + "central_longitude = 0.0\n", + "\n", + "ds = open_mfdataset(datafiles, decode_times=False)\n", + "primalVar = ds[primalVarName].isel(Time=0).values\n", + "dualVar = ds[dualVarName].isel(Time=0).values\n", + "\n", + "# Fetch lat and lon coordinates for the primal and dual mesh.\n", + "lonCell = ds[\"lonCell\"].values * 180.0 / math.pi\n", + "latCell = ds[\"latCell\"].values * 180.0 / math.pi\n", + "lonCell = ((lonCell - 180.0) % 360.0) - 180.0\n", + "\n", + "lonVertex = ds[\"lonVertex\"].values * 180.0 / math.pi\n", + "latVertex = ds[\"latVertex\"].values * 180.0 / math.pi\n", + "lonVertex = ((lonVertex - 180.0) % 360.0) - 180.0" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "### Generate triangular mesh using MPAS connectivity information" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "#### Primal mesh\n", + "\n", + "In this example, we use the MPAS `cellsOnVertex` auxilliary variable, which defines mesh connectivity for the MPAS grid. Specifically, this variable tells us the cell indices contained by each cell." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + ":::{important}\n", + "The benefits of this are twofold: \n", + "1. We're using the actual mesh description from the MPAS output file\n", + "2. For large grid this is *much* faster than synthesizing the connectivity information (e.g. by triangulating them with, for example, Delaunay triangulation).\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "So, first let's: \n", + "- Get the triangle indices (i.e. MPAS connectivity),\n", + "- Make sure it is all in counter-clockwise order,\n", + "- \"Unzip\" the mesh along the longitude (for now, assuming the central longitude from the map projection is 0.0, i.e. cutting the mesh where longitude wraps around from -180.0 to 180.0)\n", + "- Project vertices from geographic to PCS coordinates\n", + "- Create Holoviews TriMesh" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{note}\n", + "Indexing in MPAS starts from 1, not zero :-(\n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "tris = ds.cellsOnVertex.values - 1\n", + "\n", + "tris = orderCCW(lonCell, latCell, tris)\n", + "\n", + "tris = unzipMesh(lonCell, tris, 90.0)\n", + "\n", + "projection = ccrs.Robinson(central_longitude=central_longitude)\n", + "xPCS, yPCS, _ = projection.transform_points(ccrs.PlateCarree(), lonCell, latCell).T\n", + "\n", + "trimesh = createHVTriMesh(xPCS, yPCS, tris, primalVar, n_workers=n_workers)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Dual mesh\n", + "\n", + "In this example though, we use the MPAS `verticesOnCell` and `nEdgesOnCell` auxilliary variables, which defines mesh connectivity for the MPAS dual grid. A lot of details in the following code is similarv to those in the primal mesh's case except `triangulatePoly()` wheere we decompose each cell into triangles since for the dual mesh, the data are located on triangle centers, which correspond to cell (polygon) vertices." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "verticesOnCell = ds.verticesOnCell.values - 1\n", + "nEdgesOnCell = ds.nEdgesOnCell.values\n", + "\n", + "tris_dual = triangulatePoly(verticesOnCell, nEdgesOnCell)\n", + "\n", + "tris_dual = unzipMesh(lonVertex, tris_dual, 90.0)\n", + "\n", + "projection = ccrs.Robinson(central_longitude=central_longitude)\n", + "xPCS_dual, yPCS_dual, _ = projection.transform_points(\n", + " ccrs.PlateCarree(), lonVertex, latVertex\n", + ").T\n", + "\n", + "trimesh_dual = createHVTriMesh(\n", + " xPCS_dual, yPCS_dual, tris_dual, dualVar, n_workers=n_workers\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interactive HoloViz Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{note} \n", + "Since the emphasis of these plots is the interactivity, we will not provide much detail on the features/arguments other than interactivity-related ones, but such details to some extent can still be found in the code comments.\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "HoloViz tools provide essential interactivity features such as panning, zooming, \n", + "hovering, clickable/selectable legends, etc. on a plot." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{note} \n", + "The first step to achieve interactivity with HoloViz tools is to choose `bokeh` as \n", + "the backend: \n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "hv.extension(\"bokeh\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's define keyword arguments that will be commonly used for all of the plots and their `opts` throughout. For that, the plotting arguments are `aggregator=mean` to use the average of node values for shading in the rasterization and `precompute=True` to cache the data internally; the `opts` arguments are nothing more than opting in to showing the colorbar and using a colormap called \"coolwarm\":" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_kwargs = {\n", + " \"aggregator\": \"mean\",\n", + " \"precompute\": True,\n", + "}\n", + "\n", + "opts_kwargs = {\n", + " \"colorbar\": True,\n", + " \"cmap\": \"coolwarm\",\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Holoviews' options system: `opts`\n", + "\n", + "HoloViz packages provide high-level visualization functions that allow the customization of plot features through optional arguments. Furthermore, the HoloViews options system allows customizing a plot from title to size, legend, axes, and more. This system works by calling the `.opts()` method through a HoloViews plot object. We will set the `opts` defaults below and then customize plots throughout the notebook via `.opts()`. Keep an eye on them!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "opts.defaults(\n", + " opts.Image(frame_width=600, data_aspect=1),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pan & zoom\n", + "\n", + "Bokeh-rendered HoloViz plots come with default pan and zoom (box & wheel), tools as well as save and reset options in the very right. Just look at the following plot for those features, which uses Datashader's rasterization method to visualize the data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "rasterized = hds_rasterize(trimesh, **plot_kwargs)\n", + "rasterized.opts(**opts_kwargs) * gf.coastline(projection=projection)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Hover\n", + "\n", + "Use `tools=['hover']` can be used to get the plot to display data values while hovering over the plot (there are other tools that can be selected, too, but those will not be covered in this notebook). See the below plot (let's use the dual mesh data this time just for the sake of assortment of visualizations) and hover you cursor over it:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "rasterized = hds_rasterize(trimesh_dual, **plot_kwargs)\n", + "rasterized.opts(tools=[\"hover\"], clim=(-4e-4, 4e-4), **opts_kwargs) * gf.coastline(\n", + " projection=projection\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Dynamic plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can decide whether to return a dynamic plot that sends updates on widget and zoom/pan \n", + "events or whether all the data should be embedded by using the boolean `dynamic` argument. \n", + "The default value of this argument is `True`, so all the previously rendered plots were \n", + "already dynamic. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{warning}\n", + "If you are looking at this through the Pythia Cookbook gallery webpage (i.e. this is a rendered notebook as html), you cannot experience how the \"dynamic\" option works. To do that, you have to execute the notebook live (e.g. through the Binder link we provide for the cookbook or by executing it in your local).\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's look at a non-dynamic plot (i.e. `dynamic=False`):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rasterized = hds_rasterize(trimesh, **plot_kwargs, dynamic=False)\n", + "rasterized.opts(**opts_kwargs) * gf.coastline(projection=projection)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{tip}\n", + "This plot does not look immediately different than the prior ones where `dynamic=True`, right? Yes, in the global view, but try zooming into the plot (again, if you are executing notebook rather than looking at the rendered one), and see how dynamic vs. non-dynamic plots look like. See the below comparison for instance:\n", + ":::" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Non-Dynamic-Zoomed | Dynamic-Zoomed\n", + ":-------------------------:|:-------------------------:\n", + "![non-dynamic-zoomed](images/hv-non-dynamic-zoomed.png) | ![dynamic-zoomed](images/hv-dynamic-zoomed.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Further interactivity considerations\n", + "\n", + "Even though we will not cover in the scope of this cookbook, you may want to have a look at the following HoloViz resources for further interactivity in your plots:\n", + "\n", + "- [Mandelbrot](https://holoviews.org/gallery/apps/bokeh/mandelbrot.html)\n", + "- [Player](https://holoviews.org/reference/apps/bokeh/player.html)\n", + "- [HoloMap](https://holoviews.org/reference/containers/bokeh/HoloMap.html)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "HoloViz technologies can be used in order to create interactive plots when Bokeh is used as the backend. Even though a really large dataset (e.g. km-scale) is not \n", + "showcased in this notebook, HoloViz packages are reasonably performant with visualization of such data, too.\n", + "\n", + "### What's next?\n", + "\n", + "The end of this notebook wraps up this cookbook as well. Thanks for reading!\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Resources and references\n", + "\n", + "- [MPAS Mesh Specification](https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf) for those interested" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.6" + }, + "widgets": { + "application/vnd.jupyter.widget-state+json": { + "state": {}, + "version_major": 2, + "version_minor": 0 + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/mpas-datashader.ipynb b/notebooks/mpas-datashader.ipynb deleted file mode 100644 index d0b4613..0000000 --- a/notebooks/mpas-datashader.ipynb +++ /dev/null @@ -1,452 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "tags": [] - }, - "source": [ - "# MPAS with Datashader and Geoviews" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Overview\n", - "\n", - "This example of interactively plotting unstructured grid MPAS data with Datashader and Geoviews demonstrates making use of the MPAS file's connectivity information to render data on the native grid, and also\n", - "avoid costly Delaunay triangulation that is required if the MPAS connectivity information is not used, rendering data that is sampled on both the 'primal' and 'dual' MPAS mesh, using geoviews/holoviews for interactive plotting in a Jupyter Notebook. The plotting is interactive in the sense that you can pan and zoom the data. Doing so will reveal greater and greater data fidelity, and sing Datashader to perform background rendering in place of Matplotlib. Unlike Matplotlib, Datashaderwas designed for performance with large data sets.\n", - "\n", - "1. Utility Functions\n", - "2. Loading Data\n", - "3. Using MPAS's cell connectivity array to plot primal mesh data\n", - "4. Synthesizing triangles from points using Delaunay triangulation\n", - "5. Using MPAS's cell connectivity array to plot dual mesh data" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Prerequisites\n", - "\n", - "| Concepts | Importance | Notes |\n", - "| --- | --- | --- |\n", - "| []() | Necessary | |\n", - "\n", - "- **Time to learn**: X minutes" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "---" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "import cartopy.crs as ccrs\n", - "import numpy as np\n", - "import pandas as pd\n", - "\n", - "import math as math\n", - "\n", - "import geocat.datafiles as gdf # Only for reading-in datasets\n", - "\n", - "from xarray import open_mfdataset\n", - "\n", - "from numba import jit\n", - "\n", - "import dask.dataframe as dd\n", - "\n", - "import holoviews as hv\n", - "from holoviews import opts\n", - "\n", - "from holoviews.operation.datashader import rasterize as hds_rasterize \n", - "#import geoviews.feature as gf # only needed for coastlines\n", - "\n", - "hv.extension(\"bokeh\",\"matplotlib\")\n", - "\n", - "opts.defaults(\n", - " opts.Image(frame_width=600, data_aspect=1),\n", - " opts.RGB(frame_width=600, data_aspect=1))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Utility functions\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# This funtion splits a global mesh along longitude\n", - "#\n", - "# Examine the X coordinates of each triangle in 'tris'. Return an array of 'tris' where only those triangles\n", - "# with legs whose length is less than 't' are returned. \n", - "# \n", - "def unzipMesh(x,tris,t):\n", - " return tris[(np.abs((x[tris[:,0]])-(x[tris[:,1]])) < t) & (np.abs((x[tris[:,0]])-(x[tris[:,2]])) < t)]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Compute the signed area of a triangle\n", - "#\n", - "def triArea(x,y,tris):\n", - " return ((x[tris[:,1]]-x[tris[:,0]]) * (y[tris[:,2]]-y[tris[:,0]])) - ((x[tris[:,2]]-x[tris[:,0]]) * (y[tris[:,1]]-y[tris[:,0]]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Reorder triangles as necessary so they all have counter clockwise winding order. CCW is what Datashader and MPL\n", - "# require.\n", - "#\n", - "def orderCCW(x,y,tris):\n", - " tris[triArea(x,y,tris)<0.0,:] = tris[triArea(x,y,tris)<0.0,::-1]\n", - " return(tris)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Create a Holoviews Triangle Mesh suitable for rendering with Datashader\n", - "#\n", - "# This function returns a Holoviews TriMesh that is created from a list of coordinates, 'x' and 'y',\n", - "# an array of triangle indices that addressess the coordinates in 'x' and 'y', and a data variable 'var'. The\n", - "# data variable's values will annotate the triangle vertices\n", - "#\n", - "\n", - "def createHVTriMesh(x,y,triangle_indices, var,n_workers=1):\n", - " # Declare verts array\n", - " verts = np.column_stack([x, y, var])\n", - "\n", - "\n", - " # Convert to pandas\n", - " verts_df = pd.DataFrame(verts, columns=['x', 'y', 'z'])\n", - " tris_df = pd.DataFrame(triangle_indices, columns=['v0', 'v1', 'v2'])\n", - "\n", - " # Convert to dask\n", - " verts_ddf = dd.from_pandas(verts_df, npartitions=n_workers)\n", - " tris_ddf = dd.from_pandas(tris_df, npartitions=n_workers)\n", - "\n", - " # Declare HoloViews element\n", - " tri_nodes = hv.Nodes(verts_ddf, ['x', 'y', 'index'], ['z'])\n", - " trimesh = hv.TriMesh((tris_ddf, tri_nodes))\n", - " return(trimesh)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Triangulate MPAS primary mesh:\n", - "#\n", - "# Triangulate each polygon in a heterogenous mesh of n-gons by connecting\n", - "# each internal polygon vertex to the first vertex. Uses the MPAS\n", - "# auxilliary variables verticesOnCell, and nEdgesOnCell.\n", - "#\n", - "# The function is decorated with Numba's just-in-time compiler so that it is translated into\n", - "# optimized machine code for better peformance\n", - "#\n", - "\n", - "@jit(nopython=True)\n", - "def triangulatePoly(verticesOnCell, nEdgesOnCell):\n", - "\n", - " # Calculate the number of triangles. nEdgesOnCell gives the number of vertices for each cell (polygon)\n", - " # The number of triangles per polygon is the number of vertices minus 2.\n", - " #\n", - " nTriangles = np.sum(nEdgesOnCell - 2)\n", - "\n", - " triangles = np.ones((nTriangles, 3), dtype=np.int64)\n", - " nCells = verticesOnCell.shape[0]\n", - " triIndex = 0\n", - " for j in range(nCells):\n", - " for i in range(nEdgesOnCell[j]-2):\n", - " triangles[triIndex][0] = verticesOnCell[j][0]\n", - " triangles[triIndex][1] = verticesOnCell[j][i+1]\n", - " triangles[triIndex][2] = verticesOnCell[j][i+2]\n", - " triIndex += 1\n", - "\n", - " return triangles" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Load data and coordinates\n", - "\n", - "The global data sets used in this example are from the same experiment, but run at several resolutions from\n", - "30km to 3.75km. Due to their size, the higher resolution data sets are only distributed with two variables\n", - "in them:\n", - "+ relhum_200hPa: Relative humidity vertically interpolated to 200 hPa\n", - "+ vorticity_200hPa: Relative vorticity vertically interpolated to 200 hPa\n", - "\n", - "The dyamond_1 data set is available in several resolutions, ranging from 30 km to 3.75 km.\n", - "\n", - "Currently, the 30-km resolution dataset used in this example is available at geocat-datafiles.\n", - "However, the other resolutions of these data are stored on Glade because of their size.\n", - "\n", - "The relhum_200hPa is computed on the MPAS 'primal' mesh, while the vorticity_200hPa is computed on the MPAS\n", - "'dual' mesh. Note that data may also be sampled on the edges of the primal mesh. This example does not\n", - "include/cover edge-centered data.\n", - "\n", - "These data are courtesy of NCAR's Falko Judt, and were produced as part of the DYAMOND initiative: \n", - " http://dx.doi.org/10.1186/s40645-019-0304-z.\n" - ] - }, - { - "cell_type": "raw", - "metadata": { - "tags": [] - }, - "source": [ - "# Load data\n", - "\n", - "datafiles = (gdf.get(\"netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/diag.2016-08-20_00.00.00_subset.nc\"),\n", - " gdf.get(\"netcdf_files/MPAS/FalkoJudt/dyamond_1/30km/x1.655362.grid_subset.nc\") )\n", - "\n", - "primalVarName = 'relhum_200hPa'\n", - "dualVarName = 'vorticity_200hPa'\n", - "central_longitude = 0.0\n", - "\n", - "ds = open_mfdataset(datafiles, decode_times=False)\n", - "primalVar = ds[primalVarName].isel(Time=0).values\n", - "dualVar = ds[dualVarName].isel(Time=0).values\n", - "\n", - "# Fetch lat and lon coordinates for the primal and dual mesh.\n", - "lonCell = ds['lonCell'].values * 180.0 / math.pi\n", - "latCell = ds['latCell'].values * 180.0 / math.pi\n", - "lonCell = ((lonCell - 180.0) % 360.0) - 180.0\n", - "\n", - "lonVertex = ds['lonVertex'].values * 180.0 / math.pi\n", - "latVertex = ds['latVertex'].values * 180.0 / math.pi\n", - "lonVertex = ((lonVertex - 180.0) % 360.0) - 180.0" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "tags": [] - }, - "source": [ - "## Using MPAS's cell connectivity array to plot primal mesh data\n", - "\n", - "In this example we use the MPAS `cellsOnVertex` auxilliary variable, which defines mesh connectivity for the MPAS grid.\n", - "Specifically, this variable tells us the cell IDs for each cell that contains each vertex.\n", - "\n", - "The benefits of this are twofold: 1. We're using the actual mesh description from the MPAS output file; and 2. \n", - "For large grid this is *much* faster than synthesizing the connectivity information as is done\n", - "in the next example\n" - ] - }, - { - "cell_type": "raw", - "metadata": { - "tags": [] - }, - "source": [ - "# Get triangle indices for each vertex in the MPAS file. Note, indexing in MPAS starts from 1, not zero :-(\n", - "#\n", - "tris = ds.cellsOnVertex.values - 1\n", - "\n", - "# The MPAS connectivity array unforunately does not seem to guarantee consistent clockwise winding order, which\n", - "# is required by Datashader (and Matplotlib)\n", - "#\n", - "tris = orderCCW(lonCell,latCell,tris)\n", - "\n", - "# Lastly, we need to \"unzip\" the mesh along a constant line of longitude so that when we project to PCS coordinates\n", - "# cells don't wrap around from east to west. The function below does the job, but it assumes that the \n", - "# central_longitude from the map projection is 0.0. I.e. it will cut the mesh where longitude \n", - "# wraps around from -180.0 to 180.0. We'll need to generalize this\n", - "#\n", - "tris = unzipMesh(lonCell,tris,90.0)\n", - "\n", - "\n", - "# Project verts from geographic to PCS coordinates\n", - "#\n", - "projection = ccrs.Robinson(central_longitude=central_longitude)\n", - "xPCS, yPCS, _ = projection.transform_points(ccrs.PlateCarree(), lonCell, latCell).T\n", - "\n", - "\n", - "trimesh = createHVTriMesh(xPCS,yPCS,tris, primalVar,n_workers=n_workers)" - ] - }, - { - "cell_type": "raw", - "metadata": { - "tags": [] - }, - "source": [ - "# Use precompute so it caches the data internally\n", - "rasterized = hds_rasterize(trimesh, aggregator='mean', precompute=True)\n", - "rasterized.opts(tools=['hover'], colorbar=True, cmap='coolwarm') * gf.coastline(projection=projection)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Synthesizing triangles from points using Delaunay triangulation\n", - "\n", - "In this second example we do not use the triangle connectivity information stored in the MPAS file. Instead we\n", - "use Delaunay triangulation to artifically create a triangle mesh. The benefit of this approach is that we do not\n", - "need the MPAS cellsOnVertex variable if it is not available. Also, since the triangulation algorithm is run on the \n", - "coordinates after they are projected to meters we do not have to worry about wraparound. The downside is that for\n", - "high-resolution data Delaunay triangulation is prohibitively expensive. The highest resolution data set included\n", - "in this notebook (3.75km) will not triangulate in a reasonable amount of time, if at all \n" - ] - }, - { - "cell_type": "raw", - "metadata": {}, - "source": [ - "# Use Delaunay triangulation to generate the triangle connectivity. Note, it's important that the coordinate \n", - "# arrays already be in PCS coordinates (not lat-lon) for the triangulation to perform optimally\n", - "#\n", - "\n", - "from matplotlib.tri import Triangulation\n", - "\n", - "tris = Triangulation(xPCS,yPCS).triangles\n", - "\n", - "trimesh = createHVTriMesh(xPCS,yPCS,tris, primalVar,n_workers=n_workers)" - ] - }, - { - "cell_type": "raw", - "metadata": {}, - "source": [ - "# Use precompute so it caches the data internally\n", - "rasterized = hds_rasterize(trimesh, aggregator='mean', precompute=True)\n", - "rasterized.opts(tools=['hover'], colorbar=True, cmap='coolwarm') * gf.coastline(projection=projection)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Using MPAS's cell connectivity array to plot dual mesh data\n", - "\n", - "In this example we use the MPAS `verticesOnCell` and `nEdgesOnCell` auxilliary variables, which defines mesh connectivity for the\n", - "MPAS dual grid.\n", - "\n", - "As with the first example using the MPAS primal grid, data on the dual grid could be plotted by first\n", - "triangulating them with, for example, Delaunay triangulation. But using grid's native connectivity information \n", - "is faster.\n" - ] - }, - { - "cell_type": "raw", - "metadata": {}, - "source": [ - "verticesOnCell = ds.verticesOnCell.values - 1\n", - "nEdgesOnCell = ds.nEdgesOnCell.values\n", - "\n", - "# For the dual mesh the data are located on triangle centers, which correspond to cell (polygon) vertices. Here\n", - "# we decompose each cell into triangles\n", - "#\n", - "tris = triangulatePoly(verticesOnCell, nEdgesOnCell)\n", - "\n", - "tris = unzipMesh(lonVertex,tris,90.0)\n", - "\n", - "# Project verts from geographic to PCS coordinates\n", - "#\n", - "projection = ccrs.Robinson(central_longitude=central_longitude)\n", - "xPCS, yPCS, _ = projection.transform_points(ccrs.PlateCarree(), lonVertex, latVertex).T\n", - "\n", - "trimesh = createHVTriMesh(xPCS,yPCS,tris, dualVar,n_workers=n_workers)" - ] - }, - { - "cell_type": "raw", - "metadata": {}, - "source": [ - "rasterized = hds_rasterize(trimesh, aggregator='mean', precompute=True)\n", - "rasterized.opts(tools=['hover'], colorbar=True, cmap='coolwarm') * gf.coastline(projection=projection)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Summary\n", - "\n", - "### What's next?\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Resources and references\n", - "\n", - "- []()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.12" - }, - "widgets": { - "application/vnd.jupyter.widget-state+json": { - "state": {}, - "version_major": 2, - "version_minor": 0 - } - } - }, - "nbformat": 4, - "nbformat_minor": 4 -}