In [1]:

    

%run basics
%matplotlib


    

    




Using matplotlib backend: Qt4Agg

In [42]:

    

site_name_list = ["AdelaideRiver","Calperum","CapeTribulation","CowBay","CumberlandPlains",
                  "DalyPasture","DalyUncleared","DryRiver","Emerald","FoggDam","Gingin",
                  "GreatWesternWoodlands","HowardSprings","Otway","RedDirtMelonFarm","RiggsCreek",
                  "RobsonCreek","Samford","SturtPlains","Tumbarumba","Whroo","WombatStateForest",
                  "Yanco"]


    

In [2]:

    

site_name_list = ["Calperum"]


    

In [46]:

    

results = numpy.array([])
for site_name in site_name_list:
    print "Processing site: ",site_name
    in_name = "/mnt/OzFlux/Sites/"+site_name+"/Data/Portal/"+site_name+"_L3.nc"
    if not os.path.exists(in_name):
        print " File not found: "+in_name
        continue
    ds = qcio.nc_read_series(in_name)
    ts = int(ds.globalattributes["time_step"])
    ntsInDay = int(float(24.0*60.0/float(ts)))
    nts_threshold = int(ntsInDay*0.8+0.5)
    ldt = ds.series["DateTime"]["Data"]
    for year in range(ldt[0].year,ldt[-1].year+1):
    #for year in [2015]:
        start_date = datetime.datetime(year,1,1)+datetime.timedelta(minutes=ts)
        si = qcutils.GetDateIndex(ldt,str(start_date),ts=ts,match="startnextday",default=0)
        end_date = datetime.datetime(year+1,1,1)
        ei = qcutils.GetDateIndex(ldt,str(end_date),ts=ts,match="endpreviousday",default=len(ldt)-1)
        dt = numpy.array(ldt[si:ei+1])
        nDays = int(float(len(dt))/ntsInDay)
        if nDays<250:
            print "Not enough days in year ",year,nDays
            continue
        # get the data
        Fa,f,a = qcutils.GetSeriesasMA(ds,"Fa",si=si,ei=ei)
        Fe,f,a = qcutils.GetSeriesasMA(ds,"Fe",si=si,ei=ei)
        Fh,f,a = qcutils.GetSeriesasMA(ds,"Fh",si=si,ei=ei)
        mask = numpy.ma.mask_or(numpy.ma.getmaskarray(Fa),numpy.ma.getmaskarray(Fe))
        mask = numpy.ma.mask_or(mask,numpy.ma.getmaskarray(Fh))
        Fa = numpy.ma.array(Fa,mask=mask)
        Fe = numpy.ma.array(Fe,mask=mask)
        Fh = numpy.ma.array(Fh,mask=mask)
        # reshape into 2D array of (days,timesteps)
        dt_daily = dt.reshape(nDays,ntsInDay)
        Fa_daily = Fa.reshape(nDays,ntsInDay)
        Fe_daily = Fe.reshape(nDays,ntsInDay)
        Fh_daily = Fh.reshape(nDays,ntsInDay)
        # get the daily averages and the number of non-masked points in each day
        Fa_daily_avg = numpy.ma.average(Fa_daily,axis=1)
        Fa_daily_count = numpy.ma.count(Fa_daily,axis=1)
        Fe_daily_avg = numpy.ma.average(Fe_daily,axis=1)
        Fe_daily_count = numpy.ma.count(Fe_daily,axis=1)
        Fh_daily_avg = numpy.ma.average(Fh_daily,axis=1)
        Fh_daily_count = numpy.ma.count(Fh_daily,axis=1)
        # only use days where more than 80% of the data is present
        condition = (Fa_daily_count<nts_threshold)|(Fe_daily_count<nts_threshold)|(Fh_daily_count<nts_threshold)
        Fa_daily_avg = numpy.ma.masked_where(condition,Fa_daily_avg)
        Fe_daily_avg = numpy.ma.masked_where(condition,Fe_daily_avg)
        Fh_daily_avg = numpy.ma.masked_where(condition,Fh_daily_avg)
        Fa_daily_avg = numpy.ma.compressed(Fa_daily_avg)
        Fe_daily_avg = numpy.ma.compressed(Fe_daily_avg)
        Fh_daily_avg = numpy.ma.compressed(Fh_daily_avg)
        FepFh = Fe_daily_avg + Fh_daily_avg
        if len(Fa_daily_avg)<2:
            print "Not enough good data for",site_name,year
            continue
        x = numpy.array(Fa_daily_avg,copy=True)
        y = numpy.array(FepFh,copy=True)
        x = x[:,numpy.newaxis]
        a, _, _, _ = numpy.linalg.lstsq(x, y)
        pr = scipy.stats.stats.pearsonr(Fa_daily_avg,FepFh)
        print site_name,year,a[0],pr[0]
        results = numpy.append(results,a[0])


    

    




Processing site:  AdelaideRiver
Not enough days in year  2007 75
AdelaideRiver 2008 0.998253211457 0.954920650843
Not enough days in year  2009 143
Processing site:  Calperum
Not enough days in year  2010 153
Calperum 2011 0.886211740219 0.991108563456
Calperum 2012 0.918026731725 0.992434829507
Calperum 2013 0.915439974595 0.989566040668
Calperum 2014 0.884973398339 0.990591701153
Calperum 2015 0.875843886491 0.988307898064
Processing site:  CapeTribulation
Not enough good data for CapeTribulation 2013
Not enough good data for CapeTribulation 2014
Not enough days in year  2015 245
Processing site:  CowBay
Not enough good data for CowBay 2009
Not enough good data for CowBay 2010
CowBay 2011 0.929385885001 0.835468434272
CowBay 2012 0.919783669279 0.894046934041
CowBay 2013 0.995664645354 0.690791998986
CowBay 2014 0.891102634422 0.747830889456
Not enough days in year  2015 131
Processing site:  CumberlandPlains
Not enough days in year  2012 74
CumberlandPlains 2013 0.780353628324 0.948494638038
CumberlandPlains 2014 0.710517869226 0.969683876873
CumberlandPlains 2015 0.77761108941 0.975118471942
Processing site:  DalyPasture
DalyPasture 2008 0.902095561745 0.90010938832
DalyPasture 2009 0.94110968103 0.858207360511
DalyPasture 2010 0.810983531565 0.856413467968
DalyPasture 2011 0.812939369902 0.876201408825
DalyPasture 2012 0.858889565031 0.958334499856
DalyPasture 2013 0.86306861537 0.931604265142
Processing site:  DalyUncleared
DalyUncleared 2008 0.997401567233 0.900781415033
DalyUncleared 2009 1.00588498972 0.525207272459
DalyUncleared 2010 0.890144624656 0.919303050796
DalyUncleared 2011 0.856232058521 0.881939322541
DalyUncleared 2012 0.896612498489 0.960097762782
DalyUncleared 2013 0.861783887373 0.926091180184
DalyUncleared 2014 0.902231157424 0.936631273813
Not enough days in year  2015 216
Processing site:  DryRiver
Not enough days in year  2008 122
DryRiver 2009 0.833852882194 0.910841242168
DryRiver 2010 0.795158289554 0.890069570051
DryRiver 2011 0.783756757169 0.830583297214
DryRiver 2012 0.844732816263 0.92294816612
DryRiver 2013 0.829041039577 0.92679909387
DryRiver 2014 0.84623075596 0.923552962269
Not enough days in year  2015 215
Processing site:  Emerald
Not enough days in year  2011 205
Emerald 2012 0.716992778073 0.920550302074
Emerald 2013 0.723053101928 0.928373130058
Processing site:  FoggDam
FoggDam 2006 0.849564740825 0.691493217885
FoggDam 2007 0.743934562331 0.953302930004
FoggDam 2008 0.739498073098 0.716591700319
Processing site:  Gingin
Not enough days in year  2011 79
Gingin 2012 0.94397241883 0.985554147935
Gingin 2013 0.913320322246 0.986703764546
Gingin 2014 0.883535623099 0.983134394007
Gingin 2015 0.866100101494 0.981421879281
Processing site:  GreatWesternWoodlands
GreatWesternWoodlands 2013 0.723690504729 0.949015114419
GreatWesternWoodlands 2014 0.760979837408 0.96531717337
GreatWesternWoodlands 2015 0.792463421326 0.947815312891
Processing site:  HowardSprings
Not enough good data for HowardSprings 2002
HowardSprings 2003 0.758149871246 0.547640393063
Not enough good data for HowardSprings 2004
Not enough good data for HowardSprings 2005
HowardSprings 2006 0.867947718057 0.814621307371
HowardSprings 2007 0.950830128075 0.90562135144
HowardSprings 2008 0.960591771252 0.915273192476
Not enough good data for HowardSprings 2009
HowardSprings 2010 1.03471045564 0.918420325163
HowardSprings 2011 0.971905936713 0.892802260535
HowardSprings 2012 0.978188838237 0.921272827953
HowardSprings 2013 0.968341253051 0.939342746816
HowardSprings 2014 0.918841236705 0.951596818537
Not enough days in year  2015 245
Processing site:  Otway
Not enough days in year  2007 143
Otway 2008 0.865354996181 0.951564741243
Otway 2009 0.879474612095 0.949745284186
Otway 2010 0.900089745627 0.949785863357
Processing site:  RedDirtMelonFarm
Not enough days in year  2011 99
RedDirtMelonFarm 2012 0.889472420004 0.925599545315
Not enough days in year  2013 201
Processing site:  RiggsCreek
RiggsCreek 2011 0.716048502807 0.970021259588
RiggsCreek 2012 0.79584938276 0.947969692242
RiggsCreek 2013 0.882123430086 0.984782182325
RiggsCreek 2014 0.839255233829 0.981088247368
Not enough days in year  2015 186
Processing site:  RobsonCreek
Not enough days in year  2013 152
RobsonCreek 2014 0.708986130754 0.692742234737
Not enough days in year  2015 131
Processing site:  Samford
Not enough days in year  2010 212
Samford 2011 0.737226428293 0.992060124175
Samford 2012 0.734524947993 0.982876810125
Samford 2013 0.722667196108 0.976972387378
Samford 2014 0.716362097953 0.978403486732
Not enough days in year  2015 180
Processing site:  SturtPlains
Not enough days in year  2008 125
SturtPlains 2009 1.03318285151 0.924838524349
SturtPlains 2010 1.02982405882 0.926863484704
SturtPlains 2011 1.03117738572 0.896488152657
SturtPlains 2012 0.961422550423 0.896931738793
SturtPlains 2013 0.966973307854 0.921122487111
SturtPlains 2014 0.932040740245 0.889814668448
Not enough days in year  2015 245
Processing site:  Tumbarumba
Not enough good data for Tumbarumba 2002
Tumbarumba 2003 1.25294429243 0.343656498236
Not enough good data for Tumbarumba 2004
Not enough good data for Tumbarumba 2005
Not enough good data for Tumbarumba 2006
Not enough good data for Tumbarumba 2007
Not enough good data for Tumbarumba 2008
Not enough good data for Tumbarumba 2009
Not enough good data for Tumbarumba 2010
Not enough good data for Tumbarumba 2011
Not enough good data for Tumbarumba 2012
Not enough good data for Tumbarumba 2013
Not enough good data for Tumbarumba 2014
Processing site:  Whroo
Not enough days in year  2011 30
Whroo 2012 0.897847004447 0.982304435084
Whroo 2013 0.903700648637 0.977837267698
Whroo 2014 0.908223026775 0.976654362005
Not enough days in year  2015 245
Processing site:  WombatStateForest
WombatStateForest 2010 0.715834777233 0.959078508411
WombatStateForest 2011 0.69060175623 0.967828756037
WombatStateForest 2012 0.77063347491 0.975703603846
WombatStateForest 2013 0.844281783408 0.978020409856
WombatStateForest 2014 0.866010391683 0.966289705077
Not enough days in year  2015 62
Processing site:  Yanco
Not enough days in year  2012 216
Yanco 2013 0.732402441515 0.942788168894
Yanco 2014 0.719534050378 0.94580178016
Not enough days in year  2015 173

In [47]:

    

print len(results)


    

    




83

In [55]:

    

bins = [0.70,0.75,0.80,0.85,0.90,0.95,1.00,1.05]
bin_labels = [">0.70",">0.75",">0.80",">0.85",">0.90",">0.95",">1.00",">1.05"]
#bins = [0.70,0.80,0.90,1.00,1.10]
#bin_labels = [">0.70",">0.80",">0.90",">1.00",">1.10"]
hist,edges = numpy.histogram(results,bins=bins)
hist = hist.astype(float)/float(numpy.sum(hist))
bar_width = 0.25
fig,axs = plt.subplots(1,figsize=(6,6))
index = numpy.arange(len(hist))
axs.bar(index,hist)
axs.set_xticks(index+bar_width)
axs.set_xticklabels(bin_labels)
axs.set_ylabel("Frequency",fontsize=16)
axs.set_xlabel("SEB ratio",fontsize=16)
plt.show()


    

In [ ]: