Use qhull for delaunay triangulation

Author: ianthomas23

CHANGELOG

@@ -1,6 +1,11 @@
+2013-11-28 Added qhull extension module to perform Delaunay triangulation more
+           robustly than before.  It is used by tri.Triangulation (and hence
+           all pyplot.tri* methods) and mlab.griddata.  Deprecated
+           matplotlib.delaunay module. - IMT
+
 2013-10-27 Added get_rlabel_position and set_rlabel_position methods to
-	   PolarAxes to control angular position of radial tick labels.
-	   
+           PolarAxes to control angular position of radial tick labels.
+
 2013-10-06 Add stride-based functions to mlab for easy creation of 2D arrays
            with less memory.
 
@@ -27,16 +32,16 @@
            and matplotlib.units.Registry.
 
 2013-06-26 Refactored the axes module: the axes module is now a folder,
-	   containing the following submodule:
-	      - _subplots.py, containing all the subplots helper methods
-	      - _base.py, containing several private methods and a new
-		_AxesBase class. This _AxesBase class contains all the methods
-		that are not directly linked to plots of the "old" Axes
-	      - _axes.py contains the Axes class. This class now inherits from
-		_AxesBase: it contains all "plotting" methods and labelling
-		methods.
-	    This refactoring should not affect the API. Only private methods
-	    are not importable from the axes module anymore.
+           containing the following submodule:
+              - _subplots.py, containing all the subplots helper methods
+              - _base.py, containing several private methods and a new
+                _AxesBase class. This _AxesBase class contains all the methods
+                that are not directly linked to plots of the "old" Axes
+              - _axes.py contains the Axes class. This class now inherits from
+                _AxesBase: it contains all "plotting" methods and labelling
+                methods.
+           This refactoring should not affect the API. Only private methods
+           are not importable from the axes module anymore.
 
 2013-05-18 Added support for arbitrary rasterization resolutions to the
            SVG backend. Previously the resolution was hard coded to 72
@@ -52,22 +57,21 @@
 
 2013-04-25 Changed all instances of:
 
-	   from matplotlib import MatplotlibDeprecationWarning as mplDeprecation
-	   to:
-
-	   from cbook import mplDeprecation
+           from matplotlib import MatplotlibDeprecationWarning as mplDeprecation
+           to:
 
-	   and removed the import into the matplotlib namespace in __init__.py
-	   Thomas Caswell
+           from cbook import mplDeprecation
 
+           and removed the import into the matplotlib namespace in __init__.py
+           Thomas Caswell
 
 2013-04-15 Added 'axes.xmargin' and 'axes.ymargin' to rpParams to set default
-	   margins on auto-scaleing. - TAC
+           margins on auto-scaleing. - TAC
 
 2013-04-16 Added patheffect support for Line2D objects.  -JJL
 
 2013-03-19 Added support for passing `linestyle` kwarg to `step` so all `plot`
-      	   kwargs are passed to the underlying `plot` call.  -TAC
+           kwargs are passed to the underlying `plot` call.  -TAC
 
 2013-02-25 Added classes CubicTriInterpolator, UniformTriRefiner, TriAnalyzer
            to matplotlib.tri module. - GBy

MANIFEST.in

@@ -14,7 +14,5 @@ recursive-include LICENSE *
 recursive-include examples *
 recursive-include doc *
 recursive-include src  *.cpp *.c *.h *.m
-recursive-include CXX  *.cxx *.hxx *.c *.h
-recursive-include agg24 *
 recursive-include lib *
-recursive-include ttconv *.cpp *.h
+recursive-include extern *

examples/pylab_examples/griddata_demo.py

@@ -6,36 +6,22 @@
 #npts = int(raw_input('enter # of random points to plot:'))
 seed(0)
 npts = 200
-x = uniform(-2,2,npts)
-y = uniform(-2,2,npts)
-z = x*np.exp(-x**2-y**2)
+x = uniform(-2, 2, npts)
+y = uniform(-2, 2, npts)
+z = x*np.exp(-x**2 - y**2)
 # define grid.
-xi = np.linspace(-2.1,2.1,100)
-yi = np.linspace(-2.1,2.1,200)
+xi = np.linspace(-2.1, 2.1, 100)
+yi = np.linspace(-2.1, 2.1, 200)
 # grid the data.
-zi = griddata(x,y,z,xi,yi,interp='linear')
+zi = griddata(x, y, z, xi, yi, interp='linear')
 # contour the gridded data, plotting dots at the nonuniform data points.
-CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
-CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.rainbow, 
+CS = plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
+CS = plt.contourf(xi, yi, zi, 15, cmap=plt.cm.rainbow,
                   vmax=abs(zi).max(), vmin=-abs(zi).max())
-plt.colorbar() # draw colorbar
+plt.colorbar()  # draw colorbar
 # plot data points.
-plt.scatter(x,y,marker='o',c='b',s=5,zorder=10)
-plt.xlim(-2,2)    
-plt.ylim(-2,2)     
+plt.scatter(x, y, marker='o', c='b', s=5, zorder=10)
+plt.xlim(-2, 2)
+plt.ylim(-2, 2)
 plt.title('griddata test (%d points)' % npts)
 plt.show()
-
-# test case that scikits.delaunay fails on, but natgrid passes..
-#data = np.array([[-1, -1], [-1, 0], [-1, 1],
-#                 [ 0, -1], [ 0, 0], [ 0, 1],
-#                 [ 1, -1 - np.finfo(np.float_).eps], [ 1, 0], [ 1, 1],
-#                ])
-#x = data[:,0]
-#y = data[:,1]
-#z = x*np.exp(-x**2-y**2)
-## define grid.
-#xi = np.linspace(-1.1,1.1,100)
-#yi = np.linspace(-1.1,1.1,100)
-## grid the data.
-#zi = griddata(x,y,z,xi,yi)

lib/matplotlib/delaunay/triangulate.py

@@ -11,6 +11,11 @@
 
 from matplotlib._delaunay import delaunay
 from .interpolate import LinearInterpolator, NNInterpolator
+from matplotlib.cbook import warn_deprecated
+warn_deprecated('1.4',
+                name='matplotlib.delaunay',
+                alternative='matplotlib.tri.Triangulation',
+                obj_type='module')
 
 __all__ = ['Triangulation', 'DuplicatePointWarning']
 

lib/matplotlib/mlab.py

@@ -3356,124 +3356,131 @@ def newfunc(val, mask, mval):
     if opened:
         fh.close()
 
-def griddata(x,y,z,xi,yi,interp='nn'):
-    """
-    ``zi = griddata(x,y,z,xi,yi)`` fits a surface of the form *z* =
-    *f*(*x*, *y*) to the data in the (usually) nonuniformly spaced
-    vectors (*x*, *y*, *z*).  :func:`griddata` interpolates this
-    surface at the points specified by (*xi*, *yi*) to produce
-    *zi*. *xi* and *yi* must describe a regular grid, can be either 1D
-    or 2D, but must be monotonically increasing.
-
-    A masked array is returned if any grid points are outside convex
-    hull defined by input data (no extrapolation is done).
-
-    If interp keyword is set to '`nn`' (default),
-    uses natural neighbor interpolation based on Delaunay
-    triangulation.  By default, this algorithm is provided by the
-    :mod:`matplotlib.delaunay` package, written by Robert Kern.  The
-    triangulation algorithm in this package is known to fail on some
-    nearly pathological cases. For this reason, a separate toolkit
-    (:mod:`mpl_tookits.natgrid`) has been created that provides a more
-    robust algorithm fof triangulation and interpolation.  This
-    toolkit is based on the NCAR natgrid library, which contains code
-    that is not redistributable under a BSD-compatible license.  When
-    installed, this function will use the :mod:`mpl_toolkits.natgrid`
-    algorithm, otherwise it will use the built-in
-    :mod:`matplotlib.delaunay` package.
-
-    If the interp keyword is set to '`linear`', then linear interpolation
-    is used instead of natural neighbor. In this case, the output grid
-    is assumed to be regular with a constant grid spacing in both the x and
-    y directions. For regular grids with nonconstant grid spacing, you
-    must use natural neighbor interpolation.  Linear interpolation is only valid if
-    :mod:`matplotlib.delaunay` package is used - :mod:`mpl_tookits.natgrid`
-    only provides natural neighbor interpolation.
-
-    The natgrid matplotlib toolkit can be downloaded from
-    http://sourceforge.net/project/showfiles.php?group_id=80706&package_id=142792
-    """
-    try:
-        from mpl_toolkits.natgrid import _natgrid, __version__
-        _use_natgrid = True
-    except ImportError:
-        import matplotlib.delaunay as delaunay
-        from matplotlib.delaunay import  __version__
-        _use_natgrid = False
-    if not griddata._reported:
-        if _use_natgrid:
-            verbose.report('using natgrid version %s' % __version__)
-        else:
-            verbose.report('using delaunay version %s' % __version__)
-        griddata._reported = True
+
+def griddata(x, y, z, xi, yi, interp='nn'):
+    """Interpolates from a nonuniformly spaced grid to some other
+    grid.
+
+    Fits a surface of the form z = f(`x`, `y`) to the data in the
+    (usually) nonuniformly spaced vectors (`x`, `y`, `z`), then
+    interpolates this surface at the points specified by
+    (`xi`, `yi`) to produce `zi`.
+
+    Parameters
+    ----------
+    x, y, z : 1d array_like
+        Coordinates of grid points to interpolate from.
+    xi, yi : 1d or 2d array_like
+        Coordinates of grid points to interpolate to.
+    interp : string key from {'nn', 'linear'}
+        Interpolation algorithm, either 'nn' for natural neighbor, or
+        'linear' for linear interpolation.
+
+    Returns
+    -------
+    2d float array
+        Array of values interpolated at (`xi`, `yi`) points.  Array
+        will be masked is any of (`xi`, `yi`) are outside the convex
+        hull of (`x`, `y`).
+
+    Notes
+    -----
+    If `interp` is 'nn' (the default), uses natural neighbor
+    interpolation based on Delaunay triangulation.  This option is
+    only available if the mpl_toolkits.natgrid module is installed.
+    This can be downloaded from https://github.com/matplotlib/natgrid.
+    The (`xi`, `yi`) grid must be regular and monotonically increasing
+    in this case.
+
+    If `interp` is 'linear', linear interpolation is used via
+    matplotlib.tri.LinearTriInterpolator.
+
+    Instead of using `griddata`, more flexible functionality and other
+    interpolation options are available using a
+    matplotlib.tri.Triangulation and a matplotlib.tri.TriInterpolator.
+    """
+    # Check input arguments.
+    x = np.asanyarray(x, dtype=np.float64)
+    y = np.asanyarray(y, dtype=np.float64)
+    z = np.asanyarray(z, dtype=np.float64)
+    if x.shape != y.shape or x.shape != z.shape or x.ndim != 1:
+        raise ValueError("x, y and z must be equal-length 1-D arrays")
+
+    xi = np.asanyarray(xi, dtype=np.float64)
+    yi = np.asanyarray(yi, dtype=np.float64)
     if xi.ndim != yi.ndim:
-        raise TypeError("inputs xi and yi must have same number of dimensions (1 or 2)")
-    if xi.ndim != 1 and xi.ndim != 2:
-        raise TypeError("inputs xi and yi must be 1D or 2D.")
-    if not len(x)==len(y)==len(z):
-        raise TypeError("inputs x,y,z must all be 1D arrays of the same length")
-    # remove masked points.
-    if hasattr(z,'mask'):
-        # make sure mask is not a scalar boolean array.
-        if z.mask.ndim:
-            x = x.compress(z.mask == False)
-            y = y.compress(z.mask == False)
-            z = z.compressed()
-    if _use_natgrid: # use natgrid toolkit if available.
-        if interp != 'nn':
-            raise ValueError("only natural neighor interpolation"
-            " allowed when using natgrid toolkit in griddata.")
+        raise ValueError("xi and yi must be arrays with the same number of "
+                         "dimensions (1 or 2)")
+    if xi.ndim == 2 and xi.shape != yi.shape:
+        raise ValueError("if xi and yi are 2D arrays, they must have the same "
+                         "shape")
+    if xi.ndim == 1:
+        xi, yi = np.meshgrid(xi, yi)
+
+    if interp == 'nn':
+        use_nn_interpolation = True
+    elif interp == 'linear':
+        use_nn_interpolation = False
+    else:
+        raise ValueError("interp keyword must be one of 'linear' (for linear "
+                         "interpolation) or 'nn' (for natural neighbor "
+                         "interpolation).  Default is 'nn'.")
+
+    # Remove masked points.
+    mask = np.ma.getmask(z)
+    if not (mask is np.ma.nomask):
+        x = x.compress(~mask)
+        y = y.compress(~mask)
+        z = z.compressed()
+
+    if use_nn_interpolation:
+        try:
+            from mpl_toolkits.natgrid import _natgrid
+        except ImportError:
+            raise RuntimeError("To use interp='nn' (Natural Neighbor "
+                "interpolation) in griddata, natgrid must be installed.  "
+                "Either install it from http://sourceforge.net/projects/"
+                "matplotlib/files/matplotlib-toolkits, or use interp='linear' "
+                "instead.")
+
         if xi.ndim == 2:
-            xi = xi[0,:]
-            yi = yi[:,0]
-        # override default natgrid internal parameters.
-        _natgrid.seti('ext',0)
-        _natgrid.setr('nul',np.nan)
-        # cast input arrays to doubles (this makes a copy)
-        x = x.astype(np.float)
-        y = y.astype(np.float)
-        z = z.astype(np.float)
-        xo = xi.astype(np.float)
-        yo = yi.astype(np.float)
-        if min(xo[1:]-xo[0:-1]) < 0 or min(yo[1:]-yo[0:-1]) < 0:
-            raise ValueError('output grid defined by xi,yi must be monotone increasing')
-        # allocate array for output (buffer will be overwritten by nagridd)
-        zo = np.empty((yo.shape[0],xo.shape[0]), np.float)
-        _natgrid.natgridd(x,y,z,xo,yo,zo)
-    else: # use Robert Kern's delaunay package from scikits (default)
-        if xi.ndim != yi.ndim:
-            raise TypeError("inputs xi and yi must have same number of dimensions (1 or 2)")
-        if xi.ndim != 1 and xi.ndim != 2:
-            raise TypeError("inputs xi and yi must be 1D or 2D.")
-        if xi.ndim == 1:
-            xi,yi = np.meshgrid(xi,yi)
-        # triangulate data
-        tri = delaunay.Triangulation(x,y)
-        # interpolate data
-        if interp == 'nn':
-            interp = tri.nn_interpolator(z)
-            zo = interp(xi,yi)
-        elif interp == 'linear':
-            # make sure grid has constant dx, dy
-            dx = xi[0,1:]-xi[0,0:-1]
-            dy = yi[1:,0]-yi[0:-1,0]
-            epsx = np.finfo(xi.dtype).resolution
-            epsy = np.finfo(yi.dtype).resolution
-            if dx.max()-dx.min() > epsx or dy.max()-dy.min() > epsy:
-                raise ValueError("output grid must have constant spacing"
-                                 " when using interp='linear'")
-            interp = tri.linear_interpolator(z)
-            zo = interp[yi.min():yi.max():complex(0,yi.shape[0]),
-                        xi.min():xi.max():complex(0,xi.shape[1])]
-        else:
-            raise ValueError("interp keyword must be one of"
-            " 'linear' (for linear interpolation) or 'nn'"
-            " (for natural neighbor interpolation). Default is 'nn'.")
-    # mask points on grid outside convex hull of input data.
-    if np.any(np.isnan(zo)):
-        zo = np.ma.masked_where(np.isnan(zo),zo)
-    return zo
-griddata._reported = False
+            # natgrid expects 1D xi and yi arrays.
+            xi = xi[0, :]
+            yi = yi[:, 0]
+
+        # Override default natgrid internal parameters.
+        _natgrid.seti('ext', 0)
+        _natgrid.setr('nul', np.nan)
+
+        if np.min(np.diff(xi)) < 0 or np.min(np.diff(yi)) < 0:
+            raise ValueError("Output grid defined by xi,yi must be monotone "
+                             "increasing")
+
+        # Allocate array for output (buffer will be overwritten by natgridd)
+        zi = np.empty((yi.shape[0], xi.shape[0]), np.float64)
+
+        # Natgrid requires each array to be contiguous rather than e.g. a view
+        # that is a non-contiguous slice of another array.  Use numpy.require
+        # to deal with this, which will copy if necessary.
+        x = np.require(x, requirements=['C'])
+        y = np.require(y, requirements=['C'])
+        z = np.require(z, requirements=['C'])
+        xi = np.require(xi, requirements=['C'])
+        yi = np.require(yi, requirements=['C'])
+        _natgrid.natgridd(x, y, z, xi, yi, zi)
+
+        # Mask points on grid outside convex hull of input data.
+        if np.any(np.isnan(zi)):
+            zi = np.ma.masked_where(np.isnan(zi), zi)
+        return zi
+    else:
+        # Linear interpolation performed using a matplotlib.tri.Triangulation
+        # and a matplotlib.tri.LinearTriInterpolator.
+        from .tri import Triangulation, LinearTriInterpolator
+        triang = Triangulation(x, y)
+        interpolator = LinearTriInterpolator(triang, z)
+        return interpolator(xi, yi)
+
 
 ##################################################
 # Linear interpolation algorithms