In [1]:

    

%run basics


    

In [2]:

    

ncFullName=qcio.get_filename_dialog(path="../../Sites")


    

In [3]:

    

netCDF4.default_encoding = 'latin-1'
ds = qcio.DataStructure()
# check to see if the requested file exists, return empty ds if it doesn't
if not qcutils.file_exists(ncFullName,mode="quiet"):
    log.error(' netCDF file '+ncFullName+' not found')
    raise Exception("nc_read_series: file not found")
# file probably exists, so let's read it
ncFile = netCDF4.Dataset(ncFullName,'r')
# now deal with the global attributes
gattrlist = ncFile.ncattrs()
if len(gattrlist)!=0:
    for gattr in gattrlist:
        ds.globalattributes[gattr] = getattr(ncFile,gattr)
    if "time_step" in ds.globalattributes: c.ts = ds.globalattributes["time_step"]
# get a list of the variables in the netCDF file (not their QC flags)
varlist = [x for x in ncFile.variables.keys() if "_QCFlag" not in x]
for ThisOne in varlist:
    # skip variables that do not have time as a dimension
    dimlist = [x.lower() for x in ncFile.variables[ThisOne].dimensions]
    if "time" not in dimlist: continue
    # create the series in the data structure
    ds.series[unicode(ThisOne)] = {}
    # get the data and the QC flag
    data,flag,attr = qcio.nc_read_var(ncFile,ThisOne)
    ds.series[ThisOne]["Data"] = data
    ds.series[ThisOne]["Flag"] = flag
    ds.series[ThisOne]["Attr"] = attr
ncFile.close()
# make sure all values of -9999 have non-zero QC flag
qcutils.CheckQCFlags(ds)
# put a list of Python datetime objects into ds
qcutils.get_datetimefromymdhms(ds)


    

In [11]:

    

nRecs = int(ds.globalattributes["nc_nrecs"])
ts = int(ds.globalattributes["time_step"])
ldt = ds.series["DateTime"]["Data"]
nc_time = netCDF4.date2num(ldt,"days since 1800-01-01 00:00:00.0",calendar="gregorian")
flag = numpy.zeros(nRecs,dtype=numpy.int32)
attr = qcutils.MakeAttributeDictionary(long_name="time",standard_name="time",units="days since 1800-01-01 00:00:00.0",
                                       calendar="gregorian")
qcutils.CreateSeries(ds,"time",nc_time,Flag=flag,Attr=attr)


    

In [32]:

    

def rounddt(dt,ts=30):
    dt += datetime.timedelta(minutes=int(ts/2))
    dt -= datetime.timedelta(minutes=dt.minute % int(ts),seconds=dt.second,microseconds=dt.microsecond)
    return dt


    

In [45]:

    

dt = numpy.array([(ldt[i]-ldt[i-1]).total_seconds() for i in range(1,len(ldt))])
idx = numpy.where(dt!=ts*60)[0]
if len(idx)!=0:
    print "Non-integral time steps found",len(idx),"times out of",str(nRecs)
    print "Maximum time step was",numpy.max(dt),"seconds, minimum time step was",numpy.min(dt)
    dt_diffs = numpy.array([(ldt[i]-rounddt(ldt[i],ts=ts)).total_seconds() for i in range(1,len(ldt))])
    print "Maximum drift was",numpy.max(dt_diffs),"seconds, minimum drift was",numpy.min(dt_diffs),"seconds"
    ans = raw_input("Do you want to [Q]uit, [I]nterploate or [R]ound? ")
    if ans.lower()=="q":
        print "Quiting ..."
        #sys.exit()
    if ans.lower()=="i": print "Interpolate"
    if ans.lower()=="r":
        print "Rounding to the nearest time step"
        ldt_rounded = [rounddt(dt,ts=ts) for dt in ldt]
        rdt = numpy.array([(ldt_rounded[i]-ldt_rounded[i-1]).total_seconds() for i in range(1,len(ldt))])
        print "Maximum time step is now",numpy.max(rdt),"seconds, minimum time step is now",numpy.min(rdt)


    

    




Non-integral time steps found 80 times out of 17520
Maximum time step was 1801.0 seconds, minimum time step was 1799.0
Maximum drift was 50.0 seconds, minimum drift was 0.0 seconds
Do you want to [Q]uit, [I]nterploate or [R]ound? r
Rounding to the nearest time step
Maximum time step is now 1800.0 seconds, minimum time step is now 1800.0

In [8]:

    

fig = plt.figure()
plt.plot(ldt[1:],dt,'b.')
plt.show()


    

In [34]:

    

# round the start and end datetimes to the nearest integral time steps
print ldt[0],ldt[-1]
sd = rounddt(ldt[0],ts=ts)
ed = rounddt(ldt[-1],ts=ts)
print sd,ed


    

    




2010-01-01 00:30:00 2011-01-01 00:00:50
2010-01-01 00:30:00 2011-01-01 00:00:00

In [ ]: