3D surface plot using matplotlib in python

Question

I have three separate 1-D arrays that contain data points for the surface of an ellipsoid. The data points were generated through a C code and the points were stored in a .CSV file.

I want to plot a 3D surface from these data points using the plot_surface function in matplotlib. I started with plotting a 3D scatter plot as a sanity check to ensure that the data points I have belong to an ellipsoid and not some random shape.

Then I tried the plot_surface() function in matplotlib and I tried various methods that have been discussed so far on SO and you can see them in the code given towards the end. I'll list two of the outputs, since my low rep doesn't allow me to post all my results.

1. The post titled "surface plots in matplotlib", which resulted in the output: This method is implemented in my code under the comment method-4. I also implemented something similar to it, but without the unique() function and it gave me a similar output again. The last part is implemented under method-3 in my code.
2. The post titled "plotting 3-tuple data points on a surface", which resulted in the output: This part is implemented under method-2 in the code.

As you can see, none of the methods really worked for me and although the scatter plot confirms that the data points indeed belong to an ellipsoid, the surface plots given an erroneous result.

Am I doing something wrong? If not, then can you suggest an alternative method to correctly plot the 3D surface for my case? The data I am using can be found at the following link: https://drive.google.com/file/d/0BwTffmdLhwB3b0JOMXdHYzFTSGc/view?usp=sharing

I am using python2.7, OS: ubuntu-14.04. I am new to python so it'll be great if you can also provide explanation along with the solution. Thank you very much.

'''

Copyright (c) 2016 Abhishek Agrawal ([email protected])
Distributed under the MIT License.
See accompanying file LICENSE.md or copy at http://opensource.org/licenses/MIT
'''

# Set up modules and packages
# I/O
import csv
from pprint import pprint

# Numerical
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
import math

# 3D visualization special package
import mayavi
from mayavi import mlab

# System
import sys
import time
from tqdm import tqdm

print ""
print "---------------------------------------------------------------------------------"
print "                                 NAOS                                            "
print "                                                                                 "
print "         Copyright (c) 2016, A. Agrawal ([email protected])        "
print "---------------------------------------------------------------------------------"
print ""

# Start timer.
start_time = time.time( )

# Get plotting packages
import matplotlib
import matplotlib.colors
import matplotlib.axes
import matplotlib.lines as mlines
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from matplotlib import rcParams
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import axes3d
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.tri as tri

# Operations
# Read data in csv file. data returned as a panda series.
data = pd.read_csv( '../data/ellipsoidSurfacePoints.csv' )

# Plot 3D surface of the ellipsoid
fig = plt.figure()
ax = fig.gca( projection = '3d' )
ax.set_xlabel('x [km]')
ax.set_ylabel('y [km]')
ax.set_zlabel('z [km]')
ax.ticklabel_format(style='sci', axis='both', scilimits=(0,0))

x = data['X'].values
# print x[1:100]
y = data['Y'].values
# print y
z = data['Z'].values
# print z
r = np.sqrt( x**2 + y**2 + z**2 )

# **************** trisurf, scatter and wireframe ************ #
# triang = tri.Triangulation( x, y )
# ax.plot_trisurf( x, y, z, triangles=triang.triangles, cmap=cm.jet, linewidth=0.1 )
# ax.scatter( x, y, z )
# ax.plot_wireframe( x, y, z )
# plt.show()

# **************** Method - 1 ******************************** #
# pts = mayavi.mlab.points3d( x, y, z, z )
# mesh = mayavi.mlab.pipeline.delaunay2d( pts )
# pts.remove( )
# surf = mayavi.mlab.pipeline.surface( mesh )
# mayavi.mlab.show( )

# **************** Method - 2 ******************************** #
# x1 = np.linspace( x.min(), x.max() )
# y1 = np.linspace( y.min(), y.max() )
# xx, yy = np.meshgrid( x1, y1 )
# zz = griddata( ( x, y ), z, ( x1, y1 ), method='cubic' )
# ax.plot_surface( xx, yy, zz, rstride=5, cstride=5, cmap=cm.jet, linewidth=0.1, antialiased=False )
# plt.show()

# **************** Method - 3 ******************************** #
x1 = np.linspace( x.min(), x.max() )
y1 = np.linspace( y.min(), y.max() )
xx, yy = np.meshgrid( x1, y1 )
zz = griddata( ( x, y ), z, ( xx, yy ), method='cubic' )
ax.plot_surface( xx, yy, zz, rstride=5, cstride=5, cmap=cm.jet, linewidth=0.1, antialiased=False )
plt.show()

# **************** Method - 4 ******************************** #
# x1 = np.linspace( x.min(), x.max(), len( data['X'].unique() ) )
# y1 = np.linspace( y.min(), y.max(), len( data['Y'].unique() ) )
# xx, yy = np.meshgrid( x1, y1 )
# zz = griddata( ( x, y ), z, ( xx, yy ), method='cubic' )
# ax.plot_surface( xx, yy, zz, rstride=5, cstride=5, cmap=cm.jet, linewidth=0.1, antialiased=False )
# plt.show()

# **************** Method - 5 ******************************** #
# xx, yy = np.mgrid[ min(x):max(x):100j, min(y):max(y):100j ]
# zz = griddata( ( x, y ), z, ( xx, yy ), method='cubic' )
# ax.plot_surface( xx, yy, zz, rstride=5, cstride=5, cmap=cm.jet, linewidth=0.1, antialiased=False )
# plt.show()

# Stop timer
end_time = time.time( )

# Print elapsed time
print "Script time: " + str("{:,g}".format(end_time - start_time)) + "s"

print ""
print "------------------------------------------------------------------"
print "                         Exited successfully!                     "
print "------------------------------------------------------------------"
print ""

---

Answer 1

I think what you're looking for is plot_trisurf which can handle unstructured 1D vectors of data. The code below uses your data to make the plots. The key issue with your data is that each X/Y point has 3 z values which plot_trisurf cannot handle automatically (though there are ways round this, see questions about plotting spheres for example). I have dealt with this problem by splitting the replicated X/Y points into separate dataframes.

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

fileDir = "data.csv"

data = pd.read_csv(fileDir, engine = 'c', float_precision = 'round_trip', dtype=np.float64)

dataTop = data.drop_duplicates(subset=['x', 'y'], keep='first', inplace=False)
XTop = dataTop['x']
YTop = dataTop['y']
ZTop = dataTop['z']

dataMid = data.drop_duplicates(subset=['x', 'y'], keep=False, inplace=False)
XMid = dataMid['x']
YMid = dataMid['y']
ZMid = dataMid['z']

dataBottom = data.drop_duplicates(subset=['x', 'y'], keep='last', inplace=False)
XBottom = dataBottom['x']
YBottom = dataBottom['y']
ZBottom = dataBottom['z']

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

ax.plot_trisurf(XTop, YTop, ZTop, cmap='viridis', alpha=0.5)
ax.plot_trisurf(XMid, YMid, ZMid, cmap='viridis', alpha=0.5)
ax.plot_trisurf(XBottom, YBottom, ZBottom, cmap='viridis', alpha=0.5)

plt.show()

If you can tell me more about the plot you expect to see I can try be more specific. See Matplotlib like matlab's trisurf for more info on triangulation.