feat(vtk): add to_pyvista() method (#1771)

* feat(vtk): add to_pyvista() method
* reformat test_grid.py

autotest/test_export.py

@@ -1265,6 +1265,36 @@ def test_vtk_unstructured(function_tmpdir, example_data_path):
     assert np.allclose(np.ravel(top), top2), "Field data not properly written"
 
 
+@requires_pkg("pyvista")
+def test_vtk_to_pyvista(function_tmpdir, example_data_path):
+    from autotest.test_mp7_cases import Mp7Cases
+
+    case_mf6 = Mp7Cases.mf6(function_tmpdir)
+    success, buff = case_mf6.run_model()
+    assert success, f"MP7 model ({case_mf6.name}) failed"
+
+    gwf = case_mf6.flowmodel
+    plf = PathlineFile(Path(case_mf6.model_ws) / f"{case_mf6.name}.mppth")
+    pls = plf.get_alldata()
+
+    vtk = Vtk(model=gwf, binary=True, smooth=False)
+    assert not any(vtk.to_pyvista())
+
+    vtk.add_model(gwf)
+    grid = vtk.to_pyvista()
+    assert grid.n_cells == gwf.modelgrid.nnodes
+
+    vtk.add_pathline_points(pls)
+    grid, pathlines = vtk.to_pyvista()
+    n_pts = sum([pl.shape[0] for pl in pls])
+    assert pathlines.n_points == n_pts
+    assert pathlines.n_cells == n_pts + len(pls) 
+
+    # uncomment to debug
+    # grid.plot()
+    # pathlines.plot()
+
+
 @pytest.mark.mf6
 @requires_pkg("vtk")
 def test_vtk_vertex(function_tmpdir, example_data_path):

autotest/test_grid.py

@@ -1037,7 +1037,9 @@ def test_structured_thickness(structured_grid):
     thickness = structured_grid.cell_thickness
     assert np.allclose(thickness, 5.0), "thicknesses != 5."
 
-    sat_thick = structured_grid.saturated_thickness(structured_grid.botm + 10.0)
+    sat_thick = structured_grid.saturated_thickness(
+        structured_grid.botm + 10.0
+    )
     assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."
 
     sat_thick = structured_grid.saturated_thickness(structured_grid.botm + 5.0)
@@ -1049,7 +1051,9 @@ def test_structured_thickness(structured_grid):
     sat_thick = structured_grid.saturated_thickness(structured_grid.botm)
     assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
 
-    sat_thick = structured_grid.saturated_thickness(structured_grid.botm - 100.0)
+    sat_thick = structured_grid.saturated_thickness(
+        structured_grid.botm - 100.0
+    )
     assert np.allclose(sat_thick, 0.0), "saturated thicknesses != 0."
 
 
@@ -1082,10 +1086,14 @@ def test_unstructured_thickness(unstructured_grid):
     )
     assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."
 
-    sat_thick = unstructured_grid.saturated_thickness(unstructured_grid.botm + 5.0)
+    sat_thick = unstructured_grid.saturated_thickness(
+        unstructured_grid.botm + 5.0
+    )
     assert np.allclose(sat_thick, thickness), "saturated thicknesses != 5."
 
-    sat_thick = unstructured_grid.saturated_thickness(unstructured_grid.botm + 2.5)
+    sat_thick = unstructured_grid.saturated_thickness(
+        unstructured_grid.botm + 2.5
+    )
     assert np.allclose(sat_thick, 2.5), "saturated thicknesses != 2.5"
 
     sat_thick = unstructured_grid.saturated_thickness(unstructured_grid.botm)

docs/flopy_method_dependencies.md

@@ -32,4 +32,4 @@ Additional dependencies to use optional FloPy helper methods are listed below.
 | `get_sciencebase_xml_metadata()` in `flopy.export.metadata` `acdd` class             | **defusedxml**                                                           |
 | `flopy.utils.geospatial_utils` `GeoSpatialUtil` class                                | **geojson**                                                              |
 | `flopy.utils.geospatial_utils` `GeoSpatialCollection` class                          | **geojson**                                                              |
-| `flopy.export.vtk` `Vtk` class                                                       | **vtk**                                                                  |
+| `flopy.export.vtk` `Vtk` class                                                       | **vtk**, optionally **pyvista**                                          |

flopy/export/vtk.py

@@ -1466,6 +1466,22 @@ def write(self, f: Union[str, os.PathLike], kper=None):
 
             self.pvd.write(pvdfile)
 
+    def to_pyvista(self):
+        """
+        Convert VTK object to PyVista meshes. If the VTK object contains 0
+        or multiple meshes a list of meshes is returned. Otherwise the one
+        mesh is returned alone. PyVista must be installed for this method.
+
+        Returns
+        -------
+        pyvista.DataSet or list of pyvista.DataSet
+            PyVista mesh or list of meshes
+        """
+        pv = import_optional_dependency("pyvista")
+        grids = [self.vtk_grid, self.vtk_polygons, self.vtk_pathlines]
+        meshes = [pv.wrap(grid) for grid in grids if grid is not None]
+        return meshes[0] if len(meshes) == 1 else meshes
+
     def __create_transient_vtk_path(self, path, kper):
         """
         Method to set naming convention for transient vtk file series