error using annotate_contour with exodus II data

[ahzeeshan]

Bug report

Bug summary
Cannot plot contours using data from exodus II file.

Code for reproduction

Exodus file is gold.e from http://use.yt/upload/a113f8e3

ds = yt.load("/Users/zeeshan/Downloads/ExodusII/gold.e")
print(ds.field_list)
slc = yt.SlicePlot(ds, "z", ('all','forced'))
slc.annotate_contour(("all", "forced"))

Actual outcome

[('all', 'forced'), ('connect1', 'forced')]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/miniconda3/envs/py39/lib/python3.9/site-packages/yt/visualization/plot_window.py:1297, in PWViewerMPL.run_callbacks(self)
   1296 try:
-> 1297     callback(cbw)
   1298 except (NotImplementedError, YTDataTypeUnsupported):

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/yt/visualization/plot_modifications.py:955, in ContourCallback.__call__(self, plot)
    954     triangulation = Triangulation(x, y)
--> 955     zi = LinearTriInterpolator(triangulation, z)(xi, yi)
    956 elif plot._type_name == "OffAxisProjection":

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/matplotlib/tri/_triinterpolate.py:258, in LinearTriInterpolator.__init__(self, triangulation, z, trifinder)
    257 def __init__(self, triangulation, z, trifinder=None):
--> 258     super().__init__(triangulation, z, trifinder)
    260     # Store plane coefficients for fast interpolation calculations.

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/matplotlib/tri/_triinterpolate.py:39, in TriInterpolator.__init__(self, triangulation, z, trifinder)
     38 if self._z.shape != self._triangulation.x.shape:
---> 39     raise ValueError("z array must have same length as triangulation x"
     40                      " and y arrays")
     42 _api.check_isinstance((TriFinder, None), trifinder=trifinder)

ValueError: z array must have same length as triangulation x and y arrays

The above exception was the direct cause of the following exception:

YTPlotCallbackError                       Traceback (most recent call last)
File ~/miniconda3/envs/py39/lib/python3.9/site-packages/IPython/core/formatters.py:344, in BaseFormatter.__call__(self, obj)
    342     method = get_real_method(obj, self.print_method)
    343     if method is not None:
--> 344         return method()
    345     return None
    346 else:

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/yt/visualization/_commons.py:169, in validate_plot.<locals>.newfunc(self, *args, **kwargs)
    165     self._recreate_profile()
    166 if not self._plot_valid:
    167     # it is the responsibility of _setup_plots to
    168     # call plot.run_callbacks()
--> 169     self._setup_plots()
    170 retv = f(self, *args, **kwargs)
    171 return retv

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/yt/visualization/plot_window.py:1240, in PWViewerMPL._setup_plots(self)
   1237         self.plots[f]._toggle_axes(draw_axes, draw_frame)
   1239 self._set_font_properties()
-> 1240 self.run_callbacks()
   1242 if self._flip_horizontal or self._flip_vertical:
   1243     # some callbacks (e.g., streamlines) fail when applied to a
   1244     # flipped axis, so flip only at the end.
   1245     for f in field_list:

File ~/miniconda3/envs/py39/lib/python3.9/site-packages/yt/visualization/plot_window.py:1301, in PWViewerMPL.run_callbacks(self)
   1299         raise
   1300     except Exception as e:
-> 1301         raise YTPlotCallbackError(callback._type_name) from e
   1302 for key in self.frb.keys():
   1303     if key not in keys:

YTPlotCallbackError: annotate_contour callback failed

Expected outcome
Contours on top of the color plot.

Version Information

• Operating System: MacOS Ventura
• Python Version: 3.9.18
• yt version: 4.1.4 installed using conda

[neutrinoceros]

hi @ahzeeshan , thank you for reporting this

from what I can grasp, what's happening internally is that the "all", "forced" field has a shape (1000, 4), where x/y image coordinate fields have (1000,)

yt/yt/visualization/plot_modifications.py

Line 1000 in 8a43695

z = data[self.field][wI]

I am not familiar with exodusII data or unstructured mesh, so I don't know off-hand if this is to be expected. Is it ?

[ahzeeshan]

Thanks @neutrinoceros. This shape mismatch error seems related to #1916 and #1918. However, I am not sure how to modify my code. I am also not very familiar with exodusII data but my guess is the 4 comes from the 4 nodes corresponding to each element of the QUAD4 mesh.

[neutrinoceros]

thanks for digging up those relevant issues/PRs, this is helpful.
I think this would benefit from @matthewturk's expertise with unstructured meshes ?

[matthewturk]

@neutrinoceros @ahzeeshan Thanks for the ping! I absolutely would like to look at this. I really appreciate the relevant PRs and the thoughtful method of supplying relevant test data. I will attempt to dig in this afternoon or tomorrow!

[matthewturk]

Sorry for taking so long, but I've figure out the first bit that's an issue. The current setup is feeding in the vertex data, which is of shape (N, 4) as the "z" array. The x and y are of shape N, however. This is because the slice functionality for unstructured mesh data utilizes all four values to compute the local points (i.e., doing intra-cell computation). I believe that this is tractable by utilizing a different pixelizer on the yt side. I am going to continue to investigate.

[matthewturk]

After digging in, I'm not entirely certain I know the specific right thing to do. I had forgotten the specifics of how this worked, but now that I am reminding myself, I see the issue. It does not pixelize -- instead, it feeds the raw data locations in to matplotlib, and then uses that to generate the results. The main upshot of this is that for cell-centered data that occupies more than one pixel per voxel, it prevents really gross "let's trace the boxes with our contours" situations. It's slower than our pixelizer, but that is what it is.

So, some options:

1. We could pixelize the mesh elements onto the plane and run the contours over that.
2. We could unravel the x, y positions of the nodal data and feed them in. (I am leaning toward this one.)
3. We could use mean of the nodal elements as the input to the triangulator. (This is definitely the easiest, and probably would give reasonable results since we're only looking at contours rather than raw data values.)

I'll give it some more thought and experiment a bit.

[matthewturk]

I should note, item 2 on that list is made more complicated in 3D.

[ahzeeshan]

Can confirm I used option 3 as a temporary solution and it gave satisfactory results. However, the element value may not always be the mean of the node values. The form of interpolation depends on the shape functions used.