pyvista / pyvista-support

[moved] Please head over to the Discussions tab of the PyVista repository
https://github.com/pyvista/pyvista/discussions
59 stars 4 forks source link

Get Individual Contours #505

Closed adamgranthendry closed 2 years ago

adamgranthendry commented 2 years ago

Looking at the Contouring example, is there a way that I can get the list of contours at a specific height/elevation so I can analyze their values individually?

i.e. Say if at a height "z = 10" there were 3 isolines (i.e. "contours"), can I look at each contour separately and analyze their values individually?

My data is "4D". At each (x, y, z) point on my surface (which looks like the Rolling Hills demo data surface), I have unique scalar values. I want to slice my data at different "z"'s and sum the scalar_array values on each separate contour.

adamgranthendry commented 2 years ago

For instance, what if my data was this and I wanted to sum the values on each individual contour:

    mesh = examples.load_random_hills()

    contours = mesh.contour(scalars="Elevation")

    rng = np.random.default_rng(seed=12345)
    mesh.point_arrays["Data"] = rng.integers(low=0, high=10, size=mesh.points.shape[0])
    mesh.set_active_scalars("Data")

    p = pv.Plotter()
    p.add_mesh(mesh, opacity=0.85)
    p.add_mesh(contours, color="white", line_width=5)
    p.show()

image

adamgranthendry commented 2 years ago

@akaszynski and/or @adeak If you could assist me with this, I'd appreciate it! Thank you in advance!

adeak commented 2 years ago

I probably can, but this is my last time responding to such pings this fast :P

I don't believe you can directly query the results of a contour filter like that. If you look at contour.lines you'll see that it's completely made of 2-length line segments. You could check each segment's starting point (contours.lines[1::3]) and check the corresponding value in contours['Elevation'] to see which isovalue a given line segment corresponds to.

But your question might be an XY problem. Aren't you actually looking for slice_along_axis()? Calling mesh.slice_along_axis(axis='z', n=10) gives you a MultiBlock for 10 "isovalues", each having a distinct 'Elevation' scalar value. Then you can do your post-processing on each block in the MultiBlock.

adamgranthendry commented 2 years ago

@adeak You're a genius!! Yes! I am looking for slice_along_axis(). I didn't know that existed!

I hope I can repay your help in full sometime soon! Thank you!!

adamgranthendry commented 2 years ago

@adeak Thank you for your help! Here is a full example for those who may want to do the same thing:

Code

import numpy as np
import pyvista as pv
from matplotlib import pyplot as plt
from pyvista import examples

fig_levels, ax_levels = plt.subplots(2, 5)
fig_comps, ax_comps = plt.subplots(2, 5)

if __name__ == "__main__":
    p = pv.Plotter(shape=(2, 2))

    mesh_linear = examples.load_random_hills()
    mesh_random = mesh_linear.copy()

    # Linear Data
    mesh_linear.point_arrays["DataLinear"] = np.linspace(
        0, 10, mesh_linear.points.shape[0]
    )

    # Random Data
    rng = np.random.default_rng(seed=12345)
    mesh_random.point_arrays["DataRandom"] = rng.integers(
        low=0, high=10, size=mesh_random.points.shape[0]
    )

    # Add Plots
    p.subplot(0, 1)
    mesh_linear.set_active_scalars("DataLinear")
    p.add_mesh(mesh_linear)

    p.subplot(1, 1)
    mesh_random.set_active_scalars("DataRandom")
    p.add_mesh(mesh_random)

    # Create Averages Data Arrays
    mesh_linear.point_arrays["DataLinearContourAverages"] = np.zeros_like(
        mesh_linear["DataLinear"]
    )
    mesh_random.point_arrays["DataRandomContourAverages"] = np.zeros_like(
        mesh_random["DataRandom"]
    )

    # Slice Data Along Z-Azis
    n_slices = 10

    mesh_linear.set_active_scalars("Elevation")
    mesh_random.set_active_scalars("Elevation")

    contours_linear = mesh_linear.slice_along_axis(n=n_slices, axis="z")
    contours_random = mesh_random.slice_along_axis(n=n_slices, axis="z")

    # Average Linear Data for Connected Component in Contour and Plot
    for ndx, contour in enumerate(contours_linear):
        connectivity = contour.connectivity(largest=False)

        labels = connectivity.point_arrays["RegionId"]
        num_labels = len(np.unique(connectivity.point_arrays["RegionId"]))

        for id_ in range(num_labels):
            data = connectivity["DataLinear"][labels == id_]
            connectivity["DataLinearContourAverages"][labels == id_] = np.mean(data)

        p.subplot(0, 0)
        connectivity.set_active_scalars("DataLinearContourAverages")
        p.add_mesh(connectivity)

    # Average Random Data for Connected Component in Contour and Plot
    for ndx, contour in enumerate(contours_random):
        connectivity = contour.connectivity(largest=False)

        labels = connectivity.point_arrays["RegionId"]
        num_labels = len(np.unique(connectivity.point_arrays["RegionId"]))

        for id_ in range(num_labels):
            data = connectivity["DataRandom"][labels == id_]
            connectivity["DataRandomContourAverages"][labels == id_] = np.mean(data)

        p.subplot(1, 0)
        connectivity.set_active_scalars("DataRandomContourAverages")
        p.add_mesh(connectivity)

    mesh_linear.set_active_scalars("DataLinear")
    mesh_random.set_active_scalars("DataRandom")

    p.show()

Output

image