fatiando / harmonica

Forward modeling, inversion, and processing gravity and magnetic data
https://www.fatiando.org/harmonica
BSD 3-Clause "New" or "Revised" License
211 stars 69 forks source link

pyvista plot zero volumn #392

Closed brianchen123 closed 1 year ago

brianchen123 commented 1 year ago

Description of the problem:

Full code that generated the error

import numpy as np
import verde as vd
import harmonica as hm

# Build coordinates
# -----------------
# Build horizontal coordinates
easting = np.linspace(0, 300, 301)
northing = np.linspace(-100, 100, 101)
easting, northing = np.meshgrid(easting, northing)

# Define the top and bottom surfaces of the thin block
upward = easting * np.sin(-30) - 5
surface_top = easting * np.sin(np.radians(-30)) 
#shift the top surface by 50 elements (100*1 =100 m)
for i, row in enumerate(surface_top):
    surface_top[i, :] = np.concatenate((np.zeros(100), row[:-100]))
#shift the bottom surface by 30 elements (80*1 =80 m)
surface_bottom = easting * np.sin(np.radians(-30)) 
for i, row in enumerate(surface_bottom):
    surface_bottom[i, :] = np.concatenate((np.zeros(80), row[:-80]))

# Define density values (kg/m3)
# ---------------------
density_a, density_b, density_c = 2500, 2000, 2700
# Create density array for the prism layer
density_bottom = density_c * np.ones_like(surface_top)
density_block = density_b * np.ones_like(surface_top)
density_top = density_a * np.ones_like(surface_top)

# Create prism layer model
# ------------------------
model_block = hm.prism_layer(
    coordinates=(easting, northing),
    surface=surface_top,
    reference=surface_bottom,
    properties={"density": density_block},
)
model_bottom = hm.prism_layer(
    coordinates=(easting, northing),
    surface=surface_bottom ,
    reference= -200 * np.ones_like(surface_bottom),
    properties={"density": density_bottom},
)
model_top = hm.prism_layer(
    coordinates=(easting, northing),
    surface= 0 * np.ones_like(surface_top) ,
    reference= surface_top ,
    properties={"density": density_top},
)

# Plot using pyvista
# ------------------
# Try to import pyvista
try:
    import pyvista as pv
except ImportError:
    pv = None

# Plot only if pyvista is installed
if pv is not None:
    plotter = pv.Plotter()
    plotter.add_mesh(model_block.prism_layer.to_pyvista(), scalars="density")
    plotter.add_mesh(model_top.prism_layer.to_pyvista(), scalars="density")
    plotter.add_mesh(model_bottom.prism_layer.to_pyvista(), scalars="density")
    plotter.show_axes()
    plotter.show() 

Full error message

PASTE ERROR MESSAGE HERE

System information

Output of conda list

PASTE OUTPUT OF CONDA LIST HERE