Contour plots with Excel and Matplotlib – 2 

The procedures discussed in the previous post required the contour data to be arranged on a regular rectangular grid, with the data points listed in a 2D array. As far as I know there is no built in alternative in Excel, but Matplotlib provides functions that allow contour plots to be generated from data at random locations. Examples of two alternative methods are given at: Contour plot of irregularly spaced data

I have modified the code to plot the data in Excel, and to return two separate plots:

@xl_func
def contour_examples(out=1):
    
    np.random.seed(19680801)
    npts = 200
    ngridx = 100
    ngridy = 200
    x = np.random.uniform(-2, 2, npts)
    y = np.random.uniform(-2, 2, npts)
    z = x * np.exp(-x**2 - y**2)

    # fig, (ax1, ax2) = plt.subplots(nrows=2)
    if out == 1:
        fig, ax1 = plt.subplots()

        # -----------------------
        # Interpolation on a grid
        # -----------------------
        # A contour plot of irregularly spaced data coordinates
        # via interpolation on a grid.

        # Create grid values first.
        xi = np.linspace(-2.1, 2.1, ngridx)
        yi = np.linspace(-2.1, 2.1, ngridy)

        # Linearly interpolate the data (x, y) on a grid defined by (xi, yi).
        triang = tri.Triangulation(x, y)
        interpolator = tri.LinearTriInterpolator(triang, z)
        Xi, Yi = np.meshgrid(xi, yi)
        zi = interpolator(Xi, Yi)

        # Note that scipy.interpolate provides means to interpolate data on a grid
        # as well. The following would be an alternative to the four lines above:
        # from scipy.interpolate import griddata
        # zi = griddata((x, y), z, (xi[None, :], yi[:, None]), method='linear')

        ax1.contour(xi, yi, zi, levels=14, linewidths=0.5, colors='k')
        cntr1 = ax1.contourf(xi, yi, zi, levels=14, cmap="RdBu_r")

        fig.colorbar(cntr1, ax=ax1)
        ax1.plot(x, y, 'ko', ms=3)
        ax1.set(xlim=(-2, 2), ylim=(-2, 2))
        ax1.set_title('grid and contour (%d points, %d grid points)' %
                    (npts, ngridx * ngridy))
        pyxll.plot(fig)
        return 'Contour example1'
    else:

    # ----------
    # Tricontour
    # ----------
    # Directly supply the unordered, irregularly spaced coordinates
    # to tricontour.
        fig, ax2 = plt.subplots()

        ax2.tricontour(x, y, z, levels=14, linewidths=0.5, colors='k')
        cntr2 = ax2.tricontourf(x, y, z, levels=14, cmap="RdBu_r")

        fig.colorbar(cntr2, ax=ax2)
        ax2.plot(x, y, 'ko', ms=3)
        ax2.set(xlim=(-2, 2), ylim=(-2, 2))
        ax2.set_title('tricontour (%d points)' % npts)

        plt.subplots_adjust(hspace=0.5)
        # plt.show()
        pyxll.plot(fig)
        return 'Contour example2'

The resulting plots are identical to those from the original code, except that the axes are not plotted to scale:

Interpolate on a grid
Use Matplotlib tricontour

The tricontour option has been used with the data from the simple FEA analysis used in the previous post, allowing the input data to include all the mid-side nodes. The code has also been modified to allow plotting to scale, and use of a Matplotlib colormap. The resulting plot is shown alongside the original Strand7 plot:

In the revised code the data is input as three 1D arrays, the plot width and height may be specified, or set height = 0 to plot to scale, and the number of contour levels, and any colormap may be specified:

@xl_func
@xl_arg('x', 'numpy_array', ndim = 1)
@xl_arg('y', 'numpy_array', ndim = 1)
@xl_arg('z', 'numpy_array', ndim = 1)
@xl_arg('margins', 'numpy_array', ndim = 1)
@xl_arg('levels', 'int')
@xl_arg('title', 'str')
def contour_plot(x, y, z, width = 10., height = 0., ms = 0, title = '', margins = [], levels = 8, cmap = 'RdBu_r'):
    # ----------
    # Tricontour
    # ----------
    # Directly supply the unordered, irregularly spaced coordinates
    # to tricontour.
    xmin = np.min(x)
    xmax = np.max(x)
    ymin = np.min(y)
    ymax = np.max(y)
    if height == 0: height = (width-1) * (ymax-ymin)/(xmax-xmin)
    fig, ax = plt.subplots(figsize=(width, height), dpi=72)
    if margins == []:
        plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1)
    else:
        left = 0.1 if np.isnan(margins[0]) else margins[0] 
        right = 0.9 if np.isnan(margins[1]) else margins[1]
        bottom = 0.1  if np.isnan(margins[2]) else margins[2]
        top = 0.9 if np.isnan(margins[3]) else margins[3]
        plt.subplots_adjust(left=left, right=right, bottom=bottom, top=top)
   
    ax.tricontour(x, y, z, levels=levels, linewidths=0.0, colors='k')
    cntr2 = ax.tricontourf(x, y, z, levels=levels, cmap=cmap)

    fig.colorbar(cntr2, ax=ax)
    ax.plot(x, y, 'ko', ms=ms)
    xmin = np.min(x)
    xmax = np.max(x)
    ymin = np.min(y)
    ymax = np.max(y)
    ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax))
    ax.set_title(title)

    pyxll.plot(fig)
    return 'Contour plot'

A more complex example illustrates a problem with the tricontour function when plotting over an area with concave segments. The image below shows vertical deflections from a Strand7 analysis of a retaining wall with piled foundations:

The Matplotlib plot with the same data inserts an additional triangular area in front of the wall:

Limiting the plot to the region to the right of the wall face avoids this problem, and generates a plot closely matching the Strand7 output:

It seems that it is nor straightforward to create a contour plot on a shape with a concave boundary. The best information I could find was:

However, this may not solve all problems, particular where the problem is defined on an irregular domain. Also, in the case where the domain has one or more concave areas, the delaunay triangulation may result in generating spurious triangles exterior to the domain. In such cases, these rogue triangles have to be removed from the triangulation in order to achieve the correct surface representation. For these situations, the user may have to explicitly include the delaunay triangulation calculation so that these triangles can be removed programmatically. Under these circumstances, the following code could replace the previous plot code:

Response from Graham G at: https://stackoverflow.com/questions/9170838/surface-plots-in-matplotlib

See the Stackoverflow discussion for more details and suggested code.

This entry was posted in Arrays, Charts, Excel, Finite Element Analysis, Link to Python, Newton, NumPy and SciPy, PyXLL, Strand7, UDFs and tagged , , , , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.