#!/usr/bin/env python
# obs.py

"""
Author : peters 

Revision History:
File created on 28 Jul 2010.

"""
import os
import sys
import logging
        
import datetime as dtm
from string import strip
from numpy import array, ndarray, mean, logical_and
sys.path.append(os.getcwd())
sys.path.append('../../')

identifier = 'CarbonTracker CO2 and C13 mole fractions'
version = '0.0'

from da.baseclasses.obs import Observations
import da.tools.io4 as io
import da.tools.rc as rc
################### Begin Class ObsPackObservations ###################

class ObsPackObservations(Observations):
    """ an object that holds data + methods and attributes needed to manipulate mole fraction values """

    def setup(self, dacycle):

        self.startdate = dacycle['time.sample.start']
        self.enddate = dacycle['time.sample.end']

        op_id = dacycle.dasystem['obspack.input.id']
        op_dir = dacycle.dasystem['obspack.input.dir']

        op_c13_id = dacycle.dasystem['obspack.c13.input.id']
        op_c13_dir = dacycle.dasystem['obspack.c13.input.dir']

        if not os.path.exists(op_dir):
            msg = 'Could not find  the required ObsPack distribution (%s) ' % op_dir
            logging.error(msg)
            raise IOError, msg
        else:
            self.obspack_dir = op_dir
            self.obspack_id = op_id
            self.obspackc13_dir = op_c13_dir
            self.obspackc13_id = op_c13_id
        self.datalist = []
        self.datalist_c13 = [] #test 
        self.datalist_co2 = [] #test 


    def getvalues_co2(self, name, constructor=array):

        result = constructor([getattr(o, name) for o in self.datalist_co2])
        if isinstance(result, ndarray):
            return result.squeeze()
        else:
            return result


    def getvalues_c13(self, name, constructor=array):

        result = constructor([getattr(o, name) for o in self.datalist_c13])
        if isinstance(result, ndarray):
            return result.squeeze()
        else:
            return result


    def add_observations(self):
        """ Returns a MoleFractionList holding individual MoleFractionSample objects for all obs in a file
      
            The ObsPack mole fraction files are provided as time series per site with all dates in sequence. 
            We will loop over all site files in the ObsPackage, and subset each to our needs
            
        """

        # Step 1: Read list of available site files in package

        infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,))
        f = open(infile, 'r')
        lines = f.readlines()
        f.close()

        ncfilelist = []
        for line in lines:
            if line.startswith('#'): continue # header

            items = line.split()
            #ncfile, lab , start_date, stop_date, data_comparison = items[0:5]
            ncfile, lab , start_date, stop_date, data_comparison= line[:105].split()


            ncfilelist += [ncfile]

        logging.debug("ObsPack dataset info read, proceeding with %d netcdf files" % len(ncfilelist))

        for ncfile in ncfilelist:

            infile = os.path.join(self.obspack_dir, 'data', 'nc', ncfile + '.nc')
            ncf = io.ct_read(infile, 'read')
            idates = ncf.get_variable('time_components')
            dates = array([dtm.datetime(*d) for d in idates])

            subselect = logical_and(dates >= self.startdate , dates <= self.enddate).nonzero()[0]

            dates = dates.take(subselect, axis=0)
            
            obspacknum = ncf.get_variable('obspack_num').take(subselect)  # or should we propagate obs_num which is not unique across datasets??
            obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0)
            obspackid = [s.tostring().lower() for s in obspackid]
            obspackid = map(strip, obspackid)
            datasetname = ncfile  # use full name of dataset to propagate for clarity
            lats = ncf.get_variable('latitude').take(subselect, axis=0)
            lons = ncf.get_variable('longitude').take(subselect, axis=0)
            alts = ncf.get_variable('altitude').take(subselect, axis=0)
            obs = ncf.get_variable('value').take(subselect, axis=0)
            species = ncf.get_attribute('dataset_parameter')
            flags = ncf.get_variable('obs_flag').take(subselect, axis=0)
            ncf.close()

            for n in range(len(dates)): 
                self.datalist.append(MoleFractionSample(obspacknum[n], dates[n], datasetname, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], obspackid[n], species, 1, 0.0, infile)) 
                self.datalist_co2.append(MoleFractionSample(obspacknum[n], dates[n], datasetname, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], obspackid[n], species, 1, 0.0, infile))

            logging.debug("Added %d observations from file (%s) to the Data list" % (len(dates), ncfile)) 

        logging.info("CO2 observations list now holds %d values" % len(self.datalist_co2))
        logging.info("Total observations list now holds %d values" % len(self.datalist))

    def add_simulations(self, filename, silent=False):
        """ Adds model simulated values to the mole fraction objects """


        if not os.path.exists(filename):
            msg = "Sample output filename for observations could not be found : %s" % filename 
            logging.error(msg)
            logging.error("Did the sampling step succeed?")
            logging.error("...exiting")
            raise IOError, msg

        ncf = io.ct_read(filename, method='read')
        ids = ncf.get_variable('obs_num')
        simulated = ncf.get_variable('flask')
        obs_species = self.getvalues('species')
        dim_ens = len(simulated[0,:])/2
	simulated = simulated[:,:dim_ens] 	
        ncf.close()
        logging.info("Successfully read data from model sample file (%s)" % filename)

        obs_ids = self.getvalues_co2('id').tolist()
        ids = map(int, ids)

        missing_samples = []

        for idx, val,obs_specie in zip(ids, simulated, obs_species): 
            if idx in obs_ids and obs_specie == 'co2':
                index = obs_ids.index(idx)
                logging.debug("add_simulations %d, %d, %s" %(index,idx, mean(array(val))))
                self.datalist[index].simulated = val  # in mol/mol
            else:     
                missing_samples.append(idx)

        if not silent and missing_samples != []:
            logging.warning('Model samples were found that did not match any ID in the observation list. Skipping them...')

        logging.debug("Added %d simulated values to the Data list" % (len(ids) - len(missing_samples)))

    def write_sample_coords(self, obsinputfile):
        """ 
            Write the information needed by the observation operator to a file. Return the filename that was written for later use

        """

        f = io.CT_CDF(obsinputfile, method='create')
        logging.debug('Creating new observations file for ObservationOperator (%s)' % obsinputfile)

        dimid = f.add_dim('obs', len(self.datalist))
        dim200char = f.add_dim('string_of200chars', 200)
        dim10char = f.add_dim('string_of10chars', 10)
        dimcalcomp = f.add_dim('calendar_components', 6)

        if len(self.datalist) == 0:
            f.close()
            #return obsinputfile

        for key, value in self.site_move.iteritems():
            msg = "Site is moved by %3.2f degrees latitude and %3.2f degrees longitude" % value 
            f.add_attribute(key, msg)

        data = self.getvalues('id')

        savedict = io.std_savedict.copy() 
        savedict['name'] = "obs_num"
        savedict['dtype'] = "int"
        savedict['long_name'] = "Unique_Dataset_observation_index_number"
        savedict['units'] = ""
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
        f.add_data(savedict)

        data = [[d.year, d.month, d.day, d.hour, d.minute, d.second] for d in self.getvalues('xdate') ]

        savedict = io.std_savedict.copy() 
        savedict['dtype'] = "int"
        savedict['name'] = "date_components"
        savedict['units'] = "integer components of UTC date/time"
        savedict['dims'] = dimid + dimcalcomp
        savedict['values'] = data
        savedict['missing_value'] = -9
        savedict['comment'] = "Calendar date components as integers. Times and dates are UTC." 
        savedict['order'] = "year, month, day, hour, minute, second"
        f.add_data(savedict)

        data = self.getvalues('lat')

        savedict = io.std_savedict.copy() 
        savedict['name'] = "latitude"
        savedict['units'] = "degrees_north"
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['missing_value'] = -999.9
        f.add_data(savedict)

        data = self.getvalues('lon')

        savedict = io.std_savedict.copy() 
        savedict['name'] = "longitude"
        savedict['units'] = "degrees_east"
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['missing_value'] = -999.9
        f.add_data(savedict)

        data = self.getvalues('height')

        savedict = io.std_savedict.copy() 
        savedict['name'] = "altitude"
        savedict['units'] = "meters_above_sea_level"
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['missing_value'] = -999.9
        f.add_data(savedict)

        data = self.getvalues('samplingstrategy')

        savedict = io.std_savedict.copy() 
        savedict['dtype'] = "int"
        savedict['name'] = "sampling_strategy"
        savedict['units'] = "NA"
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['missing_value'] = -9
        f.add_data(savedict)

        data = self.getvalues('evn')

        savedict = io.std_savedict.copy() 
        savedict['dtype'] = "char"
        savedict['name'] = "obs_id"
        savedict['units'] = "ObsPack datapoint identifier"
        savedict['dims'] = dimid + dim200char
        savedict['values'] = data
        savedict['missing_value'] = '!'
        f.add_data(savedict)

        f.close()

        logging.debug("Successfully wrote data to obs file")
        logging.info("Sample input file for obs operator now in place [%s]" % obsinputfile)

        

    def add_model_data_mismatch(self, filename):
        """ 
            Get the model-data mismatch values for this cycle.

                (1) Open a sites_weights file
                (2) Parse the data
                (3) Compare site list against data
                (4) Take care of double sites, etc

        """    

        if not os.path.exists(filename):
            msg = 'Could not find  the required sites.rc input file (%s) ' % filename
            logging.error(msg)
            raise IOError, msg
        else:
            self.sites_file = filename

        sites_weights = rc.read(self.sites_file)

        self.rejection_threshold = int(sites_weights['obs.rejection.threshold'])
        self.global_R_scaling = float(sites_weights['global.R.scaling'])
        self.n_site_categories = int(sites_weights['n.site.categories'])

        logging.debug('Model-data mismatch rejection threshold: %d ' % self.rejection_threshold)
        logging.warning('Model-data mismatch scaling factor     : %f ' % self.global_R_scaling)
        logging.debug('Model-data mismatch site categories    : %d ' % self.n_site_categories)
   
        cats = [k for k in sites_weights.keys() if 'site.category' in k] 

        site_categories = {}
        for key in cats:
            name, error, may_localize, may_reject = sites_weights[key].split(';')
            name = name.strip().lower()
            error = float(error)
            may_reject = ("TRUE" in may_reject.upper())
            may_localize = ("TRUE" in may_localize.upper())
            site_categories[name] = {'category': name, 'error': error, 'may_localize': may_localize, 'may_reject': may_reject}

        site_info = {}
        site_move = {}
        site_hourly = {}   # option added to include only certain hours of the day (for e.g. PAL) IvdL
        site_incalt = {} # option to increase sampling altitude for sites specified in sites and weights file 
        for key, value in sites_weights.iteritems():
            if 'co2_' in key or 'sf6' in key:  # to be fixed later, do not yet know how to parse valid keys from rc-files yet.... WP
                sitename, sitecategory = key, value
                sitename = sitename.strip()
                sitecategory = sitecategory.split()[0].strip().lower()
                site_info[sitename] = site_categories[sitecategory]
            if 'site.move' in key:
                identifier, latmove, lonmove = value.split(';')
                site_move[identifier.strip()] = (float(latmove), float(lonmove))
            if 'site.hourly' in key:
                identifier, hourfrom, hourto = value.split(';')
                site_hourly[identifier.strip()] = (int(hourfrom), int(hourto))
            if 'site.incalt' in key:
                identifier, incalt = value.split(';')
                site_incalt[identifier.strip()] = (int(incalt))
                logging.debug('test')

        for obs in self.datalist:  # loop over all available data points
            if obs.species == 'co2':
                obs.mdm = 1000.0  # default is very high model-data-mismatch, until explicitly set by script
                obs.flag = 99  # default is do-not-use , until explicitly set by script
                exclude_hourly = False # default is that hourly values are not included

                identifier = obs.code
                species, site, method, lab, datasetnr = identifier.split('_')

                if site_info.has_key(identifier):
                    if site_hourly.has_key(identifier):
                        hourf, hourt = site_hourly[identifier]
                        if int(obs.xdate.hour) >= hourf and int(obs.xdate.hour) <= hourt:
                            logging.warning("Observation in hourly dataset INCLUDED, while sampling time %s was between %s:00-%s:00"%(obs.xdate.time(),hourf,hourt))
                        else:
                            logging.warning("Observation in hourly dataset EXCLUDED, while sampling time %s was outside %s:00-%s:00"%(obs.xdate.time(),hourf,hourt))
                            exclude_hourly = True
                    if site_info[identifier]['category'] == 'do-not-use' or exclude_hourly:
                        logging.warning("Observation found (%s, %d), but not used in assimilation !!!" % (identifier, obs.id))
                        obs.mdm = site_info[identifier]['error'] * self.global_R_scaling
                        obs.may_localize = site_info[identifier]['may_localize']
                        obs.may_reject = site_info[identifier]['may_reject']
                        obs.flag = 99
                    else:
                        logging.warning("Observation found (%s, %d)" % (identifier, obs.id))
                        obs.mdm = site_info[identifier]['error'] * self.global_R_scaling
                        obs.may_localize = site_info[identifier]['may_localize']
                        obs.may_reject = site_info[identifier]['may_reject']
                        obs.flag = 0

                else:
                    logging.warning("Observation NOT found (%s, %d), please check sites.rc file (%s)  !!!" % (identifier, obs.id, self.sites_file))

                if site_move.has_key(identifier):

                    movelat, movelon = site_move[identifier]
                    obs.lat = obs.lat + movelat
                    obs.lon = obs.lon + movelon

                    logging.warning("Observation location for (%s, %d), is moved by %3.2f degrees latitude and %3.2f degrees longitude" % (identifier, obs.id, movelat, movelon))
                if site_incalt.has_key(identifier):
                    incalt = site_incalt[identifier]
                    obs.height = obs.height + incalt 

                    logging.warning("Observation location for (%s, %d), is moved by %3.2f meters in altitude" % (identifier, obs.id, incalt))




        # Add site_info dictionary to the Observation object for future use

        self.site_info = site_info
        self.site_move = site_move
        self.site_hourly = site_hourly
        self.site_incalt = site_incalt

        logging.debug("Added Model Data Mismatch to all samples ")

    def write_sample_auxiliary(self, auxoutputfile):
        """ 
            Write selected information contained in the Observation object to a file. 

        """

        f = io.CT_CDF(auxoutputfile, method='create')
        logging.debug('Creating new auxiliary file for postprocessing (%s)' % auxoutputfile)

        dimid = f.add_dim('obs', len(self.datalist))
        dim200char = f.add_dim('string_of200chars', 200)
        dim10char = f.add_dim('string_of10chars', 10)
        dimcalcomp = f.add_dim('calendar_components', 6)

        if len(self.datalist) == 0:
            f.close()
            #return outfile

        for key, value in self.site_move.iteritems():
            msg = "Site is moved by %3.2f degrees latitude and %3.2f degrees longitude" % value 
            f.add_attribute(key, msg)

        data = self.getvalues('id')

        savedict = io.std_savedict.copy() 
        savedict['name'] = "obs_num"
        savedict['dtype'] = "int"
        savedict['long_name'] = "Unique_Dataset_observation_index_number"
        savedict['units'] = ""
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
        f.add_data(savedict)

        data = [[d.year, d.month, d.day, d.hour, d.minute, d.second] for d in self.getvalues('xdate')]

        savedict = io.std_savedict.copy() 
        savedict['dtype'] = "int"
        savedict['name'] = "date_components"
        savedict['units'] = "integer components of UTC date/time"
        savedict['dims'] = dimid + dimcalcomp
        savedict['values'] = data
        savedict['missing_value'] = -9
        savedict['comment'] = "Calendar date components as integers. Times and dates are UTC." 
        savedict['order'] = "year, month, day, hour, minute, second"
        f.add_data(savedict)

        data = self.getvalues('obs')

        savedict = io.std_savedict.copy()
        savedict['name'] = "observed"
        savedict['long_name'] = "observedvalues"
        savedict['units'] = "mol mol-1"
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['comment'] = 'Observations used in optimization'
        f.add_data(savedict)

        data = self.getvalues('mdm')

        savedict = io.std_savedict.copy()
        savedict['name'] = "modeldatamismatch"
        savedict['long_name'] = "modeldatamismatch"
        savedict['units'] = "[mol mol-1]"
        savedict['dims'] = dimid
        savedict['values'] = data.tolist()
        savedict['comment'] = 'Standard deviation of mole fractions resulting from model-data mismatch'
        f.add_data(savedict)

        data = self.getvalues('simulated') 

        dimmembers = f.add_dim('members', data.shape[1])

        savedict = io.std_savedict.copy()
        savedict['name'] = "modelsamples"
        savedict['long_name'] = "modelsamples for all ensemble members"
        savedict['units'] = "mol mol-1"
        savedict['dims'] = dimid + dimmembers
        savedict['values'] = data.tolist()
        savedict['comment'] = 'simulated mole fractions based on optimized state vector'
        f.add_data(savedict)

        data = self.getvalues('fromfile') 

        savedict = io.std_savedict.copy()
        savedict['name'] = "inputfilename"
        savedict['long_name'] = "name of file where original obs data was taken from"
        savedict['dtype'] = "char"
        savedict['dims'] = dimid + dim200char
        savedict['values'] = data
        savedict['missing_value'] = '!'
        f.add_data(savedict)

        f.close()

        logging.debug("Successfully wrote data to auxiliary sample output file (%s)" % auxoutputfile)

        #return outfile

	
    def add_observations_c13(self):
        """ Returns a MoleFractionList holding individual MoleFractionSample objects for all obs in a file
      
            The ObsPack mole fraction files are provided as time series per site with all dates in sequence. 
            We will loop over all site files in the ObsPackage, and subset each to our needs
            
        """

        # Step 1: Read list of available site files in package

        infile = os.path.join(self.obspackc13_dir, 'summary', '%s_dataset_summary.txt' % (self.obspackc13_id,))
        f = open(infile, 'r')
        lines = f.readlines()
        f.close()

        ncfilelist = []
        for line in lines:
            if line.startswith('#'): continue # header

            items = line.split()
            #ncfile, lab , start_date, stop_date, data_comparison = items[0:5]
            ncfile, lab , start_date, stop_date, data_comparison= line[:105].split()


            ncfilelist += [ncfile]

        logging.debug("ObsPack dataset info read, proceeding with %d netcdf files" % len(ncfilelist))

        for ncfile in ncfilelist:

            infile = os.path.join(self.obspackc13_dir, 'data', 'nc', ncfile + '.nc')
            ncf = io.ct_read(infile, 'read')
            idates = ncf.get_variable('time_components')
            dates = array([dtm.datetime(*d) for d in idates])

            subselect = logical_and(dates >= self.startdate , dates <= self.enddate).nonzero()[0]

            dates = dates.take(subselect, axis=0)
            
            obspacknum = ncf.get_variable('obspack_num').take(subselect)  # or should we propagate obs_num which is not unique across datasets??
            obspackid = ncf.get_variable('obspack_id').take(subselect, axis=0)
            obspackid = [s.tostring().lower() for s in obspackid]
            obspackid = map(strip, obspackid)
            datasetname = ncfile  # use full name of dataset to propagate for clarity
            lats = ncf.get_variable('latitude').take(subselect, axis=0)
            lons = ncf.get_variable('longitude').take(subselect, axis=0)
            alts = ncf.get_variable('altitude').take(subselect, axis=0)
            obs = ncf.get_variable('value').take(subselect, axis=0)
            species = ncf.get_attribute('dataset_parameter')
            flags = ncf.get_variable('obs_flag').take(subselect, axis=0)
            ncf.close()

            for n in range(len(dates)): 
                self.datalist.append(MoleFractionSample(obspacknum[n], dates[n], datasetname, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], obspackid[n], species, 1, 0.0, infile))
                self.datalist_c13.append(MoleFractionSample(obspacknum[n], dates[n], datasetname, obs[n], 0.0, 0.0, 0.0, 0.0, flags[n], alts[n], lats[n], lons[n], obspackid[n], species, 1, 0.0, infile))

            logging.debug("Added %d observations from file (%s) to the Data list" % (len(dates), ncfile)) 

        logging.info("C13 observations list now holds %d values" % len(self.datalist_c13))
        logging.info("Total observations list now holds %d values" % len(self.datalist))


    def add_simulations_c13(self, filename, silent=False):
        """ Adds model simulated values to the mixing ratio objects """

        import da.tools.molefraction_to_permil as mtp
        if not os.path.exists(filename):
            msg = "Sample output filename for observations could not be found : %s" % filename 
            logging.error(msg)
            logging.error("Did the sampling step succeed?")
            logging.error("...exiting")
            raise IOError, msg

        ncf = io.ct_read(filename, method='read')
        ids = ncf.get_variable('obs_num')
        simulated = ncf.get_variable('flask')
        obs_species = self.getvalues('species')

 	simulated = mtp.molefraction_to_permil(filename,simulated)
        dim_ens = len(simulated[0,:])/2
	simulated = simulated[:,dim_ens:]
        ncf.close()
        logging.info("Successfully read data from model sample file (%s)" % filename)

        obs_ids = self.getvalues('id').tolist()
        ids = map(int, ids)

        obs_ids_co2 = self.getvalues_co2('id').tolist()

        obs_ids_c13 = self.getvalues_c13('id').tolist()

        missing_samples = []

        for idx, val,obs_specie in zip(ids, simulated, obs_species): 
            logging.debug("idx %d,%d"  %(idx,obs_ids.index(idx)))
            if idx in obs_ids_c13 and obs_specie == 'co2c13':
                index = obs_ids_c13.index(idx)
                index = index + len(obs_ids_co2)
                logging.debug("add_simulations_c13 %d, %d, %s" %(index,idx,mean(array(val))))
                self.datalist[index].simulated = val
            else:     
                missing_samples.append(idx)

        if not silent and missing_samples != []:
            logging.warning('Model samples were found that did not match any ID in the observation list. Skipping them...')
            print "missing_samples", missing_samples

        logging.debug("Added %d simulated values to the Data list" % (len(ids) - len(missing_samples)))

        for i in range(len(obs_ids)):
            logging.debug("add all simulatons %d,%s" %(i,mean(array(self.datalist[i].simulated))))



    def add_model_data_mismatch_c13(self,filename):
        """ 
            Get the model-data mismatch values for this cycle.

                (1) Open a sites_weights file
                (2) Parse the data
                (3) Compare site list against data
                (4) Take care of double sites, etc

        """

        if not os.path.exists(filename):
            msg = 'Could not find  the required sites.rc input file (%s) ' % filename
            logging.error(msg)
            raise IOError, msg
        else:
            self.sites_file = filename

        sites_weights = rc.read(self.sites_file)
        
        self.rejection_threshold = int(sites_weights['obs.rejection.threshold'])
        self.global_R_scaling = float(sites_weights['global.R.scaling'])
        self.n_site_categories = int(sites_weights['n.site.categories'])

        logging.debug('Model-data mismatch rejection threshold: %d ' % self.rejection_threshold)
        logging.warning('Model-data mismatch scaling factor     : %f ' % self.global_R_scaling)
        logging.debug('Model-data mismatch site categories    : %d ' % self.n_site_categories)
   
        cats = [k for k in sites_weights.keys() if 'site.category' in k] 
        
        site_categories = {}
        for key in cats:
            name, error, may_localize, may_reject = sites_weights[key].split(';')
            name = name.strip().lower()
            error = float(error)
            may_reject = ("TRUE" in may_reject.upper())
            may_localize = ("TRUE" in may_localize.upper())
            site_categories[name] = {'category': name, 'error': error, 'may_localize': may_localize, 'may_reject': may_reject}
        
        site_info = {}
        site_move = {}
        site_hourly = {}   # option added to include only certain hours of the day (for e.g. PAL) IvdL
        site_incalt = {} # option to increase sampling altitude for sites specified in sites and weights file 
        for key, value in sites_weights.iteritems():
            if 'co2c13_' in key or 'sf6' in key:  # to be fixed later, do not yet know how to parse valid keys from rc-files yet.... WP
                sitename, sitecategory = key, value
                sitename = sitename.strip()
                sitecategory = sitecategory.split()[0].strip().lower()
                site_info[sitename] = site_categories[sitecategory]
            if 'site.move' in key:
                identifier, latmove, lonmove = value.split(';')
                site_move[identifier.strip()] = (float(latmove), float(lonmove))
            if 'site.hourly' in key:
                identifier, hourfrom, hourto = value.split(';')
                site_hourly[identifier.strip()] = (int(hourfrom), int(hourto))
            if 'site.incalt' in key:
                identifier, incalt = value.split(';')
                site_incalt[identifier.strip()] = (int(incalt))
                logging.debug('test')

        for obs in self.datalist:  # loop over all available data points
            if obs.species =='co2c13':
                obs.mdm = 1000.0  # default is very high model-data-mismatch, until explicitly set by script
                obs.flag = 99  # default is do-not-use , until explicitly set by script
                exclude_hourly = False # default is that hourly values are not included

                identifier = obs.code
                species, site, method, lab, datasetnr = identifier.split('_')

                if site_info.has_key(identifier):
                    if site_hourly.has_key(identifier):
                        hourf, hourt = site_hourly[identifier]
                        if int(obs.xdate.hour) >= hourf and int(obs.xdate.hour) <= hourt:
                            logging.warning("Observation in hourly dataset INCLUDED, while sampling time %s was between %s:00-%s:00"%(obs.xdate.time(),hourf,hourt))
                        else:
                            logging.warning("Observation in hourly dataset EXCLUDED, while sampling time %s was outside %s:00-%s:00"%(obs.xdate.time(),hourf,hourt))
                            exclude_hourly = True
                    if site_info[identifier]['category'] == 'do-not-use' or exclude_hourly:
                        logging.warning("Observation found (%s, %d), but not used in assimilation !!!" % (identifier, obs.id))
                        obs.mdm = site_info[identifier]['error']
                        obs.may_localize = site_info[identifier]['may_localize']
                        obs.may_reject = site_info[identifier]['may_reject']
                        obs.flag = 99
                    else:
                        logging.warning("Observation found (%s, %d)" % (identifier, obs.id))
                        obs.mdm = site_info[identifier]['error']
                        obs.may_localize = site_info[identifier]['may_localize']
                        obs.may_reject = site_info[identifier]['may_reject']
                        obs.flag = 0

                else:
                    logging.warning("Observation NOT found (%s, %d), please check sites.rc file (%s)  !!!" % (identifier, obs.id, self.sites_file))

                if site_move.has_key(identifier):

                    movelat, movelon = site_move[identifier]
                    obs.lat = obs.lat + movelat
                    obs.lon = obs.lon + movelon

                    logging.warning("Observation location for (%s, %d), is moved by %3.2f degrees latitude and %3.2f degrees longitude" % (identifier, obs.id, movelat, movelon))

        # Add SiteInfo dictionary to the Observation object for future use

        self.site_info = site_info
        self.site_move = site_move
        self.site_hourly = site_hourly
        self.site_incalt = site_incalt

        logging.debug("Added Model Data Mismatch to all samples ")


################### End Class ObsPackObservations ###################



################### Begin Class MoleFractionSample ###################

class MoleFractionSample(object):
    """ 
        Holds the data that defines a mole fraction Sample in the data assimilation framework. Sor far, this includes all
        attributes listed below in the __init__ method. One can additionally make more types of data, or make new
        objects for specific projects.

    """

    def __init__(self, idx, xdate, code='XXX', obs=0.0, simulated=0.0, resid=0.0, hphr=0.0, mdm=0.0, flag=0, height=0.0, lat= -999., lon= -999., evn='0000', species='co2', samplingstrategy=1, sdev=0.0, fromfile='none.nc'):
        self.code = code.strip()      # dataset identifier, i.e., co2_lef_tower_insitu_1_99
        self.xdate = xdate             # Date of obs
        self.obs = obs               # Value observed
        self.simulated = simulated         # Value simulated by model
        self.resid = resid             # Mole fraction residuals
        self.hphr = hphr              # Mole fraction prior uncertainty from fluxes and (HPH) and model data mismatch (R)
        self.mdm = mdm               # Model data mismatch
        self.may_localize = True           # Whether sample may be localized in optimizer
        self.may_reject = True              # Whether sample may be rejected if outside threshold
        self.flag = flag              # Flag
        self.height = height            # Sample height in masl
        self.lat = lat               # Sample lat
        self.lon = lon               # Sample lon
        self.id = idx               # Obspack ID within distrution (integer), e.g., 82536
        self.evn = evn               # Obspack Number within distrution (string), e.g., obspack_co2_1_PROTOTYPE_v0.9.2_2012-07-26_99_82536
        self.sdev = sdev              # standard deviation of ensemble
        self.masl = True              # Sample is in Meters Above Sea Level
        self.mag = not self.masl     # Sample is in Meters Above Ground
        self.species = species.strip()
        self.samplingstrategy = samplingstrategy
        self.fromfile = fromfile   # netcdf filename inside ObsPack distribution, to write back later

################### End Class MoleFractionSample ###################


if __name__ == "__main__":

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

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

    logging.root.setLevel(logging.DEBUG)

    DaCycle = CycleControl(args={'rc':'../../ctdas-ei-gfed2-bcb-6x4-gridded-test.rc'})
    #DaCycle.Initialize()
    DaCycle.ParseTimes()

    DaSystem = CtDaSystem('../rc/carbontracker_geocarbon_gridded.rc')
    #DaSystem.Initialize()

    DaCycle.DaSystem = DaSystem

    obs = ObsPackObservations()
    obs.DaCycle = DaCycle

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

        DaCycle['time.sample.start'] = DaCycle['time.start']
        DaCycle['time.sample.end'] = DaCycle['time.end']

        obs.Initialize()
        #obs.Validate()
        obs.add_observations()
        obs.add_model_data_mismatch()

        print(obs.getvalues('obs'))
        print(obs.getvalues('mdm'))

        DaCycle.AdvanceCycleTimes()



