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

Terrain Following Mesh from numpy mgrid? #369

Open msadowski opened 3 years ago

msadowski commented 3 years ago

Description

Hi!

Many thanks for the great package! I'm wondering if you would have some hints on how to create a terrain model from an mgrid? I've looked at the tutorial on the terrain following mesh but for the life of me I couldn't find a way to modify my data to a format that would allow me to follow the tutorial.

Here is how my data looks like when plotted using matplotlib:

image

And here is my starter code:

import numpy as np
import pyvista as pv
import matplotlib.pyplot as plt

pc = np.load('test_file.npy')

plt.figure(figsize=(10,10))
plt.imshow(pc, origin="lower")

Would you have any tips on how to convert my datatype to DEM or find any other way to create a Terrain Following Mesh?

Example Data

Please find the zip file attached test_file.zip

msadowski commented 3 years ago

I've gave it another go, trying to separate all my data into separate axes, and then use it as vertices:

vertices = np.vstack(pc.reshape(3,-1).T)
print(vertices)
cloud = pv.PolyData(vertices)
cloud.plot(eye_dome_lighting=True)

however the results are not amusing:

image

Then I decided to take a step back, and try to work with my raw XYZ data, getting this result:

image

which looks like the data is too sparse to interpolate it.

Another thing I've tried is using delaunay_3d on my data, producing almost OK result, except the fact that the depth is not correct

image

here is interpolated depth plotted in matplotlib:

image

Any tips on what else I could try out would be highly appreciated!

akaszynski commented 3 years ago

Looks like you're going to want to have a separate scale on each axis. See https://docs.pyvista.org/plotting/plotting.html#pyvista.BasePlotter.set_scale

msadowski commented 3 years ago

Many thanks for the hint Alex! This works! (and hints to me that I probably made a mistake in unit conversion at some point.

I think the last problem I'm facing is setting custom color map. I've tried hackig it as follows using pyplot 'jet' color:

cloud = pv.PolyData(points)
surf = cloud.delaunay_3d(alpha=0.9).extract_surface().triangulate().clean()

boring_cmap = plt.cm.get_cmap("jet", 20)

plotter = pv.Plotter()
surf.point_arrays['scalars'] = np.linspace(np.nanmin(di_copy), np.nanmax(di_copy), surf.n_points)

plotter.add_mesh(surf, lighting=True, show_edges=False, cmap=boring_cmap, scalars='scalars')
plotter.set_scale(xscale=1.0, yscale=1.0, zscale=3.28084)

plotter.show()

and getting this result:

image

Would you have any tips on how to make the coloring work along the z axis? Thanks!

akaszynski commented 3 years ago

Issue seems to be:

surf.point_arrays['scalars'] = np.linspace(np.nanmin(di_copy), np.nanmax(di_copy), surf.n_points)

You're going to want to instead use:

surf.point_arrays['scalars'] = surf.points[:, 2]

That way you're simply using the "Z" coordinates of the points.

For clarity in your plot, you may want to use the following instead:

surf.point_arrays['Z Coordinates'] = surf.points[:, 2]

And then:

plotter.add_mesh(surf, lighting=True, show_edges=False, cmap=boring_cmap, scalars='Z Coordinates')