expand_mixingratios.py

#!/usr/bin/env python
# expand_fluxes.py
import sys
sys.path.append('../../')
import os
import getopt
from datetime import datetime, timedelta
from da.tools.general import CreateDirs
import numpy as np
from pylab import date2num, num2date

"""
Author: Wouter Peters (Wouter.Peters@wur.nl)

Revision History:
File created on 11 May 2012.

"""

def WriteMixingRatios(DaCycle):
    """ 
    
    Write Sample information to NetCDF files. These files are organized by site and 
    have an unlimited time axis to which data is appended each cycle.

    The needed information is obtained from the sampleinfo.nc files and the original input data files from ObsPack. 

    The steps are:

    (1) Create a directory to hold timeseries output files
    (2) Read the sampleinfo.nc file for this cycle and get a list of original files they were obtained from
    (3) For each file, copy the original data file from ObsPack (if not yet present)
    (4) Open the copied file, find the index of each observation, fill in the simulated data
    
    """
    import da.tools.io4 as io
    from string import join
    import shutil
    import logging

       
    dirname     = 'data_molefractions'
#
    dirname     = CreateDirs(os.path.join(DaCycle['dir.analysis'],dirname))
#
# Some help variables
#
    dectime0    = date2num(datetime(2000,1,1))
    dt          = DaCycle['cyclelength']
    startdate   = DaCycle['time.start'] 
    enddate     = DaCycle['time.end'] 

    msg = "DA Cycle start date is %s"   % startdate.strftime('%Y-%m-%d %H:%M')      ; logging.debug(msg)
    msg = "DA Cycle end   date is %s"   % enddate.strftime('%Y-%m-%d %H:%M')        ; logging.debug(msg)

    DaCycle['time.sample.stamp']    = "%s_%s"%(startdate.strftime("%Y%m%d%H"),enddate.strftime("%Y%m%d%H"),)

    # Step (1): Get the sample output data file for this cycle

    infile   = os.path.join(DaCycle['dir.output'],'sampleinfo_%s.nc'% DaCycle['time.sample.stamp'])

    ncf_in   = io.CT_Read(infile,'read')

    obs_num   = ncf_in.GetVariable('obs_num')
    obs_val   = ncf_in.GetVariable('observed')
    mdm       = ncf_in.GetVariable('modeldatamismatch')
    simulated = ncf_in.GetVariable('modelsamples')
    infilename= ncf_in.GetVariable('inputfilename')

    infiles   = [join(s.compressed(),'') for s in infilename]

    ncf_in.close()

    for orig_file in set(infiles):

        if not os.path.exists(orig_file):
            msg = "The original input file (%s) could not be found, continuing to next file..."%orig_file ; logging.error(msg)
            continue


        copy_file = os.path.join(dirname,os.path.split(orig_file)[-1])
        if not os.path.exists(copy_file):
            shutil.copy(orig_file,copy_file)
            msg = "Copied a new original file (%s) to the analysis directory"%orig_file ; logging.debug(msg)

            ncf_out = io.CT_CDF(copy_file,'write')

            # Modify the attributes of the file to reflect added data from CTDAS properly
     

            ncf_out.Caution  = '==================================================================================='
            try:
                ncf_out.History += '\nOriginal observation file modified by user %s on %s\n' % (os.environ['USER'], datetime.today().strftime('%F'),)
            except:
                ncf_out.History = '\nOriginal observation file modified by user %s on %s\n' % (os.environ['USER'], datetime.today().strftime('%F'),)
            ncf_out.CTDAS_info = 'Simulated values added from a CTDAS run by %s on %s\n' % (os.environ['USER'], datetime.today().strftime('%F'),)\
                               + '\nCTDAS was run on platform %s' % (os.environ['HOST'],)\
                               + '\nCTDAS job directory was %s' % (DaCycle['dir.da_run'],)\
                               + '\nCTDAS Da System was %s' % (DaCycle['da.system'],)\
                               + '\nCTDAS Da ObsOperator was %s' % (DaCycle['da.obsoperator'],)
            ncf_out.CTDAS_startdate = DaCycle['time.start'].strftime('%F')
            ncf_out.CTDAS_enddate   = DaCycle['time.finish'].strftime("%F")
            ncf_out.original_file = orig_file


            # get nobs dimension

            if ncf_out.dimensions.has_key('id'): 
                dimidob     = ncf_out.dimensions['id']
                dimid       = ('id',)
            elif ncf_out.dimensions.has_key('obs'): 
                dimidob     = ncf_out.dimensions['obs']
                dimid       = ('obs',)

            if dimidob.isunlimited:
                nobs = ncf_out.inq_unlimlen()
            else:
                nobs = len(dimid)


            # add nmembers dimension

            dimmembersob = ncf_out.createDimension('nmembers',size=simulated.shape[1])
            dimmembers   = ('nmembers',)
            nmembers     = len(dimmembers)

            # Create empty arrays

            savedict                = io.std_savedict.copy()
            savedict['name']        = "modeldatamismatch"
            savedict['long_name']   = "modeldatamismatch"
            savedict['units']       = "[mol mol-1]^2"
            savedict['dims']        = dimid
            savedict['comment']     = 'Variance of mole fractions resulting from model-data mismatch'
            status                  = ncf_out.AddVariable(savedict)


            savedict                = io.std_savedict.copy()
            savedict['name']        = "modelsamplesmean"
            savedict['long_name']   = "mean modelsamples"
            savedict['units']       = "mol mol-1"
            savedict['dims']        = dimid
            savedict['comment']     = 'simulated mixing ratios based on optimized state vector'
            status                  = ncf_out.AddVariable(savedict)

            savedict                = io.std_savedict.copy()
            savedict['name']        = "modelsamplesstandarddeviation"
            savedict['long_name']   = "standard deviaton of modelsamples over all ensemble members"
            savedict['units']       = "mol mol-1"
            savedict['dims']        = dimid
            savedict['comment']     = 'std dev of simulated mixing ratios based on optimized state vector'
            status                  = ncf_out.AddVariable(savedict)

            savedict                = io.std_savedict.copy()
            savedict['name']        = "modelsamplesensemble"
            savedict['long_name']   = "modelsamples over all ensemble members"
            savedict['units']       = "mol mol-1"
            savedict['dims']        = dimid+dimmembers
            savedict['comment']     = 'ensemble of simulated mixing ratios based on optimized state vector'
            status                  = ncf_out.AddVariable(savedict)

        else:
            msg = "Modifying existing file (%s) in the analysis directory"%copy_file ; logging.debug(msg)

            ncf_out = io.CT_CDF(copy_file,'write')


        # Get existing file obs_nums to determine match to local obs_nums

        if ncf_out.variables.has_key('id'):
            file_obs_nums = ncf_out.GetVariable('id')
        elif ncf_out.variables.has_key('obspack_num'):
            file_obs_nums = ncf_out.GetVariable('obspack_num')

        # Get all obs_nums related to this file, determine their indices in the local arrays

        selected_obs_nums = [num for infile,num in zip(infiles,obs_num) if infile == orig_file]

        # For each index, get the data and add to the file in the proper file index location

        for num in selected_obs_nums:

            model_index = obs_num.tolist().index(num)
            file_index  = file_obs_nums.tolist().index(num)

            var = ncf_out.variables['modeldatamismatch']
            var[file_index] = mdm[model_index]

            var = ncf_out.variables['modelsamplesmean']
            var[file_index] = simulated[model_index,0]

            var = ncf_out.variables['modelsamplesstandarddeviation']
            var[file_index] = simulated[model_index,1:].std()

            var = ncf_out.variables['modelsamplesensemble']
            var[file_index] = simulated[model_index,:]

        # close the file

        status  = ncf_out.close()

    return None



if __name__ == "__main__":

    import logging
    from da.tools.initexit import CycleControl
    from da.ct.dasystem import CtDaSystem

    sys.path.append('../../')

    logging.root.setLevel(logging.DEBUG)

    DaCycle = CycleControl(args={'rc':'../../da.rc'})
    DaCycle.Initialize()
    DaCycle.ParseTimes()

    DaSystem    = CtDaSystem('../rc/carbontracker.rc')
    DaSystem.Initialize()

    DaCycle.DaSystem    = DaSystem

    while DaCycle['time.start'] < DaCycle['time.finish']:

        outfile = WriteMixingRatios(DaCycle)

        DaCycle.AdvanceCycleTimes()

    sys.exit(0)