def get_tindex(t,start_date,end_date,stride=None):
        
    tindex = []
    tindex.append(date2index(start_date,t,select='before'))
    tindex.append(date2index(end_date,t,select='after') + 1)
    if stride is None:
        tindex.append(1)
    else:
        tindex.append(stride)
    return tindex

def surf_vel(x,y,url,date_mid=datetime.datetime.utcnow(),
    uvar='u',vvar='v',isurf_layer=0,lonvar='lon',latvar='lat',
    tvar='time',hours_ave=24,lon360=False,ugrid=False,lonlat_sub=1,time_sub=1):
            
    nc=netCDF4.Dataset(url)
    lon = nc.variables[lonvar][:]-360.*lon360
    lat = nc.variables[latvar][:]
    
    if ugrid:
        lon2d=lon
        lat2d=lat
    elif lon.ndim==1:
        # ai and aj are logical arrays, True in subset region
        igood = np.where((lon>=x.min()) & (lon<=x.max()))
        jgood = np.where((lat>=y.min()) & (lat<=y.max()))
        bi=np.arange(igood[0].min(),igood[0].max(),lonlat_sub)
        bj=np.arange(jgood[0].min(),jgood[0].max(),lonlat_sub)
        [lon2d,lat2d]=np.meshgrid(lon[bi],lat[bj]) 
    elif lon.ndim==2:
        igood=np.where(((lon>=x.min())&(lon<=x.max())) & ((lat>=y.min())&(lat<=y.max())))
        bj=np.arange(igood[0].min(),igood[0].max(),lonlat_sub)
        bi=np.arange(igood[1].min(),igood[1].max(),lonlat_sub)
        lon2d=nc.variables[lonvar][bj,bi]
        lat2d=nc.variables[latvar][bj,bi]
    else:
        print 'uh oh'
        
    #desired_stop_date=datetime.datetime(2011,9,9,17,00)  # specific time (UTC)
    desired_stop_date=date_mid+datetime.timedelta(0,3600.*hours_ave/2.)  
    istop = netCDF4.date2index(desired_stop_date,nc.variables[tvar],select='nearest')   
    actual_stop_date=netCDF4.num2date(nc.variables[tvar][istop],nc.variables[tvar].units)
    actual_date_mid_est=actual_stop_date-datetime.timedelta(0,3600.*hours_ave/2.+3600.*5)
    start_date=actual_stop_date-datetime.timedelta(0,3600.*hours_ave)
    istart = netCDF4.date2index(start_date,nc.variables[tvar],select='nearest')
    print(date_mid.strftime('Requested mid-date: %I:00 %p on %B %d, %Y'))     
    print(actual_date_mid_est.strftime('Returned mid-date (EST): %I:00 %p on %B %d, %Y')) 
    print(start_date.strftime('start: %I:00 %p on %B %d, %Y'))   
    print(actual_stop_date.strftime('stop: %I:00 %p on %B %d, %Y'))
    if ugrid:
        u1=np.mean(nc.variables[uvar][istart:istop:time_sub,isurf_layer,:],axis=0)
        v1=np.mean(nc.variables[vvar][istart:istop:time_sub,isurf_layer,:],axis=0)
    else:
        print('reading u...')
        u1=np.mean(nc.variables[uvar][istart:istop:time_sub,isurf_layer,bj,bi],axis=0)
        print('reading v...')
        v1=np.mean(nc.variables[vvar][istart:istop:time_sub,isurf_layer,bj,bi],axis=0)

    xx2,yy2=np.meshgrid(x,y)
    ui=scipy.interpolate.griddata((lon2d.flatten(),lat2d.flatten()),u1.flatten(),(xx2,yy2),method='linear',fill_value=0.0)
    vi=scipy.interpolate.griddata((lon2d.flatten(),lat2d.flatten()),v1.flatten(),(xx2,yy2),method='linear',fill_value=0.0)
    ui[np.isnan(ui)]=0.0
    vi[np.isnan(vi)]=0.0

    
    return ui,vi,actual_date_mid_est

def get_nam_ts(url,vname,start=None,stop=None,j=None,i=None):
    nc = netCDF4.Dataset(url)
    ncv = nc.variables
    time_var = ncv['time']
    dtime = netCDF4.num2date(time_var[:],time_var.units)
    istart = netCDF4.date2index(start,time_var,select='nearest')
    istop = netCDF4.date2index(stop,time_var,select='nearest')
    var = ncv[vname]
    v = var[istart:istop,j,i]
    tim = dtime[istart:istop]
    return v,tim

def ndbc2df(collector, ndbc_id):
    """
    Ugly hack because `collector.raw(responseFormat="text/csv")`
    Usually times out.

    """
    from netCDF4 import MFDataset, date2index, num2date
    # FIXME: Only sea_water_temperature for now.
    if len(collector.variables) > 1:
        msg = "Expected only 1 variables to download, got {}".format
        raise ValueError(msg(collector.variables))
    if collector.variables[0] == 'sea_water_temperature':
        columns = 'sea_water_temperature (C)'
        ncvar = 'sea_surface_temperature'
        data_type = 'stdmet'
        # adcp, adcp2, cwind, dart, mmbcur, ocean, oceansites, pwind,
        # swden, tao-ctd, wlevel, z-hycom
    else:
        msg = "Do not know how to download {}".format
        raise ValueError(msg(collector.variables))

    uri = 'http://dods.ndbc.noaa.gov/thredds/dodsC/data/{}'.format(data_type)
    url = ('%s/%s/' % (uri, ndbc_id))
    urls = url_lister(url)

    filetype = "*.nc"
    file_list = [filename for filename in fnmatch.filter(urls, filetype)]
    files = [fname.split('/')[-1] for fname in file_list]
    urls = ['%s/%s/%s' % (uri, ndbc_id, fname) for fname in files]

    if not urls:
        raise Exception("Cannot find data at {!r}".format(url))
    nc = MFDataset(urls)

    kw = dict(calendar='gregorian', select='nearest')
    time_dim = nc.variables['time']
    time = num2date(time_dim[:], units=time_dim.units,
                    calendar=kw['calendar'])

    idx_start = date2index(collector.start_time.replace(tzinfo=None),
                           time_dim, **kw)
    idx_stop = date2index(collector.end_time.replace(tzinfo=None),
                          time_dim, **kw)
    if idx_start == idx_stop:
        raise Exception("No data within time range"
                        " {!r} and {!r}".format(collector.start_time,
                                                collector.end_time))
    data = nc.variables[ncvar][idx_start:idx_stop, ...].squeeze()

    time_dim = nc.variables['time']
    time = time[idx_start:idx_stop].squeeze()
    df = pd.DataFrame(data=data, index=time, columns=[columns])
    df.index.name = 'date_time'
    return df

def surf_vel_roms(x,y,url,date_mid=datetime.datetime.utcnow,hours_ave=24,tvar='ocean_time',lonlat_sub=1,time_sub=6):
    #url = 'http://testbedapps-dev.sura.org/thredds/dodsC/alldata/Shelf_Hypoxia/tamu/roms/tamu_roms.nc'

    #url='http://tds.ve.ismar.cnr.it:8080/thredds/dodsC/field2_test/run1/his'
    #####################################################################################

    nc = netCDF4.Dataset(url)
    mask = nc.variables['mask_rho'][:]
    lon_rho = nc.variables['lon_rho'][:]
    lat_rho = nc.variables['lat_rho'][:]
    anglev = nc.variables['angle'][:]

    desired_stop_date = date_mid+datetime.timedelta(0,3600.*hours_ave/2.)  # specific time (UTC)
    istop = netCDF4.date2index(desired_stop_date,nc.variables[tvar],select='nearest')   
    actual_stop_date=netCDF4.num2date(nc.variables[tvar][istop],nc.variables[tvar].units)   
    start_date=actual_stop_date-datetime.timedelta(0,3600.*hours_ave)
    istart = netCDF4.date2index(start_date,nc.variables[tvar],select='nearest') 

    uvar='u'
    vvar='v'
    isurf_layer = -1
    print('reading u...')
    u=np.mean(nc.variables[uvar][istart:istop:time_sub,isurf_layer,:,:],axis=0)
    print('reading v...')
    v=np.mean(nc.variables[vvar][istart:istop:time_sub,isurf_layer,:,:],axis=0)
    print('done reading data...')
    u = roms_utils.shrink(u, mask[1:-1, 1:-1].shape)
    v = roms_utils.shrink(v, mask[1:-1, 1:-1].shape)

    u, v = roms_utils.rot2d(u, v, anglev[1:-1, 1:-1])


    # <codecell>

    lon=lon_rho[1:-1,1:-1]
    lat=lat_rho[1:-1,1:-1]


    # <codecell>


    # <codecell>

    xx2,yy2=np.meshgrid(x,y)
    print('interpolating u to uniform grid...')
    ui=scipy.interpolate.griddata((lon.flatten(),lat.flatten()),u.flatten(),(xx2,yy2),method='linear',fill_value=0.0)
    print('interpolating v to uniform grid...')
    vi=scipy.interpolate.griddata((lon.flatten(),lat.flatten()),v.flatten(),(xx2,yy2),method='linear',fill_value=0.0)
    ui[np.isnan(ui)]=0.0
    vi[np.isnan(vi)]=0.0

    
    return ui,vi

def mean_precip(nc,bbox=None,start=None,stop=None):
    lon=nc.variables['lon'][:]
    lat=nc.variables['lat'][:]
    tindex0=netCDF4.date2index(start,nc.variables['time'],select='nearest')
    tindex1=netCDF4.date2index(stop,nc.variables['time'],select='nearest')
    bi=(lon>=box[0])&(lon<=box[2])
    bj=(lat>=box[1])&(lat<=box[3])
    p=nc.variables['precip_mean'][tindex0:tindex1,bj,bi]
    latmin=np.min(lat[bj])
    p=np.mean(p,axis=0)
    lon=lon[bi]
    lat=lat[bj]
    return p,lon,lat

    def test_select_dummy(self):
        nutime = self.TestTime(datetime(1950, 1, 1), 366, 24, "hours since 1400-01-01", "standard")

        dates = [datetime(1950, 1, 2, 6), datetime(1950, 1, 3), datetime(1950, 1, 3, 18)]

        t = date2index(dates, nutime, select="before")
        assert_equal(t, [1, 2, 2])

        t = date2index(dates, nutime, select="after")
        assert_equal(t, [2, 2, 3])

        t = date2index(dates, nutime, select="nearest")
        assert_equal(t, [1, 2, 3])

    def test_nonuniform(self):
        """Check that the fallback mechanism works. """
        nutime = self.TestTime(datetime(1950, 1, 1), 366, 24, "hours since 1900-01-01", "standard")

        # Let's remove the second entry, so that the computed stride is not
        # representative and the bisection method is needed.
        nutime._data = nutime._data[numpy.r_[0, slice(2, 200)]]

        t = date2index(datetime(1950, 2, 1), nutime)
        assert_equal(t, 30)

        t = date2index([datetime(1950, 2, 1), datetime(1950, 2, 3)], nutime)
        assert_equal(t, [30, 32])

    def test_select_dummy(self):
        nutime = self.TestTime(datetime(1950, 1, 1), 366, 24,
                               'hours since 1400-01-01', 'standard')

        dates = [datetime(1950, 1, 2, 6), datetime(
            1950, 1, 3), datetime(1950, 1, 3, 18)]

        t = date2index(dates, nutime, select='before')
        assert_equal(t, [1, 2, 2])

        t = date2index(dates, nutime, select='after')
        assert_equal(t, [2, 2, 3])

        t = date2index(dates, nutime, select='nearest')
        assert_equal(t, [1, 2, 3])

def getFVCOM_bottom_tempsalt_netcdf(lati,loni,starttime,endtime,layer,vname):#vname='temp'or'salinity'
        '''
        Function written by Yacheng Wang
        generates model data as a DataFrame
        according to time and local position
        different from getFVCOM_bottom_temp:
        this function only return time-temp dataframe and ues netcdf4
        getFVCOM_bottom_temp return depth and temp
        '''
        urlfvcom = 'http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3'
        nc = netCDF4.Dataset(urlfvcom)
        nc.variables
        lat = nc.variables['lat'][:]
        lon = nc.variables['lon'][:]
        times = nc.variables['time']
        jd = netCDF4.num2date(times[:],times.units)
        var = nc.variables[vname]
        inode = nearlonlat(lon,lat,loni,lati)
        modindex=netCDF4.date2index([starttime.replace(tzinfo=None),endtime.replace(tzinfo=None)],times,select='nearest')
        modtso = pd.DataFrame()
        # CHUNK THROUGH 'XDAYS' AT A TIME SINCE IT WAS HANGING UP OTHERWISE
        xdays=100
        for k in range(0,(endtime-starttime).days*24,xdays*24):
          print 'Generating a dataframe of model data requested from '+str(k)+' to ',str(k+xdays*24)
          #modtso=pd.DataFrame(var[modindex[0]:modindex[1],layer,inode],index=jd[modindex[0]:modindex[1]])
          modtso1=pd.DataFrame(var[modindex[0]+k:modindex[0]+k+xdays*24,layer,inode],index=jd[modindex[0]+k:modindex[0]+k+xdays*24])
          modtso=pd.concat([modtso,modtso1])
        return modtso

def main(O3File, MeteoMask):
#    
#    

    #If a file with data on ozone profiles already contains information on the 
    #height of the tropopause - just updated it, otherwise - create a new variable.
    FO3 = ncdf.Dataset(O3File,'r+')
    FMeteo = ncdf.MFDataset(MeteoMask)
    
    if not 'HTropo' in FO3.variables:
        var=FO3.createVariable('HTropo','float32',('Time',), zlib=True, complevel=9, fill_value=np.nan)
        var.units = 'm.'
        var.description = 'Height of the tropopause.'
    else:
        var=FO3.variables['HTropo']
    
    TimeO3 = FO3.variables['Time']
    TimeMet = ncdf.MFTime(FMeteo.variables['Time'])
    
    dtTimeO3 = ncdf.num2date(TimeO3, TimeO3.units, TimeO3.calendar)
    idxMet = ncdf.date2index(dtTimeO3, TimeMet, calendar=TimeMet.calendar, select='nearest')
    HTropo = FMeteo.variables['HTropo'][idxMet,0]
    
    var[...] = HTropo
    
    
    
    FO3.close()
    FMeteo.close()
    return 0

    def get_time_variable_overlap(self, dates):
        """Figure out if a new date array has a overlap with the already existing time
        variable.

        Return the index of the existing time variable where the new dates
        should be located.

        At the moment this only handles cases where all dates are new or none
        are new.

        Parameters
        ----------
        dates: list
            list of datetime objects


        Returns
        -------
        indexes: np.ndarray
           Array of indexes that overlap

        """
        timevar = self.dataset.variables[self.time_var]
        if timevar.size == 0:
            indexes = np.array([0])
        else:
            try:
                indexes = netCDF4.date2index(
                    dates, timevar)
            except ValueError:
                indexes = np.array([timevar.size])

        return indexes

def getFVCOM_bottom_tempsaltvel_netcdf(lati,loni,starttime,endtime,layer,vname):#vname='temp'or'salinity'
        '''
        Function written by Yacheng Wang
        generates model data as a DataFrame
        according to time and local position
        different from getFVCOM_bottom_temp:
        this function only return time-temp dataframe and ues netcdf4
        getFVCOM_bottom_temp return depth and temp
        '''
        urlfvcom = 'http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3'
        nc = netCDF4.Dataset(urlfvcom)
        nc.variables
        times = nc.variables['time']
        jd = netCDF4.num2date(times[:],times.units)
        if vname=='vel':
            u = nc.variables['u']
            v = nc.variables['v']
            lat = nc.variables['latc'][:]
            lon = nc.variables['lonc'][:]
        else:
            var=nc.variables[vname]
            lat = nc.variables['lat'][:]
            lon = nc.variables['lon'][:]
        inode = nearlonlat(lon,lat,loni,lati)
        modindex=netCDF4.date2index([starttime.replace(tzinfo=None),endtime.replace(tzinfo=None)],times,select='nearest')
        print modindex
        #print [u[modindex[0]:modindex[1],layer,inode][0],v[modindex[0]:modindex[1],layer,inode][0]]
        if vname=='vel':
            modtso=pd.DataFrame(np.array([u[modindex[0]:modindex[1],layer,inode],v[modindex[0]:modindex[1],layer,inode]]).T,index=jd[modindex[0]:modindex[1]])
            #modtso=pd.DataFrame(np.array([u[modindex[0],layer,inode],v[modindex[0],layer,inode]]).T,index=jd[modindex[0]])
        else:
            modtso=pd.DataFrame(var[modindex[0]:modindex[1],layer,inode],index=jd[modindex[0]:modindex[1]])
        return modtso

    def runTest(self):
        # Get the real dates
        # skip this until cftime pull request #55 is in a released
        # version (1.0.1?). Otherwise, fix for issue #808 breaks this
        if parse_version(cftime.__version__) >= parse_version('1.0.1'):
            dates = []
            for file in self.files:
                f = Dataset(file)
                t = f.variables['time']
                dates.extend(num2date(t[:], t.units, t.calendar))
                f.close()

        # Compare with the MF dates
        f = MFDataset(self.files,check=True)
        t = f.variables['time']
        mfdates = num2date(t[:], t.units, t.calendar)

        T = MFTime(t)
        assert_equal(len(T), len(t))
        assert_equal(T.shape, t.shape)
        assert_equal(T.dimensions, t.dimensions)
        assert_equal(T.typecode(), t.typecode())
        # skip this until cftime pull request #55 is in a released
        # version (1.0.1?). Otherwise, fix for issue #808 breaks this
        if parse_version(cftime.__version__) >= parse_version('1.0.1'):
            assert_array_equal(num2date(T[:], T.units, T.calendar), dates)
        assert_equal(date2index(datetime.datetime(1980, 1, 2), T), 366)
        f.close()

 def test_singletime(self):
     # issue 215 test (date2index with time variable length == 1)
     f = Dataset(self.file)
     time2 = f.variables['time2']
     result_index = date2index(self.first_timestamp, time2, select="exact")
     assert_equal(result_index, 0)
     f.close()

def get_roms(url, time_slice, n=3):
    url = parse_url(url)
    with Dataset(url) as nc:
        ncv = nc.variables
        time = ncv['ocean_time']
        tidx = date2index(time_slice, time, select='nearest')
        time = num2date(time[tidx], time.units, time.calendar)

        mask = ncv['mask_rho'][:]
        lon_rho = ncv['lon_rho'][:]
        lat_rho = ncv['lat_rho'][:]
        anglev = ncv['angle'][:]

        u = ncv['u'][tidx, -1, ...]
        v = ncv['v'][tidx, -1, ...]

        u = shrink(u, mask[1:-1, 1:-1].shape)
        v = shrink(v, mask[1:-1, 1:-1].shape)

        u, v = rot2d(u, v, anglev[1:-1, 1:-1])

        lon = lon_rho[1:-1, 1:-1]
        lat = lat_rho[1:-1, 1:-1]

        u, v = u[::n, ::n], v[::n, ::n]
        lon, lat = lon[::n, ::n], lat[::n, ::n]

        u = ma.masked_invalid(u)
        v = ma.masked_invalid(v)
    return dict(lon=lon, lat=lat, u=u, v=v, time=time)

def netcdf2PCRobj(ncFile,varName,dateInput):
    # EHS (04 APR 2013): To convert netCDF (tss) file to PCR file.
    # The cloneMap is globally defined (outside this method).
    
    # Get netCDF file and variable name:
    f = nc.Dataset(ncFile)
    varName = str(varName)

    # date
    date = dateInput
    if isinstance(date, str) == True: date = \
                    datetime.datetime.strptime(str(date),'%Y-%m-%d') 
    date = datetime.datetime(date.year,date.month,date.day)
    
    # time index (in the netCDF file)
    nctime = f.variables['time']  # A netCDF time variable object.
    idx = nc.date2index(date, nctime, calendar=nctime.calendar, \
                                                 select='exact') 
    
    # convert to PCR object and close f
    outPCR = pcr.numpy2pcr(pcr.Scalar,(f.variables[varName][idx].data), \
                             float(f.variables[varName]._FillValue))
    f.close(); f = None ; del f
    # PCRaster object
    return (outPCR)

def get_ncfiles_catalog(station_id, jd_start, jd_stop):
    station_name = station_id.split(":")[-1]
    uri = 'http://dods.ndbc.noaa.gov/thredds/dodsC/data/stdmet'
    url = ('%s/%s/' % (uri, station_name))
    urls = _url_lister(url)
    filetype = "*.nc"
    file_list = [filename for filename in fnmatch.filter(urls, filetype)]
    files = [fname.split('/')[-1] for fname in file_list]
    urls = ['%s/%s/%s' % (uri, station_name, fname) for fname in files]

    try:
        nc = MFDataset(urls)

        time_dim = nc.variables['time']
        calendar = 'gregorian'
        idx_start = date2index(jd_start, time_dim, calendar=calendar,
                               select='nearest')
        idx_stop = date2index(jd_stop, time_dim, calendar=calendar,
                              select='nearest')

        dir_dim = np.array(nc.variables['wind_dir'][idx_start:idx_stop, 0, 0], dtype=float)
        speed_dim = np.array(nc.variables['wind_spd'][idx_start:idx_stop, 0, 0])
        # Replace fill values with NaN
        speed_dim[speed_dim == nc.variables['wind_spd']._FillValue] = np.nan

        if dir_dim.ndim != 1:
            dir_dim = dir_dim[:, 0]
            speed_dim = speed_dim[:, 0]
        time_dim = nc.variables['time']
        dates = num2date(time_dim[idx_start:idx_stop],
                         units=time_dim.units,
                         calendar='gregorian').squeeze()
        mask = np.isfinite(speed_dim)

        data = dict()
        data['wind_speed (m/s)'] = speed_dim[mask]
        data['wind_from_direction (degree)'] = dir_dim[mask]
        time = dates[mask]

        # columns = ['wind_speed (m/s)',
        #            'wind_from_direction (degree)']
        df = DataFrame(data=data, index=time)
        return df
    except Exception as e:
        print str(e)
        df = DataFrame()
        return df

 def test_failure(self):
     nutime = self.TestTime(datetime(1950, 1, 1), 366, 24, "hours since 1900-01-01", "standard")
     try:
         t = date2index(datetime(1949, 2, 1), nutime)
     except ValueError:
         pass
     else:
         raise ValueError("This test should have failed.")

 def test_issue444(self):
     # make sure integer overflow not causing error in
     # calculation of nearest index when sum of adjacent
     # time values won't fit in 32 bits.
     ntimes = 20
     f = Dataset(self.file, 'r')
     query_time = datetime(2037, 1, 3, 21, 12)
     index = date2index(query_time, f.variables['time3'], select='nearest')
     assert(index == 11)
     f.close()