0

I have an array of temperature, and an array of 3D points such that temperature[n] is the temperatureat point[n]. How can I make a contour plot of that data (looking at 2D planes)?

I thought on creating a function extract_plane in which a parametric equation of a plane would be passed as an argument, then the points in it and the corresponding temperatures returned (that is what I need help with).

As an example:

import numpy as np
import matplotlib.pyplot as plt

points = np.array([[0., 0., 0.],
                      [1., 0., 0.],
                      [1., 1., 0.],
                      [0., 1., 0.],
                      [0., 1., 1.],
                      [0., 0., 1.],
                      [1., 0., 1.],
                      [1., 1., 1.]])
temperature = np.array([0, 0, 0, 0, 1, 1, 1, 1.])

I need help creating the following function. As it is, it only extracts the points laying in the plane z=0.

def extract_plane(points, temperature, equation):
    """
    Given a set of 3D points, and their corresponding temperatures,
    extracts a set of points that are in the plane defined by equation
    along their temperatures.

    Parameters
    ----------
    points : ndarray (3D)
        The set of points.
    temperature : ndarray (1D)
        The temperatures at the points such that temperature[n] is 
        the temperature in points[n].
    equation : ????
        The equation that defines a 2D plane (cross-section) where
        the temperature is wanted to be plotted.

    Returns
    -------
    coord : ndarray (1D)
        The set of points that are in the plane defined by equation.
    temp : ndarray (1D)
        The set of temperatures in which temp[n] coresponds to coord[n].
    """
    temp = []
    coord = []
    # plane z=0
    for n in range(points.shape[0]):
        if (points[n,2] == 0.):
            temp += [temperature[n]]
            coord += [points[n]]

    temp = np.array(temp)
    coord = np.array(coord)
    return coord, temp

And use griddata found in this cookbook to reshape temp so it can be plotted:

# griddata.py - 2010-07-11 ccampo
def griddata(x, y, z, binsize=0.01, retbin=True, retloc=True):
    """
    Place unevenly spaced 2D data on a grid by 2D binning (nearest
    neighbor interpolation).

    Parameters
    ----------
    x : ndarray (1D)
        The idependent data x-axis of the grid.
    y : ndarray (1D)
        The idependent data y-axis of the grid.
    z : ndarray (1D)
        The dependent data in the form z = f(x,y).
    binsize : scalar, optional
        The full width and height of each bin on the grid.  If each
        bin is a cube, then this is the x and y dimension.  This is
        the step in both directions, x and y. Defaults to 0.01.
    retbin : boolean, optional
        Function returns `bins` variable (see below for description)
        if set to True.  Defaults to True.
    retloc : boolean, optional
        Function returns `wherebins` variable (see below for description)
        if set to True.  Defaults to True.

    Returns
    -------
    grid : ndarray (2D)
        The evenly gridded data.  The value of each cell is the median
        value of the contents of the bin.
    bins : ndarray (2D)
        A grid the same shape as `grid`, except the value of each cell
        is the number of points in that bin.  Returns only if
        `retbin` is set to True.
    wherebin : list (2D)
        A 2D list the same shape as `grid` and `bins` where each cell
        contains the indicies of `z` which contain the values stored
        in the particular bin.

    Revisions
    ---------
    2010-07-11  ccampo  Initial version
    """
    # get extrema values.
    xmin, xmax = x.min(), x.max()
    ymin, ymax = y.min(), y.max()

    # make coordinate arrays.
    xi      = np.arange(xmin, xmax+binsize, binsize)
    yi      = np.arange(ymin, ymax+binsize, binsize)
    xi, yi = np.meshgrid(xi,yi)

    # make the grid.
    grid           = np.zeros(xi.shape, dtype=x.dtype)
    nrow, ncol = grid.shape
    if retbin: bins = np.copy(grid)

    # create list in same shape as grid to store indices
    if retloc:
        wherebin = np.copy(grid)
        wherebin = wherebin.tolist()

    # fill in the grid.
    for row in range(nrow):
        for col in range(ncol):
            xc = xi[row, col]    # x coordinate.
            yc = yi[row, col]    # y coordinate.

            # find the position that xc and yc correspond to.
            posx = np.abs(x - xc)
            posy = np.abs(y - yc)
            ibin = np.logical_and(posx < binsize/2., posy < binsize/2.)
            ind  = np.where(ibin == True)[0]

            # fill the bin.
            bin = z[ibin]
            if retloc: wherebin[row][col] = ind
            if retbin: bins[row, col] = bin.size
            if bin.size != 0:
                binval         = np.median(bin)
                grid[row, col] = binval
            else:
                grid[row, col] = np.nan   # fill empty bins with nans.

    # return the grid
    if retbin:
        if retloc:
            return grid, bins, wherebin
        else:
            return grid, bins
    else:
        if retloc:
            return grid, wherebin
        else:
            return grid

Then plot:

coord, temp = extract_plane(points, temperature, None)
x = coord[:,0]
y = coord[:,1]
g = griddata(x, y, temp, 1., False, False)
plt.contourf(g)
  • Take a look at [`scipy.interpolate.griddata`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html). It will take your data and put it on a 2D grid just like you need. You may have to fiddle with the input just a bit though. – Mad Physicist Oct 07 '17 at 16:25
  • I think you first need to become clear about what you actually want to show. You have 4 different points in your data and each has two different values. That would be the edges of a cube in 3D, but if you want to plot it in 2D, it's not really clear what you would like to show. – ImportanceOfBeingErnest Oct 07 '17 at 16:37

1 Answers1

3

The functions from the question look a lot more complicated than it needs to be. Using scipy.interpolate.griddata allows to interpolate values on a grid.

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt

points = np.array([[0., 0., 0.],
                      [1., 0., 0.],
                      [1., 1., 0.],
                      [0., 1., 0.],
                      [0., 1., 1.],
                      [0., 0., 1.],
                      [1., 0., 1.],
                      [1., 1., 1.]])
temperature = np.array([0, 0, 0, 0, 1, 1, 1, 1.])

grid_y, grid_x = np.mgrid[0:1:25j, 0:1:25j]
# equation along which to interpolate
equation = lambda x,y : 0.8*(1-x)
grid_z = equation(grid_x, grid_y)

interp = griddata(points, temperature, (grid_x, grid_y, grid_z), method='linear')

plt.subplot(121)
#plt.contourf(grid_x,grid_y, interp,  origin='lower',vmin=0,vmax=1)
plt.imshow(interp,  origin='lower',vmin=0,vmax=1)
plt.title('temperature along 0.8*(1-x)')
plt.xlabel("x")
plt.ylabel("y")

from mpl_toolkits.mplot3d import Axes3D
ax = plt.subplot(122, projection=Axes3D.name)
ax.scatter(points[:,0], points[:,1], points[:,2], c=temperature)
ax.set_zlim(-.1,1.1)
ax.plot_surface(grid_x,grid_y,grid_z, facecolors=plt.cm.viridis(interp),
                linewidth=0, antialiased=False, shade=False)
ax.set_xlabel("x")
ax.set_ylabel("y")
plt.show()

enter image description here

For equation = lambda x,y : x*(y**.5):

enter image description here

Of course using contourf is equally possible, plt.contourf(grid_x,grid_y, interp, origin='lower',vmin=0,vmax=1):

enter image description here

ImportanceOfBeingErnest
  • 321,279
  • 53
  • 665
  • 712
  • I'm trying to do the same thing but interpolated data is always nan. I guess it's because it tries to interpolate points outside of convex hull and i didn't define fill_value. But why all of them are nan? I defined grid like this grid_y, grid_x = np.mgrid[min_y:max_y:25j, min_x:max_x:25j] – ivan.zd Jan 10 '18 at 13:07