Changeset 57 for CedaCdoScripts


Ignore:
Timestamp:
23/09/13 17:14:25 (6 years ago)
Author:
mjuckes
Message:

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

Location:
CedaCdoScripts/Misc
Files:
3 added
1 edited

Legend:

Unmodified
Added
Removed
  • CedaCdoScripts/Misc/cdo_basic.py

    r54 r57  
    88debug = True 
    99cleanUpAsUGo = True 
    10 args = getopt.getopt( sys.argv[1:], 'i:m:e:r:b:x:o:', ['onlyVars=','scratch=','omitVars='] ) 
     10args = getopt.getopt( sys.argv[1:], 'i:m:e:r:b:x:o:', ['onlyVars=','scratch=','omitVars=','version=','log='] ) 
    1111ad = { 'i':'inst', 'm':'model', 'e':'ensem', 'r':'realm', 'b':'batch', 'x':'expt', 'o':'outDir' } 
    1212 
    13 knownops = ['am','bm','mclim5','cdd','tr'] 
     13knownops = ['am','sm','bm','mclim5','cdd','tr', 'fd'] 
    1414optag = args[-1][0] 
    1515assert optag in knownops, 'First argument [%s] should be operation in: %s' % (optag, str(knmownops)) 
     
    2222fullYears = True 
    2323 
    24 print args 
    2524edict = {} 
    2625for a in args[0]: 
    2726  edict[a[0][1:]] = a[1]  
     27 
     28def writeat( fh, tag, str,log): 
     29  if len(str) + len(tag) + 2 > 1024 : 
     30    log.warn( 'TRUNCATING global attribute %s: %s' % (tag,str) ) 
     31    strl = 1022 - len(tag) 
     32    str = str[:strl] 
     33  fh.write( '%s %s\n' % (tag,str) ) 
    2834 
    2935fullYearNeeded = False 
     
    5056##batch = 'b011_' 
    5157##onlyVars = ['tas'] 
     58project = 'cmip5' 
    5259ensem = edict['e'] 
    5360inst = edict['i'] 
     
    5663experiment = edict['x'] 
    5764batch = edict['b'] + '_' 
     65processing_batch = edict['b'] 
    5866opdir = edict['o'] 
    5967onlyVars = None 
     
    6472## qcDir = '%s%s/%s/%s/%s/' % (qcBaseDir,product,inst, model, expt) 
    6573 
    66 tstring = '%s%2.2i%2.2i-%2.2i%2.2i%2.2i' % time.gmtime()[0:6] 
    6774vstring = '%s%2.2i%2.2i' % time.gmtime()[0:3] 
    6875doctag = optag 
    6976if optag in ['cdd']: 
    7077  doctag = 'precip_indices' 
    71 if optag in ['tr']: 
     78if optag in ['tr','fd']: 
    7279  doctag = 'temperature_indices' 
    7380doc_filename = 'ceda_cmip5_processing_%s.pdf' % doctag 
    74 log_filename = 'ceda_cmip5_processing_log_%s_%s.txt' % (optag,tstring) 
     81tstring = '%s%2.2i%2.2i-%2.2i%2.2i%2.2i' % time.gmtime()[0:6] 
     82log_filename = edict.get( '-log', 'ceda_cmip5_processing_log_%s_%s.txt' % (optag,tstring) ) 
    7583batching='byFile' 
    7684ofreq = 'yr' 
     
    7987  batching = '5yr' 
    8088  monthsPerBatch = 60 
    81 elif optag in ['cdd','tr']: 
     89elif optag in ['cdd','tr','fd']: 
    8290  batching = 'yr' 
    8391  ifreq = 'day' 
     
    9098 
    9199scriptName = 'cdo_basic.py' 
    92 scriptVersion = 0.8 
     100scriptVersion = 0.9 
    93101log = logging.getLogger( __name__ ) 
    94102fh = logging.FileHandler( log_filename ) 
     
    108116   omitVars = string.split( edict['-omitVars'], ',' ) 
    109117   log.info( 'Omitting variables: %s' % edict['-omitVars'] ) 
     118if edict.has_key( '-version' ): 
     119   version = edict['-version'] 
     120   log.info( 'Processing version: %s' % edict['-version'] ) 
     121else: 
     122   version = None 
    110123## log.info( 'CDO version information:' ) 
    111124## os.popen( 'cdo -V > .cdo_version' ).readlines() 
     
    123136if not os.path.isdir(tmpdir): 
    124137  os.mkdir( tmpdir ) 
    125   print 'Creating temporary directory: %s ' % tmpdir 
     138  log.info( 'Creating temporary directory: %s ' % tmpdir ) 
    126139 
    127140if tmpdir[-1] != '/': 
     
    133146assert os.path.isfile( '%scdo_basic_lock' % tmpdir ), 'Failed to make lock file in work directory %s' % tmpdir 
    134147 
    135  
    136148plurality = 'single' 
    137149mdbase = "http://badc.nerc.ac.uk/browse/badc/cmip5/metadata/processing/"  
     
    139151rcout = open( 'cdo_stdout.txt', 'w' ) 
    140152rcerr = open( 'cdo_stderr.txt', 'w' ) 
    141 def run_cmd( cmd, comment='', strict=True ): 
     153def run_cmd( cmd, comment='', strict=True, append=None, out=None ): 
    142154    if type(cmd) == type('x'): 
    143155      args = map( lambda x: string.strip(x, "'"), string.split( cmd )) 
     
    145157      args = cmd 
    146158    log.info( cmd ) 
    147     proc = subprocess.Popen( args, stdout=rcout, stderr=rcerr, env={ 'LD_LIBRARY_PATH':'/usr/local/lib'} ) 
     159    if append != None: 
     160      fho = open( append, 'a' ) 
     161      clo = True 
     162    elif out != None: 
     163      fho = open( out, 'w' ) 
     164      clo = True 
     165    else: 
     166      fho = rcout 
     167      clo = False 
     168    proc = subprocess.Popen( args, stdout=fho, stderr=rcerr, env={ 'LD_LIBRARY_PATH':'/usr/local/lib'} ) 
    148169    std = proc.communicate() 
    149170    rc = proc.returncode 
     
    156177      elif rc != 0: 
    157178        print 'Failed (%s) [rc=%s] executing %s' % (comment,rc,str(cmd)) 
     179    if clo: 
     180      fho.close() 
    158181    return rc 
    159182 
    160183dir_tmpl_s = '/disks/fuse/cmip5/output1/MOHC/HadGEM2-ES/%(expt)s/%(ifreq)s/atmos/Amon/' 
     184dir_tmpl_s = '/badc/cmip5/data/cmip5/output1/MOHC/HadGEM2-ES/%(expt)s/%(ifreq)s/atmos/Amon/' 
    161185if ifreq == 'mon': 
    162186  mip = {'atmos':'Amon', 'ocean':'Omon', 'land':'Lmon' }[realm] 
     
    166190  raise 'Should not be here' 
    167191 
    168 dir_tmpl_r = '/disks/fuse/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/' 
    169 dir_tmpl_e = '/disks/fuse/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/latest/' 
     192dir_tmpl_r = '/badc/cmip5/data/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/' 
     193dir_tmpl_e = '/badc/cmip5/data/cmip5/output1/%(inst)s/%(model)s/%(expt)s/%(ifreq)s/%(realm)s/%(mip)s/%(ensem)s/latest/' 
    170194dir_tmpl = dir_tmpl_e + '%(var)s/' 
    171195bldir_tmpl = '/home/users/prototype/cmip_processing/cmip5_bl/output1_%(inst)s_%(model)s_%(expt)s/' 
     
    180204  if len(b2) == 0: 
    181205    return None 
    182   return b2[-1] 
    183  
    184 def get_version_dir(): 
     206  var = string.split( b2[-1], '_' )[0] 
     207  return var 
     208 
     209def get_version_dir(version=None): 
    185210 
    186211  pdict = { 'inst':inst, 'model':model, 'expt':expt, 'ensem':ensem, 'ifreq':ifreq, 'realm':realm, 'mip':mip } 
    187212  rdir = dir_tmpl_r % pdict 
     213  if version != None: 
     214    this_v = rdir + version 
     215    assert os.path.isdir( this_v ), 'Version directory %s not found' % this_v 
     216    return this_v 
     217 
    188218  vlist = glob.glob( rdir + 'v*' ) 
    189219  assert len(vlist) > 0, 'No version directories found in %s' % rdir 
     
    200230  ## edir = dir_tmpl % { 'inst':inst, 'model':model, 'expt':expt, 'var':var, 'ensem':ensem, 'ifreq':ifreq, 'realm':realm, 'mip':mip } 
    201231  assert os.path.isdir( edir), 'Error in get_expt_fl: %s is not a directory' % edir 
     232  log.info( 'edir:: %s' % edir ) 
    202233  fl = glob.glob( edir + '*.nc' ) 
    203234  assert len(fl) > 0, 'No netcdf files found in %s' % edir 
     
    211242    if os.path.isfile( bldir + blfile ): 
    212243      blvars = map( get_var_from_bl_record, open( bldir + blfile ).readlines() ) 
    213   if var in blvars: 
    214     print 'skipping ',var 
    215     log.info( 'Skipping black-listed ADS expt=%s, var=%s' % (expt,var) ) 
    216     return fm 
     244      if var in blvars: 
     245        log.info( 'Skipping black-listed ADS expt=%s, var=%s' % (expt,var) ) 
     246        return fm 
     247      else: 
     248        log.info( 'Var %s not in black list [%s] %s' % (var, bldir + blfile, str(blvars)) ) 
     249    else: 
     250      log.info( 'Black list file %s not present' % (bldir+blfile) ) 
     251  else: 
     252    log.info( 'Black list directory %s not present' % bldir ) 
    217253  for f in fl: 
    218254    fn = string.split(f, '/')[-1] 
     
    268304   onlyVars = ['pr'] 
    269305   omitVars = None 
    270 if optag in ['tr']: 
     306if optag in ['tr','fd']: 
    271307   onlyVars = ['tasmin'] 
    272308   omitVars = None 
     
    274310for expt in expts: 
    275311  ##dir = dir_tmpl_e % locals() 
    276   dir = get_version_dir() 
     312  dir = get_version_dir( version=version) 
    277313  vars = glob.glob( '%s/*' % dir ) 
    278   print expt, dir, vars 
    279314  thisvarl = map( lambda x: string.split(x,'/')[-1], vars ) 
    280315  thisvarl.sort() 
     
    312347lon1,lon2,lat1,lat2 = (60.,120.,30.,60.) 
    313348cdoop, ofreq,opName, proc_description = { 'am':('divdpy -yearsum -muldpm', 'yr','Annual mean','Annual mean of monthly mean data'), \ 
     349                        'sm':('divdpy -seassum -muldpm', 'yr','Seasonal mean','Annual seasonal of monthly mean data'), \ 
    314350                        'bm':('fldmean -sellonlatbox,%s,%s,%s,%s' % (lon1,lon2,lat1,lat2),ifreq,'Box mean','Spatial mean over a latitude-longitude box'), \ 
    315351                        'mclim5':('ymonavg', '5yr','5-year monthly climatology','Monthly climatology for 5 year period: 5 year mean for each month'), 
    316352                        'cdd':('eca_cdd','yr','Continuous dry days','Maximum of continuous dry days in a year'), 
    317                         'tr':('eca_tr','yr','Tropical nights','Count of tropical nights per year')  }[optag] 
     353                        'tr':('eca_tr','yr','Tropical nights','Count of tropical nights per year'), 
     354                        'fd':('eca_fd','yr','Frost days','Count of frost days (Tmin below zero) per year') }[optag] 
    318355## opdir = '%s%s' % (batch,optag) 
    319356if go and not os.path.isdir( opdir ): 
     
    327364k=-1 
    328365print 'INFO0080: ',k, str(k), str( [1,2,3] ) 
    329 for expt,var,fm in ye: 
     366##for expt,var,fm in ye: 
     367def run_evf( expt,var,fm ): 
    330368    kk = 0 
    331     dir = get_version_dir() 
     369    dir = get_version_dir( version=version) 
    332370    ##edir = dir_tmpl % locals() 
    333371    edir = '%s/%s/' % (dir,var) 
     
    367405        ifiles[key] += '%s,' % file[0] 
    368406 
     407    startstart = None 
    369408    nbfb = 0 
    370409    nbfbb = 0 
     
    395434      nodir = False 
    396435      if nbfb > 0: 
    397         print 'FILE DATE BOUNDARIES DO NOT MATCH ASSUMPTIONS: %s, %s' % (expt,var) 
    398         print 'merging files' 
    399         print nbfb, nbfbb 
    400436        log.info( 'File date boundaries do not match end of year (%s,%s): merging: %s, %s' % (nbfb, nbfbb, expt,var) ) 
    401437        if nbfbb > 0: 
     
    441477            cmd1 = None 
    442478            cmd = '/usr/local/2/bin/cdo -s %s %s/%s %s' % (cdoop, tmpdir, '.tmp1.nc', ofile) 
    443           print cmd 
     479          log.info( 'cdo command: %s\n' % cmd ) 
    444480          ## raise 'bad file boundary' 
    445481            
     
    471507          nf = len(fm) 
    472508          if smem[0][0] != 1: 
    473             print 'Omitting first file: %s' % str(fm[0]) 
     509            log.info(  'Omitting first file: %s' % str(fm[0]) ) 
    474510            i0 = 1 
    475511            nf = len(fm)-1 
     
    514550          ofiles[key] += ' %s' % (ofile) 
    515551      if nbfb > 0: 
    516           print 'FILE DATE BOUNDARIES DO NOT MATCH ASSUMPTIONS: %s, %s' % (expt,var) 
     552          log.info( 'FILE DATE BOUNDARIES DO NOT MATCH ASSUMPTIONS: %s, %s' % (expt,var) ) 
    517553    elif batching in ['5yr','yr']: 
    518554      startstart = fm[0][1] 
     
    566602 
    567603        for k in range(len(tfm)-1): 
    568             print 'INFO0100: ',k, str(k), str( [1,2,3] ) 
    569             print 'INFO0101: ',k, tfm[k], tfm[k+1] 
     604            log.info( 'INFO0100: ',k, str(k), str( [1,2,3] ) ) 
     605            log.info( 'INFO0101: ',k, tfm[k], tfm[k+1] ) 
    570606            a = str(tfm[k]) 
    571607            b = str(tfm[k+1]) 
     
    587623            if debug: 
    588624              log.info( 'Trying to log file maximum value' ) 
    589               cmd0 = 'echo %s%s >> /home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' % (edir,tfm[0][0]) 
    590               cmd = '/usr/local/2/bin/cdo info %s 2>&1 >> /home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' % ofile 
    591               rc = run_cmd( cmdo ) 
    592               rc = run_cmd( cmd ) 
     625              cmdo = 'echo %s%s' % (edir,tfm[0][0]) 
     626              cmd = '/usr/local/2/bin/cdo info %s' % ofile 
     627              rc = run_cmd( cmdo, append='/home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' ) 
     628              rc = run_cmd( cmd, append='/home/users/prototype/cmip_processing/wrk/cdd1/debug.txt' ) 
    593629          else: 
    594630            print 'single: %s' % cmd 
     
    635671        else: 
    636672          raise 'Code missing' 
    637  
    638 attfile = '%s/temp_attribute_file.txt' % tmpdir 
    639 log.info( 'Merging files' ) 
    640 if go: 
    641   for expt in expts: 
     673    return tid, crd, table_id, realization, initialization_method, physics_version, startstart 
     674 
     675 
     676for expt,var,fm in ye: 
     677    tid, crd, table_id, realization, initialization_method, physics_version, startstart = run_evf( expt,var,fm ) 
     678 
     679def proc_expt(expt,version,optag): 
     680    attfile = '%s/temp_attribute_file.txt' % tmpdir 
    642681    varl = varld[expt] 
    643682    ## for var in varld[expt]: 
    644683    for var in varl: 
     684      if optag in ['am', 'sm']: 
     685        if var == 'tasmax': 
     686          cmval = "time: maximum within days time: mean over days" 
     687        elif var == 'tasmin': 
     688          cmval = "time: minimum within days time: mean over days" 
     689        else: 
     690          cmval = "time: mean" 
     691      else: 
     692        assert True, 'Need to set cell_methods for this operator' 
     693 
    645694      key = '%s_%s' % (expt,var) 
    646695      if key in badkeys: 
    647         print 'Skipping %s because or earlier error' % key 
     696        log.info( 'Skipping %s because or earlier error' % key ) 
    648697      elif string.strip( ofiles[key] ) == '': 
    649         print 'Skipping %s because no output files registered' % key 
     698        log.info( 'Skipping %s because no output files registered' % key  ) 
    650699        log.warning( 'WARN004: unexplained missing output file for key %s' % (key) ) 
    651700      else: 
    652         print 'Trying %s' % key 
    653         dir = get_version_dir() 
     701        log.info( 'Trying %s' % key ) 
     702        dir = get_version_dir( version=version) 
    654703    ##    vdir = dir_tmpl % locals() 
    655704        vdir = '%s/%s/' % (dir,var) 
     
    667716        for k in range(36): 
    668717          tid += cx[random.randint(0,len(cx)-1 )] 
     718        this_mip  = string.split( table_id )[1] 
     719        version = string.split( string.strip( dir, '/' ), '/' )[-1] 
     720        esgf_dsid = '%s.%s.%s.%s.%s.%s.%s.%s.%s.%s' % (project,product,inst,expt,model,ifreq,realm,this_mip,ensem,version) 
     721##cmip5.output1.MIROC.MIROC5.sstClim.monClim.atmos.Amon.r1i1p1.v20111129.xml 
    669722        oo = open( attfile, 'w' ) 
    670723        kkey = '%s_%s' % (expt,var) 
     
    681734        oo.write( '%s %s\n' % ('operation_qualifiers ', ' ' ) ) 
    682735        oo.write( '%s %s\n' % ('input_frequency ', ifreq ) ) 
     736        oo.write( '%s %s\n' % ('input_esgf_dsid ', esgf_dsid ) ) 
    683737        oo.write( '%s %s\n' % ('model_id ', model ) ) 
    684738        oo.write( '%s %s\n' % ('experiment_id ', expt ) ) 
     
    694748        oo.write( '%s %s\n' % ('input_dir ', vdir ) ) 
    695749        oo.write( '%s %s\n' % ('plurality ', plurality ) ) 
     750        oo.write( '%s %s\n' % ('proc_code_repository_url ', 'http://proj.badc.rl.ac.uk/svn/exarch/CedaCdoScripts' ) 
     751        oo.write( '%s %s\n' % ('proc_code_home_url ', 'http://proj.badc.rl.ac.uk/exarch/wiki/PackageRBPS' ) 
    696752        oo.write( '%s %s\n' % ('proc_method_url ', '%s/documentation/%s/%s' % (mdbase, plurality, doc_filename)  ) ) 
    697753        oo.write( '%s %s\n' % ('proc_log_url ', '%s/logfiles/%s/%s/%s' % (mdbase, plurality, optag, log_filename)  ) ) 
     
    715771            attr += ',' + dat 
    716772          attr += ']' 
    717           oo.write( '%s %s\n' % ('input_files ', attr) ) 
     773          writeat(oo, 'input_files ', attr, log ) 
    718774        else: 
    719           oo.write( '%s %s\n' % ('input_files ', ifiles[kkey] ) ) 
     775          writeat(oo, 'input_files ', ifiles[kkey], log ) 
    720776 
    721777        if len( tids[kkey] ) > 511: 
     
    726782        oo.write( '%s %s\n' % ('input_creation_date ', crds[kkey]) ) 
    727783        oo.write( '%s %s\n' % ('tracking_id ', tid) ) 
    728         oo.write( '%s %s; %s\n' % ('history ', hists[kkey], '%s: %s: %s' % (creation_date, 'CDO script', proc_description)) ) 
     784        writeat( oo, 'history ', '%s; %s' % (string.replace(hists[kkey],'\n',';;'), '%s: %s: %s' % (creation_date, 'CDO script', proc_description)), log ) 
    729785        oo.close() 
    730786        tmpofilename = tmpdir + '.cdo_basic_temp.nc' 
     
    745801          assert os.path.isfile( string.strip(ofiles[kkey]) ), 'File "%s" not found [2], kkey = %s' % (ofiles[kkey], kkey) 
    746802          cmd = '/usr/local/2/bin/cdo -s setgatts,%s %s %s/%s' % (attfile,ofiles[kkey],opdir,ofilename) 
    747           print cmd 
     803          log.info( 'Command: %s\n' %  cmd ) 
    748804          rc = run_cmd( cmd ) 
    749805 
    750 log.info( 'Script completed' ) 
    751 rcout.close() 
    752 rcerr.close() 
    753 fh.close() 
    754 for k in eeproc.keys(): 
    755   sh[procid][k] = eeproc[k] 
    756 print sh.keys() 
    757 sh.close() 
     806        ###cmd = '/usr/local/2/bin/ncatted -a cell_methods,%s,o,c,%s %s' % (var,cmval,ofilename) 
     807        cmd = ['/usr/local/2/bin/ncatted','-a','cell_methods,%s,o,c,%s' % (var,cmval),'%s/%s' % (opdir,ofilename)] 
     808        rc = run_cmd( cmd ) 
     809 
     810def close_all(): 
     811  log.info( 'Script completed' ) 
     812  fh.close() 
     813  rcout.close() 
     814  rcerr.close() 
     815  for k in eeproc.keys(): 
     816    sh[procid][k] = eeproc[k] 
     817  sh.close() 
     818 
     819if go: 
     820  log.info( 'Merging files' ) 
     821  for expt in expts: 
     822    proc_expt(expt,version,optag) 
     823 
     824close_all() 
    758825 
    759826os.popen( 'rm %scdo_basic_lock' % tmpdir ).readlines() 
Note: See TracChangeset for help on using the changeset viewer.