This new analysis script writes a summary of the observations that were...

This new analysis script writes a summary of the observations that were sampled to a html file. Note that it currently only works with an ObsPack observation object,
because of the file structure assumed on the input side.

da/analysis/summarize_obs.py

+#!/usr/bin/env python
+import sys
+sys.path.append('../../')
+import os
+import numpy as np
+import string
+import datetime as dt
+import logging
+import da.tools.io4 as io
+from da.tools.general import CreateDirs
+
+fontsize=10
+
+def nice_lat(cls):
+  #
+  # Convert latitude from decimal to cardinal
+  #
+  if cls > 0:
+     h = 'N'
+  else:
+     h = 'S'
+  
+  dec, deg = np.math.modf(cls)
+
+  return string.strip('%2d %2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
+
+def nice_lon(cls):
+  #
+  # Convert longitude from decimal to cardinal
+  #
+  if cls > 0:
+     h = 'E'
+  else:
+     h = 'W'
+  
+  dec, deg = np.math.modf(cls)
+
+  return string.strip('%3d %2d\'%s' % (abs(deg), round(abs(60 * dec), 0), h))
+
+def nice_alt(cls):
+  #
+  # Reformat elevation or altitude
+  #
+  return string.strip('%10.1f masl' % round(cls, -1))
+
+
+def SummarizeObs(DaCycle,printfmt='html'):
+    """***************************************************************************************
+    Call example:
+
+    python summarize_obs.py 
+
+    Option printfmt    : [tex,scr,html] print summary table in latex, terminal, or html format 
+
+    Other options are all those needed to create a DaCycle object
+
+    OR:
+
+    call directly from a python script as:
+
+    q=SummarizeObs(DaCycle,printfmt='html')
+
+    ***************************************************************************************"""
+
+    sumdir=os.path.join(DaCycle['dir.analysis'],'summary')
+    if not os.path.exists(sumdir):
+        logging.info( "Creating new directory "+sumdir )
+        os.makedirs(sumdir)
+
+    mrdir=os.path.join(DaCycle['dir.analysis'],'data_molefractions')
+    if not os.path.exists(mrdir):
+        logging.error( "Input directory does not exist (%s), exiting... "%mrdir )
+        return None
+
+    mrfiles = os.listdir(mrdir)
+    infiles = [os.path.join(mrdir,f) for f in mrfiles if f.endswith('.nc')]
+
+    if printfmt == 'tex': 
+        print '\\begin{tabular*}{\\textheight}{l l l l r r r r}'
+        print 'Code &  Name & Lat, Lon, Elev & Lab &  N (flagged) & $\\sqrt{R}$  &Inn \\XS &Bias\\\\'
+        print '\hline\\\\ \n\multicolumn{8}{ c }{Semi-Continuous Surface Samples}\\\\[3pt] '
+        fmt= '%8s  & '+' %55s  & '+'%20s &'+'%6s &'+' %4d (%d)  & '+' %5.2f  & '+' %5.2f & '+'%+5.2f  \\\\'
+    elif printfmt == 'html':
+        tablehead = \
+              "<TR>\n <TH> Site code </TH> \
+                   <TH> Sampling Type </TH> \
+                   <TH> Lab. </TH> \
+                   <TH> Country </TH> \
+                   <TH> Lat, Lon, Elev. (m ASL) </TH> \
+                   <TH> No. Obs. Avail. </TH>  \
+                   <TH> &#8730;R (&mu;mol mol<sup>-1</sup>) </TH>  \
+                   <TH> &#8730;HPH  (&mu;mol mol<sup>-1</sup>) </TH> \
+                   <TH> H(x)-y (&mu;mol mol<sup>-1</sup>) </TH> \
+                   <TH> H(x)-y (JJAS) (&mu;mol mol<sup>-1</sup>) </TH>  \
+                   <TH> H(x)-y (NDJFMA) (&mu;mol mol<sup>-1</sup>) </TH> \n \
+               </TR>\n"
+
+        fmt=  """<TR> \n \
+                <TD><a href='javascript:LoadCO2Tseries("%s")'>%s </a></TD>\
+                <TD>%s</TD>\
+                <TD>%s</TD>\
+                <TD>%40s</TD>\
+                <TD>%s</TD>\
+                <TD>%d</TD>\
+                <TD>%+5.2f</TD>\
+                <TD>%+5.2f</TD>\
+                <TD>%+5.2f&plusmn;%5.2f</TD>\
+                <TD>%+5.2f&plusmn;%5.2f</TD>\
+                <TD>%+5.2f&plusmn;%5.2f</TD>\n \
+               </TR>\n"""
+    elif printfmt == 'scr':
+        print 'Code   Site     NObs flagged  R  Inn X2'
+        fmt= '%8s '+' %55s  %s %s'+' %4d '+' %4d '+' %5.2f '+' %5.2f'
+
+    table=[]
+    lons=[]
+    lats=[]
+    names=[]
+    nobs=[]
+    for infile in infiles:
+
+            #logging.debug( infile )
+            f         = io.CT_CDF(infile,'read')
+            date      = f.GetVariable('time')
+            obs       = f.GetVariable('value')*1e6
+            mdm       = f.GetVariable('modeldatamismatch')*1e6
+            simulated = f.GetVariable('modelsamplesmean')*1e6
+            simulated_std = f.GetVariable('modelsamplesstandarddeviation')*1e6
+
+            pydates = [dt.datetime(1970,1,1)+dt.timedelta(seconds=int(d)) for d in date]
+
+            summer = [i for i,d in enumerate(pydates) if d.month in [6,7,8,9] ]
+            winter = [i for i,d in enumerate(pydates) if d.month in [11,12,1,2,3,4] ]
+
+            diff=((simulated-obs).mean())
+            diffsummer=((simulated-obs).take(summer).mean())
+            diffwinter=((simulated-obs).take(winter).mean())
+            diffstd=((simulated-obs).std())
+            diffsummerstd=((simulated-obs).take(summer).std())
+            diffwinterstd=((simulated-obs).take(winter).std())
+            longsitestring=f.site_name+', '+f.site_country
+            location=nice_lat(f.site_latitude)+', '+ nice_lon(f.site_longitude)+', '+nice_alt(f.site_elevation)
+
+            if printfmt == 'html':
+                ss=(f.site_code.upper(),
+                    f.site_code.upper(),
+                    f.dataset_project,
+                    f.lab_abbr,
+                    f.site_country,
+                    location,
+                    len(mdm.compressed()),
+                    mdm.mean(),
+                    np.sqrt((simulated_std**2).mean()),
+                    diff,diffstd,
+                    diffsummer,diffsummerstd,
+                    diffwinter,diffwinterstd)
+            
+            table.append(ss)
+            lons.append(f.site_longitude)
+            lats.append(f.site_latitude)
+            names.append(f.site_code)
+            nobs.append(len(mdm.compressed()))
+            f.close()
+
+    if printfmt == 'tex':
+        saveas=os.path.join(sumdir,'site_table.tex')
+        f=open(saveas,'w')
+    elif printfmt == 'html':
+        saveas=os.path.join(sumdir,'site_table.html')
+        f=open(saveas,'w')
+        txt = "<meta http-equiv='content-type' content='text/html;charset=utf-8' />\n"
+        f.write(txt)
+        txt="<table border=1 cellpadding=2 cellspacing=2 width='100%' bgcolor='#EEEEEE'>\n"
+        f.write(txt)
+
+    f.write(tablehead)
+
+    for i,ss in enumerate(table):
+
+        f.write(fmt%ss)
+        if (i+1)%15 == 0:
+            f.write(tablehead)
+
+    if printfmt == 'tex': 
+        f.write( '\cline{2-8}\\\\' )
+        f.write( '\hline \\\\')
+        f.write( '\end{tabular*}')
+    else:
+        txt="\n</table>"
+        f.write(txt)
+    f.close()
+
+    logging.info("File written with summary: %s" % saveas)
+
+# main body if called as script
+
+if __name__ == '__main__':    # started as script
+
+    from da.tools.initexit import CycleControl
+    from da.ct.dasystem import CtDaSystem 
+    from da.ct.statevector import CtStateVector 
+
+    sys.path.append('../../')
+
+    logging.root.setLevel(logging.DEBUG)
+
+    DaCycle = CycleControl(args={'rc':'../../ctdas-od-gfed2-glb6x4-obspack.rc'})
+    DaCycle.Initialize()
+    DaCycle.ParseTimes()
+
+    DaSystem    = CtDaSystem('../rc/carbontracker_ct09_op.rc')
+    DaSystem.Initialize()
+
+    DaCycle.DaSystem    = DaSystem
+
+    q=SummarizeObs(DaCycle)
+
+    sys.exit(0)
+
+