42

I m trying to plot simultaneously a plane and some points in 3D with Matplotlib. I have no errors just the point will not appear. I can plot at different times some points and planes but never at same time. The part of the code looks like :

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

point  = np.array([1, 2, 3])
normal = np.array([1, 1, 2])

point2 = np.array([10, 50, 50])

# a plane is a*x+b*y+c*z+d=0
# [a,b,c] is the normal. Thus, we have to calculate
# d and we're set
d = -point.dot(normal)

# create x,y
xx, yy = np.meshgrid(range(10), range(10))

# calculate corresponding z
z = (-normal[0] * xx - normal[1] * yy - d) * 1. /normal[2]

# plot the surface
plt3d = plt.figure().gca(projection='3d')
plt3d.plot_surface(xx, yy, z, alpha=0.2)


#and i would like to plot this point : 
ax.scatter(point2[0] , point2[1] , point2[2],  color='green')

plt.show()
Trenton McKinney
  • 56,955
  • 33
  • 144
  • 158
user3601754
  • 3,792
  • 11
  • 43
  • 77
  • @AndreyRubshtein how is that related? The question your link to is about `matlab`, this is about `matplotlib` – tmdavison Aug 16 '17 at 11:12
  • Hey @tom, it has an explanation about the math formula, the one that is called 'Calculate corresponding z' in the original question – Andrey Rubshtein Aug 16 '17 at 12:41
  • Ok sure, although that's not really what this question and answers are about – tmdavison Aug 16 '17 at 13:12
  • 1
    Just to mention it: this method will fail for normal vectors which have a z-component equal to zero! So, for example, `normal = np.array([1,0,0])` will throw `RuntimeWarning: divide by zero encountered in true_divide` and no plane is drawn! – Douglas James Bock Jul 01 '21 at 13:26

3 Answers3

28

Just to add to @suever's answer, you there's no reason why you can't create the Axes and then plot both the surface and the scatter points on it. Then there's no need to use ax.hold():

# Create the figure
fig = plt.figure()

# Add an axes
ax = fig.add_subplot(111,projection='3d')

# plot the surface
ax.plot_surface(xx, yy, z, alpha=0.2)

# and plot the point 
ax.scatter(point2[0] , point2[1] , point2[2],  color='green')
tmdavison
  • 64,360
  • 12
  • 187
  • 165
21

You will need to tell the axes that you want new plots to add to the current plots on the axes rather than overwriting them. To do this, you will need to use axes.hold(True)

# plot the surface
plt3d = plt.figure().gca(projection='3d')
plt3d.plot_surface(xx, yy, z, alpha=0.2)

# Ensure that the next plot doesn't overwrite the first plot
ax = plt.gca()
ax.hold(True)

ax.scatter(points2[0], point2[1], point2[2], color='green')

enter image description here

UPDATE

As @tcaswell pointed out in the comments, they are considering discontinuing support for hold. As a result, a better approach may be to use the axes directly to add more plots as in @tom's answer.

Community
  • 1
  • 1
Suever
  • 64,497
  • 14
  • 82
  • 101
  • Do you actually use `hold` day-to-day? We have been talking about deprecating it and I would like to talk to someone who actually uses it.... – tacaswell Mar 18 '16 at 02:08
  • @tcaswell I do tend to, but I think it's mainly due to the habits I developed from using MATLAB. That being said, I actually do like the alternate approach mentioned by [tom](http://stackoverflow.com/a/36061408/670206) but I never think to use it! – Suever Mar 18 '16 at 02:38
  • Also, I think most things default to `hold==True`. If you find a use case where you _can not_ just remove it please make an issue on GH – tacaswell Mar 18 '16 at 02:40
  • AttributeError: 'Axes3D' object has no attribute 'hold' – Mona Jalal Jan 26 '21 at 21:09
8

Just elaborating a little more on the math part (and how it works), may be useful for someone, the equation of a plane with unit normal vector n containing a point a is shown as follows:

enter image description here

so that the equation of the plane here is x + y + 2*z = 9 and the following code can simply be used to plot the given plane:

# create the figure
fig = plt.figure()

# add axes
ax = fig.add_subplot(111,projection='3d')

xx, yy = np.meshgrid(range(10), range(10))
z = (9 - xx - yy) / 2 

# plot the plane
ax.plot_surface(xx, yy, z, alpha=0.5)

plt.show()

enter image description here

Plotting the point is straightforward with scatter()

Sandipan Dey
  • 21,482
  • 2
  • 51
  • 63