Resolve example issues with pyvista 0.42.2

bjlittle commented 1 year ago

🐛 Bug Report

Only noticed as part of manual post-testing of geovista 0.4.0 in a new environment with the latest version of pyvista 0.42.2, that there are issues with some projection examples, namely:

`from_2d__orca_moll`

`from_unstructured__fesom_fouc`

`from_unstructured__icon_eqc`

`from_unstructured__icosahedral_poly`

`from_unstructured__tri_hammer`

How to Reproduce

geovista examples --run from_2d__orca_moll
geovista examples --run from_unstructured__fesom_fouc
geovista examples --run from_unstructured__icon_eqc
geovista examples --run from_unstructured__icosahedral_poly
geovista examples --run from_unstructured__tri_hammer

Environment

Please share the report generated by either of the following CLI commands:

scooby --report geovista

or

python -c "import geovista; print(geovista.Report())"

----------------------------------------------------------------------------------------
  Date: Tue Sep 19 00:36:48 2023 BST

                 OS : Linux
             CPU(s) : 4
            Machine : x86_64
       Architecture : 64bit
        Environment : Python
         GPU Vendor : VMware, Inc.
       GPU Renderer : llvmpipe (LLVM 7.0, 256 bits)
        GPU Version : 3.3 (Core Profile) Mesa 18.3.4

  Python 3.11.5 | packaged by conda-forge | (main, Aug 27 2023, 03:34:09) [GCC 12.3.0]

            cartopy : 0.22.0
              click : 8.1.7
click-default-group : 1.2.4
            cmocean : v3.0.3
           colorcet : 3.0.1
           geovista : 0.4.0
         matplotlib : 3.7.2
            netcdf4 : 1.6.4
              numpy : 1.26.0
       platformdirs : 3.10.0
              pooch : v1.7.0
           pykdtree : 1.3.7.post0
             pyproj : 3.6.0
            pyvista : 0.42.2
             scooby : 0.7.2
                vtk : 9.2.6
        fastparquet : 2023.8.0
             pandas : 2.1.0
----------------------------------------------------------------------------------------

bjlittle commented 1 year ago

Temporarily pinned back geovista from using pyvista 0.42.2 for now, see https://github.com/conda-forge/geovista-feedstock/pull/13

tkoyama010 commented 1 year ago

I have confirmed that this bug can be resolved by reverting https://github.com/pyvista/pyvista/pull/4853. However, it needs to be added because of another bug. I am trying to reproduce the occurrence of the error in a simple PyVista example and to find an appropriate solution to this problem.

bjlittle commented 1 year ago

Thanks @tkoyama010

Brilliant detective work!

Let me know if there is anything I can do this end to help.

I'm pretty keen to get the examples running within CI with image testing to automate and catch these sort of issues earlier 👍

tkoyama010 commented 1 year ago

I'm pretty keen to get the examples running within CI with image testing to automate and catch these sorts of issues earlier 👍

Quite right. https://github.com/pyvista/pyvista/pull/4853 bug could have been prevented before release if GeoVista had implemented it, since PyVista tests GeoVista with CI. (So #429 is important.)

I also think that this issue is a good example of how to report to the PyVista project that was identified in GeoVista.

I have now used Python's traceback module to see which part of GeoVista is calling in the problem line.

--- a/pyvista/core/utilities/geometric_objects.py
+++ b/pyvista/core/utilities/geometric_objects.py
@@ -29,6 +29,8 @@ from .arrays import _coerce_pointslike_arg
 from .helpers import wrap
 from .misc import check_valid_vector

+import traceback
+
 NORMALS = {
     'x': [1, 0, 0],
     'y': [0, 1, 0],
@@ -57,6 +59,7 @@ def translate(surf, center=(0.0, 0.0, 0.0), direction=(1.0, 0.0, 0.0)):
     """
     normx = np.array(direction) / np.linalg.norm(direction)
     # assume temporary normy to calculate normz
+    traceback.print_stack()
     norm_y_temp = [0.0, 1.0, 0.0]
     normz = np.cross(normx, norm_y_temp)
     if np.array_equal(normz, (0.0, 0.0, 0.0)):

It turns out that the problematic part is being called in the Plane function.

  File "geovista/src/geovista/examples/from_2d__orca_moll.py", line 48, in <module>
    main()
  File "geovista/src/geovista/examples/from_2d__orca_moll.py", line 43, in main
    plotter.add_mesh(mesh)
  File "pyvista/venv/Lib/site-packages/geovista/geoplotter.py", line 591, in add_mesh
    sliced_mesh = slice_mesh(mesh, rtol=rtol, atol=atol)
  File "pyvista/venv/Lib/site-packages/geovista/core.py", line 898, in slice_mesh
    else slice_cells(mesh, antimeridian=True, rtol=rtol, atol=atol)
  File "pyvista/venv/Lib/site-packages/geovista/core.py", line 629, in slice_cells
    remeshed, remeshed_west, remeshed_east = remesh(
  File "pyvista/venv/Lib/site-packages/geovista/filters.py", line 114, in remesh
    poly1 = pv.Plane(
  File "pyvista/pyvista/core/utilities/geometric_objects.py", line 462, in Plane
    translate(surf, center, direction)
  File "pyvista/pyvista/core/utilities/geometric_objects.py", line 62, in translate
    traceback.print_stack()

Next, I plan to check what parameters are used to execute the Palne function.

tkoyama010 commented 1 year ago

@bjlittle Hi. I have a question. I know that I can modify the Plane function in PyVista to work properly with the following lines, but first I need to understand the remesh function it is used with. Can you please elaborate on what the purpose of this function is and why we need the Plane function for that? https://github.com/bjlittle/geovista/blob/933fe2679aa6059d6c1f79ba4f187345460e8217/src/geovista/filters.py#L114-L121

bjlittle commented 1 year ago

@tkoyama010 The geovista.filters.remesh function is fundamental to the projection pipeline within geovista.

Essentially, it is provided with only the cells from a parent mesh that require to be "ripped" as they cross a meridian of choice i.e., the antimeridian on the Earth.

The remesh function uses the pyvista.Plane to represent the plane that passes through the meridian that "rips" or "cuts" the mesh. The cells of the mesh traversing the meridian are sliced through the point of intersection with the meridian plane, and triangulated, thus breaking the cells apart and their connectivity.

Does this help?

If I get time, I'll also provide you with a visualisation of what's happening here, which will be way, way clearer than words.

tkoyama010 commented 1 year ago

@bjlittle Thanks, I think we will be able to communicate better if we explain to each other in a figure as you say. I will be prepared to explain in a diagram how the Plane function works with respect to the direction argument.

tkoyama010 commented 1 year ago

We are trying to solve this problem in the https://github.com/pyvista/pyvista/pull/5188, but it will take some more time to consider. We have therefore given up on including it in the next version.

bjlittle commented 1 year ago

@tkoyama010 Thanks for letting me know, much appreciated.

I'll see if I can generate some visualization to help clarify what's happening to give some context.

banesullivan commented 11 months ago

Tracking the issue upstream in https://github.com/pyvista/pyvista/issues/5405

user27182 commented 11 months ago

I believe adding something like poly1.rotate_y(90, inplace=True) after instantiating the plane in remesh (before the call to rotate_z) may fix this.

https://github.com/bjlittle/geovista/blob/f7c84b31185f23d8cd9cb741eb5cc343b3fa5ba8/src/geovista/filters.py#L119

More generally, I believe the fix for this is to rotate the plane 90 degrees wherever direction=(0,1,0) is used internally in GeoVista.

This will effectively 'undo' the fix from pyvista/pyvista#4853 and restore the plane orientation behaviour for the y-direction case.

This fix is similar to what is proposed in pyvista/pyvista#5435. However, implementing the changes in GeoVista instead of PyVista would remove the need to re-introduce a bug into PyVista's methods.