Changeset 57 for CedaCdoScripts

Timestamp:

23/09/13 17:14:25 (7 years ago)

Author:

mjuckes

Message:

some updates to cdo_basic.py, and adding other bits of code

Location:

CedaCdoScripts/Misc

Files:

• cdo_basic.py (modified) (35 diffs)
• derived_bundle.py (added)
• run_cdo_mf.py (added)
• var_mf.txt (added)

CedaCdoScripts/Misc/cdo_basic.py

@@ -8 +8 @@
 debug = True
 cleanUpAsUGo = True
-args = getopt.getopt( sys.argv[1:], 'i:m:e:r:b:x:o:', ['onlyVars=','scratch=','omitVars='] )
+args = getopt.getopt( sys.argv[1:], 'i:m:e:r:b:x:o:', ['onlyVars=','scratch=','omitVars=','version=','log='] )
 ad = { 'i':'inst', 'm':'model', 'e':'ensem', 'r':'realm', 'b':'batch', 'x':'expt', 'o':'outDir' }
 
-knownops = ['am','bm','mclim5','cdd','tr']
+knownops = ['am','sm','bm','mclim5','cdd','tr', 'fd']
 optag = args[-1][0]
 assert optag in knownops, 'First argument [%s] should be operation in: %s' % (optag, str(knmownops))
@@ -22 +22 @@
 fullYears = True
 
-print args
 edict = {}
 for a in args[0]:
   edict[a[0][1:]] = a[1] 
+
+def writeat( fh, tag, str,log):
+  if len(str) + len(tag) + 2 > 1024 :
+    log.warn( 'TRUNCATING global attribute %s: %s' % (tag,str) )
+    strl = 1022 - len(tag)
+    str = str[:strl]
+  fh.write( '%s %s\n' % (tag,str) )
 
 fullYearNeeded = False
@@ -50 +56 @@
 ##batch = 'b011_'
 ##onlyVars = ['tas']
+project = 'cmip5'
 ensem = edict['e']
 inst = edict['i']
@@ -56 +63 @@
 experiment = edict['x']
 batch = edict['b'] + '_'
+processing_batch = edict['b']
 opdir = edict['o']
 onlyVars = None
@@ -64 +72 @@
 ## qcDir = '%s%s/%s/%s/%s/' % (qcBaseDir,product,inst, model, expt)
 
-tstring = '%s%2.2i%2.2i-%2.2i%2.2i%2.2i' % time.gmtime()[0:6]
 vstring = '%s%2.2i%2.2i' % time.gmtime()[0:3]
 doctag = optag
 if optag in ['cdd']:
   doctag = 'precip_indices'
-if optag in ['tr']:
+if optag in ['tr','fd']:
   doctag = 'temperature_indices'
 doc_filename = 'ceda_cmip5_processing_%s.pdf' % doctag
-log_filename = 'ceda_cmip5_processing_log_%s_%s.txt' % (optag,tstring)
+tstring = '%s%2.2i%2.2i-%2.2i%2.2i%2.2i' % time.gmtime()[0:6]
+log_filename = edict.get( '-log', 'ceda_cmip5_processing_log_%s_%s.txt' % (optag,tstring) )
 batching='byFile'
 ofreq = 'yr'
@@ -79 +87 @@
   batching = '5yr'
   monthsPerBatch = 60
-elif optag in ['cdd','tr']:
+elif optag in ['cdd','tr','fd']:
   batching = 'yr'
   ifreq = 'day'
@@ -90 +98 @@
 
 scriptName = 'cdo_basic.py'
-scriptVersion = 0.8
+scriptVersion = 0.9
 log = logging.getLogger( __name__ )
 fh = logging.FileHandler( log_filename )
@@ -108 +116 @@
    omitVars = string.split( edict['-omitVars'], ',' )
    log.info( 'Omitting variables: %s' % edict['-omitVars'] )
+if edict.has_key( '-version' ):
+   version = edict['-version']
+   log.info( 'Processing version: %s' % edict['-version'] )
+else:
+   version = None
 ## log.info( 'CDO version information:' )
 ## os.popen( 'cdo -V > .cdo_version' ).readlines()
@@ -123 +136 @@
 if not os.path.isdir(tmpdir):
   os.mkdir( tmpdir )
-  print 'Creating temporary directory: %s ' % tmpdir
+  log.info( 'Creating temporary directory: %s ' % tmpdir )
 
 if tmpdir[-1] != '/':
@@ -133 +146 @@
 assert os.path.isfile( '%scdo_basic_lock' % tmpdir ), 'Failed to make lock file in work directory %s' % tmpdir
 
-
 plurality = 'single'
 mdbase = "http://badc.nerc.ac.uk/browse/badc/cmip5/metadata/processing/" 
@@ -139 +151 @@
 rcout = open( 'cdo_stdout.txt', 'w' )
 rcerr = open( 'cdo_stderr.txt', 'w' )
-def run_cmd( cmd, comment='', strict=True ):
+def run_cmd( cmd, comment='', strict=True, append=None, out=None  ):
     if type(cmd) == type('x'):
       args = map( lambda x: string.strip(x, "'"), string.split( cmd ))
@@ -145 +157 @@
       args = cmd
     log.info( cmd )
-    proc = subprocess.Popen( args, stdout=rcout, stderr=rcerr, env={ 'LD_LIBRARY_PATH':'/usr/local/lib'} )
+    if append != None:
+      fho = open( append, 'a' )
+      clo = True
+    elif out != None:
+      fho = open( out, 'w' )
+      clo = True
+    else:
+      fho = rcout
+      clo = False
+    proc = subprocess.Popen( args, stdout=fho, stderr=rcerr, env={ 'LD_LIBRARY_PATH':'/usr/local/lib'} )
     std = proc.communicate()
     rc = proc.returncode
@@ -156 +177 @@
       elif rc != 0:
         print 'Failed (%s) [rc=%s] executing %s' % (comment,rc,str(cmd))
+    if clo:
+      fho.close()
     return rc
 
 dir_tmpl_s = '/disks/fuse/cmip5/output1/MOHC/HadGEM2-ES/%(expt)s/%(ifreq)s/atmos/Amon/'
+dir_tmpl_s = '/badc/cmip5/data/cmip5/output1/MOHC/HadGEM2-ES/%(expt)s/%(ifreq)s/atmos/Amon/'
 if ifreq == 'mon':
   mip = {'atmos':'Amon', 'ocean':'Omon', 'land':'Lmon' }[realm]
@@ -166 +190 @@
   raise 'Should not be here'
 
-dir_tmpl_r = '/disks/fuse/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/'
-dir_tmpl_e = '/disks/fuse/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/latest/'
+dir_tmpl_r = '/badc/cmip5/data/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/'
+dir_tmpl_e = '/badc/cmip5/data/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/latest/'
 dir_tmpl = dir_tmpl_e + '%(var)s/'
 bldir_tmpl = '/home/users/prototype/cmip_processing/cmip5_bl/output1_%(inst)s_%(model)s_%(expt)s/'
@@ -180 +204 @@
   if len(b2) == 0:
     return None
-  return b2[-1]
-
-def get_version_dir():
+  var = string.split( b2[-1], '_' )[0]
+  return var
+
+def get_version_dir(version=None):
 
   pdict = { 'inst':inst, 'model':model, 'expt':expt, 'ensem':ensem, 'ifreq':ifreq, 'realm':realm, 'mip':mip }
   rdir = dir_tmpl_r % pdict
+  if version != None:
+    this_v = rdir + version
+    assert os.path.isdir( this_v ), 'Version directory %s not found' % this_v
+    return this_v
+
   vlist = glob.glob( rdir + 'v*' )
   assert len(vlist) > 0, 'No version directories found in %s' % rdir
@@ -200 +230 @@
   ## edir = dir_tmpl % { 'inst':inst, 'model':model, 'expt':expt, 'var':var, 'ensem':ensem, 'ifreq':ifreq, 'realm':realm, 'mip':mip }
   assert os.path.isdir( edir), 'Error in get_expt_fl: %s is not a directory' % edir
+  log.info( 'edir:: %s' % edir )
   fl = glob.glob( edir + '*.nc' )
   assert len(fl) > 0, 'No netcdf files found in %s' % edir
@@ -211 +242 @@
     if os.path.isfile( bldir + blfile ):
       blvars = map( get_var_from_bl_record, open( bldir + blfile ).readlines() )
-  if var in blvars:
-    print 'skipping ',var
-    log.info( 'Skipping black-listed ADS expt=%s, var=%s' % (expt,var) )
-    return fm
+      if var in blvars:
+        log.info( 'Skipping black-listed ADS expt=%s, var=%s' % (expt,var) )
+        return fm
+      else:
+        log.info( 'Var %s not in black list [%s] %s' % (var, bldir + blfile, str(blvars)) )
+    else:
+      log.info( 'Black list file %s not present' % (bldir+blfile) )
+  else:
+    log.info( 'Black list directory %s not present' % bldir )
   for f in fl:
     fn = string.split(f, '/')[-1]
@@ -268 +304 @@
    onlyVars = ['pr']
    omitVars = None
-if optag in ['tr']:
+if optag in ['tr','fd']:
    onlyVars = ['tasmin']
    omitVars = None
@@ -274 +310 @@
 for expt in expts:
   ##dir = dir_tmpl_e % locals()
-  dir = get_version_dir()
+  dir = get_version_dir( version=version)
   vars = glob.glob( '%s/*' % dir )
-  print expt, dir, vars
   thisvarl = map( lambda x: string.split(x,'/')[-1], vars )
   thisvarl.sort()
@@ -312 +347 @@
 lon1,lon2,lat1,lat2 = (60.,120.,30.,60.)
 cdoop, ofreq,opName, proc_description = { 'am':('divdpy -yearsum -muldpm', 'yr','Annual mean','Annual mean of monthly mean data'), \
+                        'sm':('divdpy -seassum -muldpm', 'yr','Seasonal mean','Annual seasonal of monthly mean data'), \
                         'bm':('fldmean -sellonlatbox,%s,%s,%s,%s' % (lon1,lon2,lat1,lat2),ifreq,'Box mean','Spatial mean over a latitude-longitude box'), \
                         'mclim5':('ymonavg', '5yr','5-year monthly climatology','Monthly climatology for 5 year period: 5 year mean for each month'),
                         'cdd':('eca_cdd','yr','Continuous dry days','Maximum of continuous dry days in a year'),
-                        'tr':('eca_tr','yr','Tropical nights','Count of tropical nights per year')  }[optag]
+                        'tr':('eca_tr','yr','Tropical nights','Count of tropical nights per year'),
+                        'fd':('eca_fd','yr','Frost days','Count of frost days (Tmin below zero) per year') }[optag]
 ## opdir = '%s%s' % (batch,optag)
 if go and not os.path.isdir( opdir ):
@@ -327 +364 @@
 k=-1
 print 'INFO0080: ',k, str(k), str( [1,2,3] )
-for expt,var,fm in ye:
+##for expt,var,fm in ye:
+def run_evf( expt,var,fm ):
     kk = 0
-    dir = get_version_dir()
+    dir = get_version_dir( version=version)
     ##edir = dir_tmpl % locals()
     edir = '%s/%s/' % (dir,var)
@@ -367 +405 @@
         ifiles[key] += '%s,' % file[0]
 
+    startstart = None
     nbfb = 0
     nbfbb = 0
@@ -395 +434 @@
       nodir = False
       if nbfb > 0:
-        print 'FILE DATE BOUNDARIES DO NOT MATCH ASSUMPTIONS: %s, %s' % (expt,var)
-        print 'merging files'
-        print nbfb, nbfbb
         log.info( 'File date boundaries do not match end of year (%s,%s): merging: %s, %s' % (nbfb, nbfbb, expt,var) )
         if nbfbb > 0:
@@ -441 +477 @@
             cmd1 = None
             cmd = '/usr/local/2/bin/cdo -s %s %s/%s %s' % (cdoop, tmpdir, '.tmp1.nc', ofile)
-          print cmd
+          log.info( 'cdo command: %s\n' % cmd )
           ## raise 'bad file boundary'
            
@@ -471 +507 @@
           nf = len(fm)
           if smem[0][0] != 1:
-            print 'Omitting first file: %s' % str(fm[0])
+            log.info(  'Omitting first file: %s' % str(fm[0]) )
             i0 = 1
             nf = len(fm)-1
@@ -514 +550 @@
           ofiles[key] += ' %s' % (ofile)
       if nbfb > 0:
-          print 'FILE DATE BOUNDARIES DO NOT MATCH ASSUMPTIONS: %s, %s' % (expt,var)
+          log.info( 'FILE DATE BOUNDARIES DO NOT MATCH ASSUMPTIONS: %s, %s' % (expt,var) )
     elif batching in ['5yr','yr']:
       startstart = fm[0][1]
@@ -566 +602 @@
 
         for k in range(len(tfm)-1):
-            print 'INFO0100: ',k, str(k), str( [1,2,3] )
-            print 'INFO0101: ',k, tfm[k], tfm[k+1]
+            log.info( 'INFO0100: ',k, str(k), str( [1,2,3] ) )
+            log.info( 'INFO0101: ',k, tfm[k], tfm[k+1] )
             a = str(tfm[k])
             b = str(tfm[k+1])
@@ -587 +623 @@
             if debug:
               log.info( 'Trying to log file maximum value' )
-              cmd0 = 'echo %s%s >> /home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' % (edir,tfm[0][0])
-              cmd = '/usr/local/2/bin/cdo info %s 2>&1 >> /home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' % ofile
-              rc = run_cmd( cmdo )
-              rc = run_cmd( cmd )
+              cmdo = 'echo %s%s' % (edir,tfm[0][0])
+              cmd = '/usr/local/2/bin/cdo info %s' % ofile
+              rc = run_cmd( cmdo, append='/home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' )
+              rc = run_cmd( cmd, append='/home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' )
           else:
             print 'single: %s' % cmd
@@ -635 +671 @@
         else:
           raise 'Code missing'
-
-attfile = '%s/temp_attribute_file.txt' % tmpdir
-log.info( 'Merging files' )
-if go:
-  for expt in expts:
+    return tid, crd, table_id, realization, initialization_method, physics_version, startstart
+
+
+for expt,var,fm in ye:
+    tid, crd, table_id, realization, initialization_method, physics_version, startstart = run_evf( expt,var,fm )
+
+def proc_expt(expt,version,optag):
+    attfile = '%s/temp_attribute_file.txt' % tmpdir
     varl = varld[expt]
     ## for var in varld[expt]:
     for var in varl:
+      if optag in ['am', 'sm']:
+        if var == 'tasmax':
+          cmval = "time: maximum within days time: mean over days"
+        elif var == 'tasmin':
+          cmval = "time: minimum within days time: mean over days"
+        else:
+          cmval = "time: mean"
+      else:
+        assert True, 'Need to set cell_methods for this operator'
+
       key = '%s_%s' % (expt,var)
       if key in badkeys:
-        print 'Skipping %s because or earlier error' % key
+        log.info( 'Skipping %s because or earlier error' % key )
       elif string.strip( ofiles[key] ) == '':
-        print 'Skipping %s because no output files registered' % key
+        log.info( 'Skipping %s because no output files registered' % key  )
         log.warning( 'WARN004: unexplained missing output file for key %s' % (key) )
       else:
-        print 'Trying %s' % key
-        dir = get_version_dir()
+        log.info( 'Trying %s' % key )
+        dir = get_version_dir( version=version)
     ##    vdir = dir_tmpl % locals()
         vdir = '%s/%s/' % (dir,var)
@@ -667 +716 @@
         for k in range(36):
           tid += cx[random.randint(0,len(cx)-1 )]
+        this_mip  = string.split( table_id )[1]
+        version = string.split( string.strip( dir, '/' ), '/' )[-1]
+        esgf_dsid = '%s.%s.%s.%s.%s.%s.%s.%s.%s.%s' % (project,product,inst,expt,model,ifreq,realm,this_mip,ensem,version)
+##cmip5.output1.MIROC.MIROC5.sstClim.monClim.atmos.Amon.r1i1p1.v20111129.xml
         oo = open( attfile, 'w' )
         kkey = '%s_%s' % (expt,var)
@@ -681 +734 @@
         oo.write( '%s %s\n' % ('operation_qualifiers ', ' ' ) )
         oo.write( '%s %s\n' % ('input_frequency ', ifreq ) )
+        oo.write( '%s %s\n' % ('input_esgf_dsid ', esgf_dsid ) )
         oo.write( '%s %s\n' % ('model_id ', model ) )
         oo.write( '%s %s\n' % ('experiment_id ', expt ) )
@@ -694 +748 @@
         oo.write( '%s %s\n' % ('input_dir ', vdir ) )
         oo.write( '%s %s\n' % ('plurality ', plurality ) )
+        oo.write( '%s %s\n' % ('proc_code_repository_url ', 'http://proj.badc.rl.ac.uk/svn/exarch/CedaCdoScripts' )
+        oo.write( '%s %s\n' % ('proc_code_home_url ', 'http://proj.badc.rl.ac.uk/exarch/wiki/PackageRBPS' )
         oo.write( '%s %s\n' % ('proc_method_url ', '%s/documentation/%s/%s' % (mdbase, plurality, doc_filename)  ) )
         oo.write( '%s %s\n' % ('proc_log_url ', '%s/logfiles/%s/%s/%s' % (mdbase, plurality, optag, log_filename)  ) )
@@ -715 +771 @@
             attr += ',' + dat
           attr += ']'
-          oo.write( '%s %s\n' % ('input_files ', attr) )
+          writeat(oo, 'input_files ', attr, log )
         else:
-          oo.write( '%s %s\n' % ('input_files ', ifiles[kkey] ) )
+          writeat(oo, 'input_files ', ifiles[kkey], log )
 
         if len( tids[kkey] ) > 511:
@@ -726 +782 @@
         oo.write( '%s %s\n' % ('input_creation_date ', crds[kkey]) )
         oo.write( '%s %s\n' % ('tracking_id ', tid) )
-        oo.write( '%s %s; %s\n' % ('history ', hists[kkey], '%s: %s: %s' % (creation_date, 'CDO script', proc_description)) )
+        writeat( oo, 'history ', '%s; %s' % (string.replace(hists[kkey],'\n',';;'), '%s: %s: %s' % (creation_date, 'CDO script', proc_description)), log )
         oo.close()
         tmpofilename = tmpdir + '.cdo_basic_temp.nc'
@@ -745 +801 @@
           assert os.path.isfile( string.strip(ofiles[kkey]) ), 'File "%s" not found [2], kkey = %s' % (ofiles[kkey], kkey)
           cmd = '/usr/local/2/bin/cdo -s setgatts,%s %s %s/%s' % (attfile,ofiles[kkey],opdir,ofilename)
-          print cmd
+          log.info( 'Command: %s\n' %  cmd )
           rc = run_cmd( cmd )
 
-log.info( 'Script completed' )
-rcout.close()
-rcerr.close()
-fh.close()
-for k in eeproc.keys():
-  sh[procid][k] = eeproc[k]
-print sh.keys()
-sh.close()
+        ###cmd = '/usr/local/2/bin/ncatted -a cell_methods,%s,o,c,%s %s' % (var,cmval,ofilename)
+        cmd = ['/usr/local/2/bin/ncatted','-a','cell_methods,%s,o,c,%s' % (var,cmval),'%s/%s' % (opdir,ofilename)]
+        rc = run_cmd( cmd )
+
+def close_all():
+  log.info( 'Script completed' )
+  fh.close()
+  rcout.close()
+  rcerr.close()
+  for k in eeproc.keys():
+    sh[procid][k] = eeproc[k]
+  sh.close()
+
+if go:
+  log.info( 'Merging files' )
+  for expt in expts:
+    proc_expt(expt,version,optag)
+
+close_all()
 
 os.popen( 'rm %scdo_basic_lock' % tmpdir ).readlines()