How to use the sunpy.Grid function in sunpy

interpLinear method in sunpy needs to be called layer by layer
        """
        # Put the input data back into a 2D array (Nz,Nc)
        self.datatmp[self.mask3d]=data

        # Loop through layer by layer
        dataout = np.zeros_like(X)
        for kk in range(0,self.Nkmax):
            ind = operator.and_(k == kk,self.cellind!=-1)
            if any(ind):
                dataout[ind] = self.interpLinear(self.datatmp[kk,:],X[ind],Y[ind],self.cellind[ind],k=kk)

        return dataout

class interp2Dmesh(GridSearch,Grid):
    """
    2D interpolation class for an unstructured grid
    
    Uses knowledge of the mesh to efficiently update the velocities based on 
    the particle location.
    """
    def __init__(self,x,y,cells,nfaces,method='nearest',grdfile=None):
        self.method=method
        # Initialise the trisearch array
        GridSearch.__init__(self,x,y,cells,nfaces=nfaces,force_inside=True)

        #if self.method == 'linear':
        #    #Grid.__init__(self,grdfile)
	    #self.datatmp = np.zeros(mask.shape,dtype=np.double)

f=mlab.gcf()
        f.scene.background = (0.,0.,0.)
        d = mlab.pipeline.add_dataset(ug)
        h=mlab.pipeline.surface(d,colormap='gist_earth')
        mlab.colorbar(object=h,orientation='vertical')
        mlab.view(0,0)
        
        outfile = self.suntanspath+'/depths.png'
        f.scene.save(outfile)      
        print 'Figure saved to %s.'%outfile
        
        #mlab.show()


class AverageDepth(Grid):
    """
    Returns the average of all depths inside each cells
    """
    # Projection conversion info for input data
    convert2utm=False
    CS='NAD83'
    utmzone=15
    isnorth=True
    vdatum = 'MSL'

    def __init__(self,suntanspath,**kwargs):
        
        self.__dict__.update(kwargs)
        Grid.__init__(self,suntanspath)
        
        # Initialise the trisearch object

def GridParticles(grdfile,dx,dy,nz,xypoly=None,splitvec=1):
    """
    Returns the locations of particles on a regular grid inside of suntans grid

    Inputs:
    	grdfile - netcdf filename containing the suntans grid
	dx,dy - resolution in x and y component respectively.
	nz - number of particles in the vertical. Particles are arranged in sigma layers.
	xypoly - [optional] coordinates of an additional bounding polygon [Nx2 array]
    """
    
    # Load the suntans grid
    sun = Grid(grdfile)

    # Load a trisearch object    
    tri = GridSearch(sun.xp,sun.yp,sun.cells,nfaces=sun.nfaces,verbose=False)
    
    if xypoly == None:
        xlims = [sun.xlims[0],sun.xlims[1]]
        ylims = [sun.ylims[0],sun.ylims[1]]
    else:
        xlims = [xypoly[:,0].min(),xypoly[:,0].max()]
        ylims = [xypoly[:,1].min(),xypoly[:,1].max()]

    # Construct a 2D mesh of particles
    x = np.arange(xlims[0],xlims[1],dx)
    y = np.arange(ylims[0],ylims[1],dy)
    
    X,Y = np.meshgrid(x,y)

vname = 'Mesh2_sea_surface_elevation'
dimensions = ('nMesh2_face','nMesh2_data_time')
attributes = {
		'long_name': "sea surface elevation" ,\
        'units':'m',\
        'coordinates': "Mesh2_face_x Mesh2_face_y" \
        }

dtype = 'f8'		
untrim_ugrid.update({vname:{'dimensions':dimensions,'attributes':attributes,'dtype':dtype,'zlib':True,'complevel':1,'fill_value':fillval}})





class UNTRIMGrid(Grid):
    """
    UnTRIM grid class for loading from a netcdf file
    """ 
    def __init__(self,ncfile):
        self.ncfile=ncfile
        Grid.__init__(self,ncfile,gridvars=untrim_gridvars,griddims=untrim_griddims)

class UNTRIMSpatial(Spatial):
    """
    Class for handling UnTRIM netcdf hydrodynamic output
    """
    _FillValue = -9999
    def __init__(self,ncfile,**kwargs):
        Spatial.__init__(self,ncfile,gridvars=untrim_gridvars,griddims=untrim_griddims,**kwargs)

        # Make sure the number of faces array is correct

@author: mrayson
"""

from sunpy import Grid
from netCDF4 import Dataset
import getopt, sys, time
import numpy as np
from multiprocessing import Pool
import pdb

# Globals
gridvars=['suntans_mesh','cells','face','edges','neigh','grad','nfaces','mnptr','eptr','xv','yv','xp','yp','xe','ye','normal','n1','n2','df','dg','def','Ac','dz','z_r','z_w','Nk','Nke','dv','time']	

nowritevars = ['face','neigh','grad'] # These need special attention

class JoinSuntans(Grid):
    """
    Class for joining suntans NetCDF file
    """
    def __init__(self,suntanspath,basename,numprocs,outvars=None):
        tic=time.clock()
        
        print '########################################################'
        print '     Initializing python SUNTANS NetCDF joining script...'
        print '########################################################'
        
        # Step 1) Read the main grid
        #print 'Loading suntans grid points...'
        Grid.__init__(self,suntanspath)

        self.numprocs = numprocs
        self.suntanspath=suntanspath

def __init__(self,suntanspath,timeinfo,**kwargs):
        """
        Initialise boundary path.
        
        To load data directly from an existing file:
            Boundary('boundary_ncfile.nc',0)
        """
        self.__dict__.update(**kwargs)

        if os.path.isdir(suntanspath) or self.loadfromnc:
        
            self.suntanspath = suntanspath
            self.timeinfo = timeinfo
            
            self.grd = sunpy.Grid(suntanspath)
            
            self._loadBoundary()
            
            self.getTime()
            
            # Initialise the output arrays
            self.initArrays()
            
        else:
            print 'Loading boundary data from a NetCDF file...'
            self.infile = suntanspath
            print '\t%s'%self.infile
            self._loadBoundaryNC()

            self.tsec = self.ncTime()

def __call__(self,suntanspath,depthmax=0.0,scalefac=-1.0, interpnodes=True):
        
        self.suntanspath=suntanspath
        
        # Initialise the interpolation points
        print 'Loading suntans grid points...'
        self.grd = sunpy.Grid(self.suntanspath)


        if interpnodes:
            print 'Interpolating depths onto nodes and taking min...'
            self.xy = np.column_stack((self.grd.xp,self.grd.yp))


        else:
            print 'Interpolating depths straight to cell centres...'
            self.xy = np.column_stack((self.grd.xv,self.grd.yv))


        # Initialise the Interpolation class
        print 'Building interpolant class...'
        self.F = interpXYZ(self.indata.XY,self.xy,method=self.interpmethod,NNear=self.NNear,\
                p=self.p,varmodel=self.varmodel,nugget=self.nugget,sill=self.sill,vrange=self.vrange)