refactored/Makefile

+.DEFAULT_GOAL: clean
+
+.PHONY: help clean
+
+help:
+	@echo " "
+	@echo " Usage:"
+	@echo "   make clean          # remove temporary (log) files"
+	@echo " "
+
+clean: 
+
+# -------------------
+	/bin/rm -f jb.[0-9]*.{jb,rc,log}
+	/bin/rm -f */*/jb.[0-9]*.{rc,jb,log}
+	/bin/rm -f *.pyc
+	/bin/rm -f */*/*.pyc
+	/bin/rm -f das.[0-9]*.log
+	/bin/rm -f *.[eo]*
+	/bin/rm -f *.p[eo]*
+
+
+

refactored/README.push-to-svn.txt

+How to push to maunaloa svn
+===========================
+
+You can push remotec changes to the maunaloa SVN server with the
+following steps (the first 3 are only needed the first time):
+
+Create your branch
+------------------
+
+* Add a specialized branch. The following adds a branch named remotec/refactored
+  $ svn copy https://USER@maunaloa.wur.nl/subversion/das/ct/trunk \
+             https://USER@maunaloa.wur.nl/subversion/das/ct/branches/remotec/refactored \
+      -m "Creating a branch for upstreaming from Mercurial to Subversion"
+
+Safety checks
+-------------
+
+* Make sure you have Mercurial and hgsubversion installed. If not, get them via 
+  $ easy_install --prefix ~ -U mercurial # might need adjustments to PYTHONPATH and PATH 
+  $ hg clone http://bitbucket.org/durin42/hgsubversion/ ~/hgsubversion
+
+* Make sure hgsubversion is active. To test: hg help svn. If not:
+  $ echo "[extensions]" >> ~/.hgrc
+  $ echo "hgsubversion=~/hgsubversion/hgsubversion" >> ~/.hgrc
+
+* Make sure you are in the bridge-repo. 
+  $ pwd # the path should end in ct-svn-bridge
+  Note: Never develop in the bridge-repo! (to avoid confusion)
+
+* Check that the push path is correct:
+  $ hg paths # should show default = svn+https://USER@maunaloa.wur.nl/subversion/das/ct
+
+* Make sure that svn can access the maunaloa repo
+  $ svn ls https://maunaloa.wur.nl/subversion/das/ct/branches/remotec/refactored
+  You might have to accept the host as safe or to provide the correct authentication.
+
+* Ensure that you are on the remotec/refactored branch:
+  $ hg up remotec/refactored
+
+Actual transfer
+---------------
+
+* Check your local repo for new changes
+  $ hg in ../CT # shows <rev> IDs you can use later
+
+* Pull the changes you want to get into svn from a local repo
+  $ hg pull ../CT
+
+* Graft all changes you want to contribute:
+  $ hg graft [rev] [rev] # rev could be 645 or b7a92bbb0e0b. Best use the second>. 
+  You need to specify every rev individually, 
+  but you can graft multiple revs in one command.
+
+* Check what you would push:
+  $ hg outgoing
+
+* Push the changes:
+  $ hg push
+  This *might* show some unrelated pulled revisions 
+  and *should* show your new revisions as pulled, 
+  along with paths to backup bundles (which you should not need).
+
+Done.
+
+Note: hgsubversion pushes your changes to the subversion server, then pulls them *as new revisions* and finally removes the local versions. (reason: subversion can change some metadata - like timeflags - so that it is not possible to transfer the changes exactly).
+
+Note: You should **NEVER EVER MERGE INTO REMOTEC/REFACTORED**. Hgsubversion would not be able to represent the merge correctly (because subversion does not know real merges), so it bails out when you try to push it. That could mean having to get out the mq knife and strip the merges from the history. Just use graft. Is is as good as merge, but additionally it can reorder changesets.
+
+More advanced refactoring
+-------------------------
+
+If you want to refactor your history to make it nicer to read and review for subversion users, you can use a middle branch. I use `remotec` for that.
+
+* Switch to the middle branch:  
+  $ hg up remotec
+
+* Graft all changesets you want to include, just like before.
+
+* Switch to the refactored branch:  
+  $ hg up remotec/refactored
+
+* Gather all revisions you want to put into one single refactored commit (check them via `hg log`). To do that, call the following for each revision:
+  $ hg export REV | hg import --force --no-commit -
+
+* Do that single commit:  
+  $ hg commit -m "Comment about the refactored commit"  
+  Note: To only commit parts of the changes, you can use the record extension. Info on extensions:  
+  $ hg help extensions
+  With it you can select hunks of changes one by one:
+  $ hg record -m "Comment about the partial commit you want to do"
+
+* Check outgoing and push, just like you would in the simple case.
+
+To make the gathering of multiple commits easier for me, I added the following aliases in my ~/.hgrc:
+
+    [alias]
+    getrevs = !$HG export $@ | $HG import --force --no-commit - ; for i in $@; do hg log --template "{node|short}: {desc}\n" -r $i >> .hg/upcoming-commitmessage.txt; done
+    finrevs = !echo $@ > .hg/commitmessage-tmp.txt ; if test "#@" != "" ; then echo >> .hg/commitmessage-tmp.txt; fi; cat .hg/upcoming-commitmessage.txt >> .hg/commitmessage-tmp.txt ; $HG ci -l .hg/commitmessage-tmp.txt ; mv .hg/upcoming-commitmessage.txt .hg/upcoming-commitmessage.txt.bak-$(date +%s) ; rm .hg/commitmessage-tmp.txt
+    checkrevs = !cat .hg/upcoming-commitmessage.txt
+
+With those I gather multiple commits into one refactored commit with the following commands:
+
+    hg getrevs REV REV REV … # gather all revisions
+    hg checkrevs # Check the included revs and their part of the commit message
+    hg finrevs whatever I want to say about the commit.

refactored/carbontracker.rc

+!!! Info for the CarbonTracker data assimilation system
+
+! datadir         : /home/gosat/software/REMOTEC-MOD/CT/data/ct09
+datadir         : /home/gosat/software/REMOTEC-MOD/CT/data/ctdas_2012
+! obs.input.dir   : /home/gosat/software/REMOTEC-MOD/CT/data/ct08/obsnc/with_fillvalue
+obs.input.dir   : /home/gosat/software/REMOTEC-MOD/CT/data/ctdas_2012/obsnc/with_fillvalue
+obs.input.fname : obs_forecast_remotec.nc
+ocn.covariance  : ${datadir}/covariances/ocean_oif/oif_p3_era40.dpco2.2000.01.hdf 
+bio.covariance  : ${datadir}/covariances/olson/covariance_bio_olson19.nc
+deltaco2.prefix : oif_p3_era40.dpco2
+regtype         : olson19_oif30
+nparameters     : 240
+random.seed     : 4385
+regions.dir     : /home/gosat/software/REMOTEC-MOD/CT/data/sharedaux
+regionsfile     : ${regions.dir}/regions.nc
+
+! Include a naming scheme for the variables
+
+#include da/rc/NamingScheme.wp_Mar2011.rc 
+
+! Info on the sites file used
+
+obs.sites.rc        : ${datadir}/sites_and_weights_co2.ct10.rc
+

refactored/da.rc

+! Info on the data assimilation cycle
+
+time.restart        : False
+time.start          : 2009-01-02 00:00:00
+time.finish         : 2009-01-31 00:00:00
+! length of one unit in days (typically 7)
+time.cycle          : 7
+! number of units in state vector (5 in the example above, just like in carbontracker)
+time.nlag           : 5
+! cycle 1, nlag 3 means: 3 days in the state vector.
+! cycle 7 nlag 5 means: 7*5 days in the state vector, advance by 7 days per step.
+dir.da_run          : /home/gosat/software/REMOTEC-MOD/work/CTDAS
+
+
+! Info on the DA system used
+
+da.system           : CarbonTracker
+da.system.rc        : /home/gosat/software/REMOTEC-MOD/CT/carbontracker.rc
+
+! Info on the forward model to be used
+
+da.obsoperator         : TM5
+da.obsoperator.rc      : /home/gosat/software/REMOTEC-MOD/TM5/tm5-ctdas.rc
+da.optimizer.nmembers  : 20
+!da.optimizer.nmembers  : 150

refactored/da/analysis/expand_fluxes.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@noaa.gov)
+
+Revision History:
+File created on 21 Ocotber 2008.
+
+"""
+
+def proceed_dialog(txt, yes=['y','yes'], all=['a','all', 'yes-to-all']):
+    """ function to ask whether to proceed or not """
+    response=raw_input(txt)
+    if response.lower() in yes:
+       return 1
+    if response.lower() in all:
+       return 2
+    return 0
+
+def _getpriordata(DaCycle, StateVector, startdate, enddate, nlag, cyclelength):
+    """Read gridmean and gridvariance from the oldest available nlag.
+
+    :returns: gridmean, gridvariance, used nlag."""
+    import logging
+    tried = []
+    if enddate - timedelta(cyclelength) != startdate:
+        logging.warn("the last cycle is shorter than the cycle length. Take the results with a grain of salt.")
+    for n in range(nlag,0,-1):
+        # TODO: This breaks, when enddate - startdate is not a multiple of the cyclelength.
+        # TODO: enddate -> startdate, timedelta(cyclelength*n) -> timedelta(cyclelength*(n-1))
+        priordate       = startdate - timedelta(cyclelength*(n-1))
+        savedir         = DaCycle['dir.output'].replace(startdate.strftime('%Y%m%d'),priordate.strftime('%Y%m%d') )
+        filename        = os.path.join(savedir,'savestate.nc')
+        tried.append(filename)
+        if os.path.exists(filename):
+            dummy                   = StateVector.ReadFromFile(filename)
+            gridmean,gridvariance   = StateVector.StateToGrid(lag=n)
+            
+            msg = 'Read prior dataset from file %s, sds %d: '%(filename,n) ; logging.debug(msg)
+            return gridmean, gridvariance, n
+    raise ValueError("No savestate.nc found in the possible paths:" + str(tried) + ". Likely the timeframe is no multiple of the cycle length: " + str(cyclelength) + ". startdate: " + str(startdate) + " enddate: " + str(enddate))
+
+
+def SaveWeeklyAvg1x1Data(DaCycle, StateVector):
+    """
+        Function creates a NetCDF file with output on 1x1 degree grid. It uses the flux data written by the 
+        :class:`~da.baseclasses.obsoperator.ObsOperator.py`, and multiplies these with the mapped parameters and
+        variance (not covariance!) from the :class:`~da.baseclasses.statevector.StateVector`.
+        
+           :param DaCycle: a :class:`~da.tools.initexit.CycleControl` object
+           :param StateVector: a :class:`~da.baseclasses.statevector.StateVector`
+           :rtype: None
+    """
+
+    import da.tools.io4 as io
+    import logging
+#
+    dirname     = 'data_flux1x1_weekly'
+#
+    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'] 
+    nlag        = StateVector.nlag
+
+    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)
+
+#
+# Create or open NetCDF output file
+#
+    saveas      = os.path.join(dirname,'flux_1x1.nc')
+    ncf         = io.CT_CDF(saveas,'write')
+
+#
+# Create dimensions and lat/lon grid
+#
+    dimgrid     = ncf.AddLatLonDim()
+    dimdate     = ncf.AddDateDim()
+#
+# set title and tell GMT that we are using "pixel registration"
+#
+    setattr(ncf,'Title','CarbonTracker fluxes')
+    setattr(ncf,'node_offset',1)
+#
+# skip dataset if already in file
+#
+    ncfdate = date2num(startdate)-dectime0+dt.days/2.0
+    skip    = ncf.has_date(ncfdate)
+    if skip:
+        msg = 'Skipping writing of data for date %s : already present in file %s' % (startdate.strftime('%Y-%m-%d'),saveas) ; logging.warning(msg)
+    else:
+        
+#
+# if not, process this cycle. Start by getting flux input data from CTDAS
+#
+        filename = os.path.join(DaCycle['dir.output'],'flux1x1_%s_%s.nc'%(startdate.strftime('%Y%m%d%H'),enddate.strftime('%Y%m%d%H'),) )
+
+        file                = io.CT_Read(filename,'read')
+        bio                 = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.bio.flux']))
+        ocean               = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.ocean.flux']))
+        fire                = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.fires.flux']))
+        fossil              = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.fossil.flux']))
+        #mapped_parameters   = np.array(file.GetVariable(DaCycle.DaSystem['final.param.mean.1x1']))
+        file.close()
+
+        next=ncf.inq_unlimlen()[0]
+
+
+# Start adding datasets from here on, both prior and posterior datasets for bio and ocn
+
+        for prior in [True, False]:
+#
+# Now fill the StateVector with the prior values for this time step. Note that the prior value for this time step
+# occurred nlag time steps ago, so we make a shift in the output directory, but only if we are more than nlag cycle away from the start date..
+#
+
+            if prior:
+                qual_short='prior'
+                gridmean, gridvariance, choicelag = _getpriordata(DaCycle, StateVector, startdate, enddate, nlag, cyclelength=dt.days)
+            else:
+                qual_short='opt'
+                savedir                 = DaCycle['dir.output']
+                filename                = os.path.join(savedir,'savestate.nc')
+                dummy                   = StateVector.ReadFromFile(filename)
+                gridmean,gridvariance   = StateVector.StateToGrid(lag=nlag)
+
+                msg = 'Read posterior dataset from file %s, sds %d: '%(filename,nlag) ; logging.debug(msg)
+#
+# if prior, do not multiply fluxes with parameters, otherwise do
+#
+            biomapped=bio*gridmean 
+            oceanmapped=ocean*gridmean 
+            biovarmapped=bio*gridvariance
+            oceanvarmapped=ocean*gridvariance
+
+#
+#
+#  For each dataset, get the standard definitions from the module mysettings, add values, dimensions, and unlimited count, then write
+#
+            savedict=ncf.StandardVar(varname='bio_flux_'+qual_short)
+            savedict['values']=biomapped.tolist()
+            savedict['dims']=dimdate+dimgrid
+            savedict['count']=next
+            ncf.AddData(savedict)
+#
+            savedict=ncf.StandardVar(varname='ocn_flux_'+qual_short)
+            savedict['values']=oceanmapped.tolist()
+            savedict['dims']=dimdate+dimgrid
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+            savedict=ncf.StandardVar(varname='bio_flux_%s_cov'%qual_short)
+            savedict['values']=biovarmapped.tolist()
+            savedict['dims']=dimdate+dimgrid
+            savedict['count']=next
+            ncf.AddData(savedict)
+#
+            savedict=ncf.StandardVar(varname='ocn_flux_%s_cov'%qual_short)
+            savedict['values']=oceanvarmapped.tolist()
+            savedict['dims']=dimdate+dimgrid
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+        # End prior/posterior block
+
+        savedict=ncf.StandardVar(varname='fire_flux_imp')
+        savedict['values']=fire.tolist()
+        savedict['dims']=dimdate+dimgrid
+        savedict['count']=next
+        ncf.AddData(savedict)
+#
+        savedict=ncf.StandardVar(varname='fossil_flux_imp')
+        savedict['values']=fossil.tolist()
+        savedict['dims']=dimdate+dimgrid
+        savedict['count']=next
+        ncf.AddData(savedict)
+#
+        savedict=ncf.StandardVar(varname='date')
+        savedict['values']=date2num(startdate)-dectime0+dt.days/2.0
+        savedict['dims']=dimdate
+        savedict['count']=next
+        ncf.AddData(savedict)
+
+        sys.stdout.write('.')
+        sys.stdout.flush()
+#
+#   Done, close the new NetCDF file
+#
+    ncf.close()
+#
+#   Return the full name of the NetCDF file so it can be processed by the next routine
+#
+    msg="Gridded weekly average fluxes now written" ; logging.info(msg)
+
+    return saveas
+
+def SaveWeeklyAvgStateData(DaCycle, StateVector):
+    """
+        Function creates a NetCDF file with output for all parameters. It uses the flux data written by the 
+        :class:`~da.baseclasses.obsoperator.ObsOperator.py`, and multiplies these with the mapped parameters and
+        variance (not covariance!) from the :class:`~da.baseclasses.statevector.StateVector`.
+        
+           :param DaCycle: a :class:`~da.tools.initexit.CycleControl` object
+           :param StateVector: a :class:`~da.baseclasses.statevector.StateVector`
+           :rtype: None
+    """
+
+    import da.tools.io4 as io
+    import logging
+    from da.analysis.tools_regions import globarea
+    
+#
+    dirname     = 'data_state_weekly'
+#
+    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'] 
+    nlag        = StateVector.nlag
+
+    area        = globarea()
+    vectorarea  = StateVector.GridToVector(griddata=area,method='sum')
+
+    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)
+
+#
+# Create or open NetCDF output file
+#
+    saveas      = os.path.join(dirname,'statefluxes.nc')
+    ncf         = io.CT_CDF(saveas,'write')
+
+#
+# Create dimensions and lat/lon grid
+#
+    dimregs     = ncf.AddDim('nparameters',StateVector.nparams)
+    dimmembers  = ncf.AddDim('nmembers',StateVector.nmembers)
+    dimdate     = ncf.AddDateDim()
+#
+# set title and tell GMT that we are using "pixel registration"
+#
+    setattr(ncf,'Title','CarbonTracker fluxes')
+    setattr(ncf,'node_offset',1)
+#
+# skip dataset if already in file
+#
+    ncfdate = date2num(startdate)-dectime0+dt.days/2.0
+    skip    = ncf.has_date(ncfdate)
+    if skip:
+        msg = 'Skipping writing of data for date %s : already present in file %s' % (startdate.strftime('%Y-%m-%d'),saveas) ; logging.warning(msg)
+    else:
+
+        next=ncf.inq_unlimlen()[0]
+
+#
+# if not, process this cycle. Start by getting flux input data from CTDAS
+#
+        filename = os.path.join(DaCycle['dir.output'],'flux1x1_%s_%s.nc'%(startdate.strftime('%Y%m%d%H'),enddate.strftime('%Y%m%d%H'),) )
+
+        file                = io.CT_Read(filename,'read')
+        bio                 = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.bio.flux']))
+        ocean               = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.ocean.flux']))
+        fire                = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.fires.flux']))
+        fossil              = np.array(file.GetVariable(DaCycle.DaSystem['background.co2.fossil.flux']))
+        #mapped_parameters   = np.array(file.GetVariable(DaCycle.DaSystem['final.param.mean.1x1']))
+        file.close()
+
+        next=ncf.inq_unlimlen()[0]
+
+        vectorbio           = StateVector.GridToVector(griddata=bio*area,method='sum')
+        vectorocn           = StateVector.GridToVector(griddata=ocean*area,method='sum')
+        vectorfire          = StateVector.GridToVector(griddata=fire*area,method='sum')
+        vectorfossil        = StateVector.GridToVector(griddata=fossil*area,method='sum')
+
+
+# Start adding datasets from here on, both prior and posterior datasets for bio and ocn
+
+        for prior in [True, False]:
+#
+# Now fill the StateVector with the prior values for this time step. Note that the prior value for this time step
+# occurred nlag time steps ago, so we make a shift in the output directory, but only if we are more than nlag cycle away from the start date..
+#
+
+            if prior:
+                qual_short='prior'
+                for n in range(nlag,0,-1):
+                    priordate       = enddate - timedelta(dt.days*n)
+                    savedir         = DaCycle['dir.output'].replace(startdate.strftime('%Y%m%d'),priordate.strftime('%Y%m%d') )
+                    filename        = os.path.join(savedir,'savestate.nc')
+                    if os.path.exists(filename):
+                        dummy                   = StateVector.ReadFromFile(filename)
+                        choicelag               = n
+
+                        msg = 'Read prior dataset from file %s, sds %d: '%(filename,n) ; logging.debug(msg)
+                        break
+            else:
+                qual_short='opt'
+                savedir                 = DaCycle['dir.output']
+                filename                = os.path.join(savedir,'savestate.nc')
+                dummy                   = StateVector.ReadFromFile(filename)
+                choicelag               = nlag
+
+                msg = 'Read posterior dataset from file %s, sds %d: '%(filename,nlag) ; logging.debug(msg)
+#
+# if prior, do not multiply fluxes with parameters, otherwise do
+#
+            data                = StateVector.EnsembleMembers[choicelag-1][0].ParameterValues*vectorbio # units of mole region-1 s-1
+
+            savedict            = ncf.StandardVar(varname='bio_flux_%s'%qual_short)
+            savedict['values']  = data
+            savedict['dims']    = dimdate+dimregs
+            savedict['count']   = next
+            ncf.AddData(savedict)
+
+            members             = StateVector.EnsembleMembers[choicelag-1] 
+            deviations          = np.array([mem.ParameterValues*data for mem in members])
+
+            savedict=ncf.StandardVar(varname='bio_flux_%s_cov'%qual_short)
+
+            savedict['values']=deviations.tolist()
+            savedict['dims']=dimdate+dimmembers+dimregs
+            savedict['comment']="This is the matrix square root, use (M x M^T)/(nmembers-1) to make covariance"
+            savedict['units']="mol region-1 s-1"
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+            savedict=ncf.StandardVar('unknown')
+            savedict['name'] = 'bio_flux_%s_std'%qual_short
+            savedict['long_name'] = 'Biosphere flux standard deviation, %s'%qual_short
+            savedict['values']=deviations.std(axis=0)
+            savedict['dims']=dimdate+dimregs
+            savedict['comment']="This is the standard deviation on each parameter"
+            savedict['units']="mol region-1 s-1"
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+            data                = StateVector.EnsembleMembers[choicelag-1][0].ParameterValues*vectorocn # units of mole region-1 s-1
+
+            savedict=ncf.StandardVar(varname='ocn_flux_%s'%qual_short)
+            savedict['values']=data
+            savedict['dims']=dimdate+dimregs
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+            deviations          = np.array([mem.ParameterValues*data for mem in members])
+
+            savedict=ncf.StandardVar(varname='ocn_flux_%s_cov'%qual_short)
+            savedict['values']=deviations.tolist()
+            savedict['dims']=dimdate+dimmembers+dimregs
+            savedict['comment']="This is the matrix square root, use (M x M^T)/(nmembers-1) to make covariance"
+            savedict['units']="mol region-1 s-1"
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+            savedict=ncf.StandardVar('unknown')
+            savedict['name'] = 'ocn_flux_%s_std'%qual_short
+            savedict['long_name'] = 'Ocean flux standard deviation, %s'%qual_short
+            savedict['values']=deviations.std(axis=0)
+            savedict['dims']=dimdate+dimregs
+            savedict['comment']="This is the standard deviation on each parameter"
+            savedict['units']="mol region-1 s-1"
+            savedict['count']=next
+            ncf.AddData(savedict)
+
+        data                = vectorfire
+
+        savedict=ncf.StandardVar(varname='fire_flux_imp')
+        savedict['values']=data
+        savedict['dims']=dimdate+dimregs
+        savedict['count']=next
+        ncf.AddData(savedict)
+
+        data                = vectorfossil
+
+        savedict=ncf.StandardVar(varname='fossil_flux_imp')
+        savedict['values']=data
+        savedict['dims']=dimdate+dimregs
+        savedict['count']=next
+        ncf.AddData(savedict)
+
+        savedict=ncf.StandardVar(varname='date')
+        savedict['values']  = ncfdate
+        savedict['dims']    = dimdate
+        savedict['count']   = next
+        ncf.AddData(savedict)
+
+        sys.stdout.write('.')
+        sys.stdout.flush()
+#
+#   Done, close the new NetCDF file
+#
+    ncf.close()
+#
+#   Return the full name of the NetCDF file so it can be processed by the next routine
+#
+    msg="Vector weekly average fluxes now written" ; logging.info(msg)
+
+    return saveas
+
+
+def SaveWeeklyAvgTCData(DaCycle, StateVector):
+    """
+        Function creates a NetCDF file with output on TransCom regions. It uses the flux input from the 
+        function `SaveWeeklyAvg1x1Data` to create fluxes of length `nparameters`, which are then projected
+        onto TC regions using the internal methods from :class:`~da.baseclasses.statevector.StateVector`.
+        
+           :param DaCycle: a :class:`~da.tools.initexit.CycleControl` object
+           :param StateVector: a :class:`~da.baseclasses.statevector.StateVector`
+           :rtype: None
+
+        This function only read the prior fluxes from the flux_1x1.nc files created before, because we want to convolve 
+        these with the parameters in the StateVector. This creates posterior fluxes, and the posterior covariance for the complete
+        StateVector in units of mol/box/s which we then turn into TC fluxes and covariances.
+    """
+    from da.analysis.tools_regions import globarea, StateToGrid
+    from da.analysis.tools_transcom import StateToTranscom, StateCovToTranscom, transcommask
+    import da.tools.io4 as io
+    import logging
+#
+    dirname     = 'data_tc_weekly'
+#
+    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'] 
+    ncfdate     = date2num(startdate)-dectime0+dt.days/2.0
+
+    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)
+
+    # Write/Create NetCDF output file
+    #
+    saveas          = os.path.join(dirname,'tcfluxes.nc')
+    ncf             = io.CT_CDF(saveas,'write')
+    dimdate         = ncf.AddDateDim()
+    dimidateformat  = ncf.AddDateDimFormat()
+    dimregs         = ncf.AddRegionDim(type='tc')
+#
+# set title and tell GMT that we are using "pixel registration"
+#
+    setattr(ncf,'Title','CarbonTracker TransCom fluxes')
+    setattr(ncf,'node_offset',1)
+    #
+
+    skip             = ncf.has_date(ncfdate)
+    if skip:
+        msg = 'Skipping writing of data for date %s : already present in file %s' % (startdate.strftime('%Y-%m-%d'),saveas) ; logging.warning(msg)
+    else:
+
+        # Get input data
+
+        area=globarea()
+
+        infile=os.path.join(DaCycle['dir.analysis'],'data_state_weekly','statefluxes.nc')
+        if not os.path.exists(infile):
+            msg="Needed input file (%s) does not exist yet, please create file first, returning..."%infile ; logging.error(msg)
+            return None
+
+        ncf_in      = io.CT_Read(infile,'read')
+
+        # Transform data one by one
+
+        # Get the date variable, and find index corresponding to the DaCycle date
+
+        try:
+            dates       = ncf_in.variables['date'][:]
+        except KeyError:
+            msg         = "The variable date cannot be found in the requested input file (%s) " % infile ; logging.error(msg)
+            msg         = "Please make sure you create gridded fluxes before making TC fluxes " ; logging.error(msg)
+            raise KeyError
+
+        try:
+            index       = dates.tolist().index(ncfdate)
+        except ValueError:
+            msg         = "The requested cycle date is not yet available in file %s "% infile ; logging.error(msg)
+            msg         = "Please make sure you create state based fluxes before making TC fluxes " ; logging.error(msg)
+            raise ValueError
+
+        # First add the date for this cycle to the file, this grows the unlimited dimension
+
+        savedict            = ncf.StandardVar(varname='date')
+        savedict['values']  = ncfdate
+        savedict['dims']    = dimdate
+        savedict['count']   = index
+        dummy               = ncf.AddData(savedict)
+
+        # Now convert other variables that were inside the flux_1x1 file
+
+        vardict     = ncf_in.variables
+        for vname, vprop in vardict.iteritems():
+
+            data    = ncf_in.GetVariable(vname)[index]
+
+            savedict            = ncf.StandardVar(varname=vname)
+
+            if vname=='latitude': continue
+            elif vname=='longitude': continue
+            elif vname=='date': continue
+            elif vname=='idate': continue
+            elif 'cov' in vname:
+                
+                #AB? Only works with recent version of numpy, not with 1.4.1 installed on the suns
+                #cov                 = data.transpose().dot(data) 
+                cov                 = np.dot(data.transpose(),data)
+                tcdata              = StateVector.VectorToTC(vectordata=cov,cov=True) # vector to TC
+
+                savedict['units']   = '[mol/region/s]**2'
+                savedict['dims']    = dimdate+dimregs+dimregs
+                
+            else:
+
+                tcdata              = StateVector.VectorToTC(vectordata=data) # vector to TC
+
+                savedict['dims']    = dimdate+dimregs
+                savedict['units']   = 'mol/region/s'
+
+            savedict['values']  = tcdata
+            ncf.AddData(savedict)
+
+        ncf_in.close()
+    ncf.close()
+
+    msg="TransCom weekly average fluxes now written" ; logging.info(msg)
+
+    return saveas
+
+def SaveTCDataExt(rundat):
+    """ Function SaveTCDataExt saves surface flux data to NetCDF files for extended TransCom regions
+        
+        *** Inputs ***
+        rundat : a RunInfo object
+
+        *** Outputs ***
+        NetCDF file containing n-hourly global surface fluxes per TransCom region
+
+        *** Example ***
+        ./expand_savestate project=enkf_release sd=20000101 ed=20010101 """
+    from da.analysis.tools_transcom import StateCovToGrid, transcommask, ExtendedTCRegions
+    import da.tools.io4 as io
+
+    infile=os.path.join(rundat.outputdir,'data_tc_weekly','tcfluxes.nc')
+    if not os.path.exists(infile):
+        print "Needed input file (%s) does not exist yet, please create weekly tc flux files first, returning..."%infile
+        return None
+
+    # Create NetCDF output file
+    #
+    saveas          = os.path.join(rundat.outputdir,'data_tc_weekly','tc_extfluxes.nc')
+    ncf             = io.CT_CDF(saveas,'create')
+    dimdate         = ncf.AddDateDim()
+    dimidateformat  = ncf.AddDateDimFormat()
+    dimregs         = ncf.AddRegionDim(type='tc_ext')
+
+
+    ncf_in      = io.CT_Read(infile,'read')
+    vardict     = ncf_in.variables
+    for vname, vprop in vardict.iteritems():
+
+        data    = ncf_in.GetVariable(vname)[:]
+
+        if vname=='latitude': continue
+        elif vname=='longitude': continue
+        elif vname=='date': dims=dimdate
+        elif vname=='idate': dims=dimdate+dimidateformat
+
+        if vname not in ['date','idate']:
+
+                if 'cov' in vname: 
+                    alldata=[] 
+                    for dd in data:
+                        dd=StateCovToGrid(dd*area,transcommask,reverse=True)
+                        alldata.append(dd)
+                    data=array(alldata) 
+                    dims=dimdate+dimregs+dimregs
+                else:
+                    data=ExtendedTCRegions(data)
+                    dims=dimdate+dimregs
+
+        print vname,data.shape
+        savedict            = ncf.StandardVar(varname=vname)
+        savedict['values']  = data.tolist()
+        savedict['dims']    = dims
+        savedict['units']   = 'mol/region/s'
+        savedict['count']   = 0
+        ncf.AddData(savedict,nsets=data.shape[0])
+         
+    ncf.close()
+
+    return saveas
+
+def SaveEcoDataExt(rundat):
+    """ Function SaveEcoDataExt saves surface flux data to NetCDF files for extended ecoregions
+        
+        *** Inputs ***
+        rundat : a RunInfo object
+
+        *** Outputs ***
+        NetCDF file containing n-hourly global surface fluxes per TransCom region
+
+        *** Example ***
+        ./expand_savestate project=enkf_release sd=20000101 ed=20010101 """
+    from ecotools import xte
+    from da.analysis.tools_transcom import LookUpName
+    import pycdf as CDF
+
+    infile=os.path.join(rundat.outputdir,'data_eco_weekly','ecofluxes.nc')
+    if not os.path.exists(infile):
+        print "Needed input file (%s) does not exist yet, please create weekly ecoflux files first, returning..."%infile
+        return None
+
+    # Create NetCDF output file
+    #
+    saveas=os.path.join(rundat.outputdir,'data_eco_weekly','eco_extfluxes.nc')
+    ncf=CT_CDF(saveas,'create')
+    dimdate=ncf.AddDateDim()
+    dimidateformat=ncf.AddDateDimFormat()
+    dimregs         = ncf.AddRegionDim(type='tc_ext')
+
+
+    ncf_in=CDF.CDF(infile)
+    vardict=ncf_in.variables()
+    for vname, vprop in vardict.iteritems():
+        dims=()
+        data=array(ncf_in.var(vname)[:])
+        atts=ncf_in.var(vname).attributes()
+        if vname not in ['date','idate']:
+                if 'cov' in vname: 
+                    data=xte(data[:,0:rundat.n_land,0:rundat.n_land],cov=True)
+                    dims=dimdate+dimregs+dimregs
+                else: 
+                    data=xte(data[:,0:rundat.n_land])
+                    dims=dimdate+dimregs
+        elif vname=='date': dims=dimdate
+        elif vname=='idate': dims=dimdate+dimidateformat
+        else:
+             print 'Dataset with unknown dimensions encountered in file: %s'%vname
+
+        savedict=ncf.StandardVar(varname=vname)
+        savedict['units']=atts['units']
+        savedict['values']=data.tolist()
+        savedict['dims']=dims
+        ncf.AddData(savedict,nsets=data.shape[0])
+         
+    ncf.close()
+
+    return saveas
+
+def SaveTimeAvgData(rundat,infile,avg='monthly'):
+    """ Function saves time mean surface flux data to NetCDF files
+        
+        *** Inputs ***
+        rundat : a RunInfo object
+
+        *** Outputs ***
+        daily NetCDF file containing 1-hourly global surface fluxes at 1x1 degree
+
+        *** Example ***
+        ./expand_savestate project=enkf_release sd=20000101 ed=20010101 """
+
+    import da.analysis.tools_time as timetools
+    import da.tools.io4 as io
+
+    dirname,filename    = os.path.split(infile)
+    dirname             = CreateDirs(os.path.join(rundat.outputdir,dirname.replace('weekly',avg) ) )
+
+    dectime0            = date2num(datetime(2000,1,1))
+
+# Create NetCDF output file
+#
+    saveas              = os.path.join(dirname,filename)
+    ncf                 = io.CT_CDF(saveas,'create')
+    dimdate             = ncf.AddDateDim()
+#
+# Open input file specified from the command line
+#
+    if not os.path.exists(infile):
+        print "Needed input file (%s) not found. Please create this first:"%infile
+        print "returning..."
+        return None 
+    else:
+        pass
+
+    file        = io.CT_Read(infile,'read')
+    datasets    = file.variables.keys()
+    date        = file.GetVariable('date')
+    globatts    = file.ncattrs()
+
+    time        = [datetime(2000,1,1)+timedelta(days=d) for d in date]
+
+#
+# Add global attributes, sorted 
+#
+
+    keys        = globatts
+    keys.sort()
+    for k in keys:
+        setattr(ncf,k,k)
+
+# loop over datasets in infile, skip idate and date as we will make new time axis for the averaged data
+
+    for sds in datasets:
+        if sds in ['idate','date'] : continue
+        print '['
+
+# get original data
+
+        data        = array(file.GetVariable(sds))
+        varatts     = file.variables[sds].ncattrs()
+        vardims     = file.variables[sds].dimensions
+#
+# Depending on dims of input dataset, create dims for output dataset. Note that we add the new dimdate now.
+#
+
+        dims=()
+        for d in vardims:
+            if 'latitude' in d: 
+                dimgrid     = ncf.AddLatLonDim()
+                dims        += (dimgrid[0],)
+            if 'longitude' in d: 
+                dimgrid     = ncf.AddLatLonDim()
+                dims        += (dimgrid[1],)
+            for type in ['eco','eco_ext','tc','tc_ext','olson']:
+                if 'regions_%s' % type == d: 
+                    dimregs     = ncf.AddRegionDim(type)
+                    dims        += dimregs
+#
+# If the variable name from the infile is not defined in the standard, simply copy attributes from the old to the new sds
+# Do not copy over the grid information
+#
+
+        if ncf.StandardVar(sds)['name'] == 'unknown':
+
+                savedict                    = ncf.StandardVar(sds)
+                savedict['name']            = sds
+                savedict['values']          = data.tolist()
+                savedict['dims']            = dims
+                savedict['units']           = file.variables[sds].units
+                savedict['long_name']       = file.variables[sds].long_name
+                savedict['comment']         = file.variables[sds].comment
+                savedict['standard_name']   = file.variables[sds].standard_name
+                savedict['count']           =0
+                ncf.AddData(savedict)
+                sys.stdout.write('.')
+                sys.stdout.flush()
+
+        else:
+
+            if sds in ['latitude','longitude']: continue
+
+            if avg == 'monthly':
+                time_avg, data_avg      = timetools.monthly_avg(time,data)
+            elif avg == 'seasonal':
+                time_avg, data_avg      = timetools.season_avg(time,data)
+            elif avg == 'yearly':
+                time_avg, data_avg      = timetools.yearly_avg(time,data)
+            elif avg == 'longterm':
+                time_avg, data_avg      = timetools.longterm_avg(time,data)
+                time_avg                = [time_avg]
+                data_avg                = [data_avg]
+            else:
+                raise ValueError,'Averaging (%s) does not exist' % avg
+
+            count=-1
+            for dd,data in zip(time_avg,data_avg):
+                count=count+1
+
+                if not ncf.has_date(date2num(dd)-dectime0): 
+
+                    savedict=ncf.StandardVar('date')
+                    savedict['values']=date2num(dd)-dectime0
+                    savedict['dims']=dimdate
+                    savedict['count']=count
+                    ncf.AddData(savedict)
+
+                savedict=ncf.StandardVar(sds)
+                savedict['values']=data.tolist()
+                savedict['units']=file.variables[sds].units
+                if 'cov' in sds: 
+                    savedict['dims']=dimdate+dims+dims
+                else:
+                    savedict['dims']=dimdate+dims
+                savedict['count']=count
+                ncf.AddData(savedict)
+
+                sys.stdout.write('.')
+                sys.stdout.flush()
+        print ']'
+
+
+# end NetCDF file access
+    file.close()
+    ncf.close()
+
+if __name__ == "__main__":
+
+    import logging
+    from da.tools.initexit import CycleControl
+    from da.ct.dasystem import CtDaSystem 
+    from da.ctgridded.statevector import CtGriddedStateVector 
+
+    sys.path.append('../../')
+
+    logging.root.setLevel(logging.DEBUG)
+
+    DaCycle = CycleControl(args={'rc':'../../dagridded.rc'})
+    DaCycle.Initialize()
+    DaCycle.ParseTimes()
+
+    DaSystem    = CtDaSystem('../rc/carbontrackergridded.rc')
+    DaSystem.Initialize()
+
+    DaCycle.DaSystem    = DaSystem
+
+    StateVector = CtGriddedStateVector(DaCycle)
+    dummy       = StateVector.Initialize()
+
+
+    no, yes, all = range(3)
+
+    proceed = no
+
+    # first, parse results from the inverse part of the run 
+    if proceed != all:
+        proceed = proceed_dialog('Create average 1x1 flux data files? [y|yes, n|no, a|all|yes-to-all ]')
+    if proceed != no:
+        while DaCycle['time.start'] < DaCycle['time.finish']:
+
+            savedas=SaveWeeklyAvg1x1Data(DaCycle, StateVector)
+            savedas=SaveWeeklyAvgStateData(DaCycle, StateVector)
+            savedas=SaveWeeklyAvgTCData(DaCycle, StateVector)
+
+            DaCycle.AdvanceCycleTimes()
+
+        sys.exit(2)
+
+        a=SaveTimeAvgData(rundat,savedas,avg='monthly')
+        a=SaveTimeAvgData(rundat,savedas,avg='yearly')
+        a=SaveTimeAvgData(rundat,savedas,avg='longterm')
+
+    if proceed != all:
+        proceed = proceed_dialog('Create average tc flux data files? [y|yes, n|no, a|all|yes-to-all] ')
+    if proceed != no:
+        a=SaveTimeAvgData(rundat,savedas,avg='monthly')
+        a=SaveTimeAvgData(rundat,savedas,avg='seasonal')
+        a=SaveTimeAvgData(rundat,savedas,avg='yearly')
+        a=SaveTimeAvgData(rundat,savedas,avg='longterm')
+    if proceed != all:
+        proceed = proceed_dialog('Create average extended tc flux data files? [y|yes, n|no, a|all|yes-to-all] ')
+    if proceed != no:
+        savedas=SaveTCDataExt(rundat)
+        a=SaveTimeAvgData(rundat,savedas,avg='monthly')
+        a=SaveTimeAvgData(rundat,savedas,avg='seasonal')
+        a=SaveTimeAvgData(rundat,savedas,avg='yearly')
+        a=SaveTimeAvgData(rundat,savedas,avg='longterm')
+    sys.exit(0)

refactored/da/analysis/olson19_region_names

+ 1 CONIFOR  "Conifer Forest"
+ 2 BRODLFOR "Broadleaf Forest: temperate, subtropical drought"
+ 3 MIXEDFOR "Mixed Forest: conifer/broadleaf; so. Boreal"
+ 4 GRASSHRB "Grassland +/- Shrub or Tree"
+ 5 TROPICFR "Tropical/subtr. Forest: montane, seasonal, rainforest "
+ 6 SCRUBWDS "Scrub +/- Woodland &/or fields (evergreen/decid.) "
+ 7 SEMIDTUN "Semidesert shrub/steppe; Tundra (polar, alpine) "
+ 8 FLDWDSAV "Field/Woods complex &/or Savanna, tallgrass "
+ 9 NORTAIGA "Northern Boreal Taiga woodland/tundra"
+10 FORFDREV "Forest/Field; Dry Evergreen broadleaf woods "
+11 WETLAND  "Wetlands: mires (bog/fen); marsh/swamp +/- mangrove "
+12 DESERTS  "Desert: bare/alpine; sandy; salt/soda "
+13 SHRBTRST "Shrub/Tree: succulent/thorn "
+14 CROPSETL "Crop/Settlement (irrigated or not) "
+15 CONIFRFC "Conifer snowy Rainforest, coastal      "
+27 WTNDMHTH "Wooded Tundra Margin; Heath/moorland "
+19 MANGROVE "Mangrove/wet forest/thicket + tidal flat "
+25 ICEPLDE  "Ice and Polar Desert"
+19 WATER    "Water bodies"

refactored/da/analysis/runinfo.py

+#!/usr/bin/env python
+# runinfo.py
+
+"""
+Author : peters 
+
+Revision History:
+File created on 01 Mar 2011.
+
+"""
+
+import mysettings 
+import commands
+import os
+import sys
+from datetime import datetime, timedelta
+
+
+class RunInfo(object):
+    """ 
+    This is a CarbonTracker  run info object
+    Initialize it with the following keywords:
+   
+    MANDATORY: 
+    project  = Name of project, should correspond to a directory [string]
+    sd       = start date of project output [datetime object or integer yyyymmdd]
+    ed       = end date of project output [datetime object or integer yyyymmdd]
+    OPTIONAL:
+    basedir  = base directory of model output on your machine [string, default='~/Modeling/download']
+    outputdir= output directory for figures and summary files [string, default='~/Modeling/SEATA']
+
+    Tip: change the std_rundat dictionary in module filtertools to set new defaults for all values!
+    """
+    def __init__(self,infodict={},rcfile=None):
+        """ Initialize from the name of an rcfile, or by passing the dictionary directly """
+        import da.tools.rc as rc
+        import da.tools.io4 as io
+        import copy
+
+
+        if rcfile:
+            infodict    = rc.read(rcfile)
+
+        self.basedir    = infodict['dir.da_run']
+        self.proj       = os.path.split(self.basedir)[-1]
+        self.projdir    = self.basedir
+
+        self.dir_scheme = mysettings.ct_dir_scheme
+        
+        self.inputdir   = os.path.join(self.basedir,'output')
+        self.outputdir  = os.path.join(self.basedir,'analysis')
+
+        if not os.path.exists(self.outputdir): 
+            os.makedirs(self.outputdir)
+            print "Creating new output directory "+self.outputdir
+ 
+        sd=infodict['time.start']
+        ed=infodict['time.finish']
+        dt=infodict['time.cycle']
+
+        self.dt = timedelta(days=int(dt))
+
+        if not isinstance(sd,datetime):
+            self.sd = datetime.strptime(sd,'%Y-%m-%d %H:%M:%S')
+        if not isinstance(ed,datetime):
+            self.ed = datetime.strptime(ed,'%Y-%m-%d %H:%M:%S')
+
+        dd             = copy.deepcopy(self.sd)
+        self.inputfiles= []
+        self.weeks     = []
+        self.inputdict = {}
+
+        while dd < self.ed:
+            filename = os.path.join(self.inputdir,'%s'%dd.strftime('%Y%m%d'),'savestate.nc')
+            if os.path.exists(filename):
+                self.inputfiles.append(filename)
+                self.weeks.append(dd)
+                self.inputdict[dd] = filename
+            else:
+                break
+                    
+            dd = dd+self.dt
+
+        self.ed     = dd
+        self.nweeks = len(self.weeks)
+
+        # run parameters from file
+        ncf             = io.CT_Read(self.inputfiles[0],'read')
+        self.nlag       = len(ncf.dimensions['nlag'])
+        self.nmembers   = len(ncf.dimensions['nmembers'])
+        self.nparameters  = len(ncf.dimensions['nparameters'])
+        ncf.close()
+
+        self.namingscheme = 'wp_Mar2011'
+        filename          = os.path.join('NamingScheme.'+self.namingscheme+'.rc')
+
+        try: 
+            rcinfo=rc.read(filename)
+        except IOError: 
+            print '%s was specified as naming scheme, but no matching %s rc-file was found, exiting...'%(NamingScheme,filename,)
+            sys.exit(1)
+        except:
+            print 'Unknown error reading rc-file: %s, exiting...'%(filename,) 
+            sys.exit(1)
+
+
+        self.namedict=rcinfo
+    
+    def __str__(self):
+        return 'project              : '+self.projdir+\
+               '\nstart date           : '+str(self.sd)+\
+               '\nend date             : '+str(self.ed)+\
+               '\ndelta date           : '+str(self.dt)+\
+               '\nnweeks               : '+str(self.nweeks)+\
+               '\nnparameters          : '+str(self.nparameters)+\
+               '\nnmembers             : '+str(self.nmembers)+\
+               '\nnlag                 : '+str(self.nlag)+\
+               '\nDA output dir        : '+self.inputdir+\
+               '\nanalysis output dir  : '+self.outputdir+\
+               '\nnaming scheme        : '+self.namingscheme
+
+def get_rundat_settings(args):
+    """ create settings dict for rundat from scripts arguments """
+
+    settings=std_rundat.copy()  # start with defaults for rundat
+
+    for items in args:
+        items=items.lower()
+        k, v = items.split('=')
+        if settings.has_key(k):
+            if k in std_constructors:
+                settings[k] = std_constructors[k](v)
+            else:
+                settings[k] = v
+        else: raise IOError,'Parameter unknown:%s'%(v,)
+
+    return settings
+
+
+
+
+
+if __name__ == "__main__":
+    import getopt
+    import sys
+    from string import *
+
+    sys.path.append('../../')
+
+# Parse keywords
+
+    rundat=RunInfo(rcfile='../../da.rc')
+
+    print rundat
+
+    sys.exit(0)

refactored/da/analysis/t3_region_names

+ 1 T-NABR "North American Boreal"
+ 2 T-NATM "North American Temperate"
+ 3 T-SATR "South American Tropical"
+ 4 T-SATM "South American Temperate"
+ 5 T-NAFR "Northern Africa"
+ 6 T-SAFR "Southern Africa"
+ 7 T-EUBR "Eurasia Boreal"
+ 8 T-EUTM "Eurasia Temperate"
+ 9 T-ASTR "Tropical Asia"
+10 T-AUST "Australia"
+11 T-EURO "Europe"
+12 O-NPTM "North Pacific Temperate"
+13 O-WPTR "West Pacific Tropical"
+14 O-EPTR "East Pacific Tropical"
+15 O-SPTM "South Pacific Temperate"
+16 O-NOCN "Northern Ocean"
+17 O-NATM "North Atlantic Temperate"
+18 O-SATR "Atlantic Tropical"
+19 O-SATM "South Atlantic Temperate"
+20 O-SOCN "Southern Ocean"
+21 O-INTR "Indian Tropical"
+22 O-INTM "Indian Temperate"
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

refactored/da/analysis/tools_regions.py

+#!/usr/bin/env python
+import commands
+import os
+import sys
+from da.analysis.tools_transcom import *
+
+# Aggregated olson ecosystem regions for CT Europe
+
+aggregates=[\
+("Forests/Wooded" , (1, 2, 3, 5, 8, 10,) ) ,\
+("Grass/Shrubs" , (4, 6, 13) ) ,\
+("Crops/Agriculture", (14,) ) ,\
+("Tundra/Taiga" , (7,9,16) ) ,\
+("Coastal/Wet" , (11,15,17) ) ,\
+("Ice/Water/Deserts" , (12,18,19) ) \
+]
+
+agg_dict={}
+for k,v in aggregates:
+    agg_dict[k]=v
+
+ext_econams=[a for a,b in aggregates]
+ext_ecocomps=[b for a,b in aggregates]
+
+eco19_to_ecosums = zeros((19, 6),float)
+for i,k in enumerate(ext_ecocomps):
+    indices = [x-1 for x in k]
+    eco19_to_ecosums[:,i].put(indices,1.0)
+
+##### END OF REGION DEFINITIONS
+
+def StateToGrid(values,regionmap,reverse=False,avg=False):
+    """ 
+    This method converts parameters from a CarbonTracker StateVector object to a gridded map of linear multiplication values. These
+    can subsequently be used in the transport model code to multiply/manipulate fluxes
+
+    """
+    import numpy as np
+    nregions  = regionmap.max()
+
+    # dictionary for region <-> map conversions
+
+    regs={}
+    for r in np.arange(1,nregions+1):
+         sel=(regionmap.flat == r).nonzero()
+         if len(sel[0])>0: regs[r]=sel
+
+    regionselect=regs
+
+    if reverse:
+
+        """ project 1x1 degree map onto ecoregions """
+
+        result=np.zeros(nregions,float)
+        for k,v in regionselect.iteritems():
+            if avg: 
+                result[k-1]=values.ravel().take(v).mean()
+            else : 
+                result[k-1]=values.ravel().take(v).sum()
+        return result
+
+    else:
+
+        """ project ecoregion properties onto 1x1 degree map """
+
+        result=np.zeros((180,360,),float)
+        for k,v in regionselect.iteritems():
+            result.put(v,values[k-1])
+
+        return result
+
+def globarea(im=360,jm=180,silent=True):
+    """ Function calculates the surface area according to TM5 definitions"""
+    import numpy as np
+
+    radius      = 6.371e6  # the earth radius in meters
+    deg2rad     = np.pi/180.
+    g           = 9.80665 
+
+    dxx         = 360.0/im*deg2rad 
+    dyy         = 180.0/jm*deg2rad 
+    lat         = np.arange(-90*deg2rad,90*deg2rad,dyy)
+    dxy         = dxx*(np.sin(lat+dyy)-np.sin(lat))*radius**2
+    area        = np.resize(np.repeat(dxy,im,axis=0) ,[jm,im])
+    if not silent:
+        print 'total area of field = ',np.sum(area.flat)
+        print 'total earth area    = ',4*np.pi*radius**2
+    return area
+

refactored/da/analysis/tools_time.py

+#! /usr/bin/env python
+from datetime import datetime, timedelta
+import calendar
+from pylab import floor
+from matplotlib.dates import date2num, num2date
+from numpy import array, zeros, newaxis
+
+def Fromdatetime(date):
+    dt=date.timetuple()
+    return datetime(*dt[0:6])
+
+def increase_time(dd,**kwargs):
+    """ Function increases the time by specified amount"""
+    return dd+timedelta(**kwargs)
+
+def chardate(dd,cut=8):
+    return dd.strftime('%Y%m%d%H%M')[0:cut]
+
+def timegen(sd,ed,dt):
+    dd=[]
+    while sd <= ed:
+        dd.append(Fromdatetime(sd))
+        sd=sd+dt
+    return dd
+
+def itau2datetime(itau,iyear0):
+    """ Function returns a datetime object from TM5s itau times"""
+    date0=datetime(iyear0,1,1,0,0,0)
+    if len(itau)==1:
+           itau=[itau] 
+    for time in itau:
+            sec=time%60
+            min=(time/60)%60
+            hrs=(time/3600)%24
+            day=(time/86400) 
+            dt=timedelta(days=day,hours=hrs,minutes=min,seconds=sec)
+            yield date0+dt 
+
+def date2dec(date):
+    """ Function converts datetime object to a Decimal number time  (e.g., 1991.875 such as in IDL, CCG) """
+    if not isinstance(date,list): date=[date]
+
+    newdate=[]
+    for dd in date:
+        Days0=date2num(datetime(dd.year,1,1))
+        if calendar.isleap(dd.year):
+            DaysPerYear=366.
+        else:
+            DaysPerYear=365.
+        DayFrac=date2num(dd)
+        newdate.append(dd.year+(DayFrac-Days0)/DaysPerYear)
+    if len(newdate) == 1: return newdate[0] 
+    return newdate
+
+def dec2date(dectime):
+    """ Function converts decimal time from year fraction (e.g., 1991.875 such as in IDL, CCG) to a python datetime object """
+    dt=num2date(dec2num(dectime)).timetuple()
+    return datetime(*dt[0:7])
+
+def dec2num(dectime):
+    """ Function converts decimal time from year fraction (e.g., 1991.875 such as in IDL, CCG) to a python decimal numtime """
+    if not isinstance(dectime,list): dectime=[dectime]
+
+    newdectime=[]
+    for dd in dectime:
+        yr=floor(dd)
+        Days0=date2num(datetime(int(yr),1,1))
+        if calendar.isleap(yr):
+            DaysPerYear=366.
+        else:
+            DaysPerYear=365.
+        DayFrac=(dd-yr)*DaysPerYear
+        newdectime.append(Days0+DayFrac)
+    if len(newdectime) == 1: return newdectime[0] 
+    return newdectime
+
+def num2dec(numtime):
+    """ Function converts python decimal numtime to an IDL decimal time """
+    res=date2dec(num2mydate(numtime))
+    return res
+
+def num2mydate(num):
+    """ Function converts decimal time from year fraction (e.g., 1991.875 such as in IDL, CCG) to a python datetime object """
+    dt=num2date(num).timetuple()
+    return datetime(*dt[0:7])
+
+def monthgen(sd,ed):
+    """ Generate sequence of datetime objects spaced by one month"""
+    from pylab import arange
+    if ed<sd: 
+        raise ValueError,'start date exceeds end date'
+        sys.exit(2)
+    dates=[]
+    for year in arange(sd.year,ed.year+2):
+        for month in arange(1,13):
+            date=datetime(year,month,1)
+            if date > ed: return dates
+            else: dates.append(date)
+
+def nextmonth(dd):
+    """ Find next 1st of the month following the date dd"""
+
+    if dd.month == 12: 
+        cc=dd.replace(year=dd.year+1)
+        ee=cc.replace(month=1)
+    else:
+        ee=dd.replace(month=dd.month+1)
+    ff = ee.replace(day=1)
+    return ff
+
+def in_interval(start,stop,times_in):
+    """ returns a list of fractions in time interval """
+    from numpy import logical_and, arange
+    from pylab import drange, num2date, date2num
+    from datetime import timedelta
+    import copy
+
+    times=copy.copy(times_in)
+
+    interval=times[1]-times[0]
+    times.append(times[-1]+interval)  # extend by one interval
+    times_filled=[times[0]+timedelta(days=d) for d in range((times[-1]-times[0]).days)]
+
+    b=[]
+    in_int=0.0
+    for t in times_filled:   # loop over days
+        if t in times[1:]:   # if new interval starts
+            b.append(in_int) # add previous aggregate to output
+            in_int=0.0         # reset counter
+        in_int+=int(logical_and(t >= start,t < stop))  # count if in interval [start,stop >
+    b.append(in_int)
+
+    if len(b) != len(times_in) : raise ValueError
+
+    return b
+
+def yearly_avg(time,data,sdev=False):
+    """ make monthly average from array using rundat and data"""
+
+    years=array([d.year for d in time])
+    
+    aa=[]
+    ss=[]
+    tt=[]
+    dd=time[0]
+    ed=time[-1]
+    while dd <= ed:
+        ddnext=datetime(dd.year+1,1,1)
+        weights=in_interval(dd,ddnext,time)
+        if len(weights) > 1:
+            weights=array(weights)
+            if weights.sum() > 0.0:
+                weights=weights/weights.sum()
+            else:
+                weights=weights
+
+            if weights.shape[0] != data.shape[0]:
+                raise ValueError,'yearly_avg has wrongly shaped weights (%d) for data of (%d)' %(weights.shape[0],data.shape[0])
+                 
+            sel=(weights != 0.0).nonzero()[0]     
+            #print sel,array(time).take(sel),dd,ddnext
+            if data.ndim == 1:
+                avg_data=(weights.take(sel)*data.take(sel,axis=0)).sum(axis=0)
+                std_data=(weights.take(sel)*data.take(sel,axis=0)).std(axis=0)
+            elif data.ndim == 2:
+                avg_data=(weights.take(sel)[:,newaxis]*data.take(sel,axis=0)).sum(axis=0).squeeze()
+                std_data=(weights.take(sel)[:,newaxis]*data.take(sel,axis=0)).std(axis=0).squeeze()
+            elif data.ndim == 3:
+                avg_data=(weights.take(sel)[:,newaxis,newaxis]*data.take(sel,axis=0)).sum(axis=0).squeeze()
+                std_data=(weights.take(sel)[:,newaxis,newaxis]*data.take(sel,axis=0)).std(axis=0).squeeze()
+            else:
+                raise ValueError,'yearly_avg takes 1, 2, or 3d arrays only'
+        elif len(weights)==1:    
+            avg_data=data[0]
+            std_data=0.0
+        else:
+            continue # next year
+
+        aa.append(avg_data)
+        ss.append(std_data)
+        tt.append(datetime(dd.year,6,15))
+
+        dd=ddnext
+
+    aa=array(aa).squeeze()
+    ss=array(ss).squeeze()
+    time=tt
+    if len(tt) == 1: 
+        aa=aa.reshape(1,*aa.shape)
+        ss=ss.reshape(1,*ss.shape)
+    if sdev: return time,aa, ss
+    else : return time,aa
+
+def monthly_avg(time,data,sdev=False):
+    """ make monthly average from array using rundat and data"""
+
+    years=array([d.year for d in time])
+    months=array([d.month for d in time])
+    
+    mm=[]
+    ss=[]
+    tt=[]
+    dd=time[0]
+    ed=time[-1]
+
+    while dd <= ed:
+        ddnext=nextmonth(dd)
+        weights=in_interval(dd,ddnext,time)
+        if len(weights) > 1:
+            weights=array(weights)
+            if weights.sum() > 0.0:
+                weights=weights/weights.sum()
+            else:
+                weights=weights
+
+            if weights.shape[0] != data.shape[0]:
+                raise ValueError,'yearly_avg has wrongly shaped weights (%d) for data of (%d)' %(weights.shape[0],data.shape[0])
+            
+            sel=(weights != 0.0).nonzero()[0]     
+            #print sel,array(time).take(sel),dd,nextmonth(dd)
+            if data.ndim == 1:
+                avg_data=(weights.take(sel)*data.take(sel,axis=0)).sum(axis=0)
+                std_data=(weights.take(sel)*data.take(sel,axis=0)).std(axis=0)
+            elif data.ndim == 2:
+                avg_data=(weights.take(sel)[:,newaxis]*data.take(sel,axis=0)).sum(axis=0).squeeze()
+                std_data=(weights.take(sel)[:,newaxis]*data.take(sel,axis=0)).std(axis=0).squeeze()
+            elif data.ndim == 3:
+                avg_data=(weights.take(sel)[:,newaxis,newaxis]*data.take(sel,axis=0)).sum(axis=0).squeeze()
+                std_data=(weights.take(sel)[:,newaxis,newaxis]*data.take(sel,axis=0)).std(axis=0).squeeze()
+            else:
+                raise ValueError,'monthly_avg takes 1, 2, or 3d arrays only'
+        elif len(weights)==1:    
+            avg_data=data[0]
+            std_data=0.0
+        else:
+            continue # next month
+
+        mm.append(avg_data)
+        ss.append(std_data)
+        tt.append(datetime(dd.year,dd.month,15))
+
+        dd=ddnext
+
+
+    mm=array(mm).squeeze()
+    ss=array(ss).squeeze()
+    time=tt
+
+    if len(tt) == 1: 
+        mm=mm.reshape(-1,*mm.shape)
+        ss=ss.reshape(-1,*ss.shape)
+
+    if sdev: return time,mm, ss
+    else : return time,mm
+
+def season_avg(time,data,sdev=False):
+    """ make season average from array using rundat and data"""
+
+    seas = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]]
+
+    mm=[]
+    ss=[]
+    tt=[]
+    dd=time[0]
+    ed=time[-1]
+
+    while dd <= ed:
+        ddmid=nextmonth(dd)
+        ddnext=nextmonth(nextmonth(nextmonth(dd)))
+        weights=in_interval(dd,ddnext,time)
+        if len(weights) > 1:
+            weights=array(weights)
+            if weights.sum() > 0.0:
+                weights=weights/weights.sum()
+            else:
+                weights=weights
+
+            if weights.shape[0] != data.shape[0]:
+                raise ValueError,'yearly_avg has wrongly shaped weights (%d) for data of (%d)' %(weights.shape[0],data.shape[0])
+            
+            sel=(weights != 0.0).nonzero()[0]     
+            #print sel,array(time).take(sel),dd,nextmonth(dd)
+            if data.ndim == 1:
+                avg_data=(weights.take(sel)*data.take(sel,axis=0)).sum(axis=0)
+                std_data=(weights.take(sel)*data.take(sel,axis=0)).std(axis=0)
+            elif data.ndim == 2:
+                avg_data=(weights.take(sel)[:,newaxis]*data.take(sel,axis=0)).sum(axis=0).squeeze()
+                std_data=(weights.take(sel)[:,newaxis]*data.take(sel,axis=0)).std(axis=0).squeeze()
+            elif data.ndim == 3:
+                avg_data=(weights.take(sel)[:,newaxis,newaxis]*data.take(sel,axis=0)).sum(axis=0).squeeze()
+                std_data=(weights.take(sel)[:,newaxis,newaxis]*data.take(sel,axis=0)).std(axis=0).squeeze()
+            else:
+                raise ValueError,'season_avg takes 1, 2, or 3d arrays only'
+        elif len(weights)==1:    
+            avg_data=data[0]
+            std_data=0.0
+        else:
+            continue # next month
+
+        mm.append(avg_data)
+        ss.append(std_data)
+        tt.append(datetime(ddmid.year,ddmid.month,15))
+
+        dd=ddnext
+
+
+    mm=array(mm).squeeze()
+    ss=array(ss).squeeze()
+    time=tt
+
+    if len(tt) == 1: 
+        mm=mm.reshape(-1,*mm.shape)
+        ss=ss.reshape(-1,*ss.shape)
+
+    if sdev: return time,mm, ss
+    else : return time,mm
+
+def longterm_avg(time,data):
+    """ Create long term mean """
+
+    time_avg = num2date(date2num(time).mean())
+    data_avg = data.mean(axis=0)
+
+    return time_avg, data_avg
+
+
+
+if __name__ == '__main__':
+    #print monthgen(datetime(2000,1,1),datetime(2006,5,1))
+    dd=datetime(2002,3,1)
+    print nextmonth(dd),dd
+
+        
+         

refactored/da/analysis/tools_transcom.py

+#!/usr/bin/env python
+
+import os
+import sys
+sys.path.append('../../')
+rootdir     = os.getcwd().split('da/')[0]
+analysisdir = os.path.join(rootdir,'da/analysis')
+
+from string import join, split
+from numpy import array, identity, zeros, arange
+import da.tools.io4 as io
+
+# Get masks of different region definitions
+
+matrix_file     = os.path.join(analysisdir,'regions.nc')
+cdf_temp        = io.CT_CDF(matrix_file,'read')
+transcommask    = cdf_temp.GetVariable('transcom_regions')
+olson240mask    = cdf_temp.GetVariable('regions')
+olsonmask       = cdf_temp.GetVariable('land_ecosystems')
+oifmask         = cdf_temp.GetVariable('ocean_regions')
+dummy           = cdf_temp.close()
+
+# Names and short names of TransCom regions
+
+transshort  = []
+transnams   = []
+transland   = []
+temp        = open(os.path.join(analysisdir,'t3_region_names'),'r').readlines()
+for line in temp:
+    items   = line.split()
+    if items:
+        num,abbr,name=(items[0],items[1],join(items[2:]),)
+        transnams.append(name.strip('"'))
+        transshort.append(abbr)
+        if abbr.startswith('T'): transland.append(name.strip('"'))
+
+transnams.append("Non-optimized")
+transshort.append("I-NNOP")
+
+# Names and short names of Olson regions
+
+olsonnams   = []
+olsonshort  = []
+temp        = open(os.path.join(analysisdir,'olson19_region_names'),'r').readlines()
+for line in temp:
+    items   = line.split()
+    if items:
+        num,abbr,name=(items[0],items[1],join(items[2:]),)
+        olsonnams.append(name.strip('"'))
+        olsonshort.append(abbr)
+
+ext_transnams   = []
+ext_transshort  = []
+ext_transcomps  = []
+
+# Get names of aggregated regions for post aggregation
+
+matrix_file     = os.path.join(analysisdir,'postagg_definitions.nc')
+cdf_temp        = io.CT_CDF(matrix_file,'read')
+xform           = cdf_temp.GetVariable('xform')
+keys            = cdf_temp.ncattrs()
+
+keys.sort()
+for k in keys:
+    if 'longname' in k: 
+        ext_transnams.append(getattr(cdf_temp,k) )
+    if 'shortname' in k: 
+        ext_transshort.append(getattr(cdf_temp,k) )
+    if 'component' in k: 
+        ext_transcomps.append(map(int,getattr(cdf_temp,k).split(',')) )
+
+dummy           = cdf_temp.close()
+
+
+# Names of the ocean inversion flux regions, to go along with oifmask
+
+oifnams=['(1) NOCN Arctic Ocean',\
+         '(2) NAH  North Atlantic (49 - 76N)',\
+         '(3) NAM  North Atlantic (36 - 49N)',\
+         '(4) NAL  North Atlantic (18 - 36N)',\
+         '(5) NAT  North Atlantic ( 0 - 18N)',\
+         '(6) SAT  South Atlantic ( 0 - 18S)',\
+         '(7) SAL  South Atlantic (18 - 31S)',\
+         '(8) SAM  South Atlantic (31 - 44S)',\
+         '(9) SAH  South Atlantic (44 - 58S)',\
+         '(10) SOCN Southern Ocean (S of 58S)',\
+         '(11) NPHW North Pacific (N of 49N, W of 195E)',\
+         '(12) NPHE North Pacific (N of 36N, E of 195E)',\
+         '(13) NPK  North Pacific (Kuroshio Extension)',\
+         '(14) NPLW North Pacific (18N - K.Ext, W of 195E)',\
+         '(15) NPLE North Pacific (18 - 36N, E of 195E)',\
+         '(16) NPTW North Pacific ( 0 - 18N, W of 199E)',\
+         '(17) NPTE North Pacific ( 0 - 18N, E of 199E)',\
+         '(18) SPTW South Pacific ( 0 - 18S, W of 199E)',\
+         '(19) SPTE South Pacific ( 0 - 18S, E of 199E)',\
+         '(20) SPLW South Pacific (18 - 31S, W of 233E)',\
+         '(21) SPLE South Pacific (18 - 31S, E of 233E)',\
+         '(22) SPMW South Pacific (31 - 44S, W of 248E)',\
+         '(23) SPME South Pacific (31 - 44S, E of 248E, W of 278E)',\
+         '(24) SPMC South Pacific (31 - 44S, coastal E of 278E)',\
+         '(25) SPH  South Pacific (44 - 58S)              ',\
+         '(26) NI   North Indian',\
+         '(27) SIT  South Indian (0 - 18S)',\
+         '(28) SIL  South Indian (18 - 31S)',\
+         '(29) SIM  South Indian (31 - 44S)',\
+         '(30) SIH  South Indian (44 - 58S)']
+
+oiflocs=[ (200, 80,),\
+          (330, 55,),\
+          (330, 40,),\
+          (330, 22,),\
+          (330,  8,),\
+          (350,-12,),\
+          (350,-27,),\
+          (350,-40,),\
+          (350,-53,),\
+          (200,-70,),\
+          (178, 54,),\
+          (210, 40,),\
+          (165, 38,),\
+          (178, 25,),\
+          (215, 25,),\
+          (170, 8,),\
+          (230, 8,),\
+          (175,-10,),\
+          (240,-10,),\
+          (195,-27,),\
+          (265,-27,),\
+          (195,-40,),\
+          (262,-40,),\
+          (283,-40,),\
+          (220,-53,),\
+          (68,  8,),\
+          (75, -10,),\
+          (75, -27,),\
+          (75,-40,),\
+          (75, -53,)]
+
+
+translocs=[ (-177,0), \
+            (-92,53,),\
+            (-108,34,),\
+            (-66,4,),\
+            (-50,-17,),\
+            (15,17,),\
+            (26,-12,),\
+            (84,63,),\
+            (103,30,),\
+            (115,0,),\
+            (132,-25,),\
+            (9,50,),\
+            (-174,46,),\
+            (136,6,),\
+            (-108,6,),\
+            (-123,-15,),\
+            (-32,58,),\
+            (-32,38,),\
+            (-32,0,),\
+            (-32,-38,),\
+            (-14,-65,),\
+            (68,2,)]
+
+#olsonshort=[str(name.split()[1:2]).join('  ') for name in olsonnams]
+old_olsonshort=[join(split(name,' ')[1:2],' ') for name in olsonnams]
+
+olsonlabs=['Conifer Forest','Broadleaf Forest','Mixed Forest','Grass/Shrub','Tropical Forest','Scrub/Woods','Semitundra','Fields/Woods/\nSavanna',\
+     'Northern Taiga','Forest/Field','Wetland','Deserts','Shrub/Tree/\nSuc ','Crops','Conifer\n Snowy/Coastal',\
+     'Wooded tundra','Mangrove','Ice and \nPolar desert','Water'] 
+
+ecmwfnams=[ ' 1 CRPSMF Crops, mixed farming',\
+            ' 2 SHGRSS Short Grass',\
+            ' 3 EVNDLF Evergreen Needleleaf',\
+            ' 4 DECNDLF Deciduous Needleleaf',\
+            ' 5 EVBRDLF Evergreen Broadleaf',\
+            ' 6 DECBRLF Deciduous Broadleaf',\
+            ' 7 TLGRSS Tall Grass',\
+            ' 8 DES Desert',\
+            ' 9 TDR Tundra',\
+            '10 IRRCR Irrigated Crops',\
+            '11 SMDES Semidesert',\
+            '12 ICE Ice Caps',\
+            '13 BGM Bogs and Marches',\
+            '14 INW Inland Water',\
+            '15 OCE Ocean',\
+            '16 EVSHRB Evergreen Shrubs',\
+            '17 DECSHR Deciduous shrubs',\
+            '18 MXFRST Mixed Forest',\
+            '19 INTFRST Interrupted Forest'] 
+
+ecmwfshort=[str(name.split()[1:2]).join('  ') for name in ecmwfnams]
+
+ecmwflabs=['Crops, mixed farming','Short Grass','Evergreen Needleleaf','Deciduous Needleleaf','Evergreen Broadleaf',\
+      'Deciduous Broadleaf','Tall Grass','Desert',\
+     'Tundra','Irrigated Crops','Semidesert','Ice Caps','Bogs and Marches','Inland Water','Ocean',\
+     'Evergreen Shrubs','Deciduous shrubs','Mixed Forest','Interrupted Forest'] 
+
+a =array([\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,1  ,1  ,1  ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,0  ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,0  ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,1  ,1  ,1  ,1  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  1 ,1 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,1 ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,1  ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,1  ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,1  ,1  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,0  ,0  ,0  ,0  ,1 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,1  ,0  ,0  ,0 ,\
+  0 ,0 ,0 ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,0  ,1  ,1  ,0])
+O30to11     = a.reshape(11,30).transpose()
+
+O11to11     = identity(11)
+
+ntcland     = 11 # TC standard
+ntcocean    = 11 # TC standard
+
+Ols259_to_TC23  = zeros((259,23),float)
+Ols240_to_TC23  = zeros((240,23),float)
+Ols221_to_TC23  = zeros((221,23),float)
+for i in arange(ntcland): 
+    Ols259_to_TC23[i*19:(i+1)*19,i]=1.0
+
+Ols259_to_TC23[190:228,10]      = 1.0  # Europe
+Ols259_to_TC23[228:258,11:22]   = O30to11
+for i in arange(ntcland): Ols240_to_TC23[i*19:(i+1)*19,i]=1.0
+Ols240_to_TC23[209:239,11:22]   = O30to11
+for i in arange(ntcland): Ols221_to_TC23[i*19:(i+1)*19,i]=1.0
+Ols221_to_TC23[209:220:,11:22]  = O11to11
+
+Ols221_to_TC23[220,22]=1.0
+Ols240_to_TC23[239,22]=1.0
+Ols259_to_TC23[258,22]=1.0
+
+
+ntcland     = 11 # TC standard
+ntcocean    = 11 # TC standard
+
+ExtendedTCRegionsFile='postagg_definitions.nc'
+
+def StateToTranscom(runinfo,x):
+    """ convert to transcom shaped regions"""
+    from numpy import dot
+
+    try:
+        nparams     = runinfo.nparameters 
+    except:
+        nparams     = runinfo
+    if nparams==240: 
+        M=Ols240_to_TC23
+    elif nparams==221: 
+        M=Ols221_to_TC23
+    elif nparams==259: 
+        M=Ols259_to_TC23
+    else:
+        raise ValueError('Do not know how to convert %s regions to 23 transcom regions'%(nparams,))
+
+    return dot(array(x).squeeze(),M)
+
+def StateCovToTranscom(runinfo,p):
+    """ convert to transcom shaped regions"""
+    try:
+        nparams     = runinfo.nparameters 
+    except:
+        nparams     = runinfo
+    if nparams==240: 
+        M=Ols240_to_TC23
+    elif nparams==221: 
+        M=Ols221_to_TC23
+    elif nparams==259: 
+        M=Ols259_to_TC23
+    else:
+        raise ValueError('Do not know how to convert %s regions to 23 transcom regions'%(nparams,))
+
+    return transpose(dot(dot(transpose(M),p),M),(1,0,2))
+
+def ExtendedTCRegions(data,cov=False):
+    """ convert to extended transcom shaped regions"""
+    import da.tools.io4 as io 
+    from numpy import dot
+
+    nparams     = data.shape[-1]
+    if nparams != 23 :
+        raise ValueError('Do not know how to convert %s regions to 37 extended transcom regions'%(nparams,))
+    M           = xform
+
+    if not cov:
+        return dot(array(data).squeeze(),M)
+    else:
+        return transpose(dot(transpose(M),dot(array(data).squeeze(),M)),(1,0,2))
+
+def cov2corr(A):
+    b=1./sqrt(A.diagonal())
+    return A*dot(b[:,newaxis],b[newaxis,:])
+
+def map_to_tc(data):
+    """ function projects 1x1 degree map onto TransCom regions by adding gridboxes over larger areas """
+    from hdf2field import Sds2field 
+    import cPickle
+    import os
+    from plottools import rebin
+
+    transcommapfile= 'tc_land11_oif30.hdf'
+    transcomconversionfile= 'map_to_tc.pickle'
+    try:
+        regionselect = cPickle.load(open(transcomconversionfile,'rb'))
+    except:
+        # read map from NetCDF
+        print '[in map_to_tc() in tctools.py:] '+ \
+              'Creating conversion map and pickle file for future quick use, patience please...'
+        map = Sds2field(transcommapfile,'tc_region')
+
+        # create dictionary for region <-> map conversions based on 1x1 map
+        
+        regs={}
+        nregions=map.max()
+        for r in arange(1,nregions+1):
+             sel=(map.flat == r).nonzero()
+             if len(sel[0])>0: regs[r]=sel
+        regionselect=regs
+        dummy = cPickle.dump(regionselect,open(transcomconversionfile,'wb'),-1)
+    
+    result=zeros(len(regionselect.keys()),float)
+    for k,v in regionselect.iteritems():
+        result[k-1]=data.ravel().take(v).sum()
+    return result
+
+def LookUpName(rundat,reg=33,tc=False,eco=False, olson=False, longnames=False):
+    """ return name of region number reg """
+
+    if longnames:
+        econames=olsonnams
+    else :
+        econames=olsonshort
+
+    if tc:
+        return (transnams[reg-1],)
+    elif eco:
+        if reg > rundat.npparameters:
+            raise IOError,'Region number exceeds definitions'
+        elif reg > rundat.n_land and reg != rundat.nparameters: 
+            ret=('Ocean', oifnams[reg-rundat.n_land-1])
+        elif reg > 209 and reg <= rundat.n_land:
+            ret=('Europe', econames[(reg-1)%19]+"_East")
+        elif reg == rundat.nparameters:
+            ret=(transnams[-1])
+        else:
+            ret=(transnams[(reg-1)/19], econames[(reg-1)%19])
+        return ret
+    elif olson:
+        return (econames[(reg-1)%19],)
+
+if __name__ == '__main__':
+    print transnams
+    print transshort
+    print ext_transnams
+    print ext_transshort
+    print olsonnams
+    print olsonshort
+    print ext_transcomps
+    
+

refactored/da/baseclasses/dasystem.py

+#!/usr/bin/env python
+# control.py
+
+"""
+.. module:: dasystem
+.. moduleauthor:: Wouter Peters 
+
+Revision History:
+File created on 26 Aug 2010.
+
+The DaSystem class is found in the module :mod:`dasystem`, or in a specific implementation under the da/ source tree. It is derived from the standard python :class:`dictionary` object. 
+
+It describes the details of the data assimilation system used (i.e., CarbonTracker, or CT Methane, or ....) ::
+
+    datadir         : /Volumes/Storage/CO2/carbontracker/input/ct08/   ! The directory where input data is found
+    obs.input.dir   : ${datadir}/obsnc/with_fillvalue                  ! the observation input dir
+    obs.input.fname : obs_forecast.nc                                  ! the observation input file
+    ocn.covariance  : ${datadir}/oif_p3_era40.dpco2.2000.01.hdf        ! the ocean flux covariance file
+    bio.covariance  : ${datadir}/covariance_bio_olson19.nc             ! the biosphere flux covariance file
+    deltaco2.prefix : oif_p3_era40.dpco2                               ! the type of ocean product used
+    regtype         : olson19_oif30                                    ! the ecoregion definitions
+    nparameters     : 240                                              ! the number of parameters to solve for
+    random.seed     : 4385                                             ! the random seed for the first cycle
+    regionsfile     : transcom_olson19_oif30.hdf                       ! the ecoregion defintion mask file
+
+    ! Info on the sites file used
+
+    obs.sites.rc        : ${datadir}/sites_and_weights_co2.ct10.rc     ! the weights in the covariance matric of each obs
+
+The full baseclass description:
+
+.. autoclass:: da.baseclasses.dasystem.DaSystem
+   :members:
+
+"""
+
+import os
+import sys
+import logging
+import datetime
+
+################### Begin Class DaSystem ###################
+
+class DaSystem(dict):
+    """ 
+    Information on the data assimilation system used. This is normally an rc-file with settings.
+    """
+
+    def __init__(self,rcfilename):
+        """
+        Initialization occurs from passed rc-file name, items in the rc-file will be added
+        to the dictionary
+        """
+
+        self.Identifier     = 'CarbonTracker CO2'    # the identifier gives the platform name
+
+        self.LoadRc(rcfilename)
+
+        msg  = 'Data Assimilation System initialized: %s'%self.Identifier ; logging.debug(msg)
+
+
+    def __str__(self):
+        """
+        String representation of a DaInfo object
+        """
+
+        msg = "==============================================================="    ; print msg
+        msg = "DA System Info rc-file is %s" % self.RcFileName                                ; print msg
+        msg = "==============================================================="    ; print msg
+
+        return ""
+
+
+    def LoadRc(self,RcFileName):
+        """ 
+        This method loads a DA System Info rc-file with settings for this simulation 
+        """
+        import da.tools.rc as rc 
+
+        for k,v in rc.read(RcFileName).iteritems():
+            self[k]             = v
+        self.RcFileName     = RcFileName
+        self.DaRcLoaded    = True
+
+        msg                 = 'DA System Info rc-file (%s) loaded successfully'%self.RcFileName ; logging.debug(msg)
+
+        return True
+
+    def Initialize(self ):
+        """ 
+        Initialize the object. 
+        """
+
+
+    def Validate(self ):
+        """ 
+        Validate the contents of the rc-file given a dictionary of required keys
+        """
+
+        needed_rc_items={}
+
+        for k,v in self.iteritems():
+            if v == 'True' : self[k] = True
+            if v == 'False': self[k] = False
+
+        for key in needed_rc_items:
+
+            if not self.has_key(key):
+                status,msg = ( False,'Missing a required value in rc-file : %s' % key)
+                logging.error(msg)
+                raise IOError,msg
+
+        status,msg = ( True,'DA System Info settings have been validated succesfully' )  ; logging.debug(msg)
+
+        return None
+################### End Class DaSystem ###################
+
+
+if __name__ == "__main__":
+    pass

refactored/da/baseclasses/obs.py

+#!/usr/bin/env python
+# obs.py
+
+"""
+.. module:: obs
+.. moduleauthor:: Wouter Peters 
+
+Revision History:
+File created on 28 Jul 2010.
+
+.. autoclass:: da.baseclasses.obs.Observation 
+   :members: Initialize, Validate, AddObs, AddSimulations, AddModelDataMismatch, WriteSampleInfo  
+
+.. autoclass:: da.baseclasses.obs.ObservationList 
+   :members: __init__
+
+"""
+import os
+import sys
+import logging
+import datetime
+from numpy import array, ndarray
+
+identifier = 'Observation baseclass'
+version    = '0.0'
+
+################### Begin Class Observation ###################
+
+class Observation(object):
+    """ 
+    The baseclass Observation is a generic object that provides a number of methods required for any type of observations used in 
+    a data assimilation system. These methods are called from the CarbonTracker pipeline. 
+
+    .. note:: Most of the actual functionality will need to be provided through a derived Observations class with the methods 
+              below overwritten. Writing your own derived class for Observations is one of the first tasks you'll likely 
+              perform when extending or modifying the CarbonTracker Data Assimilation Shell.
+
+    Upon initialization of the class, an object is created that holds no actual data, but has a placeholder attribute `self.Data` 
+    which is an empty list of type :class:`~da.baseclasses.obs.ObservationList`. An ObservationList object is created when the 
+    method :meth:`~da.baseclasses.obs.Observation.AddObs` is invoked in the pipeline. 
+
+    From the list of observations, a file is written  by method 
+    :meth:`~da.baseclasses.obs.Observation.WriteSampleInfo`
+    with the sample info needed by the 
+    :class:`~da.baseclasses.observationoperator.ObservationOperator` object. The values returned after sampling 
+    are finally added by :meth:`~da.baseclasses.obs.Observation.AddSimulations`
+
+    """ 
+
+    def __init__(self, DaCycle = None):
+        """
+        create an object with an identifier, version, and an empty ObservationList
+        """
+        self.Identifier = self.getid()
+        self.Version    = self.getversion()
+        self.Data       = ObservationList([])  # Initialize with an empty list of obs
+
+        # The following code allows the object to be initialized with a DaCycle object already present. Otherwise, it can
+        # be added at a later moment.
+
+        if DaCycle != None:
+            self.DaCycle = DaCycle
+        else:
+            self.DaCycle = {}
+
+        msg                 = 'Observation object initialized: %s'%self.Identifier ; logging.info(msg)
+
+    def getid(self):
+        """
+        Return the identifier of the Observation Object, used for logging purposed
+        """
+        return identifier
+
+    def getversion(self):
+        """
+        Return the version of the Observation Object, used for logging purposed
+        """
+        return identifier
+
+    def __str__(self):
+        """ Prints a list of Observation objects"""
+        return "This is a %s object, version %s"%(self.Identifier,self.Version)
+
+    def Initialize(self):
+        """ Perform all steps needed to start working with observational data, this can include moving data, concatenating files,
+            selecting datasets, etc.
+        """
+
+    def Validate(self):
+        """ Make sure that data needed for the ObservationOperator (such as observation input lists, or parameter files)
+            are present.
+        """
+
+    def AddObs(self):
+        """ 
+        Add actual observation data to the Observation object. This is in a form of an 
+        :class:`~da.baseclasses.obs.ObservationList` that is contained in self.Data. The 
+        list has as only requirement that it can return the observed+simulated values 
+        through the method :meth:`~da.baseclasses.obs.ObservationList.getvalues`
+
+        """
+
+    def AddSimulations(self):
+        """ Add the simulation data to the Observation object. 
+        """
+
+    def AddModelDataMismatch(self ):
+        """ 
+            Get the model-data mismatch values for this cycle.
+
+        """
+
+    def WriteSampleInfo(self ):
+        """ 
+            Write the information needed by the observation operator to a file. Return the filename that was written for later use
+        """
+
+
+
+################### End Class Observations ###################
+
+################### Begin Class ObservationList ###################
+
+class ObservationList(list):
+    """ This is a special type of list that holds observed sample objects. It has methods to extract data from such a list easily """
+
+    def getvalues(self,name,constructor=array):
+
+        result = constructor([getattr(o,name) for o in self])
+        if isinstance(result,ndarray): 
+            return result.squeeze()
+        else:
+            return result
+
+################### End Class ObservationList ###################
+