File: create-structured-surface.py

package info (click to toggle)
python-pyvista 0.44.1-11
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 159,804 kB
  • sloc: python: 72,164; sh: 118; makefile: 68
file content (156 lines) | stat: -rw-r--r-- 4,840 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
"""
.. _create_structured:

Creating a Structured Surface
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Create a StructuredGrid surface from NumPy arrays
using :class:`pyvista.StructuredGrid`.
"""

from __future__ import annotations

import numpy as np

# sphinx_gallery_thumbnail_number = 2
import pyvista as pv
from pyvista import examples

# %%
# From NumPy Meshgrid
# +++++++++++++++++++
#
# Create a simple meshgrid using NumPy

# Make data
x = np.arange(-10, 10, 0.25)
y = np.arange(-10, 10, 0.25)
x, y = np.meshgrid(x, y)
r = np.sqrt(x**2 + y**2)
z = np.sin(r)

# %%
# Now pass the NumPy meshgrid to PyVista

# Create and plot structured grid
grid = pv.StructuredGrid(x, y, z)
grid.plot()

# %%

# Plot mean curvature as well
grid.plot_curvature(clim=[-1, 1])

# %%
# Generating a structured grid is a one-liner in this module, and the points
# from the resulting surface can be accessed as a NumPy array:

grid.points


# %%
# From XYZ Points
# +++++++++++++++
#
# Quite often, you might be given a set of coordinates (XYZ points) in a simple
# tabular format where there exists some structure such that grid could be
# built between the nodes you have. A great example is found in
# `pyvista-support#16`_ where a structured grid that is rotated from the
# cartesian reference frame is given as just XYZ points. In these cases, all
# that is needed to recover the grid is the dimensions of the grid
# (`nx` by `ny` by `nz`) and that the coordinates are ordered appropriately.
#
# .. _pyvista-support#16: https://github.com/pyvista/pyvista-support/issues/16
#
# For this example, we will create a small dataset and rotate the
# coordinates such that they are not on orthogonal to cartesian reference
# frame.

rng = np.random.default_rng(seed=0)


def make_point_set():
    """Ignore the contents of this function. Just know that it returns an
    n by 3 numpy array of structured coordinates."""
    n, m = 29, 32
    x = np.linspace(-200, 200, num=n) + rng.uniform(-5, 5, size=n)
    y = np.linspace(-200, 200, num=m) + rng.uniform(-5, 5, size=m)
    xx, yy = np.meshgrid(x, y)
    A, b = 100, 100
    zz = A * np.exp(-0.5 * ((xx / b) ** 2.0 + (yy / b) ** 2.0))
    points = np.c_[xx.reshape(-1), yy.reshape(-1), zz.reshape(-1)]
    foo = pv.PolyData(points)
    foo.rotate_z(36.6, inplace=True)
    return foo.points


# Get the points as a 2D NumPy array (N by 3)
points = make_point_set()
points[0:5, :]

# %%
# Now pretend that the (n by 3) NumPy array above are coordinates that you
# have, possibly from a file with three columns of XYZ points.
#
# We simply need to recover the dimensions of the grid that these points make
# and then we can generate a :class:`pyvista.StructuredGrid` mesh.
#
# Let's preview the points to see what we are dealing with:
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 10))
plt.scatter(points[:, 0], points[:, 1], c=points[:, 2])
plt.axis("image")
plt.xlabel("X Coordinate")
plt.ylabel("Y Coordinate")
plt.show()

# %%
# In the figure above, we can see some inherit structure to the points and thus
# we could connect the points as a structured grid. All we need to know are the
# dimensions of the grid present. In this case, we know (because we made this
# dataset) the dimensions are ``[29, 32, 1]``, but you might not know the
# dimensions of your pointset. There are a few ways to figure out the
# dimensionality of structured grid including:
#
# * manually counting the nodes along the edges of the pointset
# * using a technique like principle component analysis to strip the rotation from the dataset and count the unique values along each axis for the new;y projected dataset.

# Once you've figured out your grid's dimensions, simple create the
# :class:`pyvista.StructuredGrid` as follows:

mesh = pv.StructuredGrid()
# Set the coordinates from the numpy array
mesh.points = points
# set the dimensions
mesh.dimensions = [29, 32, 1]

# and then inspect it
mesh.plot(show_edges=True, show_grid=True, cpos="xy")


# %%
# Extending a 2D StructuredGrid to 3D
# +++++++++++++++++++++++++++++++++++
#
# A 2D :class:`pyvista.StructuredGrid` mesh can be extended into a 3D mesh.
# This is highly applicable when wanting to create a terrain following mesh
# in earth science research applications.
#
# For example, we could have a :class:`pyvista.StructuredGrid` of a topography
# surface and extend that surface to a few different levels and connect each
# "level" to create the 3D terrain following mesh.
#
# Let's start with a simple example by extending the wave mesh to 3D
struct = examples.load_structured()
struct.plot(show_edges=True)

# %%
top = struct.points.copy()
bottom = struct.points.copy()
bottom[:, -1] = -10.0  # Wherever you want the plane

vol = pv.StructuredGrid()
vol.points = np.vstack((top, bottom))
vol.dimensions = [*struct.dimensions[0:2], 2]
vol.plot(show_edges=True)