3-D plot of rectangular surface using matplotlib

Question

I have made a 3D mesh. See below:

And the code below:

import numpy as np
from matplotlib.colors import LogNorm
from mpl_toolkits.mplot3d import Axes3D, axes3d
from matplotlib import cm
import matplotlib.pyplot as plt


def surface(x1, x2, y1, y2, z1, z2, xSum, ySum, zSum, res):
    
    logRes = np.log10(res) / 3
    colVal = plt.get_cmap('jet_r')

    if x1 == min(xSum):
        ax.plot_surface(([x1, x1]), ([y1, y2], [y1, y2]), ([z1, z1], [z2, z2]), linewidth=1,
                        color=colVal(logRes),
                        rstride=1, cstride=1, antialiased=False, shade=False)
    if x2 == max(xSum):
        ax.plot_surface(([x2, x2]), ([y1, y2], [y1, y2]), ([z1, z1], [z2, z2]), linewidth=1,
                        color=colVal(logRes),
                        rstride=1, cstride=1, antialiased=False, shade=False)
    if y1 == min(ySum):
        ax.plot_surface(([x1, x2], [x1, x2]), ([y1, y1]), ([z1, z1], [z2, z2]), linewidth=1,
                        color=colVal(logRes),
                        rstride=1, cstride=1, antialiased=False, shade=False)
    if y2 == max(ySum):
        ax.plot_surface(([x1, x2], [x1, x2]), ([y2, y2]), ([z1, z1], [z2, z2]), linewidth=1,
                        color=colVal(logRes),
                        rstride=1, cstride=1, antialiased=False, shade=False)
    if z1 == max(zSum):
        ax.plot_surface(([x1, x1], [x2, x2]), ([y1, y2], [y1, y2]), ([z1, z1]), linewidth=1,
                        color=colVal(logRes),
                        rstride=1, cstride=1, antialiased=False, shade=False)
    if z2 == min(zSum):
        ax.plot_surface(([x1, x1], [x2, x2]), ([y1, y2], [y1, y2]), ([z2, z2]), linewidth=1,
                        color=colVal(logRes),
                        rstride=1, cstride=1, antialiased=False, shade=False)

def wireframe(x1, x2, y1, y2, z1, z2):
    #
    if x1 == min(xSum):
        ax.plot_wireframe(([x1, x1]), ([y1, y2], [y1, y2]), ([z1, z1], [z2, z2]), linewidth=1)
    if x2 == max(xSum):
        ax.plot_wireframe(([x2, x2]), ([y1, y2], [y1, y2]), ([z1, z1], [z2, z2]), linewidth=1)

    if y1 == min(ySum):
        ax.plot_wireframe(([x1, x2], [x1, x2]), ([y1, y1]), ([z1, z1], [z2, z2]), linewidth=1)
    if y2 == max(ySum):
        ax.plot_wireframe(([x1, x2], [x1, x2]), ([y2, y2]), ([z1, z1], [z2, z2]), linewidth=1)

    if z1 == max(zSum):
        ax.plot_wireframe(([x1, x1], [x2, x2]), ([y1, y2], [y1, y2]), ([z1, z1]), linewidth=1)
    if z2 == min(zSum):
        ax.plot_wireframe(([x1, x1], [x2, x2]), ([y1, y2], [y1, y2]), ([z2, z2]), linewidth=1)

'''
Main program starts from here
'''
resMin = 1
resMax = 1000

xStart = 0
yStart = 0
zStart = 0

x = [100, 50 , 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 50, 100]
y = [100, 50 , 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 50, 100]
z = [10, 20, 30, 40, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]


z = np.multiply(z, -1)

res = 10

xSum = [xStart, np.sum(x)]
ySum = [yStart, np.sum(y)]
zSum = [zStart, np.sum(z)]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

m = cm.ScalarMappable(cmap=cm.jet_r, norm=LogNorm())
m.set_array([resMin, resMax])


cbar = plt.colorbar(m, shrink=0.8, aspect=10)
cbar.set_label('Resistivity', rotation=270)

x1 = xStart
x2 = x1 + x[0]
for i in range(0, len(x)):
    y1 = yStart
    y2 = y1 + y[0]
    for j in range(0, len(y)):
        z1 = zStart
        z2 = (z1 + z[0])
        for k in range(0, len(z)):
            if x1 == min(xSum) or x2 ==max(xSum) or y1 == min(ySum) or y2 ==max(ySum) or z2 == min(zSum) or z1 == max(zSum):
                surface(x1, x2, y1, y2, z1, z2, xSum, ySum, zSum, res)
                wireframe(x1, x2, y1, y2, z1, z2)
            if k < len(z) - 1:
                z1 = z2
                z2 = z1 + z[1 + k]
        if j < len(y)-1:
            y1 = y2
            y2 = y1 + y[1+j]
    if i < len(x)-1:
        x1 = x2
        x2 = x1 + x[1+i]

ax.set_xlabel('X')
# ax.set_xlim3d(0, 500)
ax.set_ylabel('Y')
# ax.set_ylim3d(0, 500)
ax.set_zlabel('Z')
# ax.set_zlim3d(-200, 0)

plt.show()

The purpose is I want to create a 3D model for creating a preparation or showing a result of 3D geophysical modelling & inversion. And the code is very slow when I run it. The higher number of x, y, z the slower response of my code. So it is not effective when I try to rotate it.

My question is, how to make the code effective, light, fast when I run the code?

The code makes a rectangle surface and plots it one by one. I have tried to make it light by plotting the outer side of the cube.

Answer 1

The visualisation simply doesn't make sense. If this is a voxel-like system, you need some cutaways or else you'll only ever see the surface, and the surface is meaningless. Also, my guess is that this demonstration with a constant resistivity \$\rho\$ doesn't reflect the reality of a variable \$\rho\$.

Either this was always broken, or it has since broken due to functional deprecations: your call to colorbar fails unless you pass ax=ax; and all of your calls to plot_surface and plot_wireframe fail because of the type and shape of your z: it must be an actual ndarray, and it must be two-dimensional. I was able to guess my way through adding constructors and manual broadcasts to get this to run again and it shows the same (non-helpful) cube as you did.

You shouldn't loop to construct a mesh.

Don't use the jet colour map. Its perceptual characteristics are very bad. Use a perceptually-uniform colour map like the default Viridis instead.

I suggest deleting everything and writing

import numpy as np
import matplotlib.pyplot as plt

x = [100, 50 , 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 50, 100]
y = [100, 50 , 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 50, 100]
z = [10, 20, 30, 40, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180]
xx, yy, zz = np.meshgrid(x, y, z)
res = np.broadcast_to(10, shape=zz.shape)

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
points = ax.scatter(xx, yy, zz, c=res)
cbar = fig.colorbar(mappable=points, ax=ax, label='Resistivity')
plt.show()

or, since there are only three x-planes and three y-planes, graph three interpolated coloured planes through either x or y.