6

I have reviewed both of these threads, but am still struggling to make a 3D surface plot from a numpy array of x, y, z coordinates.

My array looks like this:

>>> points
array([[ 322697.1875    , 3663966.5       ,  -30000.        ],
       [ 325054.34375   , 3663966.5       ,  -30000.        ],
       [ 325054.34375   , 3665679.5       ,  -30000.        ],
       [ 322697.1875    , 3665679.5       ,  -30000.        ],
       [ 322697.1875    , 3663966.5       ,  -27703.12304688],
       [ 325054.34375   , 3663966.5       ,  -27703.15429688],
       [ 325054.34375   , 3665679.5       ,  -27703.70703125],
       [ 322697.1875    , 3665679.5       ,  -27703.67382812]])

ax.plot_surface accepts x, y, z points so I convert the above array into separate pieces below:

x = points[:, 0]
y = points[:, 1]
z = points[:, 2]

I then put it into a meshgrid for passing into ax.plot_surface():

import numpy as np

X, Y, Z = np.meshgrid(x, y, z)

And then try to plot:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(16,10))
ax = plt.axes(projection = '3d')
ax.plot_surface(X, Y, Z, alpha=0.5)
plt.show()

When I run this I receive an error: rows, cols = Z.shape ValueError: too many values to unpack (expected 2).

I'm not sure where to go with this now, I don't expect the answer but a push in the correct direction would be great.

I would like the output to be similar in appearance to this but with my data: enter image description here

UPDATE: If I do not include z in the meshgrid, but only x and y, I get this output when I run ax.plot_surface(X, Y, z, alpha=0.5): enter image description here

This is really close, but I want all the sides to be filled in. Only one is showing as filled in. I've added the point coordinates to show the boundaries. I feel like it has something to do with the meshgrid that I'm creating. Here is the output of X, Y:

>>> X, Y = np.meshgrid(x, y)
(array([[322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ],
       [322697.1875 , 325054.34375, 325054.34375, 322697.1875 ,
        322697.1875 , 325054.34375, 325054.34375, 322697.1875 ]]), array([[3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5],
       [3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5, 3663966.5,
        3663966.5, 3663966.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5],
       [3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5, 3665679.5,
        3665679.5, 3665679.5]]))

If I just take x, y unique values I get an error thrown:

x = np.unique(x)
y = np.unique(y)

>>> x
array([322697.1875 , 325054.34375])
>>> y
array([3663966.5, 3665679.5])

X, Y = np.meshgrid(x, y)
>>> X, Y
(array([[322697.1875 , 325054.34375],
       [322697.1875 , 325054.34375]]), array([[3663966.5, 3663966.5],
       [3665679.5, 3665679.5]]))

>>> ax.plot_surface(X, Y, z, alpha=0.5)
Traceback (most recent call last):
  File "<pyshell#61>", line 1, in <module>
    ax.plot_surface(X, Y, z, alpha=0.5)
  File "/Users/NaN/anaconda/envs/py36/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py", line 1586, in plot_surface
    X, Y, Z = np.broadcast_arrays(X, Y, Z)
  File "/Users/NaN/anaconda/envs/py36/lib/python3.6/site-packages/numpy/lib/stride_tricks.py", line 259, in broadcast_arrays
    shape = _broadcast_shape(*args)
  File "/Users/NaN/anaconda/envs/py36/lib/python3.6/site-packages/numpy/lib/stride_tricks.py", line 193, in _broadcast_shape
    b = np.broadcast(*args[:32])
ValueError: shape mismatch: objects cannot be broadcast to a single shape
10
  • At what line in your code is the error appearing? I guess it comes from Z having three dimensions instead of 2. Commented Jun 6, 2019 at 9:04
  • Don't include z in your meshgrid. Commented Jun 6, 2019 at 10:13
  • @tnknepp @cripcate Thanks, I tried just meshing with x and y and plugging in z separately at ax.plot_surface but am getting a shape mismatch error as seen in the updated post Commented Jun 6, 2019 at 14:21
  • In your case X,Y, and Z each need to be a 5x5 array if you want to use a single surface plot instead of multiple as in the linked question. I might provide an answer over there later on. Commented Jun 6, 2019 at 17:47
  • @ImportanceOfBeingErnest Yes, in reality, I have many polygons to plot in one single chart, i.e. many points arrays, therefore I need to figure out a way to put this in a function or loop to plot all at once without hardwiring any values Commented Jun 6, 2019 at 17:52

1 Answer 1

1

The arrays x, y, z need to be parametrized in two dimensions. One way of doing this is to use spherical coordinates as e.g. in Plot surfaces on a cube.

The remaining task is to distill the unique coordinates from the input data. I'm assuming here that there are only 2 distinct values per dimension.

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

def get_cube():   
    phi = np.arange(1,10,2)*np.pi/4
    Phi, Theta = np.meshgrid(phi, phi)

    x = np.cos(Phi)*np.sin(Theta)
    y = np.sin(Phi)*np.sin(Theta)
    z = np.cos(Theta)/np.sqrt(2)
    return x,y,z


points = np.array([[ 322697.1875    , 3663966.5       ,  -30000. ],
                   [ 325054.34375   , 3663966.5       ,  -30000. ],
                   [ 325054.34375   , 3665679.5       ,  -30000. ],
                   [ 322697.1875    , 3665679.5       ,  -30000. ],
                   [ 322697.1875    , 3663966.5       ,  -27703.12],
                   [ 325054.34375   , 3663966.5       ,  -27703.12],
                   [ 325054.34375   , 3665679.5       ,  -27703.12],
                   [ 322697.1875    , 3665679.5       ,  -27703.12]])

ux = np.unique(points[:,0])
uy = np.unique(points[:,1])
uz = np.unique(points[:,2])

x,y,z = get_cube()
offset = lambda X, o: o[0] + (X+.5)*np.diff(o)[0]


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

ax.plot_surface(offset(x, ux), offset(y, uy), offset(z, uz))

plt.show()

enter image description here

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.