def read_ncdf(name,path='/work/gg0877/KST/tide_gaughes/ncdf',origin='2012-01-01 00:00:00'):
  ncname = '_'.join(word[0].upper()+word[1:] for word in name.split('_'))
  if path[-6:]=='hawaii':
    nc = netCDF4.MFDataset(path+'/'+ncname+'*.nc')
  else:
    nc = netCDF4.Dataset(path+'/'+ncname+'.nc')
  ncv = nc.variables
  ot = utime('seconds since '+origin)
  ut = utime(ncv['time'].units)
  time = ot.date2num(ut.num2date(ncv['time'][:]))
  print(ncv.keys())
  scale = {'cm':100.,'m':1.0,'millimeters':1000.,'mm':1000.}

  if 'elev' in ncv:
    if 'units' in ncv['elev'].ncattrs():
      scale_factor = scale[ncv['elev'].units]
    else:
      if ncv['elev'][:].std()>100.:
        scale_factor=1000.
      elif ncv['elev'][:].std()>10.:
        scale_factor=100.
      else:
        scale_factor=1.0
    elev = ncv['elev'][:].squeeze()/scale_factor
  elif 'sea_surface_height_above_reference_level' in ncv:
    if 'units' in ncv['sea_surface_height_above_reference_level'].ncattrs():
      scale_factor = scale[ncv['sea_surface_height_above_reference_level'].units]
    else:
      scale_factor = 1000.
    elev = ncv['sea_surface_height_above_reference_level'][:].squeeze()/scale_factor
  nc.close()
  return time,elev

 def check(self, max_n_ticks, num1, num2):
     locator = NetCDFTimeDateLocator(max_n_ticks=max_n_ticks,
                                     calendar=self.calendar,
                                     date_unit=self.date_unit)
     utime = netcdftime.utime(self.date_unit, self.calendar)
     return locator.compute_resolution(num1, num2, utime.num2date(num1),
                                       utime.num2date(num2))

def read_ncom(flow_data):
    print '  > NCOM;'
    filename  = data_dir + '/ncom/ncom_glb_regp01_2012-3-4-5.nc'

    ncf   = netCDF4.Dataset(filename,'r')
    ncvar = ncf.variables
    time       = ncvar['time']
    utime      = netcdftime.utime(time.units)
    dates      = utime.num2date(time[:])
    sta_name   = filename.split('/')[-1][:3]
    
    lona,lata = np.meshgrid(ncvar['lon'][:],ncvar['lat'][:])
    dist = ((lona - lon_orig)**2 + (lata - lat_orig)**2)
    [i,j] = np.where(dist==dist.min())
    
    elev  = ncvar['surf_el'][:,i,j].flatten()
    data  = pd.DataFrame(data = elev, columns = ['elev'], index = dates)    
    data  = data.dropna()
    data  = data.resample('H')  #hourly mean
    
    date_tmp =  ps.datetime64todatetime(data.index.values)
    datenum  = [ps.datetime2datenum(datei) for datei in  date_tmp]
    
    flow_data[sta_name]['name_long'] = 'Elevation from NCOM model '
    flow_data[sta_name]['lat']       = lata[i,j].item()
    flow_data[sta_name]['lon']       = lona[i,j].item()
    flow_data[sta_name]['elev']      = data['elev'].values
    flow_data[sta_name]['date']      = date_tmp
    flow_data[sta_name]['datenum_mat']   = datenum

 def units_checker(units):
     try:
         u = utime(units)
         return u.unit_string
     except ValueError, TypeError:
         msg = 'Invalid time units format - Available format is "{}"'.format(TIME_UNITS_FORMAT)
         raise ArgumentTypeError(msg)

def read_Elgar_pressure_data(flow_data):
    print '   > Elgar pressure data;'
    dir1 = data_dir + '/elgar/*_filter_abs.nc'
    flist    = glob.glob(dir1)
    flist.sort()

    for filename in flist[:]:
        ncf   = netCDF4.Dataset(filename,'r')
        ncvar = ncf.variables
        time       = ncvar['time']
        utime      = netcdftime.utime(time.units)
        dates      = utime.num2date(time[:])
        sta_name   = filename.split('/')[-1][:3]
        if 'p' in sta_name or 'q' in sta_name:
            print '    > read  > Station name: ', sta_name
    
            water_depth = ncvar['water_depth'][:]
    
            data  = pd.DataFrame(data = water_depth, columns = ['water_depth'], index = dates)    
            data  = data.dropna()
            data  = data.resample('H')  #hourly mean
            
            date_tmp =  ps.datetime64todatetime(data.index.values)
            datenum  = [ps.datetime2datenum(datei) for datei in  date_tmp]

            
            flow_data[sta_name]['name_long'] = 'Bottom pressure '+sta_name
            flow_data[sta_name]['lat']       = ncvar['lat'][:]
            flow_data[sta_name]['lon']       = ncvar['lon'][:]
            flow_data[sta_name]['elev']      = (data['water_depth'] - data['water_depth'].mean()).values
            flow_data[sta_name]['date']      = date_tmp
            flow_data[sta_name]['datenum_mat']   = datenum

  def __init__(self,ncfile='na_aviso_2012.nc'):
    self.sshname='adt_mean'
    self.lonname='lon'
    self.latname='lat'
    self.timename='time'

    nc = netCDF4.Dataset(ncfile)
    sv = nc.variables
    self.ncv = sv
    xslice=slice(None)
    yslice=slice(None)
    self.lon = sv[self.lonname][xslice]
    self.lat = sv[self.latname][yslice]
    self.lon2,self.lat2 = np.meshgrid(self.lon,self.lat)
    self.time = sv[self.timename][:]
    self.timeunits = sv[self.timename].units
    self.ut = netcdftime.utime(self.timeunits)
    self.dates = self.ut.num2date(self.time)
    self.ssh = sv[self.sshname][yslice,xslice]

    print('  build spatial tree')
    self.mask = self.ssh.mask.squeeze()
    self.water = where(self.mask==False)
    vlon2 = self.lon2[self.water]
    vlat2 = self.lat2[self.water]
    self.tree = cKDTree(zip(vlon2,vlat2))

def get_list_dates_from_nc(nc, type_dates):
    
    '''
    Returns list of dates from NetCDF dataset.
    
    :param nc: NetCDF dataset
    :type nc: netCDF4.Dataset
    :param type_dates: type of dates ('dt' for datetime objects, 'num' for float objects) 
    :type type_dates: str

    :rtype: list of datetime/float 
    
    '''
    
    var_time = nc.variables['time']
    time_units = var_time.units # str (ex.: 'days since 1850-01-01 00:00:00')
    try:
        time_calend = var_time.calendar # str (ex.: 'standard'/'gregorian'/...)
    except:
        time_calend = 'gregorian'
    
    if type_dates == 'num':
        arr_dt = var_time[:]
        list_dt = arr_dt.tolist() # numpy array -> list
        
    if type_dates == 'dt':
        t = netcdftime.utime(time_units, time_calend) 
        arr_dt = t.num2date(var_time[:]) 
        list_dt = arr_dt.tolist() # numpy array -> list
    del arr_dt
    
    return list_dt

 def load_time(self):
   tnum=self.nc.variables['time'][-self.ntimes_max:]
   tunits=self.nc.variables['time'].units
   self.time=utime(tunits,calendar='standard').num2date(tnum)
   ind0=self.nc.variables['time'].size-self.ntimes_max
   ind1=self.nc.variables['time'].size
   self.time_inds=range(ind0,ind1)

def num2date(tnum,tunits,calendar='standard'):
  from netcdftime import utime

  #return utime(tunits,calendar).num2date(tnum)
  # for numpy 1.15.3 there are problems with masked arrays without mask!
  # for some reason, removing mask or dividing by 1 solves the problem!!!
  return utime(tunits,calendar).num2date(tnum/1)

def get_list_dates(ifile, type_dates):
    
    '''
    Returns list of dates from one file.
    
    :param ifile: NetCDF file
    :type ifile: str
    :param type_dates: type of dates ('dt' for datetime objects, 'num' for float objects) 
    :type type_dates: str

    :rtype: list of datetime/float 
    
    '''
    
    nc = Dataset(ifile, 'r')
    var_time = nc.variables['time']
    time_units = var_time.units # str (ex.: 'days since 1850-01-01 00:00:00')
    time_calend = var_time.calendar # str (ex.: 'standard'/'gregorian'/...)
    
    if type_dates == 'num':
        arr_dt = var_time[:]
        list_dt = arr_dt.tolist() # numpy array -> list
        
    if type_dates == 'dt':
        t = utime(time_units, time_calend) # <netcdftime.utime instance at 0xecae18>
        arr_dt = t.num2date(var_time[:]) # arr_dt: numpy array of dates datetime; var_time[:]: time values (ex.: [49323.5, 49353, 49382.5, ...])
        list_dt = arr_dt.tolist() # numpy array -> list
    del arr_dt
    
    nc.close()
    
    return list_dt

def _fix_netcdftime_datetime(dt, unit):
    # Fix netcdftime bug: netcdftime.datetime(2000, 03, 31).dayofyr
    # Happens to work with datetime objects too.
    fix_unit = "days since {}".format(str(dt))
    fix_utime = netcdftime.utime(fix_unit, calendar=unit.calendar)
    dt = fix_utime.num2date(0)
    return dt

def num2date(num, calend, units):
    '''
    type num: float date
    '''   
    t = utime(units, calend) 
    dt = t.num2date(num) 
    return dt

def date2num(dt, calend, units):
    '''
    type dt: datetime object
    '''
    t = utime(units, calend)
    dt_num = t.date2num(dt)
    return dt_num

def read_d3d_time(netcdf_vars):
    """
    Read time from netcdf file and return Datetime vector
    """
    utime = netcdftime.utime(netcdf_vars['time'].units)
    dates = utime.num2date(netcdf_vars['time'][:])
    return dates

 def extract_netcdf_time(cls, bandName, calendar):
     """Convert netcdf time to datetime using appropriate library"""
     from netcdftime import utime
     epoch, units = cls.extract_epoch_units(bandName)
     cdftime = utime(units, calendar)
     timestamps = cdftime.num2date([epoch])
     return timestamps[0]

 def setUp(self):
     self.cdftime_mixed = utime("hours since 0001-01-01 00:00:00")
     self.cdftime_mixed_tz = utime("hours since 0001-01-01 00:00:00 -06:00")
     self.cdftime_pg = utime("seconds since 0001-01-01 00:00:00", calendar="proleptic_gregorian")
     self.cdftime_noleap = utime("days since 1600-02-28 00:00:00", calendar="noleap")
     self.cdftime_leap = utime("days since 1600-02-29 00:00:00", calendar="all_leap")
     self.cdftime_360day = utime("days since 1600-02-30 00:00:00", calendar="360_day")
     self.cdftime_jul = utime("hours since 1000-01-01 00:00:00", calendar="julian")
     self.cdftime_iso = utime("seconds since 1970-01-01T00:00:00Z")

 def __new__(self, ncfile, name='time', units=None, calendar='standard'):
     self._nc = octant.io.Dataset(ncfile)
     data = self._nc.variables[name][:]
     data = data.view(time)
     if units == None:
         units = self._nc.variables[name].units
     data.utime = netcdftime.utime(units, calendar=calendar)
     return data

 def setUp(self):
     self.cdftime_mixed = utime('hours since 0001-01-01 00:00:00')
     self.cdftime_mixed_tz = utime('hours since 0001-01-01 00:00:00 -06:00')
     self.cdftime_pg = utime('seconds since 0001-01-01 00:00:00',
                       calendar='proleptic_gregorian')
     self.cdftime_noleap = utime('days since 1600-02-28 00:00:00',calendar='noleap')
     self.cdftime_leap = utime('days since 1600-02-29 00:00:00',calendar='all_leap')
     self.cdftime_360day = utime('days since 1600-02-30 00:00:00',calendar='360_day')
     self.cdftime_jul = utime('hours since 1000-01-01 00:00:00',calendar='julian')
     self.cdftime_iso = utime("seconds since 1970-01-01T00:00:00Z")

 def test_extract_conversion_factors_from_string(self):
     u = utime('hours since 1970-01-01 00:00:00')
     t_origin = api.Calendar(u.tzoffset).time(
         api.YMDhms(u.origin.year, u.origin.month, u.origin.day, u.origin.hour, u.origin.minute, u.origin.second))
     delta_t_dic = {'days': api.deltahours(24), 'hours': api.deltahours(1), 'minutes': api.deltaminutes(1)}
     delta_t = delta_t_dic[u.units]
     self.assertIsNotNone(u)
     self.assertEqual(delta_t, api.deltahours(1))
     self.assertEqual(t_origin, 0)

 def __new__(self, ncfile, name='time', units=None, calendar='standard'):
     raise DeprecationWarning, 'Use netCDF4.num2date instead.  Wrapper function in octant.roms.nc_time'
     self._nc = octant.io.Dataset(ncfile)
     data = self._nc.variables[name][:]
     data = data.view(time)
     if units == None:
         units = self._nc.variables[name].units
     data.utime = netcdftime.utime(units, calendar=calendar)
     return data