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

def plot2Ddata(xyz, data, vec=None, nx=100, ny=100, ax=None, mask=None, level=None, figname=None, ncontour=10, dataloc=False, contourOpts={}, clim=None):
    """
        TODOs:

            Write explanation for parameters

    """
    if ax is None:
        fig = plt.figure()
        ax = plt.subplot(111)

    xmin, xmax = xyz[:,0].min(), xyz[:,0].max()
    ymin, ymax = xyz[:,1].min(), xyz[:,1].max()
    x = np.linspace(xmin, xmax, nx)
    y = np.linspace(ymin, ymax, ny)
    X, Y = np.meshgrid(x, y)
    xy = np.c_[X.flatten(), Y.flatten()]
    if vec is None:
        F = LinearNDInterpolator(xyz[:,:2], data)
        DATA = F(xy)
        DATA = DATA.reshape(X.shape)
        cont = ax.contourf(X, Y, DATA, ncontour, **contourOpts)
        if level is not None:
            CS = ax.contour(X, Y, DATA, level, colors="k", linewidths=2)

    else:
        # Assume size of data is (N,2)
        datax = data[:,0]
        datay = data[:,1]
        Fx = LinearNDInterpolator(xyz[:,:2], datax)
        Fy = LinearNDInterpolator(xyz[:,:2], datay)
        DATAx = Fx(xy)
        DATAy = Fy(xy)
        DATA = np.sqrt(DATAx**2+DATAy**2).reshape(X.shape)
        DATAx = DATAx.reshape(X.shape)
        DATAy = DATAy.reshape(X.shape)
        cont = ax.contourf(X, Y, DATA, ncontour, **contourOpts)
        ax.streamplot(X, Y, DATAx, DATAy, color="w")
        if level is not None:
            CS = ax.contour(X, Y, DATA, level, colors="k", linewidths=2)

#         plt.clabel(CS, inline=1, fmt="%.1e")
    if dataloc:
        ax.plot(xyz[:,0], xyz[:,1], 'k.', ms=1)
    plt.gca().set_aspect('equal', adjustable='box')
    if figname:
        plt.axis("off")
        fig.savefig(figname, dpi=200)
    if level is  None:
        return True, cont, ax
    else:
        return True, cont, ax, CS

@@ -763,32 +763,143 @@
 eshape(X.shape)%0A
+        if clim is not None:%0A            DATA%5BDATA%3Cclim%5B0%5D%5D = clim%5B0%5D%0A            DATA%5BDATA%3Eclim%5B1%5D%5D = clim%5B1%5D%0A
         cont = a
@@ -938,32 +938,32 @@
  **contourOpts)%0A
-
         if level
@@ -1049,17 +1049,16 @@
 dths=2)%0A
-%0A
     else
@@ -1296,24 +1296,24 @@
 Ay = Fy(xy)%0A
-
         DATA
@@ -1351,32 +1351,143 @@
 eshape(X.shape)%0A
+        if clim is not None:%0A            DATA%5BDATA%3Cclim%5B0%5D%5D = clim%5B0%5D%0A            DATA%5BDATA%3Eclim%5B1%5D%5D = clim%5B1%5D%0A
         DATAx = 

Add clim for plotting utils