plot data stored in a file using (mayavi) mlab.contour3d?

Question

I am struggling to do a simple contour3d plot in mayavi using data that is read in from a file. The data is a regular 3d cartesian grid in the form

   x1   y1    z1    val(x1,y1,z1)  
   x1   y1    z2    val(x1,y1,z2)  
    .    .     .       .  
   x1   y1     z_K  val(x1,y1,z_K)  

   x1   y2     z1    val(x1,y2,z1)  
    .    .     .       .  
    .    .     .       .  
   x_L   y_M    z_K    val(x_L,y_M,z_K)  

(Here the first 3 values on each row give the (x,y,z) coordinates of the point while the 4th value gives the value of a scalar at that point. I can paste a minimal example file if useful)

The data file can be successfully read in using numpy.loadtxt , but how do I get from there to plotting an isosurface using (mayavi) mlab.contour3d ? I think the default output array from loadtxt is not in the right format for mlab.contour3d.

All the examples I have been able to find of mlab.contour3d generate a grid using ogrid, then plot a simple function (sin etc) of this grid. I have been able to run those examples successfully, but they don't tell me how to read in data from a file into the right format of array ready for plotting. I am sure it would help a great many newbies in a similar position if someone could give me a pointer ; plotting 3d data stored in a file generated by another program must surely be one of the most common plotting tasks scientists have to do.

Answer 1

I found the following to work well.

x,y,z,d are 1 dimensional numpy arrays from a raw input text file, where d = f(x,y,z)

# Import relevant modules
import numpy as np
import scipy.interpolate as si
from mayavi import mlab

# Generate the figure
figure = mlab.figure()

# Generate the 3D mesh
xi, yi, zi = np.meshgrid(*([np.linspace(0,1,100)] * 3))

# Interpolate d over the 3D mesh
di = si.griddata((x, y, z), d, (xi, yi, zi))

# Generate the desired figure
min_d = d.min()
max_d = d.max()
grid = mlab.pipeline.scalar_field(di)
mlab.pipeline.volume(grid, vmin=min_d, vmax=min_d + .5*(max_d-min_d))

# Other good options too
## mlab.contour3d(di, vmin=min_d, vmax=min_d + .8*(max_d-min_d))
## pts = mlab.points3d(x, y, z, d)
mlab.axes()
mlab.show()

Answer 2

Apologies to Mr E for the delayed response, I had no internet access over the weekend. I did, however, have access to a computer, and I have found an answer to my question, albeit a slightly ugly one. I am sure the following method can be improved, and, for the reasons stated above, I know this is an important basic use of 3d plotting software, so if anyone has any improvements that they can offer they will be greatly appreciated. The following explanation assumes you are using Linux, but I am sure it is straightforward to do the same thing (i.e. saving and running python files) on other OSes.

First to generate a sample data file. The real data comes from a Fortran program, but for the current test purposes I will use python to generate a sample data file.Store the following code in the file "gen-data.py" (copy-pasting into your favourite text editor and then click on "save as" etc.)

#!/usr/bin/python
import numpy as numpy

# run this program with 
#" python  gen-data.py > stream.out  "
# to generate sample data file stream.out  
x_low = -1 ; x_upp = 1 ; NX = 5
y_low = -2 ; y_upp = 2 ; NY = 3
z_low = -3 ; z_upp = 3 ; NZ = 3

x_arr = numpy.linspace(x_low,x_upp,num = NX, endpoint=True)
y_arr = numpy.linspace(y_low,y_upp,num = NY, endpoint=True)
z_arr = numpy.linspace(z_low,z_upp,num = NZ, endpoint=True)

#the following line generates data (for distance from the origin) 
# over  a structured grid
for x in x_arr:
    for y in y_arr:
        for z in z_arr:
            print x , y , z , x**2 + y**2 + z**2

Run the program using python gen-data.py > stream.out which will store a sample data file of the type described above in the datafile "stream.out". You should have a file with the values:
-1.0 -2.0 -3.0 14.0
-1.0 -2.0 0.0 5.0
-1.0 -2.0 3.0 14.0
-1.0 0.0 -3.0 10.0
-1.0 0.0 0.0 1.0
-1.0 0.0 3.0 10.0
-1.0 2.0 -3.0 14.0
-1.0 2.0 0.0 5.0
-1.0 2.0 3.0 14.0
-0.5 -2.0 -3.0 13.25
. . . .
0.5 2.0 3.0 13.25
1.0 -2.0 -3.0 14.0
1.0 -2.0 0.0 5.0
1.0 -2.0 3.0 14.0
1.0 0.0 -3.0 10.0
1.0 0.0 0.0 1.0
1.0 0.0 3.0 10.0
1.0 2.0 -3.0 14.0
1.0 2.0 0.0 5.0
1.0 2.0 3.0 14.0
Each row in the data file is of the form
x y z V(x,y,z)
where x,y,z describe the x,y,x coords of a point in space and V(x,y,z) denotes the value of a scalar at that point.

Plotting the data
Our problem is how to plot this data using mayavi - I am particularly interested in plotting isosurfaces, which can be achieved using the contour3d command. The numerous examples on the web show contour3d plotting data which is generated using the mgrid command. (There are also examples using the ogrid command , but for me mgrid was easier to understand). Strategy: If we can manipulate our data to have the same shape and from as output from the mgrid command we should be able to plot it. Analysing output from mgrid , it was clear that what was required were 3 three-dimensional numpy arrays to store x, y and z coordinates, and another three-dimensional numpy array to store the scalar values (V in the above) at these points. The following program achieves these steps. I think the program could definitely be improved: the routine fill_up_array could, I am sure, be replaced by a one-liner by someone who knows about array slicing in python, and there are probably other places than can be improved. I can't help thinking that the whole thing can probably be done in 1 or 2 lines for someone who knows what they are doing in numpy/mayavi, but the following program is , I believe, easy to understand and it works (you should see truncated spherical surfaces in the plot that appears).
Save the following file to "hope.py" and run with
python hope.py

import numpy
from mayavi import mlab

def fill_up_array(output_arr,source_arr,nx,ny,nz,posn):
# takes a slice of results from source_arr and writes to output_arr
# there is probably an easy one liner to do this ?
#TODO: add checks to ensure input is sensible
    for i in range(nx):
        for j in range(ny):
            for k in range(nz):
                output_arr[i][j][k] = source_arr[i][j][k][posn]

# these numbers have to correspond to those used in gen-data.py
NX = 5 ; NY = 3 ; NZ = 3 

NDIM = 4  # number of columns in data file, 4 for current example


#initialise arrays:
# xx will contain x coordinate of data point
# yy will contain y coordinate of data point
# zz will contain z coordinate of data point
# VV will contain sample scalar value at (x,y,z)
xx = numpy.zeros((NX,NY,NZ))
yy = numpy.zeros((NX,NY,NZ))
zz = numpy.zeros((NX,NY,NZ))
VV =  numpy.zeros((NX,NY,NZ))


#now read in values from stream.out file to these arrays
full = numpy.loadtxt("stream.out")
fy = numpy.reshape(full, (NX,NY,NZ,NDIM))


fill_up_array(xx,fy,NX,NY,NZ,0)
fill_up_array(yy,fy,NX,NY,NZ,1)
fill_up_array(zz,fy,NX,NY,NZ,2)
fill_up_array(VV,fy,NX,NY,NZ,3)


#test plot
mlab.contour3d(xx, yy, zz, VV)
mlab.show()