Plotting MPAS data

[kmpaul]

MPAS data represents one example of the use of unstructured grids. Plotting data on these grids can sometimes be done effectively with MatPlotLib's tripcolor function, or similar functions that automatically generate triangulated meshes from unstructured data.

However, plotting the MPAS data on the native grid has high value, which requires being able to plot values in "cells" that are non-triangular (hexagonal and pentagonal). The following repository demonstrates how to do this for MPAS atmospheric data in Python:

https://github.com/MiCurry/MPAS-Plotting

There are 3 things about this repository that, I believe, could be improved upon:

1. It is slow. It constructs a matplotlib.collections.PatchCollection object to represent the unstructured MPAS grid (see get_mpas_patches), and the documentation in the code suggests that this can be very slow. I believe it can probably be sped up with Numba, Cython, or any other number of tools. Maybe even smart uses of NumPy.
2. It uses basemap, which is old and no longer supported. I believe the example plotting script in the mpas_plot_pressure.py example can be reworked to use Cartopy. (And then it should probably be added as an example to GeoCAT-examples).
3. It is specific to MPAS, but I think the logic can probably be generalized to any unstructured grid.

I encourage discussion about how to generalize this code, speed it up, and how to modernize it.

CC @NCAR/geocat @NCAR/xdev

[dcherian]

holoviews can handle unstructured meshes (and already uses numba etc) so maybe @jbednar can guide us to the right function. I don't see examples with hexagonal or pentagonal cells though.

[andersy005]

This issue MPAS-Dev/MPAS-Tools#281 may be of interest

[andersy005]

Also, here's a great example of package that is using matplotlib+cartopy to visualize unstructured data following UGRID conventions.

https://psy-maps.readthedocs.io/en/latest/examples/example_ugrid.html

[jbednar]

HoloViz has only triangular unstructured meshes natively and in our current examples, but I would assume that it would be reasonable to convert other grids into trimeshes for display. We would also be happy if someone added native rendering for those other grid types to Datashader (or paid us to hire someone to do that), but that is a good bit of work, so I hope a trimesh will suffice.

[clyne]

@kmpaul what is driving the requirement for rendering on the "native" MPAS grid (as opposed to decomposing the higher sided polygons into triangles)? For reasons of simplicity and performance most 2D visualization methods operate on triangles or rectangles (even quadrilaterals can pose challenges). Even for higher order elements, subdivision in to smaller and smaller triangles generally gives acceptable results for visualization purposes.

[kmpaul]

@clyne If you are suggesting that triangulated visualization of unstructured data is sufficient, then what were you referring to this morning about matplotlib having difficulty with unstructured data? MatPlotLib already has plotting operations that pre-calculate triangulated meshes using Delaunay triangulation.

[clyne]

I would think that triangulation would be a reasonable approach for ESS grids. The one exception being when you want to plot the cell boundaries, which I think could be handled with MPL's line drawing primitives (albeit maybe not efficiently). MPL provides triangle contouring and interpolated triangle coloring, which is sufficient for contour plots (obviously) and pseudo color plots. Quiver plots are possible with triangle meshes, but locating the glyphs - typically done on a regular grid - is not easily doable without a means to sample the tri-mesh arbitrarily. Streamlines do not appear to be supported for tri-meshes.

With all of that said, I think the biggest gap in the ability to plot unstructured data is mapping the mesh from whatever form it is stored in the model output file to something that is usable by MPL. Because conventions like UGRID are not widely used (not used at all?) within the ESS community every model has to be special cased.

Lastly, are MPL's tri* plotting functions compatible with Cartopy? If not, that is another big gap.

[kmpaul]

This is great input! Thanks, @clyne! I appreciate you writing that up. Would you (or maybe @erogluorhan?) be willing to write a few GitHub issues on these needs. It seems like from your comment you have succinctly stated the necessary missing gaps in plotting data on unstructured grids:

• Plotting Cell Boundaries (i.e., plotting the mesh, which the repo above demonstrates how to do, with the stated problems associated with it)
• Quiver plots
• Streamlines
• Associated data wrangling functions to prepare the data for the above operations

And MPL's tri* plotting functions actually are compatible with Cartopy. I have an example from @allibco (Allison Baker) using tripcontour with Cartopy. Not sure how robust it is, though.

If we can split these separate needs out into new issues, then I'd be willing to close this issue.

[jbednar]

The HoloViz tools make it simple to overlay cell boundaries on top of grid cells even if the boundaries and grid are different. To try this out on some particular non-triangular mesh, you'd:

1. Create a GeoViews or HoloViews TriMesh out of your mesh by triangulation, e.g. following the Delaunay example in the docs, and visualize it without showing edges (triangles = datashade(gv.TriMesh.from_vertices(...),aggregator=ds.mean('z')).opts(filled=True)). The main docs are in HoloViews, but you should probably actually use the corresponding component from GeoViews if you want geographic projections to be handled; see e.g. https://geoviews.org/gallery/matplotlib/trimesh_uk.html.
2. Make a GeoViews or HoloViews Graph out of the initial mesh (not the TriMesh), connecting only the points that delineate the original mesh cells (boundaries = datashade(hv.Graph(...))).
3. Overlay the two, on a map or coastline if desired: gf.coastline() * triangles * boundaries

You can use either the Matplotlib or the Bokeh backends to do this, and the advantage compared to using Matplotlib directly is that you can use Datashader to render much larger meshes and render them more quickly, while easily overlaying anything that's in the same space. If you use GeoViews, it's already integrated with Cartopy (e.g. the coastlines come from Cartopy).

[kmpaul]

It would be interesting to see if we could explore whether GeoViews might reduce the complexity of some of the more complex MPL examples in the GeoCAT-examples gallery.

[erogluorhan]

@kmpaul thanks a lot for bringing this issue to the table! It indeed helped pinpoint a number of significant issues much. @clyne let me know if you are willing to put the issues discussed here into a few separate issues. @michaelavs when the time comes, having a couple examples out of these issues implemented in geocat-examples would be nice.

[clyne]

I opened issues as requested:

#24
#25
#26

and also added an issue to produce an MPAS plot example with whatever means possible: @clyne/GeoCAT-examples#221. I think we'll learn a lot once we sit down and try to generate a plot.

[kmpaul]

Then I will close this issue and direct people to those new issues. Thanks, @clyne!

[clyne]

Kevin, I'd keep it open or at least reference this issue from the others. There is a lot of helpful discuss here IMHO.

[bradyrx]

@xylar, might be good to keep an eye on this as well. FYI all, Xylar and I had been discussing this for MPAS-Ocean. I typically make visualizations of MPAS-O with ParaView since it recognizes the unstructured mesh. But that's of course very time-consuming. I'm glad you all are also thinking about a way to do this statically through python.

[xylar]

I created some python code to convert MPAS cells to triangles and plot them either with flat or gouraud shading (data at triangle centers or nodes): MPAS-Dev/MPAS-Tools#319

I also experimented a bit with datashader but didn't get as far as I would have liked. Documentation was sparse and hard to follow.

[jbednar]

For Datashader trimeshes, I'd suggest working with them via the HoloViews TriMesh element, which is a bit more convenient than Datashader's own API, and makes it simple to overlay on maps. See the example at https://examples.pyviz.org/bay_trimesh/bay_trimesh.html and the docs (also sparse and maybe hard to follow) at https://holoviews.org/reference/elements/bokeh/TriMesh.html.

[mgrover1]

I put together an example using CESM data on an unstructured grid... if someone could send over a path to MPAS data, I can modify this and submit a PR https://ncar.github.io/esds/posts/cesm-datashader/

[clyne]

@mgrover1 I'm happy to point you to some MPAS data from Falko Judt that we have permission to redistribute a subset of. It's on glade in /glade/campaign/mmm/wmr/fjudt/projects/dyamond_1. There are actually several runs of the same experiment here, but with different resolutions ranging from 120km to 3.75km. The latter are awesome candidates for Datashader because they are impractical to render with MPL due to their size.

Coincidentally, I have a notebook that @jbednar and colleagues helped me put together for a talk that demonstrates plotting the 3.7km data. It's been on my todo list to clean it up submit it to GeoCAT-examples for a while now, but I haven't been able to find the cycles. Perhaps we can merge our two efforts? One important difference with the script that the Anaconda folks helped me with is that it uses the mesh connectivity information stored in the MPAS output files. Delaunay triangulation is simply too expensive for kilometer scale resolution data (I tried it with the 3.75km data, but after my process ran for 20 minutes and gobbled up over 50GBs of memory I killed it). Using the native connectivity info the entire script runs in about a minute on my laptop. This same issue will of course come up with high res. CAM-SE data, which is why it would be good to have an example the uses the SE connectivity.

I'm starting to ramble. Let me know your thoughts.

Let me know your thoughts.

[mgrover1]

Sounds good - could you share the notebook you used to read in the connectivity information? Or do you have an idea of what variables those would be within MPAS? I also saw this notebook which was shared by @aaronspring which shows a few ways to improve the performance when calculating the triangulation and interface with xarray datasets and data arrays https://gitlab.dkrz.de/m300524/pymistral/-/blob/master/notebooks/xarray_ICON_R2B8.ipynb

[clyne]

Hi @mgrover1 , I've put my notebook under the Datashader directory here: https://github.com/clyne/Notebooks

If there is an easier way to share it let me know. I was going to simply submit a draft PR to the GeoCAT-examples repo, but the repo is currently set up with the expectation that the notebooks generate a static image, not an interactive Geoviews plot. We may need to get some guidance from @erogluorhan on how to support interactive plots in the GeoCAT-examples repo. In the interim it would be great to iterate on the examples. The MPAS data for the notebook are large (several GBs). Currently they are on Glade in /glade/p/cisl/vast/clyne/FalkoJudt/dyamond_1 . So another question for the GeoCAT team will be how to access larger data sets when needed. We have an allocation on Dash for this purpose, but I'm not sure if there is infrastructure in GeoCAT-examples for accessing it.

[erogluorhan]

Hi @clyne, your notebook does not seem to be either here: https://github.com/clyne/Notebooks or anywhere in your github account. Could you please clarify?

I was going to simply submit a draft PR to the GeoCAT-examples repo, but the repo is currently set up with the expectation that the notebooks generate a static image, not an interactive Geoviews plot.

Yes, the GeoCAT-examples gallery currently works with the expectation that the notebooks generate a static image, and that image is shown on the gallery. Even though there is one exception to this workflow, which is the matplotlib Animation example, it still works saving multiple images as a .gif file and then showing that GIF file as a slide-bar animation.

That said, we need to explore ways to show dynamically animated visualizations, which would be output of datashader notebooks. (cc'ed @anissa111 @michaelavs )

So another question for the GeoCAT team will be how to access larger data sets when needed. We have an allocation on Dash for this purpose, but I'm not sure if there is infrastructure in GeoCAT-examples for accessing it.

This is a good question! The gallery currently uses xarray.open_dataset in combination with geocat.datafiles..get(). xarray.open_dataset provides two ways for opening datasets:

1. Local files: This is what GeoCAT-examples gallery currently uses. geocat.datafiles..get() functions downloads the data set files from the geocat-datafiles repo, than xarray.open_dataset opens them.
2. OpenDAP URLs: To my knowledge, this is the only way of opening datasets through URL with xarray.open_dataset. Thus, when we got a DASH dataset of ours ready, my first attempt to access them via xarray.open_dataset was not successful (as expected).

So, we also need to explore options here for opening datasets that are on the internet (e.g. DASH URLs) and have large sizes. netcdf4 provides https access I think, but I do not know how efficient it could be for large data files. (cc'ed @anissa111 @michaelavs )

[jbednar]

the repo is currently set up with the expectation that the notebooks generate a static image, not an interactive Geoviews plot.

GeoViews and HoloViews can be set up to generate static images if that makes the gallery easier to maintain, though in that case you'll want to explain to the user how to turn that off if they do want a dynamically reconfigured image. To get a static output, you can use hv.extension('matplotlib') for the whole notebook, or if you want to keep it using Bokeh you can do hv.output(fig='png'), or to make a particular object a static PNG you can do hv.output(obj, fig='png').

[clyne]

@erogluorhan, sorry about that. The notebook is in

https://github.com/clyne/Notebooks/blob/main/Datashader/MPAS_Datashader_Trimesh.ipynb

Re the data, a couple of thoughts:

1. I believe @michaelavs may have already addressed the use of Dash for large data hosting in another example, but I could be mistaken.
2. The data are available at several resolutions from 240km (a few MBs) to 3.75km (GBs). Perhaps we put a low resolution version of the data in geocat-datafiles that would allow the notebook to handle data like other GeoCAT examples. We could then simply provide links in the comments to the high resolution data on Dask .

[clyne]

the repo is currently set up with the expectation that the notebooks generate a static image, not an interactive Geoviews plot.

GeoViews and HoloViews can be set up to generate static images if that makes the gallery easier to maintain, though in that case you'll want to explain to the user how to turn that off if they do want a dynamically reconfigured image. To get a static output, you can use hv.extension('matplotlib') for the whole notebook, or if you want to keep it using Bokeh you can to hv.output(fig='png'), or to make a particular object a static PNG you can do hv.output(obj, fig='png').

Thanks for the tip, Jim. Because of the focus on using Datashader with large data (the 3.75km data results in about a hundred data points per pixel for a typical image) it would be nice to have an interactive example that lets you drill down into the data set. But that my have to be a goal for a future example :-)

[jbednar]

Definitely! Note that there are a few different levels of interactivity here:

1. static images: No requirement for JavaScript in the browser
2. JavaScript-based interactivity: Interactive Bokeh plots with JavaScript executed in the browser for panning and zooming but no server-side rerendering
3. Python-backed interactivity: Interactive Bokeh plots that issue a call to the server to update the rendering, which is then displayed using JavaScript in the browser

It sounds like your current tooling handles case 1, while e.g. holoviews.org handles case 2, with interactive JS-based plots but no Python server behind them. But servers of type 2 won't actually show any greater resolution when the user zooms a datashader plot; e.g. see the Bokeh plot at https://examples.pyviz.org/bay_trimesh/bay_trimesh.html . The user can zoom in all they like, but the data is never re-rendered. For that you'd need a server of type 3, which is expensive to run if it's public facing, because it needs the server to run some heavy-duty computations for every user. We do that for panel.holoviz.org in some cases (e.g. https://panel.holoviz.org/gallery/demos/nyc_taxi.html), but it's not trivial to make sure such servers stay up and can handle the load of heavy usage, so I'm not sure I'd recommend that people do that just to support a typical website. More commonly you'd put a link to a Binder instance and tell people to visit that if they truly want to explore the data in its full resolution.

[clyne]

Good points! So if we want to demonstrate Datashader's level-of-detail rendering capabilities the only practical way to do this currently is to have the user download the notebook and run it on their local hw. Do i understand that correctly?

[jbednar]

That, or provide a link to a Binder instance where they can run it on remote hardware, or capture a gif of zooming in with the detail updating and show that alongside.

Ok, there's one more option, but it's probably not quite feasible yet. The relatively new package pyiodide lets you run python in the local browser, so in principle people could publish Datashader-based code that then executes on the local browser and dynamically updates. I haven't actually tried it, but it would be cool if it works!

[MiCurry]

Hi all, good to have a lot of discussion on this topic. I dove into GeoCat-Viz after we had a question on the forums and realized we could provide a bit more help. First, I really like some of the utility functions, they are quite helpful. Its nice to have a function like gvutil.set_titles_and_labels(..), as I often forget the commands to add titles to places.

Here are some plots I created using NetCDF4, Cartopy and some GeoCat-Viz utilities that plot the MPAS using tricontour functions and tripcolor that uses @xylar's mesh_to_triangles.py. The first being very similar to this GeoCat-Viz example.

Tricontour

import argparse
import sys

import numpy as np
from netCDF4 import Dataset
import cartopy.crs as ccrs

import matplotlib.pyplot as plt
import matplotlib.tri as tri

from scipy.spatial import Delaunay
import cartopy.feature as cfeature

from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter

import geocat.viz as viz
import geocat.viz.util as gvutil


def normalize_coords(lat, lon):
    lat = np.rad2deg(lat)
    lon = np.rad2deg(lon)
    lon[lon > 180.0] -= 360
    return lat, lon

parser = argparse.ArgumentParser()
parser.add_argument('var', type=str, help="Desired variable to plot")
parser.add_argument('vlevel', type=int, help="Desired vertical level")
parser.add_argument('file',
                    type=str,
                    help='''File you want to plot from''')
args = parser.parse_args()

mesh = Dataset(args.file)
varName = args.var
vlevel = args.vlevel

var = mesh.variables[varName]
field = var[:]
field = field.squeeze()

latCell = mesh.variables['latCell'][:]
lonCell = mesh.variables['lonCell'][:]
zgrid = mesh.variables['zgrid'][:]
latCell, lonCell = normalize_coords(latCell, lonCell)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
# set_global() is currently needed for global plots. It seems Cartopy (or Matplotlib) currently has
# issues with setting the extent. See https://github.com/SciTools/cartopy/issues/1345
ax.set_global()
ax.add_feature(cfeature.LAND, color='silver')

gvutil.set_axes_limits_and_ticks(ax,
                                 xlim=(-180, 180),
                                 ylim=(-90, 90),
                                 xticks=np.linspace(-180,180,13),
                                 yticks=np.linspace(-90,90,7))

gvutil.add_major_minor_ticks(ax)
gvutil.add_lat_lon_ticklabels(ax)

ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol=''))
ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol=''))

gvutil.set_titles_and_labels(ax,
                             maintitle="{0} at {1:.0f} m".format(var.name, zgrid[0,vlevel]),
                             lefttitle=var.long_name,
                             lefttitlefontsize=13,
                             righttitle=var.units,
                             righttitlefontsize=13,
                             xlabel="",
                             ylabel="")

contour = ax.tricontour(lonCell, latCell, field[:,vlevel], 
                        colors='black',
                        levels=7,
                        linewidths=0.5,
                        linestyles='-',
                        transform=ccrs.PlateCarree())

ax.clabel(contour, contour.levels, fontsize=12, fmt="%.0f")

plt.xticks(fontsize=14)
plt.yticks(fontsize=14)

plt.tight_layout()
plt.savefig('mpas_{0}_l{1}.png'.format(var.name, vlevel))

To me, this plot looks strange, and I think that is because I'm plotting the zonal winds at a height level, rather than a pressure height. Although, I didn't do so for the plot above, its quite easy to plot an MPAS field at a specific pressure level, which, perhaps might be a nice feature GeoCat-Comp could use (if it doesn't already have it?).

Tricontourf

import sys

import numpy as np
from netCDF4 import Dataset
import cartopy.crs as ccrs

import matplotlib.pyplot as plt
import matplotlib.tri as tri

from scipy.spatial import Delaunay
import cartopy.feature as cfeature

from cartopy.mpl.gridliner import LongitudeFormatter, LatitudeFormatter

import geocat.viz as viz
import geocat.viz.util as gvutil

def normalize_coords(lat, lon):
    lat = np.rad2deg(lat)
    lon = np.rad2deg(lon)
    lon[lon > 180.0] -= 360
    return lat, lon

print('Argv[1]: ',sys.argv[1])
mesh = Dataset(sys.argv[1])

ter = mesh.variables['ter']
latCell = mesh.variables['latCell'][:]
lonCell = mesh.variables['lonCell'][:]
latCell, lonCell = normalize_coords(latCell, lonCell)

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()

gvutil.set_axes_limits_and_ticks(ax,
                                 xlim=(-180, 180),
                                 ylim=(-90, 90),
                                 xticks=np.linspace(-180,180,13),
                                 yticks=np.linspace(-90,90,7))

gvutil.add_major_minor_ticks(ax, labelsize=10)
gvutil.add_lat_lon_ticklabels(ax)

gvutil.set_titles_and_labels(ax,
                             maintitle="Terrain Height",
                             lefttitle=ter.long_name,
                             lefttitlefontsize=13,
                             righttitle=ter.units,
                             righttitlefontsize=13,
                             xlabel="",
                             ylabel="")

contour = ax.tricontourf(lonCell, latCell, ter[:], 
                         levels=100,
                         cmap='terrain', 
                         transform=ccrs.PlateCarree())

fig.colorbar(contour, orientation="horizontal")

plt.savefig('mpas_terrain.png')

You'll notice the missing data in the upper left and upper right portions of the plot. I belive one can use the xr_add_cyclic_longitudes function to correct that. But I have not tried that, as that function expects an xarray DataArray, and I most often use NetCDF4-Python.

Is there any possibility that xr_add_cyclic_longitudes be refactored to use just Numpy arrays so that it can be called with whatever the user uses for IO (XArray, NetCDF4, etc.). I think as well, I could just use cartopy.util.add_cyclic_point(..)

tripcolor and mesh_to_triangles.py

Lastly, here is a plot that uses @xylar's mesh_to_triangles.py that currently lives in the MPAS-Dev/MPAS-Tools repository and Cartopy.

The mesh_to_triangles.py function definetly seems like the way to go for pcolor plots. While I've only made plots of our smaller resolution meshes, I did test generating the triangles on and it seems to handle those well.

I do think it would be nice if it only operated on Numpy arrays, rather then XArray Datasets. Perhaps it can be refactored to just use Numpy arrays, return the dsTris as Numpy arrays, and then package dsTris into either XArray Datasets in an Xarray Dataset interface function. However, it works fine for now. Here is a very simple example:

import sys
import xarray
import numpy
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation
import cartopy.crs as ccrs

from mesh_to_triangles import mesh_to_triangles

mesh = xarray.open_dataset(sys.argv[1])
meshTris = mesh_to_triangles(mesh)

ter = mesh.ter.values

terTri = ter[meshTris.triCellIndices]

terNode = (ter[meshTris.nodeCellIndices]*meshTris.nodeCellWeights).sum('nInterp')

nTriangles = meshTris.sizes['nTriangles']
tris = numpy.arange(3*nTriangles).reshape(nTriangles, 3)

lonNode = numpy.rad2deg(meshTris.lonNode.values).ravel()
latNode = numpy.rad2deg(meshTris.latNode.values).ravel()
terNode = terNode.values.ravel()

triangles = Triangulation(lonNode, latNode, tris)

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=ccrs.Orthographic(central_longitude=-85))
ax.set_global() # See Cartopy GitHub Issue #1345
ax.tripcolor(triangles, terNode, transform=ccrs.PlateCarree())
ax.coastlines()

plt.savefig('ter_tris_ortho.png')

I'll work on cleaning these up a little bit more and once I do I'll add them to my MPAS-Plotting repo.

[MiCurry]

On the topic of the MPAS-Patches code, I don't really think it is a worthwhile way of plotting data. mesh_to_triangles.py is already a lot faster, and easier to handle for graphics packages.

[clyne]

Thanks much, Miles! We'll keep you in the loop as we move forward with our unstructured grid visualization efforts.

[erogluorhan]

Definitely! Note that there are a few different levels of interactivity here:

1. static images: No requirement for JavaScript in the browser
2. JavaScript-based interactivity: Interactive Bokeh plots with JavaScript executed in the browser for panning and zooming but no server-side rerendering
3. Python-backed interactivity: Interactive Bokeh plots that issue a call to the server to update the rendering, which is then displayed using JavaScript in the browser

It sounds like your current tooling handles case 1, while e.g. holoviews.org handles case 2, with interactive JS-based plots but no Python server behind them. But servers of type 2 won't actually show any greater resolution when the user zooms a datashader plot; e.g. see the Bokeh plot at https://examples.pyviz.org/bay_trimesh/bay_trimesh.html . The user can zoom in all they like, but the data is never re-rendered. For that you'd need a server of type 3, which is expensive to run if it's public facing, because it needs the server to run some heavy-duty computations for every user. We do that for panel.holoviz.org in some cases (e.g. https://panel.holoviz.org/gallery/demos/nyc_taxi.html), but it's not trivial to make sure such servers stay up and can handle the load of heavy usage, so I'm not sure I'd recommend that people do that just to support a typical website. More commonly you'd put a link to a Binder instance and tell people to visit that if they truly want to explore the data in its full resolution.

@jbednar This was great input, and I've been exploring these for a while to include an interactive (hover, zoom, etc) Datashader example under the GeoCAT-examples gallery. You're right with your understanding that GeoCAT-examples is employing (1), while (2) could be a good way to go for interactive plots. When we say (2), I understand the use of the package nbsite; please correct me if I am wrong.

I also agree with your Binder approach that we would likely provide pre-interpreted Jupyter notebooks in the gallery but give a Binder link to the reader if they are interested in further playing with the data on their own.

Finally, I am currently reading through pyiodide as you mentioned back in 2021, but please let me know if you know some success stories or a,ny other experience with it in the meantime.

[erogluorhan]

Hi @MiCurry! In an effort to relate the existing, scattered discussions/issues to our recently awarded Project Raijin (Unstructured grid data analysis and visualization), I walked through this issue again and would like to give you a few updates in hopes that it could be helpful:

• If you haven't done so, you may find it useful to give its website a look. Uxarray is our development of a package that extends Xarray to consume/recognize/output unstructured grids (such as UGRID, SCRIP) as well as enable analysis and visualization functionalities on the native grid. Currently, you may not find as much information on the visualization part of it as the infrastructure of the project though. You are welcome to give us any feedback on our tentative Uxarray API.

• Regarding your plot examples, I think they are not yet in your repo. They would be a good fit to the GeoCAT-examples. Even though we had example(s) for using tricontour and tricontourf(e.g. NCL_station_1.py), they would still provide the readers with plotting MPAS data (we had used them for a different type of data), especially for datasets that are not too big. Please let me know if you would be willing to contribute them to the GeoCAT-examples gallery with a PR, then we can follow both them and your plotting-related improvement/other thoughts over there (otherwise, I am happy to discuss here). Just as an overall response to your few points (e.g. xr_add_cyclic_longitudes), we are preferring to use Xarray over netCDF4 both because of its advanced functionalities and to comply better with the Pangeo stack.

[jbednar]

@jbednar This was great input, and I've been exploring these for a while to include an interactive (hover, zoom, etc) Datashader example under the GeoCAT-examples gallery. You're right with your understanding that GeoCAT-examples is employing (1), while (2) could be a good way to go for interactive plots. When we say (2), I understand the use of the package nbsite; please correct me if I am wrong.

For 2, we use created nbsite to help render Jupyter notebooks into a web page, but nowadays there should be much better tools for integrating notebooks and Sphinx. I'd start with MyST-NB.

I also agree with your Binder approach that we would likely provide pre-interpreted Jupyter notebooks in the gallery but give a Binder link to the reader if they are interested in further playing with the data on their own.

Right; that's useful.

Finally, I am currently reading through pyiodide as you mentioned back in 2021, but please let me know if you know some success stories or a,ny other experience with it in the meantime.

Alas, Bokeh does not work in Pyodide, and it looks like it's going to be a lot of work to get there. So that option's off the table for now, without funding someone to work on replacing Tornado in Bokeh with something that works in Pyodide.

[erogluorhan]

@jbednar Great, thanks! Yes, I am currently trying MyST-NB for GeoCAT-examples and I did have good results with it. I also tried nbsphinx out but couldn't get it configured to run with our structure yet (Then I started playing with MyST-NB). Good to know to avoid Pyodide (at least for now).

[erogluorhan]

All, I made a draft PR at our GeoCAT-examples to add a MPAS visualization with Datashader/Geoviews/Holoviews (the plots themselves were previously implemented by @clyne by adding them as pre-executed Jupyter NBs to GeoCAT-examples is main goal here).

Please feel free to put in your feedback!

[MiCurry]

Hey @erogluorhan thanks for the info! Glad to hear Project Raijin is gaining some speed. I am no longer at NCAR and thus no longer on the MPAS project. If you're wanting an MPAS contact, please contact either Bill Skamaroc or Michael Duda.

[erogluorhan]

Okay, Miles, thanks so much!

[erogluorhan]

Since the last discussion on this, Project Raijin has made significant progress, starting the Python package UXarray for data analysis and visualization of various unstructured meshes, including MPAS. This usage example is just a starting point with UXarray, and we are heavily working on a lot of other visualization methods such as point, triangle, polygon rasterizations, etc. I think we can close this discussion here and keep it going on the UXarray/Raijin end.