Negative values are ignored in a 3D plot

Question

I have to plot a 3d function which has meaningless negative values (they should not appear in the plot). The function which has to be plot is like:

def constraint_function(x, y):
    return min(
        (1800 - 0.3 * x - 0.5 * y) / 0.4,
        (500 - 0.1 * x - 0.08 * y) / 0.12,
        (200 - 0.06 * x - 0.04 * y) / 0.05
    )

I'm calculating the function the following way:

xs = np.linspace(0, 3600, 1000)
ys = np.linspace(0, 3600, 1000)
zs = np.empty(shape=(1000, 1000))
for ix, x in enumerate(xs):
    for iy, y in enumerate(ys):
        zs[ix][iy] = constraint_function(x, y)
xs, ys = np.meshgrid(xs, ys)

The function has valid values mostly in the square [0, 3600]x[0, 3600]. The first approach I had is setting the axis limits to fit my needs:

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')
ax.azim = 20
ax.set_xlim(0, 3500)
ax.set_ylim(0, 3500)
ax.set_zlim(0, 4500)
ax.plot_surface(xs, ys, zs)

plt.show()

Which results in the following plot:

It just ignored the limits and did plot it anyway. The second approach was defining the negative values as np.nan changing the function to be as:

def constraint_function(x, y):
    temp = min(
        (1800 - 0.3 * x - 0.5 * y) / 0.4,
        (500 - 0.1 * x - 0.08 * y) / 0.12,
        (200 - 0.06 * x - 0.04 * y) / 0.05
    )
    return temp if temp >= 0 else np.nan

and setting the alpha of invalid values to zero:

plt.cm.jet.set_bad(alpha=0.0)
ax.azim = 20
ax.set_xlim(0, 3500)
ax.set_ylim(0, 3500)
ax.set_zlim(0, 4500)

ax.plot_surface(xs, ys, zs)

plt.show()

It leaves me with saw-like borders which is also something I don't want to have. Is there a way to get rid of these edges and getting a smooth line when the plot is turning negative?

Answer 1

Technically, you can skew the grid, such that the points of the grid, which would cause a zick-zack pattern are shifted such, that they lie on a line.

This is shown below.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x=np.linspace(-5,5,6)
X,Y = np.meshgrid(x,x)
Z = X+Y

X[Z==-2] = X[Z==-2]+1
Y[Z==-2] = Y[Z==-2]+1
Z[Z==-2] = 0
Z[Z<0] = np.nan

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

ax.set_zlim(0, 12)
ax.plot_surface(X, Y, Z)

plt.show()

The problem would now be to generalize this approach for arbitrary surfaces. It's sure possible but needs a bit of work.

Answer 2

First, your z-value array axes are reversed; it should be zs[iy][ix] not zs[ix][iy]. Because of this your plot is flipped left-for-right.

Second, building your z array by iterating in Python is much slower; you should instead delegate to numpy, like so:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# create axis sample
xs = np.linspace(0, 3600, 1000)
ys = np.linspace(0, 3600, 1000)

# create mesh samples
xxs, yys = np.meshgrid(xs, ys)

# create data
zzs = np.min([
    ((1800 - 0.30 * xxs - 0.50 * yys) / 0.40),
    (( 500 - 0.10 * xxs - 0.08 * yys) / 0.12),
    (( 200 - 0.06 * xxs - 0.04 * yys) / 0.05)
], axis=0)

# clip data which is below 0.0
zzs[zzs < 0.] = np.NaN

NumPy vectorized operations are many times faster.

Third, there is nothing particularly wrong with your code except the sampling resolution is too low; set it higher,

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.azim = 20
ax.set_xlim(0, 3500)
ax.set_ylim(0, 3500)
ax.set_zlim(0, 4500)
ax.plot_surface(xxs, yys, zzs, rcount=200, ccount=200)

plt.show()

produces