In [1]:
%run basics
%matplotlib
In [2]:
def round_datetime(ds,mode="nearest_timestep"):
"""
Purpose:
Round the series of Python datetimes to the nearest time based on mode
Usage:
qcutils.round_datetime(ds,mode=mode)
where;
mode = "nearest_second" rounds to the nearesy second
mode = "nearest_timestep" rounds to the nearest time step
Author: PRI
Date: February 2015
"""
# local pointer to the datetime series
ldt = ds.series["DateTime"]["Data"]
# check which rounding option has been chosen
if mode.lower()=="nearest_timestep":
# get the time step
if "time_step" in ds.globalattributes:
ts = int(ds.globalattributes["time_step"])
else:
ts = numpy.mean(qcutils.get_timestep(ds)/60)
ts = roundtobase(ts,base=30)
ds.globalattributes["time_step"] = ts
# round to the nearest time step
rldt = [rounddttots(dt,ts=ts) for dt in ldt]
elif mode.lower()=="nearest_second":
# round to the nearest second
rldt = [rounddttoseconds(dt) for dt in ldt]
else:
# unrecognised option for mode, return original datetime series
log.error(" round_datetime: unrecognised mode ("+str(mode)+")"+" ,returning original time series")
rldt = ds.series["DateTime"]["Data"]
# replace the original datetime series with the rounded one
ds.series["DateTime"]["Data"] = rldt
In [3]:
def rounddttots(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
def rounddttoseconds(dt):
dt += datetime.timedelta(seconds=0.5)
dt -= datetime.timedelta(seconds=dt.second % 1,microseconds=dt.microsecond)
return dt
def roundtobase(x,base=5):
return int(base*round(float(x)/base))
In [4]:
def read_eddypro_full(csvname):
ds = qcio.DataStructure()
csvfile = open(csvname,'rb')
csvreader = csv.reader(csvfile)
n = 0
adatetime = []
us_data_list = []
us_flag_list = []
Fh_data_list = []
Fh_flag_list = []
Fe_data_list = []
Fe_flag_list = []
Fc_data_list = []
Fc_flag_list = []
for row in csvreader:
if n==0:
header=row
elif n==1:
varlist=row
us_data_col = varlist.index('u*')
us_flag_col = varlist.index('qc_Tau')
Fh_data_col = varlist.index('H')
Fh_flag_col = varlist.index('qc_H')
Fe_data_col = varlist.index('LE')
Fe_flag_col = varlist.index('qc_LE')
Fc_data_col = varlist.index('co2_flux')
Fc_flag_col = varlist.index('qc_co2_flux')
elif n==2:
unitlist=row
else:
adatetime.append(datetime.datetime.strptime(row[1]+' '+row[2],'%Y-%m-%d %H:%M'))
us_data_list.append(float(row[us_data_col]))
us_flag_list.append(float(row[us_flag_col]))
Fh_data_list.append(float(row[Fh_data_col]))
Fh_flag_list.append(float(row[Fh_flag_col]))
Fe_data_list.append(float(row[Fe_data_col]))
Fe_flag_list.append(float(row[Fe_flag_col]))
Fc_data_list.append(float(row[Fc_data_col]))
Fc_flag_list.append(float(row[Fc_flag_col]))
n = n + 1
nRecs = len(adatetime)
ds.series['DateTime'] = {}
ds.series['DateTime']['Data'] = adatetime
round_datetime(ds,mode="nearest_timestep")
ds.series['ustar'] = {}
ds.series['ustar']['Data'] = numpy.array(us_data_list,dtype=numpy.float64)
ds.series['ustar']['Flag'] = numpy.array(us_flag_list,dtype=numpy.int32)
ds.series['Fh'] = {}
ds.series['Fh']['Data'] = numpy.array(Fh_data_list,dtype=numpy.float64)
ds.series['Fh']['Flag'] = numpy.array(Fh_flag_list,dtype=numpy.int32)
ds.series['Fe'] = {}
ds.series['Fe']['Data'] = numpy.array(Fe_data_list,dtype=numpy.float64)
ds.series['Fe']['Flag'] = numpy.array(Fe_flag_list,dtype=numpy.int32)
ds.series['Fc'] = {}
ds.series['Fc']['Data'] = numpy.array(Fc_data_list,dtype=numpy.float64)
ds.series['Fc']['Flag'] = numpy.array(Fc_flag_list,dtype=numpy.int32)
ds.globalattributes["nc_nrecs"] = nRecs
return ds
In [5]:
epname = qcio.get_filename_dialog(title='Choose an EddyPro full output file')
ofname = qcio.get_filename_dialog(title='Choose an L3 output file')
ds_ep = read_eddypro_full(epname)
ds_of = qcio.nc_read_series(ofname)
In [6]:
dt_ep = ds_ep.series['DateTime']['Data']
dt_of = ds_of.series['DateTime']['Data']
si = dt_of.index(dt_ep[0])
ei = dt_of.index(dt_ep[-1])
In [7]:
print dt_ep[0],dt_of[si]
print dt_ep[-1],dt_of[ei]
In [8]:
us_of,f,a = qcutils.GetSeriesasMA(ds_of,'ustar',si=si,ei=ei)
us_ep,f,a = qcutils.GetSeriesasMA(ds_ep,'ustar')
Fh_of,f,a = qcutils.GetSeriesasMA(ds_of,'Fh',si=si,ei=ei)
Fh_ep,f,a = qcutils.GetSeriesasMA(ds_ep,'Fh')
Fe_of,f,a = qcutils.GetSeriesasMA(ds_of,'Fe',si=si,ei=ei)
Fe_ep,f,a = qcutils.GetSeriesasMA(ds_ep,'Fe')
Fc_of,f,a = qcutils.GetSeriesasMA(ds_of,'Fc',si=si,ei=ei)
Fc_ep,f,a = qcutils.GetSeriesasMA(ds_ep,'Fc')
In [9]:
us_of.mask = numpy.ma.mask_or(us_of.mask,us_ep.mask)
us_ep.mask = numpy.ma.mask_or(us_of.mask,us_ep.mask)
Fh_of.mask = numpy.ma.mask_or(Fh_of.mask,Fh_ep.mask)
Fh_ep.mask = numpy.ma.mask_or(Fh_of.mask,Fh_ep.mask)
Fe_of.mask = numpy.ma.mask_or(Fe_of.mask,Fe_ep.mask)
Fe_ep.mask = numpy.ma.mask_or(Fe_of.mask,Fe_ep.mask)
Fc_of.mask = numpy.ma.mask_or(Fc_of.mask,Fc_ep.mask)
Fc_ep.mask = numpy.ma.mask_or(Fc_of.mask,Fc_ep.mask)
In [10]:
fig = plt.figure(1,figsize=(8,8))
qcplot.xyplot(us_ep,us_of,sub=[2,2,1],regr=1,xlabel='u*_EP (m/s)',ylabel='u*_OF (m/s)')
qcplot.xyplot(Fh_ep,Fh_of,sub=[2,2,2],regr=1,xlabel='Fh_EP (W/m2)',ylabel='Fh_OF (W/m2)')
qcplot.xyplot(Fe_ep,Fe_of,sub=[2,2,3],regr=1,xlabel='Fe_EP (W/m2)',ylabel='Fe_OF (W/m2)')
qcplot.xyplot(Fc_ep,Fc_of,sub=[2,2,4],regr=1,xlabel='Fc_EP (umol/m2/s)',ylabel='Fc_OF (umol/m2/s)')
plt.tight_layout()
plt.show()
In [19]:
fig=plt.figure()
plt.plot(dt_ep,Fh_ep,'b.')
plt.plot(dt_of[si:ei+1],Fh_of,'r+')
plt.show()
In [ ]: