Prioritization of UXarray analysis operators

[erogluorhan]

Hi folks!

This topic could be very much helpful for us to discuss prioritization to implement the analysis operators so that we can ship the most time-sensitively needed unstructured grid functions the earliest. To allow folks to upvote any individual function (which could kind of show us consensus on most commonly wanted tools?), I think it is better to list each function (with or without additional comments, description) in a separate comment.

Below are some functions as an initial list; please feel free to upvote and/or comment any of them, or add your own priorities for other functions.

In addition, please also mention folks/organizations from the community as you see fit to snow-ball the input. I am starting with an initial list from the Raijin collaborators, issue/discussion creators so far, etc:
@andrewgettelman @anissa111 @bjorn-stevens @brianpm @cecilehannay @ChrisBarker-NOAA @clyne @dcherian @falkojudt @iuryt @michaelavs @paullric @rabernat @rajeeja @rljacob @whannah1 @zarzycki

[erogluorhan]

Global and zonal means (conservative and non-conservative)

[andrewgettelman]

Global conservative would be very nice. IMO zonal means do not have to be conservative.

[falkojudt]

I'd add that from the weather side means over arbitrary areas would be useful, either defined by a polygon, or by distance from a lat/lon point.

[erogluorhan]

Cross sections

[andrewgettelman]

If we could get a zonal mean 3D field, then the cross section falls out right?

[paullric]

Cross sections can connect two arbitrary points, not just in the zonal direction.

[erogluorhan]

Regridding to lat-lon

[andrewgettelman]

Many of our tools still require re-gridding, and we may want to re-grid to lower resolution as well for storage.

[whannah1]

@andrewgettelman Shouldn't we rely on other tools for regridding though? NCL's ability to do this was always linked with ESMF, and while I used to use that alot, in retrospect I think it's much better practice to do all regridding tasks (i.e. with NCO) separate from analysis tasks.

[erogluorhan]

Gradients (derivatives)

[paullric]

Also divergence, curl, and Laplacian.

[pilotchute]

As a starting point for thinking about solutions:
A) this overlaps and extends interpolation on unstructured grids, as any interpolation scheme (other than nearest neighbor) can be thought of as creating a field around the point of interest, then pulls the value from the point in the field. Thus a numerical gradient/divergence/curl/laplacian/etc can be done using multiple samples of the field around the point of interest.
B) The math for this is well understood for regular grids, and we would be effectively creating a local regular grid around our points
C) There is at least one flaw in this method, and that is the level of abstraction away from the original dataset. The interpolation scheme will have a large effect on the calculated result so picking an interpolation scheme (or having the option to choose between several) will be an important part of creating good data. I personally would like to suggest Natural Neighbor but so far I've only found one workable voronoi calculation for a sphere, and getting the area overlaps for NN might be non-trivial.

[geoffstanley]

I advocate for using finite volume formulations to define the divergence and the curl, which make use of the divergence theorem and Stokes theorem. The result would then be consistent with finite volume numerical models, and there's no ambiguity about interpolation to produce continuous fields, as @pilotchute noted in point C above. I've been doing this with unstructured grids, at least in 2D, for neutralocean.

In 2D, let

• $m$ be the index to the central grid point where we want to calculate the divergence and curl,
• $\mathcal{N}(m)$ be the set of indices to grid points that neighbour $m$,
• $\Delta_{m,n}$ be the distance between grid points $m$ and $n$,
• $\Delta^\perp_{m,n}$ be the distance of the interface between grid points $m$ and $n$,
• $A_m$ be the area of cell $m$,
• $A_{m,n} = \Delta_{m,n} \Delta^\perp_{m,n}$ be the area associated with the interface between $m$ and $n$.*

Given a vector field $\mathbf{v}$ discretized by elements $v_{m,n}$ being the component of the local vector in the direction from $m$ to $n$, a finite volume formulation declares the discretization of $\nabla \cdot \mathbf{v}$ at cell $m$ to be $(\nabla \cdot \mathbf{v})_m = A_m^{-1} \int \nabla \cdot \mathbf{v} dA$, the average of the divergence of $\mathbf{v}$ over the area of cell $m$. The divergence theorem turns this into $A_m^{-1} \oint \mathbf{v} \cdot \hat{n} dl$ where $\hat{n}$ is the outward pointing unit normal vector and the loop integral goes around the outside of cell $m$. Given that $\mathbf{v} \cdot \hat{n}$ is $v_{m,n}$ and the corresponding $dl$ is $\Delta^\perp_{m,n}$, the loop integral is exactly turned into a sum, and we get $(\nabla \cdot \mathbf{v})_m = A_m^{-1} \sum_{n \in \mathcal{N}(m)} v_{m,n} \Delta^\perp_{m,n}$.

For the 2D curl --- $\nabla_2 \times (u, v) = \partial v / \partial x - \partial u / \partial y$ --- the finite volume formulation proceeds similarly but uses Stokes theorem. The vector field $\mathbf{v}$ is now discretized by $v_{m,n}$ being the component of the local vector in the counter-clockwise direction along the interface between $m$ and $n$. Then $(\nabla_2 \times \mathbf{v})_m = A_m^{-1} \int \nabla_2 \times \mathbf{v} dA = A_m^{-1} \oint \mathbf{v} \cdot dl = A_m^{-1} \sum_{n \in \mathcal{N}(m)} v_{m,n} \Delta^\perp_{m,n}$. Same final formula, but here $v_{m,n}$ is along rather than across the interface, so if $\mathbf{v}$ is the same vector field (e.g. velocity) as used in taking the divergence, then the 2D curl is here operating on a different grid (the vorticity grid, rather than the tracer grid).

The gradient of a scalar field (tracer) $c$ is a vector field discretized by elements $(\nabla c)_{m,n}$ being the component of $\nabla c$ in the direction from $m$ to $n$, naturally given by $(\nabla c)_{m,n} = (c_n - c_m) / \Delta_{m,n}$.

The Laplacian comes along for free given the gradient and divergence definitions above: $(\nabla^2 c)_m = A_m^{-1} \sum_{n \in \mathcal{N}(m)} (c_m - c_n) \Delta^\perp_{m,n} / \Delta_{m,n}$.

The dot product of two vector fields $\mathbf{a}$ and $\mathbf{b}$ was initially less clear to me, but we can also use a finite volume formulation for it by locally converting one of the vector fields into a gradient of a scalar field: $\mathbf{a}_{m,n} = (\chi_n - \chi_m) / \Delta_{m,n}$, where $\chi_m$ and $\chi_n$ are free parameters, and without loss of generality we can choose $\chi_m = 0$. Then $\mathbf{a} \cdot \mathbf{b} = \nabla \cdot (\chi \mathbf{b}) - \chi \nabla \cdot \mathbf{b}$, to which we can apply finite volume formulations for the divergences. The result is $(\mathbf{a} \cdot \mathbf{b})_m = (2 A_m)^{-1} \sum_{n \in \mathcal{N}(m)} a_{m,n} b_{m,n} A_{m,n}$: see Eq (63) and Appendix B of Lang et al (2023).

Most or perhaps all of this is not new, but I'm not sure off-hand where else it's written down. Hope this helps. I'm very interested in UXarray providing better numerical tools for unstructured grids! Perhaps I can eventually use it inside of neutralocean. For now, I've been implementing similar things myself over there.

* The area $A_{m,n}$ should not be calculated separately by spherical area formulas for instance; keeping it as the product of these two distances is required for the gradient and (negative) divergence operators to be adjoint; see Appendix D of Stanley et al (2021).

[iuryt]

I would add integrate and some cumulative_integrate. The latter is not yet available on xarray, but would be awesome to implement.
A nice example to have on the docs that combine derivatives and integrations could be getting some model output to compute geostrophic velocities applying thermal wind to the buoyancy field.

[paullric]

I'm confused about what you mean by "cumulative sum" vs "total sum", and "primitive discrete function". I don't think I've ever heard those terms before.

[iuryt]

Sorry, maybe I should have explained better.
Let's say we have some discrete values for the sine function.
The definite integral (total sum) from 0 to 2pi is 0, but if we calculate the indefinite integral (or inverse derivative/ cumulative sum), we get the primitive function, which is the cosine.

I am showing some example here below. I don't quite remember why I have to multiply by -1 and sum 1 (probably from the boundary condition), but this is the idea.

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

# number of points
n = 100

# sin and cosine functions
t = np.linspace(0,2*np.pi,n)
sin = xr.DataArray(np.sin(t), coords = dict(t=t))
cos = xr.DataArray(np.cos(t), coords = dict(t=t))


# create figure
fig,ax = plt.subplots()

# plot sine
sin.plot(label="sin(t)")
# plot cosine
cos.plot(lw=3,color="0.1",linestyle="--", label="$\int$ sin(t) dt = cos(t)")

# calculate delta t
dt = np.gradient(sin.t)
# cumulative integration
(-np.cumsum(sin*dt)+1).plot(color="r", label="$\sum$ sin(t) $\Delta$t")

ax.legend()

[paullric]

Understood. It seems that you could obtain the cross-section values and then integrate it in 1D using cumsum?

[iuryt]

I think so. So the idea is that the cross-section function returns a regular xarray object that is structured and then we can do whatever we usually do using xarray?

[paullric]

Yes, we can do that. Is there any standard way to represent piecewise functions? I'm envisioning that a cross-section could also be interpreted as a piecewise function, where the discontinuities in each piece occur at face edges.

[adamrher]

plotting example

I don't know if this is unrelated, but I've yet to see someone plot a proper raster plot of an unstructured grid, using python (prove me wrong!). I view this functionality as crucial, especially in the model development process. In NCL, it's as simple as feeding in the vertices of the mesh to the plot resources. See example NCL plot here for a spectral-element run.

[Edited by Rob Jacob to turn into a suggestion]

[adamrher]

ah, thx walter! yes, I do consider that a proper raster. pyNGL is supposed to mimic NCL tools in python, right? I presume it's no longer supported but the tools are still there for everyone to use?

[whannah1]

Yea, it's all in "maintenance mode" in favor of that other thing (geocat?) that is going to also extend existing python libraries (i.e. matplotlib) to add functionality from NCL. However, as I understand it, support for plotting on unstructured grid is not a high priority for them... :(

[clyne]

I don't know if this is unrelated, but I've yet to see someone plot a proper raster plot of an unstructured grid, using python (prove me wrong!). I view this functionality as crucial, especially in the model development process. In NCL, it's as simple as feeding in the vertices of the mesh to the plot resources. See example NCL plot here for a spectral-element run.

See this example that demonstrates how to interactively plot MPAS data with Datashader. Datashader was designed for scalability and is much higher performing on high resolution (e.g. storm resolving) grids than NCL or Matplotlib. The example can easily be extended to plot unstructured grids from other model outputs as well. The current example shows how to plot unstructured data stored in an Xarray object. Later, the example will be revised to demonstrate plotting with UXarray objects and should be simpler to use.

[zarzycki]

Given that PyNGL is in maintenance mode, i this capability may still be useful. I think users will want something analogous to NCL’s fill options (raster vs area vs cell) because there are times when you want an aesthetically pleasing contour field and other times you want the raster native model cells clearly shown.

[adamrher]

See this example that demonstrates how to interactively plot MPAS data with Datashader. Datashader was designed for scalability and is much higher performing on high resolution (e.g. storm resolving) grids than NCL or Matplotlib. The example can easily be extended to plot unstructured grids from other model outputs as well. The current example shows how to plot unstructured data stored in an Xarray object. Later, the example will be revised to demonstrate plotting with UXarray objects and should be simpler to use.

@clyne looking a little closer at the link you provided. Why are the MPAS raster plot examples made up of quadrilateral control volumes? E.g.,

By "proper" raster plot, I am referring to control volumes that reflect the underlying computational grid, which should be made of hexagonal elements for MPAS. Here's what a NCL's "CellFill" options looks like for MPAS (apologies for the low res image):

Can Datashader do this? Or is pyNGL still the only library that can do this?

[clyne]

@adamrher that is a good question. Datashader samples the data at a rate determined by the resolution of the output device. That's how it achieves performance with large data. However, the notebook and images in the example are pre-rendered. When you zoom in on them the resampling is performed on the image, not the data. If you download the notebook and run it on your own machine you'll see that as you zoom in more and detail should be revealed, and the blocky artifacts should go away.

[whannah1]

@clyne I don't think @adamrher is referring to the resampling of the image. The black lines that indicate the coastlines are still well resolved in the image. Adam is saying that the polygons do not appear to be the right shape for this grid, suggesting that datashader is not using the method that we need, which is to show the exact outline of each grid cell.

[falkojudt]

@clyne @whannah1 If you run the notebook and zoom in far enough, you can kinda make out the original hexagonal mesh:
As far as I understand the rendering method uses the Delauney triangles (i.e., the triangle in the figure below), which is not MPAS's computational grid (hexagonal cells in the figure):

So you have to average the grid cells that are in a Delauney triangle to get the value of the polygon that is eventually plotted rasterized = hds_rasterize(trimesh, aggregator='mean', precompute=True)

But of course this isn't the same thing as rendering each grid cell like in @adamrher's example which is probably what we want.

[whannah1]

Thanks for that explanation @falkojudt. I personally find the cell fill method much more useful than the Delauney triangle method when looking for potential grid issues, or just helping me see the scale of the model grid.

[clyne]

@whannah1 the coastlines are not part of the MPAS data. I'm guessing, but don't know for sure, that matplotlib is rendering them under the hood. Re preserving the polygon shape: virtually all low-level raster graphics packages are incapable of rendering polygons other than triangles: n-gons, where n is greater than 3 are decomposed into triangles, and the triangles are rasterized. Polygonal outlines can be drawn as collections of line segments. Datashader, to the best of my knowledge, doesn't draw lines. If drawing the mesh (outlines of the cells that make up the grid) is important, that is something we'll have to explore. Matplotlib can do it, but the performance would be unusable for larger data.

[whannah1]

@clyne The shaded triangle approach is useless to me. In the examples where the shaded cells have been informative I needed to see the exact boundary of the cell to explain something, like the checkerboard noise in E3SM. This requires the mesh information, which is straightforward to obtain. Otherwise, I'll likely need to regrid multiple data sets so that they are on a common grid, and then we're back in rectilinear world where things are easier. Personally, I would prefer if geocat didn't even bother with the delauney triangle appraoch, unless there's some use case that I'm not thinking of where it really adds value.

[ChrisBarker-NOAA]

I thnk perhaps we've been a bit distracted by visualization issues. In my mind, Uxarray should provide the tools that visualization systems required, but not the final results. To that end, the key focus should be on:

• A consistent data model (definition of the grid, etc) that all tools can rely on.
• Core analysis that takes advantage of that data model:
  • interpolation
  • basic calculus: gradients, etc.

Visualization tools can then choose to use the core data model (values on nodes, for instance) and do whatever smoothing, etc that they want, or can use the built-in analysis, like interpolation, and render that directly.

@whannah1 wrote:

"Personally, I would prefer if geocat didn't even bother with the delaunay triangle approach,"

exactly -- if a grid (e.g. traingular mesh) is already defined, then that grid should be what is used. visualization tools may want to provide utilities like delaunay triangulation for data set that don't have a mesh define, but ideally, they would create a mesh that got put into Uxarray format, and then go from there (and that's utility that could be provided by Uxarray as well, but only and an explicit option).

[clyne]

Sorry if I wasn't clear, but there is no need for a Delaunay triangulation if the mesh information is present. Examples 1 and 3 here plot the primary and dual MPAS mesh and do NOT use Delaunay triangulation, but do have to split hexgons, pentagons, etc. into triangles to render them. Though it is not done in this example, the drawing of these triangles can be done in a way so that they look like the original cells.

[brianpm]

Hi all -- @PeterHjortLauritzen was asking about this capability (i.e., visualizing the individual, correctly shaped grid cells from unstructured meshes). Did we decide on a direction? Is this in GeoCAT's hands, or is it in the uxarray roadmap?

[clyne]

Generating constant shaded cells, as in the image provided by @falkojudt above, is not as straightforward as I'd hoped. At the moment Datashader lacks the ability to draw cell-centered data. It expects node-centered data. While it is possible to render a "polygon soup" with each polygon vertex having a constant value, this is pretty inefficient. However, Datashader is continuing to evolve and we can try to work with the Anaconda folks to have this feature added (or try the polygon soup approach and see if the results are satisfactory).

In the interim, it is possible to draw wire frames of the mesh connectivity using VAPOR. An image generated with 3.75km MPAS data is attached. VAPOR is currently a desktop application. However, the VAPOR team will be releasing a Python API later this summer.

[ahijevyc]

Minimum or maximum value within an arbitrary distance of an arbitrary lon, lat point.

This is useful for tropical cyclone tracking and analysis.

For example, what is the maximum wind speed within 300 km of 88.2W, 29N?

[philipc2]

@ahijevyc

Working on this in #616

Do you have a sample grid & data file that could be used for a usage example for this type of tracking & analysis?

[ahijevyc]

Hi @philipc2
I found this example: Tropical Storm Arlene 2023 in the 2023 HWT dataset.

import pandas as pd
from pathlib import Path

idate = pd.to_datetime("20230531T00")
base_path = (
    Path("/glade/campaign/mmm/parc/schwartz/HWT2023/mpas")
    / idate.strftime("%Y%m%d%H")
    / "post/mem_4"
)
grid_path = "/glade/campaign/mmm/parc/schwartz/MPAS/15-3km_mesh/mpas_init/static.nc"
# Paths to Data Variable files
var_names = [
    d.strftime("diag.%Y-%m-%d_%H.%M.%S.nc")
    for d in pd.date_range(
        start=idate + pd.Timedelta(hours=36),
        end=idate + pd.Timedelta(hours=48),
        freq="6H",
    )
]
data_paths = [base_path / name for name in var_names]
# Tropical storm ARLENE
track = [
    (pd.to_datetime("20230601T12"), -87.5, 28.6),
    (pd.to_datetime("20230601T18"), -87.5, 28.6),
    (pd.to_datetime("20230602T00"), -87.6, 28.7),
]
uxds = uxarray.open_mfdataset(
    grid_path, data_paths, concat_dim="Time", combine="nested", use_dual=False
)

I am also working on a use case for the Jan 2024 ESDS workshop. Along the way I also added a bug issue for query_radius and a pull request for returning BallTree via face centers, not node. There's a couple other issues I haven't added yet, like BallTree.query_radius returning only one list of neighboring indices (not a list of lists) and fixing the side effect of reversing the last dimension of coords argument.

[philipc2]

Thank you for this! I've seen the issues and they are being worked on!

Let me know if you would be interested in working with me to get that notebook incorporated as a usage example, I think it would be a great addition.

[ahijevyc]

Thank you for that offer. It's great to see all the progress you've made this past year. I would love to work together on a usage case. Do you have a preferred way to collaborate, like a git repo? Currently I have a Jupyter notebook sitting in my home directory.

It would be great to show off some fixes to the latest issues. For example, I'd like to apply .isel or .where to the UxArray DataArray instead of clipping the grid file and data file individually before creating a UxArray DataArray. You already have an issue about .isel #619, so I could write up one about .where. I suspect .where requires a similar fix. I would like uxarray.DataArray.where to return an uxarray.DataArray not an ordinary xarray.DataArray.

[philipc2]

#616 implements isel for the Grid, UxDataArray and UxDataset

I haven't looked into a where method, but the .subset family of functions might handle the cases you are working with.

I can create a draft PR and a branch with your notebook & we can colaborate through that if you'd like.

[lskopintseva]

Hi,
Xarray allows pretty easy to create xarray object and following netCDF for structured data. Are the good libraries/methods which allow doing this for unstructured data? Given: only timeseries are available with lats and lons. Wish: create unstructured dataset with ability to visualize these data on unstructured grid and being able doing spatial statistic ala xarray. As I see, it boils down in creating the unstructured topology around given lats and lons...

[ChrisBarker-NOAA]

I think you are on the right place -- this project is about adding support for unstructured grids to xarray. However, I think you need to be clear on your terms, so we know what to look for:

a "unstructured grid" means that there IS a grid that defines the relationship between data points, etc -- but that the grid is not "logically rectangular" -- e.g. you can't know what a given grid cell's neighbors are simply by looking at the array indexes. But there is a well defined structure -- most commonly a triangular mesh. and this project is about defining and working with such meshes.

However, general unstructured data, like "only time series are available with lats and lons" is another story -- if you only have lat=lon points, you can't really do a lot with them without defining a relationship between them -- that is if you want to interpolate, you have to decide how to interpolate, which inherently defines a relationship between the points. For instance, one could generate a Delaunay triangulation of the points, and then you would have a traingular mesh that could be used with Uxarray.

I think that job is for other libraries, but providing easy access to it from xarray would makes sense.

Note: you can also do simple things like nearest neighbor interpolation, which may belong in Xarray -- as part of Uxarray? I don't know.

[lskopintseva]

Hi,
Are there any libraries what define this relationship and help building required topology with necessary dependencies? I am seeking for smth general in that sense. Voronoi mesh could be a generalized case. I guess, triangulation is a particular case of it? Then it is more about how to build necessary elements, indexing and avoid infinite definition of the boundaries. I would be very happy, whether such package exists already. And even more happy whether there are easy ways to create netcdfs for unstructured data.

[ChrisBarker-NOAA]

Scipy includes Delaunay triangulation:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html

But it doesn't support "constrained" Delaunay :-(

Thee are a couple wrappers around Shewcuck's traingle code (which, alas, has licencing issues)
(I actually have one too, but haven't published it due to the licencing issues)

Here's one:

https://rufat.be/triangle/

and I think there's a couple more -- I'm not sure how robust any of them are.

The gridded package has some other utilities for building the topology variables, etc. It is supposed to be able to write UGRID compliant netcdf, though it's not complete in that regard (it does write the grid itself though) -- gridded could really use some updating to use XArray, and, hopefully, the Uxarray package.

There is also some triangular mesh code in Matplotlib -- I always wished that had been pulled out into its own package, but the maintainer doesn't want the hassle of maintaining a separate package. But MPL is a pretty common dependency, so you can use parts of it easily enough.

Perhaps a new issue on this project is in order to discuss -- I'd start with specifically outlining what you want, and then the maintainers can decide if it belongs here -- and if so, you can start to fill in the pieces.

[ahijevyc]

That would be great if you could draft a PR and a branch with my notebook. Are the .subset family of functions available yet, or are they still a work in progress?

…

On Tue, Jan 2, 2024 at 4:19 PM Philip Chmielowiec ***@***.***> wrote: #616 <#616> implements isel for the Grid, UxDataArray and UxDataset I haven't looked into a where method, but the .subset family of functions might handle the cases you are working with. I can create a draft PR and a branch with your notebook & we can colaborate through that if you'd like. — Reply to this email directly, view it on GitHub <#46 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEIPI2HZZH6BA4U2TBAKBX3YMSIYPAVCNFSM5V3ZPAJKU5DIOJSWCZC7NNSXTOKENFZWG5LTONUW63SDN5WW2ZLOOQ5TOOJZG44DSMQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[philipc2]

They are in review in #616, feel free to give it a review if you'd like!

I'll get that branch set up today.