In [1]:
%run basics
%matplotlib
import pytz
In [2]:
cf = qcio.load_controlfile(path="controlfiles")
site_list = cf["Sites"]
in_path = cf["Sites"]["AdelaideRiver"]["in_filepath"]
in_name = cf["Sites"]["AdelaideRiver"]["in_filename"]
in_filename = os.path.join(in_path,in_name)
out_path = cf["Sites"]["AdelaideRiver"]["out_filepath"]
out_name = cf["Sites"]["AdelaideRiver"]["out_filename"]
out_filename = os.path.join(out_path,out_name)
interpolate = qcutils.get_keyvaluefromcf(cf,["Sites","AdelaideRiver"],"interpolate",default=True)
site_name = cf["Sites"]["AdelaideRiver"]["site_name"]
site_timezone = cf["Sites"]["AdelaideRiver"]["site_timezone"]
site_tz = pytz.timezone(site_timezone)
ds_60minutes = qcio.DataStructure()
ds_30minutes = qcio.DataStructure()
In [3]:
f = netCDF4.MFDataset(in_filename)
In [4]:
# get the local datetime
valid_date = f.variables["valid_date"][:]
valid_time = f.variables["valid_time"][:]
dt_utc_all=numpy.array([datetime.datetime.strptime(str(int(valid_date[i])*10000+int(valid_time[i])),"%Y%m%d%H%M") for i in range(0,len(valid_date))])
time_step = numpy.array([(dt_utc_all[i]-dt_utc_all[i-1]).total_seconds() for i in range(1,len(dt_utc_all))])
time_step = numpy.append(time_step,3600)
index = numpy.where(time_step!=0)[0]
dt_utc = dt_utc_all[index]
dt_utc = [x.replace(tzinfo=pytz.utc) for x in dt_utc]
dt_loc = [x.astimezone(site_tz) for x in dt_utc]
dt_loc = [x-x.dst() for x in dt_loc]
dt_loc = [x.replace(tzinfo=None) for x in dt_loc]
ds_60minutes.series["DateTime"] = {}
ds_60minutes.series["DateTime"]["Data"] = dt_loc
In [5]:
# set some global attributes
nRecs = len(dt_loc)
ds_60minutes.globalattributes["nc_nrecs"] = nRecs
ds_60minutes.globalattributes["time_step"] = 60
ds_60minutes.globalattributes["time_zone"] = site_timezone
ds_60minutes.globalattributes["site_name"] = site_name
ds_60minutes.globalattributes["xl_datemode"] = 0
# get the datetime in some different formats
qcutils.get_xldatefromdatetime(ds_60minutes)
qcutils.get_ymdhmsfromdatetime(ds_60minutes)
In [6]:
var_list= cf["Variables"].keys()
flag_60minutes = numpy.zeros(nRecs,dtype=numpy.int32)
for label in var_list:
access_name = cf["Variables"][label]["access_name"]
if access_name not in f.variables.keys():
log.error("Requested variable "+access_name+" not found in ACCESS data")
continue
attr = {}
for vattr in f.variables[access_name].ncattrs():
attr[vattr] = getattr(f.variables[access_name],vattr)
attr["missing_value"] = c.missing_value
for i in range(0,3):
for j in range(0,3):
label_ij = label+'_'+str(i)+str(j)
if len(f.variables[access_name].shape)==3:
series = f.variables[access_name][:,i,j]
elif len(f.variables[access_name].shape)==4:
series = f.variables[access_name][:,0,i,j]
else:
print "Unrecognised variable ("+label+") dimension in ACCESS file"
series = series[index]
qcutils.CreateSeries(ds_60minutes,label_ij,series,
Flag=flag_60minutes,Attr=attr)
In [7]:
# get derived quantities and adjust units
# air temperature from K to C
attr = qcutils.GetAttributeDictionary(ds_60minutes,"Ta_00")
if attr["units"] == "K":
for i in range(0,3):
for j in range(0,3):
label = "Ta_"+str(i)+str(j)
Ta,f,a = qcutils.GetSeriesasMA(ds_60minutes,label)
Ta = Ta - c.C2K
attr["units"] = "C"
qcutils.CreateSeries(ds_60minutes,label,Ta,Flag=flag_60minutes,Attr=attr)
# soil temperature from K to C
attr = qcutils.GetAttributeDictionary(ds_60minutes,"Ts_00")
if attr["units"] == "K":
for i in range(0,3):
for j in range(0,3):
label = "Ts_"+str(i)+str(j)
Ts,f,a = qcutils.GetSeriesasMA(ds_60minutes,label)
Ts = Ts - c.C2K
attr["units"] = "C"
qcutils.CreateSeries(ds_60minutes,label,Ts,Flag=flag_60minutes,Attr=attr)
# pressure from Pa to kPa
attr = qcutils.GetAttributeDictionary(ds_60minutes,"ps_00")
if attr["units"] == "Pa":
for i in range(0,3):
for j in range(0,3):
label = "ps_"+str(i)+str(j)
ps,f,a = qcutils.GetSeriesasMA(ds_60minutes,label)
ps = ps/float(1000)
attr["units"] = "kPa"
qcutils.CreateSeries(ds_60minutes,label,ps,Flag=flag_60minutes,Attr=attr)
# wind speed from components
for i in range(0,3):
for j in range(0,3):
u_label = "u_"+str(i)+str(j)
v_label = "v_"+str(i)+str(j)
Ws_label = "Ws_"+str(i)+str(j)
u,f,a = qcutils.GetSeriesasMA(ds_60minutes,u_label)
v,f,a = qcutils.GetSeriesasMA(ds_60minutes,v_label)
Ws = numpy.sqrt(u*u+v*v)
attr = qcutils.MakeAttributeDictionary(long_name="Wind speed",units="m/s",height="10m")
qcutils.CreateSeries(ds_60minutes,Ws_label,Ws,Flag=f,Attr=attr)
# wind direction from components
for i in range(0,3):
for j in range(0,3):
u_label = "u_"+str(i)+str(j)
v_label = "v_"+str(i)+str(j)
Wd_label = "Wd_"+str(i)+str(j)
u,f,a = qcutils.GetSeriesasMA(ds_60minutes,u_label)
v,f,a = qcutils.GetSeriesasMA(ds_60minutes,v_label)
Wd = float(270) - numpy.ma.arctan2(v,u)*float(180)/numpy.pi
index = numpy.ma.where(Wd>360)[0]
if len(index)>0: Wd[index] = Wd[index] - float(360)
attr = qcutils.MakeAttributeDictionary(long_name="Wind direction",units="degrees",height="10m")
qcutils.CreateSeries(ds_60minutes,Wd_label,Wd,Flag=f,Attr=attr)
# relative humidity from temperature, specific humidity and pressure
for i in range(0,3):
for j in range(0,3):
q_label = "q_"+str(i)+str(j)
Ta_label = "Ta_"+str(i)+str(j)
ps_label = "ps_"+str(i)+str(j)
RH_label = "RH_"+str(i)+str(j)
q,f,a = qcutils.GetSeriesasMA(ds_60minutes,q_label)
Ta,f,a = qcutils.GetSeriesasMA(ds_60minutes,Ta_label)
ps,f,a = qcutils.GetSeriesasMA(ds_60minutes,ps_label)
RH = mf.RHfromspecifichumidity(q, Ta, ps)
attr = qcutils.MakeAttributeDictionary(long_name='Relative humidity',units='%',standard_name='not defined')
qcutils.CreateSeries(ds_60minutes,RH_label,RH,Flag=f,Attr=attr)
# absolute humidity from temperature and relative humidity
for i in range(0,3):
for j in range(0,3):
Ta_label = "Ta_"+str(i)+str(j)
RH_label = "RH_"+str(i)+str(j)
Ah_label = "Ah_"+str(i)+str(j)
Ta,f,a = qcutils.GetSeriesasMA(ds_60minutes,Ta_label)
RH,f,a = qcutils.GetSeriesasMA(ds_60minutes,RH_label)
Ah = mf.absolutehumidityfromRH(Ta, RH)
attr = qcutils.MakeAttributeDictionary(long_name='Absolute humidity',units='g/m3',standard_name='not defined')
qcutils.CreateSeries(ds_60minutes,Ah_label,Ah,Flag=f,Attr=attr)
# soil moisture from kg/m2 to m3/m3
attr = qcutils.GetAttributeDictionary(ds_60minutes,"Sws_00")
for i in range(0,3):
for j in range(0,3):
label = "Sws_"+str(i)+str(j)
Sws,f,a = qcutils.GetSeriesasMA(ds_60minutes,label)
Sws = Sws/float(100)
attr["units"] = "frac"
qcutils.CreateSeries(ds_60minutes,label,Sws,Flag=flag_60minutes,Attr=attr)
# net radiation and upwelling short and long wave radiation
for i in range(0,3):
for j in range(0,3):
label_Fn = "Fn_"+str(i)+str(j)
label_Fsd = "Fsd_"+str(i)+str(j)
label_Fld = "Fld_"+str(i)+str(j)
label_Fsu = "Fsu_"+str(i)+str(j)
label_Flu = "Flu_"+str(i)+str(j)
label_Fn_sw = "Fn_sw_"+str(i)+str(j)
label_Fn_lw = "Fn_lw_"+str(i)+str(j)
Fsd,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fsd)
Fld,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fld)
Fn_sw,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fn_sw)
Fn_lw,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fn_lw)
Fsu = Fsd - Fn_sw
Flu = Fld - Fn_lw
Fn = (Fsd-Fsu)+(Fld-Flu)
attr = qcutils.MakeAttributeDictionary(long_name='Up-welling long wave',
standard_name='surface_upwelling_longwave_flux_in_air',
units='W/m2')
qcutils.CreateSeries(ds_60minutes,label_Flu,Flu,Flag=f,Attr=attr)
attr = qcutils.MakeAttributeDictionary(long_name='Up-welling short wave',
standard_name='surface_upwelling_shortwave_flux_in_air',
units='W/m2')
qcutils.CreateSeries(ds_60minutes,label_Fsu,Fsu,Flag=f,Attr=attr)
attr = qcutils.MakeAttributeDictionary(long_name='Calculated net radiation',
standard_name='surface_net_allwave_radiation',
units='W/m2')
qcutils.CreateSeries(ds_60minutes,label_Fn,Fn,Flag=f,Attr=attr)
# ground heat flux as residual
for i in range(0,3):
for j in range(0,3):
label_Fg = "Fg_"+str(i)+str(j)
label_Fn = "Fn_"+str(i)+str(j)
label_Fh = "Fh_"+str(i)+str(j)
label_Fe = "Fe_"+str(i)+str(j)
Fn,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fn)
Fh,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fh)
Fe,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fe)
Fg = Fn - Fh - Fe
attr = qcutils.MakeAttributeDictionary(long_name='Calculated ground heat flux',
standard_name='downward_heat_flux_in_soil',
units='W/m2')
qcutils.CreateSeries(ds_60minutes,label_Fg,Fg,Flag=f,Attr=attr)
# Available energy
for i in range(0,3):
for j in range(0,3):
label_Fg = "Fg_"+str(i)+str(j)
label_Fn = "Fn_"+str(i)+str(j)
label_Fa = "Fa_"+str(i)+str(j)
Fn,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fn)
Fg,f,a = qcutils.GetSeriesasMA(ds_60minutes,label_Fg)
Fa = Fn - Fg
attr = qcutils.MakeAttributeDictionary(long_name='Calculated available energy',
standard_name='not defined',units='W/m2')
qcutils.CreateSeries(ds_60minutes,label_Fa,Fa,Flag=f,Attr=attr)
In [8]:
ncfile = qcio.nc_open_write(out_filename)
qcio.nc_write_series(ncfile, ds_60minutes,ndims=1)
In [12]:
Fsd,f,a=qcutils.GetSeriesasMA(ds_60minutes,"Fsd_11")
fig = plt.figure()
plt.plot(dt_loc,Fsd)
plt.show()
In [9]:
Ah_00,f,a = qcutils.GetSeriesasMA(ds_60minutes,"Ah_00")
In [13]:
print valid_date[0],valid_time[0],dt_loc[0],Ah_00[0]
In [ ]: