Create numpy 2-D "meshgrid" from 3+ columns in a Pandas DataFrame

    Arguments:
      df (DataFrame): Must have 3 or 4 columns of numerical data

    Returns:
      OrderedDict: column labels from the data frame are the keys, values are 2-D matrices
        All matrices have shape NxM, where N = len(set(df.iloc[:,0])) and M = len(set(df.iloc[:,1]))

    >>> pandas_mesh(pd.DataFrame(np.arange(18).reshape(3,6),
    ...                          columns=list('ABCDEF'))).values()  # doctest: +NORMALIZE_WHITESPACE
    [array([[ 0,  6, 12],
            [ 0,  6, 12],
            [ 0,  6, 12]]),
     array([[ 1,  1,  1],
            [ 7,  7,  7],
            [13, 13, 13]]),
     array([[  2.,  nan,  nan],
            [ nan,   8.,  nan],
            [ nan,  nan,  14.]]),
     array([[  3.,  nan,  nan],
            [ nan,   9.,  nan],
            [ nan,  nan,  15.]]),
     array([[  4.,  nan,  nan],
            [ nan,  10.,  nan],
            [ nan,  nan,  16.]]),
     array([[  5.,  nan,  nan],
            [ nan,  11.,  nan],
            [ nan,  nan,  17.]])]

def pandas_mesh(df):
    """Create numpy 2-D "meshgrid" from 3+ columns in a Pandas DataFrame

    Arguments:
      df (DataFrame): Must have 3 or 4 columns of numerical data

    Returns:
      OrderedDict: column labels from the data frame are the keys, values are 2-D matrices
        All matrices have shape NxM, where N = len(set(df.iloc[:,0])) and M = len(set(df.iloc[:,1]))

    >>> pandas_mesh(pd.DataFrame(np.arange(18).reshape(3,6),
    ...                          columns=list('ABCDEF'))).values()  # doctest: +NORMALIZE_WHITESPACE
    [array([[ 0,  6, 12],
            [ 0,  6, 12],
            [ 0,  6, 12]]),
     array([[ 1,  1,  1],
            [ 7,  7,  7],
            [13, 13, 13]]),
     array([[  2.,  nan,  nan],
            [ nan,   8.,  nan],
            [ nan,  nan,  14.]]),
     array([[  3.,  nan,  nan],
            [ nan,   9.,  nan],
            [ nan,  nan,  15.]]),
     array([[  4.,  nan,  nan],
            [ nan,  10.,  nan],
            [ nan,  nan,  16.]]),
     array([[  5.,  nan,  nan],
            [ nan,  11.,  nan],
            [ nan,  nan,  17.]])]
    """
    xyz = [df[c].values for c in df.columns]
    index = pd.MultiIndex.from_tuples(zip(xyz[0], xyz[1]), names=['x', 'y'])
    # print(index)
    series = [pd.Series(values, index=index) for values in xyz[2:]]
    # print(series)
    X, Y = np.meshgrid(sorted(list(set(xyz[0]))), sorted(list(set(xyz[1]))))
    N, M = X.shape
    Zs = []
    # print(Zs)
    for k, s in enumerate(series):
        Z = np.empty(X.shape)
        Z[:] = np.nan
        for i, j in itertools.product(range(N), range(M)):
            Z[i, j] = s.get((X[i, j], Y[i, j]), np.NAN)
        Zs += [Z]
    return OrderedDict((df.columns[i], m) for i, m in enumerate([X, Y] + Zs))