pyvista / pyvista-support

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

Plotting with vertical exaggeration #8

Closed lperozzi closed 5 years ago

lperozzi commented 5 years ago

Description

I would like to plot a point dataset over a surface. I use this code:

p = pv.Plotter(notebook=False)
p.add_mesh(grav, point_size=10.0, render_points_as_spheres=True, cmap='bwr',scalar_bar_args=sargs, clim=[-50,110],stitle='Residuals')
p.add_mesh(surf, cmap='gray')
p.show()

and I obtain this:

image

Now I would like to add a vertical exageration factor of 3 to both surface and point set by using the warp_by_scalar method:

grav = grav.warp_by_scalar(scale_factor=3.0)
surf = surf.warp_by_scalar(scale_factor=3.0)

the resulting plot is:

image

It seems that the points set has not been scaled correctly. How to fix that?

Example Data

gravi.txt surf.txt

banesullivan commented 5 years ago

Hi @lperozzi! This is an excellent question and outlines a common task when dealing with regional scale geo-datasets.

Necessary Imports

import numpy as np
import pyvista as pv
import PVGeo

pv.set_plot_theme('doc')

Load the data

Looks like you have some data files in the UBC formats - PVGeo has readers for these!

# load data files using PVGeo
gravi = PVGeo.ubc.GravObsReader().apply('gravi.txt')
psurf = PVGeo.ubc.TopoReader().apply('surf.txt')

# Filter points of topo to a surface
surf = psurf.delaunay_2d()

Reproduce the first figure:

gargs = dict(point_size=5.0, render_points_as_spheres=True, 
           cmap='bwr', clim=[-50,110], stitle='Residuals')

p = pv.Plotter()
p.add_mesh(gravi, **gargs)
p.add_mesh(surf)
p.camera_position = [2,5,0.5]
p.show()

download

The warp_by_scalar filter

This probably isn't the filter of choice for this particular example. The warp_by_scalar filter distorts the surface by the active scalar on the mesh by default - the current implementation above is warping the grav mesh of scattered points by the 'Residuals' array which likely has a different dimension that the spatial reference in the z-axis.

Directly editing the Z-coordinates

It may be best to directly edit the Z-coordinates of the meshes considering they share the same coordinate datum.

Note I use a factor of 10 instead of 3 to clearly demonstrate that the exaggeration is indeed occurring:

scale_factor = 10.0
# Apply a warp
warped_gravi = gravi.copy()
warped_gravi.points[:, -1] *= scale_factor
warped_surf = surf.copy()
warped_surf.points[:, -1] *= scale_factor

p = pv.Plotter(notebook=1)
p.add_mesh(warped_gravi,**gargs)
p.add_mesh(warped_surf)
p.camera_position = [2,1,0.2]
p.show()

download

Scaling the Plotting scene

This is a somewhat buggy feature we added that has some issues resulting from the VTK backend but it's possible to scale an entire rendering scene and not actually edit the spatial reference of the meshes. This is simply changing the visualization, not the meshes:

p = pv.Plotter()
p.add_mesh(gravi, **gargs)
p.add_mesh(surf)
p.set_scale(zscale=scale_factor)
p.camera_position = [2,1,0.2]
p.show()

download

lperozzi commented 5 years ago

Hi @banesullivan , many thanks for the hints. Actually scale an entire rendering scene is exactly what I need!

banesullivan commented 5 years ago

Just an FYI, there are issues when using the scaled render view: