NetCDF data loader
        
        :parameter filename: file name
        :parameter params: a list of variables to load (default : load ALL variables).
        
        :parameter depthrange: if a depth dimension is found, subset along this dimension.
        :parameter timerange: if a time dimension is found, subset along this dimension. 

        .. note:: using :attr:`altimetry.tools.nctools.limit` allows subsetting to a given region.

        :parameter kwargs: additional arguments for subsetting along given dimensions.
        
        .. note:: You can index along any dimension by providing the name of the dimensions to subsample along. Values associated to the provided keywords should be a length 2 or 3 tuple (min,max,<step>) (cf. :func:`altimetry.data.nctools.load_ncVar`).
        
        :keyword output_is_dict: data structures are dictionnaries (eg. my_hydro_data.variable['data']). If false uses an object with attributes (eg. my_hydro_data.variable.data). 
        
        :return  {type:dict} outStr: Output data structure containing all recorded parameters as specificied by NetCDF file PARAMETER list.
        
        :author: Renaud Dussurget

def load(self, filename, params=None, force=False, depthrange=None, timerange=None, output_is_dict=True, **kwargs):
        """
        NetCDF data loader
        
        :parameter filename: file name
        :parameter params: a list of variables to load (default : load ALL variables).
        
        :parameter depthrange: if a depth dimension is found, subset along this dimension.
        :parameter timerange: if a time dimension is found, subset along this dimension. 

        .. note:: using :attr:`altimetry.tools.nctools.limit` allows subsetting to a given region.

        :parameter kwargs: additional arguments for subsetting along given dimensions.
        
        .. note:: You can index along any dimension by providing the name of the dimensions to subsample along. Values associated to the provided keywords should be a length 2 or 3 tuple (min,max,<step>) (cf. :func:`altimetry.data.nctools.load_ncVar`).
        
        :keyword output_is_dict: data structures are dictionnaries (eg. my_hydro_data.variable['data']). If false uses an object with attributes (eg. my_hydro_data.variable.data). 
        
        :return  {type:dict} outStr: Output data structure containing all recorded parameters as specificied by NetCDF file PARAMETER list.
        
        :author: Renaud Dussurget
        

        """
        
        if (params is not None) & isinstance(params,str): params=[params]
        
        #Open file
        self._filename = filename
        try:
            ncf = ncfile(self._filename, "r")
        except Exception,e:
            warn(repr(e),stacklevel=2)
            return {}
        
        #Load global attributes
        akeys = ncf.ncattrs()
        attrStr=OrderedDict()
        for A in akeys : attrStr.update({A:ncf.getncattr(A)})
     
        #Get list of recorded parameters:
        dum = ncf.variables.keys()
        nparam = np.shape(dum)[0]
    #        par_list=np.array([''.join(ncf.variables.keys()[0,i,:].compressed()) for i in np.arange(nparam)])
        par_list = np.array(['{0}'.format(v) for v in ncf.variables.keys()])
        
        #remove empty items and update nparam
        par_list = par_list.compress([len(par) != 0 for par in par_list])
        nparam = par_list.size
        
        if nparam == 0 : self.Error('File has no data ({0})'.format(self._filename))
        
        #Get dimensions
        ncdimlist = np.array(['{0}'.format(d) for d in ncf.dimensions.keys()])
        ndims = len(ncdimlist)
        dimStr = OrderedDict()
        dimStr.update({'_ndims':ndims})
    
        if ndims == 0 : self.Error('File has no dimensions ({0})'.format(self._filename))
        
        #Check for the presence of strategic dimensions
        checkedDims = np.array(['lon', 'lat', 'time', 'depth'])
        existDim = -np.ones(4,dtype=int)
        if not self.use_local_dims :
            for i,d in enumerate(ncdimlist) :
                if ( (d.lower().startswith('lon')) | (d.lower().find('longitude') != -1) ) & (d.find('LatLon') ==-1) : existDim[0]=i
                if ( (d.lower().startswith('lat')) | (d.lower().find('latitude') != -1) ) & (d.find('LatLon') ==-1): existDim[1]=i
                if (d.lower().startswith('time')) | (d.lower().startswith('date')) : existDim[2]=i
                if (d.lower().startswith('lev')) | (d.lower().startswith('dep')) : existDim[3]=i
                    
#            existDim[0] = np.where([d.lower().startswith('lon') | (d.lower().find('longitude') != -1) for d in ncdimlist])[0]
#            existDim[1] = np.where([d.lower().startswith('lat') | (d.lower().find('latitude') != -1) for d in ncdimlist])[0]
#            existDim[2] = np.where([(d.lower().startswith('time')) | (d.lower().startswith('date')) for d in ncdimlist])[0]
#            existDim[3] = np.where([(d.lower().startswith('lev')) | (d.lower().startswith('dep')) for d in ncdimlist])[0]
        
        identified = existDim > -1
        
#        checkedDims[identified]=checkedDims[identified][existDim.compress(identified).astype(int)]
#        checkedDims=checkedDims[identified]
        
#        for cn, vn in enumerate(checkedDims) : dimStr.update({cn:len(ncf.dimensions[vn])})
        
        #Update dimension structure with identified dimensions
        #Load dimensional variables
        #TODO : Add scaling here in case...
        for i,d in enumerate(existDim) : 
            if identified[i] :
                dimStr.update({ncdimlist[d]:len(ncf.dimensions[ncdimlist[d]])}) #Append dimension
                cmd =  'load_ncVar(\'' + ncdimlist[d] + '\',nc=ncf)'
                self.message(4, 'loading : {0}={1}'.format(checkedDims[i],cmd))
                locals()[checkedDims[i]]=load_ncVar(ncdimlist[d], nc=ncf,**kwargs)

        missdims=set(ncdimlist)
        missdims.difference_update(ncdimlist[existDim[identified]])
        missdims=list(missdims)

        for i,d in enumerate(missdims) :
            dimStr.update({d:len(ncf.dimensions[d])})
            if ncf.variables.has_key(d) :
                cmd = 'load_ncVar(\'' + d + '\',nc=ncf)'
                self.message(4, 'loading : {0}={1}'.format(d,cmd))
                locals()[d]=load_ncVar(d, nc=ncf,**kwargs)
                        
            #If the variable associated to the dimension do not exist, generate it
            else :
                self.message(1, '[WARNING] Netcdf file not standard - creating data for {0} dimnsion'.format(d))
                ndim=len(ncf.dimensions[d])
                cmd = '=var'
                self.message(4, 'loading : {0}={1}'.format(d,cmd))
                locals()[d]={'_dimensions':{'_ndims':1,d:ndim}, 'data':np.arange(ndim)}
        
#        #Update dimension structure with identified dimensions
#        for cn, vn in zip(*(checkedDims[identified], ncdimlist[identified])) : dimStr.update({cn:len(ncf.dimensions[vn])})
#        for vn in ncdimlist[~identified] : dimStr.update({vn:len(ncf.dimensions[vn])})
        
        
        #Load dimensional variables
        #TODO : Add scaling here in case...
#        for cn, vn in zip(*(checkedDims[identified], ncdimlist[identified])) :
##            cmd = 'self.'+cn+'=ncf.variables[\''+vn+'\'][:]'
#            cmd = cn + '=load_ncVar(\'' + vn + '\',nc=ncf)'
#            self.message(4, 'exec : ' + cmd)
#            exec(cmd)
#            
#        for vn in ncdimlist[~identified] :
##            cmd = 'self.'+vn+'=ncf.variables[\''+vn+'\'][:]'
#            cmd = vn + '=load_ncVar(\'' + vn + '\',nc=ncf)'
#            self.message(4, 'exec : ' + cmd)
#            exec(cmd)
            
        #Update dimlist with dimensions present in the object
#        dimlist = np.append(checkedDims[identified], ncdimlist[~identified])

        dimlist=ncdimlist.copy()
#        dimlist[existDim[identified]]=checkedDims[identified]
        if identified.sum() > 0 : dimlist[existDim[identified]]=checkedDims[identified]
        else : dimlist = dimlist[[]]
        
#        for d in ncdimlist[checkedDims[identified]] :
#            if not d.startswith('_') : dimlist = np.append(dimlist,d)
#        dimlist=[(d if not d.startswith('_') else None) for d in dimStr.keys()]
        
        if params is not None :
            if force : par_list = [i.upper() for i in params]
            else :par_list = list(set(params).intersection(par_list))
        else : par_list = par_list.tolist()
        
        #remove dimensional  variable
        for d in ncdimlist[existDim[identified]] :
            par_list.pop(par_list.index(d))
        
        self.message(2, 'Recorded parameters : ' + str(nparam) + ' -> ' + str(par_list))
     
        
        #Extract within limits
        if (existDim[0] > -1) & (existDim[1] > -1):
            llind, flag = in_limits(lon['data'],lat['data'], limit=self.limit)
            if isinstance(flag,tuple) :
                lon['data'] = recale(lon['data'].compress(flag[0]),degrees=True)
                lon['_dimensions'][lon['_dimensions'].keys()[1]] = flag[0].sum()
                lat['data'] = lat['data'].compress(flag[1])
                lat['_dimensions'][lat['_dimensions'].keys()[1]] = flag[1].sum()
            else :
                lon['data'] = recale(lon['data'].compress(flag),degrees=True)
                lon['_dimensions'][lon['_dimensions'].keys()[1]] = flag.sum()
                lat['data'] = lat['data'].compress(flag)
                lat['_dimensions'][lat['_dimensions'].keys()[1]] = flag.sum()
            
            locals()[ncdimlist[existDim[0]]]=lon.copy()
            locals()[ncdimlist[existDim[1]]]=lat.copy()
            
            dimStr.update({ncdimlist[existDim[0]]:len(lon['data'])})
            dimStr.update({ncdimlist[existDim[1]]:len(lat['data'])})
#            self.message(4, 'self.lon & self.lat updated')
        
        if (existDim[2] > -1):
            if (timerange is not None) : timeflag = (time['data'] >= np.min(timerange)) & (time['data'] <= np.max(timerange))
            else : timeflag = np.ones(len(time['data']), dtype=bool)
            if timeflag.sum() == 0 : self.Error('No data within specified depth range (min/max = {0}/{1})'.format(np.min(time), np.max(time)))
            time['data'] = time['data'].compress(timeflag)
            time['_dimensions'][time['_dimensions'].keys()[1]] = timeflag.sum()
            locals()[ncdimlist[existDim[2]]]=time.copy()
            dimStr.update({ncdimlist[existDim[2]]:len(time['data'])})
#            self.message(4, 'self.lon & self.lat updated')
        
        #Extract within depth range
        if (existDim[3] > -1):
            if (depthrange is not None) : depthflag = (depth['data'] >= np.min(depthrange)) & (depth['data'] <= np.max(depthrange))
            else : depthflag = np.ones(len(depth['data']), dtype=bool)
            if depthflag.sum() == 0 : self.Error('No data within specified depth range (min/max = {0}/{1})'.format(np.min(depth), np.max(depth)))
            depth['data'] = depth['data'].compress(depthflag)
            depth['_dimensions'][depth['_dimensions'].keys()[1]] = depthflag.sum()
            locals()[ncdimlist[existDim[3]]]=depth.copy()
            dimStr.update({ncdimlist[existDim[3]]:len(depth['data'])})
        
        #Create output data structure
        outStr = OrderedDict()
        outStr.update({'_dimensions':dimStr})
        outStr.update({'_attributes':attrStr})
        
        if (existDim[0] > -1) : outStr.update({ncdimlist[existDim[0]]:lon})
        if (existDim[1] > -1) : outStr.update({ncdimlist[existDim[1]]:lat})
        if (existDim[2] > -1) : outStr.update({ncdimlist[existDim[2]]:time})
        if (existDim[3] > -1) : outStr.update({ncdimlist[existDim[3]]:depth})
        
        
        #Update object with remaining variables
        for d in dimlist.compress([not outStr.has_key(f) for f in dimlist]) :
#            cmd = 'outStr.update({\''+d+'\':'+d+'[\'data\']})'
            cmd = 'outStr.update({\''+d+'\':'+d+'})'
            self.message(4, 'exec : '+cmd)
            exec(cmd)
        
        ncdimStr=outStr.copy()
        #Get dimension lengths
        shape=()
        for d in dimlist: shape += np.shape(locals()[d]['data'])
        
        ndims = np.size(shape)
     
#        #Create dimension structure
#        curDim = [str(dimname) for dimname in dimStr.keys()[1:]] #[str(dimname) for dimname in ncf.variables['LONGITUDE'].dimensions]
#        curDimval = [dimStr[dim] for dim in curDim] #[len(ncf.dimensions[dimname]) for dimname in curDim]
#   
#        outStr={'_dimensions':{'_ndims':ndims,'nbpoints':sz[0]},'lon':lon,'lat':lat,'date':date}
        
#        for d in dimlist : outStr.update({d:self.__dict__[d]})

        #Sort NCDIMLIST to match DIMLIST
#        ncdimlist[np.sort(existDim.astype(np.int)[identified])]=ncdimlist[existDim[identified].tolist()]

        #Setup kwargs with current dimensionnal properties
        for d, ncd in zip(*(dimlist,ncdimlist)):
            if not kwargs.has_key(ncd) :
                if kwargs.has_key(d) :
                    
                    kwargs.update({ncd:kwargs[d]})
                    del kwargs[d]
                else :
                    dvar=ncdimStr[d]['data']
                    if isinstance(dvar,np.ma.masked_array) : kwargs.update({ncd:(np.nanmin(dvar.data),np.nanmax(dvar.data))})
                    else : kwargs.update({ncd:(np.nanmin(dvar),np.nanmax(dvar))})
#            else :
#                outStr['NbLatitudes']['data']        
        for param in par_list :
#            dumVar = load_ncVar(param,  nc=ncf, lon=llind[0], lat=llind[1], time=np.arange(len(time['data'])).compress(timeflag),**kwargs) #Load variables
#            dumVar = load_ncVar(param,  nc=ncf, longitude=(self.limit[1],self.limit[3]), latitude=(self.limit[0],self.limit[2]), time=(self.time.min(),self.time.max()),**kwargs) #Load variables
            
            
            dumVar = load_ncVar(param,  nc=ncf, **kwargs) #Load variables
#            dimStr = dumVar['_dimensions']
            
            #update dimensions
#            curDim = [str(dimname) for dimname in dimStr.keys()[1:]] #[str(dimname) for dimname in ncf.variables['LONGITUDE'].dimensions]
#            curDimval = [dimStr[dim] for dim in curDim] #[len(ncf.dimensions[dimname]) for dimname in curDim]
            
#            curDim = dimlist[where_list(curDim, ncdimlist.tolist())] #Convert to object dimension names
#            curDim = dimlist[where_list(curDim, dimlist.tolist())] #Convert to object dimension names (???)
            
##            curDim = [str(dimname) for dimname in ncf.variables[param].dimensions]
##            curDimval = [len(ncf.dimensions[dimname]) for dimname in curDim]
#            flag = [(np.array(dimname) == outStr['_dimensions'].keys()).sum() == 0 for dimname in curDim] #find dimensions to update
#            dimUpdate = np.array(curDim).compress(flag)
#            for enum in enumerate(dimUpdate) : 
#                self.message(2, 'Appending dimensions {0}:{1} to dataStructure'.format(enum[1], np.array(curDimval).compress(flag)[enum[0]]))
#                outStr['_dimensions'].update({enum[1]:np.array(curDimval).compress(flag)[enum[0]]}) #Append new dimension
#                outStr['_dimensions']['_ndims'] += 1 #update dimension counts
            
#            cmd = 'dumStr = {\'' + param + '\':dumVar[\'data\']}'


            #Set list as variable with attributes
#            if (not output_is_dict):
#                var=dumVar.pop('data')
#                for k in dumVar.keys():
#                    setattr(var, k, dumVar.pop(k))
#                dumVar=var.copy()
            
            cmd = 'dumStr = {\'' + param + '\':dumVar}'
            self.message(4, 'exec : ' + cmd)
            exec(cmd)
            outStr.update(dumStr)
            
            #Update output dimensions with extracted dimensions
            for ddum in dumStr[param]['_dimensions'].keys()[1:] :
                if outStr['_dimensions'].get(ddum) != dumStr[param]['_dimensions'][ddum] : outStr['_dimensions'][ddum]=dumStr[param]['_dimensions'][ddum]
            
#            cmd = 'self.'+param+'='
        
        ncf.close()
        return outStr