bug: Cannot open VTK file after creating it

[sebmeng]

Describe the bug
I have a modflow 2005 file that i created using some DEM files used for thefining the lguide layers. Oustide the active domain i have -1e9 values, that's why i can't visualize it when i create the vtk file. I've tried to extract only the active cells, but is not possible, since the latest version flopy.export.vti.Vtk doesn't have the "nanvalue" option anymore.

To Reproduce
Steps to reproduce the behavior:

modelgrid = sim.modelgrid


vtkobj = flopy.export.vtk.Vtk(sim, modelgrid=modelgrid, vertical_exageration=10)


vtkobj.add_array(sim.dis.botm.array, "botm")


hk = sim.lpf.hk.array


ibound = sim.bas6.ibound.array


active_cells = ibound > 0


hk_active = np.full_like(hk, np.nan)
hk_active[active_cells] = hk[active_cells]


vtkobj.add_array(hk_active, "ActiveCells")


vtkobj.write("model.vtk")

Expected behavior
VTK file with active cells only

Screenshots

Desktop (please complete the following information):
OS: Windows 10
-Visual Studio Code v1.76.2
Thanks in advance for the help!

[jlarsen-usgs]

The .add_array() method has a masked_values= parameter that will set array values to nan. Example usage:

vtkobj.add_array(hk_active, "ActiveCells", masked_values=[1e-9, any, other, values])

[wpbonelli]

@sebmeng is it possible to share your model input files and/or the full flopy script to produce them? Then we can try to reproduce the issue.

[sebmeng]

Of course! I zipped my model folder with inputs but it is too heavy for github (1.33 GB). Myabe can i share the folder with only the modflow packages and the script that i used (Build_model_4guideLayers.ipynb), the cell of the ipynb where i used the masked values is this:

# Get the modelgrid from the model
modelgrid = sim.modelgrid

# Create a vtk object
vtkobj = flopy.export.vtk.Vtk(sim, modelgrid=modelgrid, vertical_exageration=10)

vtkobj.add_array(sim.dis.botm.array, "botm",masked_values=[-1.000000E+09])


hk = sim.lpf.hk.array


ibound = sim.bas6.ibound.array


active_cells = ibound > 0


hk_active = np.full_like(hk, np.nan)
hk_active[active_cells] = hk[active_cells]


vtkobj.add_array(hk_active, "ActiveCells", masked_values=[-1.000000E+09])


vtkobj.write("model1.vtk")

if you need the input files to run the firsts cells where the packages are created please let me know and we can find a way to share them. I really appreciate your help!

Regards, Sebastian

[sebmeng]

Uploading Gondebeek.zip…

[sebmeng]

it is still uploading :o

[vcantarella]

Try replacing the last line:

vtkobj.write("model1.vtk")

to:

vtkobj.write("model1")

I remember having a similar issue, I think flopy will decide by itself which file format to use.

[sebmeng]

Uploading Modflow.zip…

[sebmeng]

@wpbonelli did you recieve the ZIP file with the data?
I'm still struggling with this issue

[wpbonelli]

@sebmeng apologies for letting this thread drop, the zip uploads seem to have failed. Could you send the full Python script you are using to construct your model to wbonelli [at] ucar [dot] edu?

[sebmeng]

Yes! I will send it to you rn

Thanks!

[sebmeng]

Does someone had the time to see this? :(

[wpbonelli]

Apologies @sebmeng I received your email, should have a chance to take a look in the next day or two

[wpbonelli]

@sebmeng I tried plotting your grid with pyvista and it worked OK — I don't see any immediate indication of a flopy issue, but let me know if I've missed something. Details in email reply.

[sebmeng]

Thank you for your help Wes!
Just the possibility of exporting only the active domain to VTK would be awesome,

Regards,
Sebastian

[sebmeng]

I found the full workaraound to export only the active cells with pyvista and vtk:

# Load the model
modelname = 'Gondebeek_SS' 
sim = flopy.modflow.Modflow.load(modelname, model_ws='.', exe_name='MODFLOW-NWT.exe')

# Get the modelgrid from the model
modelgrid = sim.modelgrid

# Create a vtk object
vtkobj = flopy.export.vtk.Vtk(sim, modelgrid=modelgrid, vertical_exageration=10,binary=False)

# Flatten the entire 3D array of bottom elevations
botm_flattened = sim.dis.botm.array.flatten()

# Add the flattened array to the VTK object
vtkobj.add_array(botm_flattened, "botm")

# Calculate the top elevations
botm = sim.dis.botm.array
top = np.empty_like(botm)
top[0] = sim.dis.top.array
top[1:] = botm[:-1]

# Flatten the entire 3D array of top elevations
top_flattened = top.flatten()

# Add the flattened array to the VTK object
vtkobj.add_array(top_flattened, "top")

# Flatten the entire 3D array of hydraulic conductivity
hk_flattened = sim.upw.hk.array.flatten()

# Add the flattened array to the VTK object
vtkobj.add_array(hk_flattened, "hk")

# Flatten the entire 3D array of specific storage
ss_flattened = sim.upw.ss.array.flatten()

# Add the flattened array to the VTK object
vtkobj.add_array(ss_flattened, "ss")

# Flatten the entire 3D array of specific yield
sy_flattened = sim.upw.sy.array.flatten()

# Add the flattened array to the VTK object
vtkobj.add_array(sy_flattened, "sy")

# Flatten the entire 3D array of ibound
ibound_flattened = sim.bas6.ibound.array.flatten()

# Add the flattened array to the VTK object
vtkobj.add_array(ibound_flattened, "ibound")

# Write the VTK object to a file
vtkobj.write("model_SS.vtk")

# Set the backend for PyVista
pv.set_jupyter_backend('trame')

# Convert the VTK object to a PyVista mesh
mesh = vtkobj.to_pyvista()

# Apply a threshold filter to keep only cells where ibound equals 1
thresholded_mesh = mesh.threshold(0.5, scalars="ibound")

# List all available point data
print(thresholded_mesh)

# Plot the thresholded mesh with a specific property
# Replace 'property_name' with the name of the property you want to visualize
thresholded_mesh.plot(scalars='hk', cpos=[mesh.points[0], mesh.points[1], (0, 0, 0)], show_edges=False)

# Save the thresholded mesh to a VTK file
thresholded_mesh.save("Gondebeek_SS_act.vtk")

with this code, i get this:

The VTK file exported is also visible in ParaView:

I appreciate your help @wpbonelli and @vcantarella! Your comments and help were crucial for me to get this result.
Maybe in the future would be amazing to be able to export the active grid directly from the vtk module of FloPy.

Thanks!!