diff --git a/Makefile b/Makefile deleted file mode 100644 index a24e59125ac6554cae9312d8ea7f8f7bba2c503a..0000000000000000000000000000000000000000 --- a/Makefile +++ /dev/null @@ -1,29 +0,0 @@ -.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 - /bin/rm -f jb.[0-9]*.rc - /bin/rm -f jb.[0-9]*.log - /bin/rm -f *.pyc - /bin/rm -f */*/jb.[0-9]*.rc - /bin/rm -f */*/jb.[0-9]*.jb - /bin/rm -f */*/jb.[0-9]*.log - /bin/rm -f */*/*.pyc - /bin/rm -f das.[0-9]*.log - /bin/rm -f *.o* - /bin/rm -f *.po* - /bin/rm -f *.e* - /bin/rm -f *.pe* - - - diff --git a/carbontracker.rc b/carbontracker.rc deleted file mode 100644 index 87bb91b6005ccac732049b5e338f2217005fd210..0000000000000000000000000000000000000000 --- a/carbontracker.rc +++ /dev/null @@ -1,15 +0,0 @@ -!!! Info for the CarbonTracker data assimilation system - -!datadir : /lfs0/projects/co2/input/ -datadir : /data/CO2/carbontracker/ct08/ -obs.input.dir : ${datadir}/obsnc/with_fillvalue/ -obs.input.fname : obs_final.nc -ocn.covariance : ${datadir}/oif_p3_era40.dpco2.2000.01.hdf -bio.covariance : ${datadir}/covariance_bio_olson19.nc -deltaco2.prefix : oif_p3_era40.dpco2 -regtype : olson19_oif30 -nparameters : 240 -random.seed : 4385 -regionsfile : transcom_olson19_oif30.hdf - - diff --git a/carbontrackerjet.rc b/carbontrackerjet.rc deleted file mode 100644 index 97e49b83daedc0974f644655e1be32007f77cd07..0000000000000000000000000000000000000000 --- a/carbontrackerjet.rc +++ /dev/null @@ -1,13 +0,0 @@ -!!! Info for the CarbonTracker data assimilation system - -datadir : /lfs0/projects/co2/input/ct_new_2010 -obs.input.dir : ${datadir}/obsnc/with_fillvalue -obs.input.fname : obs_final.nc -ocn.covariance : ${datadir}/oif_p3_era40.dpco2.2000.01.hdf -bio.covariance : ${datadir}/covariance_bio_olson19.nc -deltaco2.prefix : oif_p3_era40.dpco2 -regtype : olson19_oif30 -nparameters : 240 -random.seed : 4385 -regionsfile : transcom_olson19_oif30.hdf - diff --git a/da.rc b/da.rc deleted file mode 100644 index 504ed648aba777fce100bef2ee76ee783ba569af..0000000000000000000000000000000000000000 --- a/da.rc +++ /dev/null @@ -1,19 +0,0 @@ -! Info on the data assimilation cycle - -time.restart : False -time.start : 2000-01-01 00:00:00 -time.finish : 2000-01-08 00:00:00 -time.cycle : 7 -time.nlag : 2 -dir.da_run : ${HOME}/tmp/test_da - -! Info on the DA system used - -da.system : CarbonTracker -da.system.rc : carbontracker.rc - -! Info on the forward model to be used - -forecast.model : TM5 -forecast.model.rc : ${HOME}/Modeling/TM5/ct_new_intel.rc -forecast.nmembers : 2 diff --git a/da/__init__.py b/da/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/da/baseclasses/__init__.py b/da/baseclasses/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/da/baseclasses/dasystem.py b/da/baseclasses/dasystem.py deleted file mode 100755 index 173556dad6a08176729fb754592b20bc25bd8ca6..0000000000000000000000000000000000000000 --- a/da/baseclasses/dasystem.py +++ /dev/null @@ -1,87 +0,0 @@ -#!/usr/bin/env python -# control.py - -""" -Author : peters - -Revision History: -File created on 26 Aug 2010. - -""" - -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 - """ - - self.LoadRc(rcfilename) - - 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.info(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 diff --git a/da/baseclasses/obs.py b/da/baseclasses/obs.py deleted file mode 100755 index 0b9af32af0b9023a7073fab0a953522b5c79f7d9..0000000000000000000000000000000000000000 --- a/da/baseclasses/obs.py +++ /dev/null @@ -1,90 +0,0 @@ -#!/usr/bin/env python -# obs.py - -""" -Author : peters - -Revision History: -File created on 28 Jul 2010. - -""" -import os -import sys -import logging -import datetime - -identifier = 'Observation baseclass' -version = '0.0' - -################### Begin Class Observation ################### - -class Observation(object): - """ an object that holds data + methods and attributes needed to manipulate observations for a DA cycle """ - - def __init__(self): - self.Identifier = self.getid() - self.Version = self.getversion() - self.Data = ObservationList([]) # Initialize with an empty list of obs - - msg = 'Observation object initialized: %s'%self.Identifier ; logging.info(msg) - - def getid(self): - return identifier - - def getversion(self): - 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 the observation data to the Observation object. This is in a form of a list that goes under the name Data. The - list has as only requirement that it can return the observed+simulated values through a method "getvalues" - """ - - def AddSimulations(self): - """ Add the simulation data to the Observation object. - """ - - def AddModelDataMismatch(self, DaCycle): - """ - Get the model-data mismatch values for this cycle. - - """ - - def WriteSampleInfo(self, DaCycle): - """ - 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 """ - from numpy import array, ndarray - - def getvalues(self,name,constructor=array): - from numpy import ndarray - result = constructor([getattr(o,name) for o in self]) - if isinstance(result,ndarray): - return result.squeeze() - else: - return result - -################### End Class MixingRatioList ################### - diff --git a/da/baseclasses/observationoperator.py b/da/baseclasses/observationoperator.py deleted file mode 100755 index 6fbcc3171f1d41d45fa72f9987a6a5498ec8b990..0000000000000000000000000000000000000000 --- a/da/baseclasses/observationoperator.py +++ /dev/null @@ -1,67 +0,0 @@ -#!/usr/bin/env python -# model.py - -""" -Author : peters - -Revision History: -File created on 30 Aug 2010. - -""" - -import os -import sys -import logging -import datetime - -identifier = 'GeneralObservationOperator' -version = '0.0' - -################### Begin Class ObservationOperator ################### -class ObservationOperator(object): - """ - - This is a class that defines an ObervationOperator. This object is used to control the sampling of - a statevector in the ensemble Kalman filter framework. The methods of this class specify which (external) code - is called to perform the sampling, and which files should be read for input and are written for output. - - The baseclasses consist mainly of empty mehtods that require an application specific application - - """ - - def __init__(self): - """ The instance of an ObservationOperator is application dependent """ - self.Identifier = self.getid() - self.Version = self.getversion() - - msg = 'Observation Operator object initialized: %s'%self.Identifier ; logging.info(msg) - - def getid(self): - return identifier - - def getversion(self): - return version - - - - def __str__(self): - return "This is a %s object, version %s"%(self.Identifier,self.Version) - - def Initialize(self,DaCycle): - """ Perform all steps necessary to start the observation operator through a simple Run() call """ - - def ValidateInput(self,DaCycle): - """ Make sure that data needed for the ObservationOperator (such as observation input lists, or parameter files) - are present. - """ - def SaveData(self): - """ Write the data that is needed for a restart or recovery of the Observation Operator to the save directory """ - - - -################### End Class ObservationOperator ################### - - - -if __name__ == "__main__": - pass diff --git a/da/baseclasses/optimizer.py b/da/baseclasses/optimizer.py deleted file mode 100755 index 639d2c7b1578e2675c0f904111b53a36db7af354..0000000000000000000000000000000000000000 --- a/da/baseclasses/optimizer.py +++ /dev/null @@ -1,259 +0,0 @@ -#!/usr/bin/env python -# optimizer.py - -""" -Author : peters - -Revision History: -File created on 28 Jul 2010. - -""" - -import os -import sys -import logging -import datetime - -identifier = 'Optimizer baseclass' -version = '0.0' - -################### Begin Class Optimizer ################### - -class Optimizer(object): - """ - This creates an instance of an optimization object. It handles the minimum least squares optimization - of the state vector given a set of sample objects. Two routines will be implemented: one where the optimization - is sequential and one where it is the equivalent matrix solution. The choice can be made based on considerations of speed - and efficiency. - """ - - def __init__(self): - self.Identifier = self.getid() - self.Version = self.getversion() - - msg = 'Optimizer object initialized: %s'%self.Identifier ; logging.info(msg) - - def getid(self): - return identifier - - def getversion(self): - return version - - def Initialize(self, dims): - - self.nlag = dims[0] - self.nmembers = dims[1] - self.nparams = dims[2] - self.nobs = dims[3] - self.localization = False - self.localization = False - - self.CreateMatrices() - - return None - - def CreateMatrices(self): - """ Create Matrix space needed in optimization routine """ - import numpy as np - - # mean state [X] - self.x = np.zeros( (self.nlag*self.nparams,), float) - # deviations from mean state [X'] - self.X_prime = np.zeros( (self.nlag*self.nparams,self.nmembers,), float) - # mean state, transported to observation space [ H(X) ] - self.Hx = np.zeros( (self.nobs,), float) - # deviations from mean state, transported to observation space [ H(X') ] - self.HX_prime = np.zeros( (self.nobs,self.nmembers), float) - # observations - self.obs = np.zeros( (self.nobs,), float) - # covariance of observations - self.R = np.zeros( (self.nobs,self.nobs,), float) - # Total covariance of fluxes and obs in units of obs [H P H^t + R] - self.HPHR = np.zeros( (self.nobs,self.nobs,), float) - # Kalman Gain matrix - self.KG = np.zeros( (self.nlag*self.nparams,self.nobs,), float) - # flags of observations - self.flags = np.zeros( (self.nobs,), int) - - def StateToMatrix(self,StateVector): - import numpy as np - - for n in range(self.nlag): - - members = StateVector.EnsembleMembers[n] - self.x[n*self.nparams:(n+1)*self.nparams] = members[0].ParameterValues[n] - self.X_prime[n*self.nparams:(n+1)*self.nparams,:] = np.transpose(np.array([m.ParameterValues for m in members])) - - self.X_prime = self.X_prime - self.x[:,np.newaxis] # make into a deviation matrix - - self.obs[:] = StateVector.EnsembleMembers[-1][0].ModelSample.Data.getvalues('obs') - self.Hx[:] = StateVector.EnsembleMembers[-1][0].ModelSample.Data.getvalues('simulated') - - for m,mem in enumerate(StateVector.EnsembleMembers[-1]): - - self.HX_prime[:,m] = mem.ModelSample.Data.getvalues('simulated') - - self.HX_prime = self.HX_prime - self.Hx[:,np.newaxis] # make a deviation matrix - - self.R[:,:] = np.identity(self.nobs) - - return None - - def MatrixToState(self,StateVector): - import numpy as np - - for n in range(self.nlag): - - members = StateVector.EnsembleMembers[n] - for m,mem in enumerate(members): - members[m].ParameterValues[:] = self.X_prime[n*self.nparams:(n+1)*self.nparams,m] + self.x[n*self.nparams:(n+1)*self.nparams] - - return None - - def SerialMinimumLeastSquares(self): - """ Make minimum least squares solution by looping over obs""" - import numpy as np - import numpy .linalg as la - - tvalue=1.97591 - for n in range(self.nobs): - - if self.flags[n] == 1: continue - - PHt = 1./(self.nmembers-1)*np.dot(self.X_prime,self.HX_prime[n,:]) - self.HPHR[n,n] = 1./(self.nmembers-1)*(self.HX_prime[n,:]*self.HX_prime[n,:]).sum()+self.R[n,n] - - self.KG[:,n] = PHt/self.HPHR[n,n] - - dummy = self.Localize(n) - - alpha = np.double(1.0)/(np.double(1.0)+np.sqrt( (self.R[n,n])/self.HPHR[n,n] ) ) - - res = self.obs[n]-self.Hx[n] - - self.x[:] = self.x + self.KG[:,n]*res - - for r in range(self.nmembers): - self.X_prime[:,r] = self.X_prime[:,r]-alpha*self.KG[:,n]*(self.HX_prime[n,r]) - -#WP !!!! Very important to first do all obervations from n=1 through the end, and only then update 1,...,n. The current observation -#WP should always be updated last because it features in the loop of the adjustments !!!! - - for m in range(n+1,self.nobs): - res = self.obs[n]-self.Hx[n] - fac = 1.0/(self.nmembers-1)*(self.HX_prime[n,:]*self.HX_prime[m,:]).sum()/self.HPHR[n,n] - self.Hx[m] = self.Hx[m] + fac * res - self.HX_prime[m,:] = self.HX_prime[m,:] - alpha*fac*self.HX_prime[n,:] - - for m in range(1,n+1): - res = self.obs[n]-self.Hx[n] - fac = 1.0/(self.nmembers-1)*(self.HX_prime[n,:]*self.HX_prime[m,:]).sum()/self.HPHR[n,n] - self.Hx[m] = self.Hx[m] + fac * res - self.HX_prime[m,:] = self.HX_prime[m,:] - alpha*fac*self.HX_prime[n,:] - - - def BulkMinimumLeastSquares(self): - """ Make minimum least squares solution by solving matrix equations""" - import numpy as np - import numpy.linalg as la - - # Create full solution, first calculate the mean of the posterior analysis - - HPH = np.dot(self.HX_prime,np.transpose(self.HX_prime))/(self.nmembers-1) # HPH = 1/N * HX' * (HX')^T - self.HPHR[:,:] = HPH+self.R # HPHR = HPH + R - HPb = np.dot(self.X_prime,np.transpose(self.HX_prime))/(self.nmembers-1) # HP = 1/N X' * (HX')^T - self.KG[:,:] = np.dot(HPb,la.inv(self.HPHR)) # K = HP/(HPH+R) - - for n in range(self.nobs): - dummy = self.Localize(n) - - self.x[:] = self.x + np.dot(self.KG,self.obs-self.Hx) # xa = xp + K (y-Hx) - - # And next make the updated ensemble deviations. Note that we calculate P by using the full equation (10) at once, and - # not in a serial update fashion as described in Whitaker and Hamill. - # For the current problem with limited N_obs this is easier, or at least more straightforward to do. - - I = np.identity(self.nlag*self.nparams) - sHPHR = la.cholesky(self.HPHR) # square root of HPH+R - part1 = np.dot(HPb,np.transpose(la.inv(sHPHR))) # HP(sqrt(HPH+R))^-1 - part2 = la.inv(sHPHR+np.sqrt(self.R)) # (sqrt(HPH+R)+sqrt(R))^-1 - Kw = np.dot(part1,part2) # K~ - self.X_prime[:,:] = np.dot(I,self.X_prime)-np.dot(Kw,self.HX_prime) # HX' = I - K~ * HX' - - P_opt = np.dot(self.X_prime,np.transpose(self.X_prime))/(self.nmembers-1) - - # Now do the adjustments of the modeled mixing ratios using the linearized ensemble. These are not strictly needed but can be used - # for diagnosis. - - part3 = np.dot(HPH,np.transpose(la.inv(sHPHR))) # HPH(sqrt(HPH+R))^-1 - Kw = np.dot(part3,part2) # K~ - self.Hx[:] = self.Hx + np.dot(np.dot(HPH,la.inv(self.HPHR)),self.obs-self.Hx) # Hx = Hx+ HPH/HPH+R (y-Hx) - self.HX_prime[:,:] = self.HX_prime-np.dot(Kw,self.HX_prime) # HX' = HX'- K~ * HX' - - msg = 'Minimum Least Squares solution was calculated, returning' ; logging.info(msg) - - return None - - def SetLocalization(self): - """ determine which localization to use """ - - self.localization = True - self.localizetype = "None" - - msg = "Current localization option is set to %s"%self.localizetype ; logging.info(msg) - - def Localize(self,n): - """ localize the Kalman Gain matrix """ - import numpy as np - - if not self.localization: return - - return - -################### End Class Optimizer ################### - - - -if __name__ == "__main__": - - sys.path.append('../../') - - import os - import sys - from da.tools.general import StartLogger - from da.tools.initexit import CycleControl - from da.ct.statevector import CtStateVector, PrepareState - from da.ct.obs import CtObservations - import numpy as np - import datetime - import da.tools.rc as rc - - opts = ['-v'] - args = {'rc':'../../da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} - - StartLogger() - DaCycle = CycleControl(opts,args) - - DaCycle.Initialize() - print DaCycle - - StateVector = PrepareState(DaCycle) - - samples = CtObservations(DaCycle.DaSystem,datetime.datetime(2005,3,5)) - dummy = samples.AddObs() - dummy = samples.AddSimulations('/Users/peters/tmp/test_da/output/20050305/samples.000.nc') - - - nobs = len(samples.Data) - dims = ( int(DaCycle['time.nlag']), - int(DaCycle['forecast.nmembers']), - int(DaCycle.DaSystem['nparameters']), - nobs, ) - - opt = CtOptimizer(dims) - - opt.StateToMatrix(StateVector) - - opt.MinimumLeastSquares() - - opt.MatrixToState(StateVector) diff --git a/da/baseclasses/platform.py b/da/baseclasses/platform.py deleted file mode 100755 index 5bcea69a0dcc1f674c61b39b100f2711a5d9dee0..0000000000000000000000000000000000000000 --- a/da/baseclasses/platform.py +++ /dev/null @@ -1,110 +0,0 @@ -#!/usr/bin/env python -# jobcontrol.py - -""" -Author : peters - -Revision History: -File created on 06 Sep 2010. - -""" - -import sys -import os -import logging -import subprocess - -std_joboptions={'jobname':'test','jobaccount':'co2','jobnodes':'nserial 1','jobshell':'/bin/sh','depends':'','jobtime':'00:30:00'} - -class PlatForm(object): - """ This specifies platform dependent options under generic object calls. A platform object is used to control and submit jobs""" - - def __init__(self): - self.Identifier = 'iPad' # the identifier gives the plaform name - self.Version = '1.0' # the platform version used - - msg1 = '%s object initialized'%self.Identifier ; logging.debug(msg1) - msg2 = '%s version: %s'%(self.Identifier,self.Version) ; logging.debug(msg2) - - - def __str__(self): - - return None - - def GetJobTemplate(self,joboptions={},block=False): - """ Return the job template for a given computing system, and fill it with options from the dictionary provided as argument""" - - template = """## \n"""+ \ - """## This is a set of dummy names, to be replaced by values from the dictionary \n"""+ \ - """## Please make your own platform specific template with your own keys and place it in a subfolder of the da package.\n """+ \ - """## \n"""+ \ - """ \n"""+ \ - """#$ jobname \n"""+ \ - """#$ jobaccount \n"""+ \ - """#$ jobnodes \n"""+ \ - """#$ jobtime \n"""+ \ - """#$ jobshell \n """ - - if 'depends' in joboptions: - template += """#$ -hold_jid depends \n""" - - # First replace from passed dictionary - for k,v in joboptions.iteritems(): - while k in template: - template = template.replace(k,v) - - # Fill remaining values with std_options - for k,v in std_joboptions.iteritems(): - while k in template: - template = template.replace(k,v) - - return template - - def GetMyID(self): - - return os.getpid() - - def WriteJob(self,DaCycle,template, jobid): - """ This method writes a jobfile to the exec dir""" - # - # Done, write jobfile - # - targetdir = os.path.join(DaCycle['dir.exec']) - jobfile = os.path.join(targetdir,'jb.%s.jb'%jobid) - f = open(jobfile,'w') - dummy = f.write(template) - dummy = f.close() - dummy = os.chmod(jobfile,477) - - msg = "A job file was created (%s)"%jobfile ; logging.debug(msg) - - return jobfile - - - def SubmitJob(self,jobfile): - """ This method submits a jobfile to the queue, and returns the queue ID """ - - cmd = ['sh',jobfile] - msg = "A new task will be started (%s)\n\n\n"%cmd ; logging.debug(msg) - process = subprocess.Popen(cmd) - jobid = subprocess.Popen(cmd).pid - dummy = subprocess.Popen(cmd).wait() - - return 0 - - def KillJob(self,jobid): - """ This method kills a running job """ - - return None - - def StatJob(self,jobid): - """ This method gets the status of a running job """ - import subprocess - - output = subprocess.Popen(['qstat',jobid], stdout=subprocess.PIPE).communicate()[0] ; logging.info(output) - - return output - - -if __name__ == "__main__": - pass diff --git a/da/baseclasses/statevector.py b/da/baseclasses/statevector.py deleted file mode 100755 index cd28062e18c999a8203aa41dde77aaaeb8a7680e..0000000000000000000000000000000000000000 --- a/da/baseclasses/statevector.py +++ /dev/null @@ -1,292 +0,0 @@ -#!/usr/bin/env python -# ct_statevector_tools.py - -""" -Author : peters - -Revision History: -File created on 28 Jul 2010. - -""" - -import os -import sys -import logging -import datetime - -identifier = 'Baseclass Statevector ' -version = '0.0' - -################### Begin Class EnsembleMember ################### - -class EnsembleMember(object): - """ - An ensemble member. This consists of - * a member number - * parameter values - * a ModelSample object to hold sampled values for this member - """ - - def __init__(self, membernumber): - self.membernumber = membernumber # the member number - self.ParameterValues = None # Parameter values of this member - self.ModelSample = None # Model Sampled Parameter values of this member - - def __str__(self): - return "%03d"%self.membernumber - - def WriteToFile(self,DaCycle): - """ Write the information needed by an external model to a netcdf file for future use """ - from da.tools.io import CT_CDF, std_savedict - import numpy as np - - outdir = DaCycle['dir.input'] - filename = os.path.join(outdir,'parameters.%03d.nc'% self.membernumber) - f = CT_CDF(filename,'create') - dimparams = f.AddParamsDim(len(self.ParameterValues)) - - data = self.ParameterValues - - savedict = std_savedict.copy() - savedict['name'] = "parametervalues" - savedict['long_name'] = "parameter_values_for_member_%d"%self.membernumber - savedict['units'] = "unitless" - savedict['dims'] = dimparams - savedict['values'] = data.tolist() - savedict['comment'] = 'These are parameter values to use for member %d'%self.membernumber - dummy = f.AddData(savedict) - - dummy = self.AddCustomFields(f,DaCycle) - - dummy = f.close() - - msg = 'Successfully wrote data from ensemble member %d to file (%s) ' % (self.membernumber,filename,) ; logging.info(msg) - - def AddCustomFields(self,filehandle,DaCycle): - """ Placeholder to add more custom fields to the output of the ensemble member""" - - return None - -################### End Class EnsembleMember ################### - -################### Begin Class StateVector ################### - -import numpy as np -class StateVector(object): - """ an object that holds data + methods and attributes needed to manipulate state vector values """ - - def __init__(self): - self.Identifier = self.getid() - self.Version = self.getversion() - - msg = 'Statevector object initialized: %s'%self.Identifier ; logging.info(msg) - - def getid(self): - return identifier - - def getversion(self): - return version - - - - def Initialize(self,dims): - - self.nlag = dims[0] - self.nmembers = dims[1] - self.nparams = dims[2] - self.isOptimized = False - - # These list objects hold the data for each time step of lag in the system. Note that the ensembles for each time step consist - # of lists of EnsembleMember objects, we define member 0 as the mean of the distribution and n=1,...,nmembers as the spread. - - self.EnsembleMembers = range(self.nlag) - - for n in range(self.nlag): - self.EnsembleMembers[n] = [] - - def __str__(self): - return "This is a base class derived state vector object" - - def MakeNewEnsemble(self,lag, covariancematrix = None): - """ Make a new ensemble, the attribute lag refers to the position in the state vector. - Note that lag=1 means an index of 0 in python, hence the notation lag-1 in the indexing below. - The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag] - - The optional DaCycle object to be passed holds info on the data assimilation system settings needed - for instance to read covariance data from file. This is not part of the base class object here, hence, this - mehtod will almost always be overwritten in the derived class for the specific application. - - The code below is just an example of what could be used - """ - - if covariancematrix == None: - covariancematrix=np.identity(self.nparams) - - # Make a cholesky decomposition of the covariance matrix - - U,s,Vh = np.linalg.svd(covariancematrix) - dof = np.sum(s)**2/sum(s**2) - C = np.linalg.cholesky(covariancematrix) - - msg = 'Cholesky decomposition has succeeded offline' ; logging.debug(msg) - msg = 'Appr. degrees of freedom in covariance matrix is %s'%(int(dof)) ; logging.info(msg) - - - # Create mean values - - NewMean = np.ones(self.nparams,float) - - # Create the first ensemble member with a deviation of 0.0 and add to list - - NewMember = self.GetNewMember(0) - NewMember.ParameterValues = np.zeros((self.nparams),float) + NewMean - dummy = self.EnsembleMembers[lag-1].append(NewMember) - - # Create members 1:nmembers and add to EnsembleMembers list - - for member in range(1,self.nmembers): - randstate = np.random.get_state() - rands = np.random.randn(self.nparams) - - NewMember = self.GetNewMember(member) - NewMember.ParameterValues = np.dot(C,rands) + NewMean - dummy = self.EnsembleMembers[lag-1].append(NewMember) - - msg = '%d new ensemble members were added to the state vector # %d'%(self.nmembers,lag) ; logging.debug(msg) - - def GetNewMember(self,memberno): - """ Return an ensemblemember object """ - - return EnsembleMember(memberno) - - - def Propagate(self): - """ - - Propagate the parameter values in the StateVector to the next cycle. This means a shift by one cycle step for all states that will - be optimized once more, and the creation of a new ensemble for the time step that just comes in for the first time (step=nlag). - In the future, this routine can incorporate a formal propagation of the statevector. - - """ - - # Remove State Vector n=1 by simply "popping" it from the list and appending a new empty list at the front. This empty list will - # hold the new ensemble for the new cycle - - dummy = self.EnsembleMembers.pop(0) - - dummy = self.EnsembleMembers.append([]) - - # And now create a new week of mean + members for n=nlag - - dummy = self.MakeNewEnsemble(self.nlag) - - msg = 'The state vector has been propagated by one cycle ' ; logging.info(msg) - - return None - - def WriteToFile(self,filename): - """ Write the StateVector object to a netcdf file for future use """ - from da.tools.io import CT_CDF - import numpy as np - - f = CT_CDF(filename,'create') - dimparams = f.AddParamsDim(self.nparams) - dimmembers = f.AddMembersDim(self.nmembers) - dimlag = f.AddLagDim(self.nlag,unlimited=True) - - for n in range(self.nlag): - - members = self.EnsembleMembers[n] - MeanState = members[0].ParameterValues - - savedict = f.StandardVar(varname='meanstate') - savedict['dims'] = dimlag+dimparams - savedict['values'] = MeanState - savedict['count'] = n - savedict['comment'] = 'this represents the mean of the ensemble' - dummy = f.AddData(savedict) - - members = self.EnsembleMembers[n] - data = np.array([m.ParameterValues for m in members])- MeanState - - savedict = f.StandardVar(varname='ensemblestate') - savedict['dims'] = dimlag+dimmembers+dimparams - savedict['values'] = data - savedict['count'] = n - savedict['comment'] = 'this represents deviations from the mean of the ensemble' - dummy = f.AddData(savedict) - - dummy = f.close() - - msg = 'Successfully wrote the State Vector to file (%s) ' % (filename,) ; logging.info(msg) - - def ReadFromFile(self,filename): - """ Read the StateVector object from a netcdf file for future use """ - import Nio - import numpy as np - - f = Nio.open_file(filename,'r') - MeanState = f.variables['statevectormean'].get_value() - EnsembleMembers = f.variables['statevectorensemble'].get_value() - dummy = f.close() - - for n in range(self.nlag): - for m in range(self.nmembers): - NewMember = self.GetNewMember(m) - NewMember.ParameterValues = EnsembleMembers[n,m,:] + MeanState[n] # add the mean to the deviations to hold the full parameter values - dummy = self.EnsembleMembers[n].append(NewMember) - - msg = 'Successfully read the State Vector from file (%s) ' % (filename,) ; logging.info(msg) - -################### End Class StateVector ################### - -if __name__ == "__main__": - - sys.path.append('../../') - - import os - import sys - from da.tools.general import StartLogger - from da.tools.initexit import CycleControl - import numpy as np - import da.tools.rc as rc - - opts = ['-v'] - args = {'rc':'../../da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} - - StartLogger() - DaCycle = CycleControl(opts,args) - - dummy = EnsembleMember(0) - print dummy - - DaCycle.Initialize() - print DaCycle - - dims = ( int(DaCycle['time.nlag']), - int(DaCycle['forecast.nmembers']), - int(DaCycle.DaSystem['nparameters']), - ) - - StateVector = StateVector(dims) - - StateVector.MakeNewEnsemble(lag=1) - StateVector.MakeNewEnsemble(lag=2) - - members = StateVector.EnsembleMembers[1] - - members[0].WriteToFile(DaCycle['dir.input']) - - StateVector.Propagate() - - savedir = DaCycle['dir.output'] - filename = os.path.join(savedir,'savestate.nc') - - StateVector.WriteToFile(filename) - - - savedir = DaCycle['dir.output'] - filename = os.path.join(savedir,'savestate.nc') - - StateVector.ReadFromFile(filename) - diff --git a/da/bin/readme_wrapper.txt b/da/bin/readme_wrapper.txt deleted file mode 100644 index f355bfec41b4dadb1fef202ed7276d1992659c4a..0000000000000000000000000000000000000000 --- a/da/bin/readme_wrapper.txt +++ /dev/null @@ -1 +0,0 @@ -mpicc tm5_mpi_wrapper.c -o tm5_mpi_wrapper diff --git a/da/bin/tm5_mpi_wrapper.c b/da/bin/tm5_mpi_wrapper.c deleted file mode 100644 index ae4f5455041863a40e9c8a5c4d5af67ec82f6f73..0000000000000000000000000000000000000000 --- a/da/bin/tm5_mpi_wrapper.c +++ /dev/null @@ -1,61 +0,0 @@ -#include <stdio.h> -#include <errno.h> -#include "mpi.h" - -#define CMDLENGTH 200 - -int safe_system (const char *command); - -int main(int argc, char *argv[]) { - - int ierr; - int myrank; - char cmd[CMDLENGTH]; - - ierr = MPI_Init(&argc, &argv); - if (argc != 2) { - fprintf(stderr, "[tm5_mpi_wrapper] Expecting 1 argument, got %d.\n",argc); - exit(-1); - } - MPI_Comm_rank(MPI_COMM_WORLD, &myrank); - snprintf(cmd,CMDLENGTH,"%s %03d",argv[1],myrank); - //snprintf(cmd,CMDLENGTH,"%s %03d >> tm5.%03d.log",argv[1],myrank,myrank);// - printf( "MPI rank %d about to execute command \"%s\".\n",myrank,cmd ); - ierr = safe_system(cmd); - if(ierr != 0) { - MPI_Abort(MPI_COMM_WORLD,ierr); - exit(ierr); - } - ierr = MPI_Finalize( ); - exit(ierr); -} - - - - - -int safe_system (const char *command) { - int pid, status; - - if (command == 0) - return 1; - pid = fork(); - if (pid == -1) // fork failed - return -1; - if (pid == 0) { // then this is the child - char *argv[4]; - argv[0] = "sh"; - argv[1] = "-c"; - argv[2] = (char*)command; - argv[3] = 0; - execv("/bin/sh", argv); - _exit(127); - } - do { - if (waitpid(pid, &status, 0) == -1) { - if (errno != EINTR) - return -1; - } else - return status; - } while(1); -} diff --git a/da/ct/__init__.py b/da/ct/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/da/ct/dasystem.py b/da/ct/dasystem.py deleted file mode 100755 index 454a695235b3b39ae1fb8ec916389948ac417dee..0000000000000000000000000000000000000000 --- a/da/ct/dasystem.py +++ /dev/null @@ -1,69 +0,0 @@ -#!/usr/bin/env python -# control.py - -""" -Author : peters - -Revision History: -File created on 26 Aug 2010. - -""" - -import os -import sys -import logging -import datetime - -################### Begin Class CtDaSystem ################### - -from da.baseclasses.dasystem import DaSystem - -class CtDaSystem(DaSystem): - """ Information on the data assimilation system used. This is normally an rc-file with settings. - """ - - def Initialize(self): - """ - Initialize the object - """ - import Nio - - mapfile = os.path.join(self['datadir'],self['regionsfile']) - ncf = Nio.open_file(mapfile,'r') - self.regionmap = ncf.variables['budget_region'].get_value() - dummy = ncf.close() - - - def Validate(self): - """ - Validate the contents of the rc-file given a dictionary of required keys - """ - - needed_rc_items = ['obs.input.dir', - 'obs.input.fname', - 'ocn.covariance', - 'nparameters', - 'bio.covariance', - 'deltaco2.prefix', - 'regtype'] - - - 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 CtDaSystem ################### - - -if __name__ == "__main__": - pass diff --git a/da/ct/io.py b/da/ct/io.py deleted file mode 100755 index 2196be26e9b71e47abe8f1b94781a6b91df91bb3..0000000000000000000000000000000000000000 --- a/da/ct/io.py +++ /dev/null @@ -1,217 +0,0 @@ -#!/usr/bin/env python -# io.py - -""" -Author : peters - -Revision History: -File created on 15 Oct 2008. -File modified for CT data assimilation system in July 2010, Wouter Peters - -""" -import standardvariables -import pycdf as CDF -import datetime as dt -from numpy import array -import os - -disclaimer = "This data belongs to the CarbonTracker project" -email = "wouter.peters@wur.nl" -url = "http://carbontracker.wur.nl" -institution = "Wageningen University and Research Center" -source = "CarbonTracker release 2.0" -conventions = "CF-1.1" -historytext = 'Created on '+dt.datetime.now().strftime('%B %d, %Y')+' by %s'%os.environ['USER'] - -std_savedict={'name':'unknown','values':[],'dims':(0,0,),'units':'','long_name':'','_FillValue':float(-999.),'comment':''} - -class CT_CDF(CDF.CDF): - """ function opens a NetCDF/HDF/GRIB file for writing of output""" - def __init__(self,filename, method='read'): - - if method not in ['read','write','create']: - raise ValueError, 'Method %s is not defined for a CarbonTracker NetCDF file object' % method - - if method == 'read': - print 'Reading from file' - super(CDF.CDF,self).__init__(filename, CDF.NC.NOWRITE) - elif method == 'write': - #print 'Adding to existing file' - super(CT_CDF,self).__init__(filename, CDF.NC.WRITE|CDF.NC.CREATE) - self.AddCTHeader() - elif method == 'create': - #print 'Creating new file' - super(CT_CDF,self).__init__(filename, CDF.NC.WRITE|CDF.NC.TRUNC|CDF.NC.CREATE) - self.AddCTHeader() - - - def AddCTHeader(self): - - self.automode() - # - setattr(self,'Institution',institution) - setattr(self,'Contact',email) - setattr(self,'URL',url) - setattr(self,'Source',source) - setattr(self,'Convention',conventions) - setattr(self,'Disclaimer',disclaimer) - setattr(self,'History',historytext) - - def AddParamsDim(self,nparams): - - self.automode() - dimparams=self.def_dim('nparameters',nparams) - - return (dimparams,) - - def AddMembersDim(self,nmembers): - - self.automode() - dimmembers=self.def_dim('nmembers',nmembers) - - return (dimmembers,) - - def AddLagDim(self,nlag,unlimited=True): - - self.automode() - if unlimited: - dimlag =self.def_dim('nlag',CDF.NC.UNLIMITED) - else: - dimlag=self.def_dim('nlag',nlag) - - return (dimlag,) - - def AddObsDim(self,nobs): - - self.automode() - dimobs=self.def_dim('nobs',nobs) - - return (dimobs,) - - def AddLatLonDim(self,istart=0,iend=360,jstart=0,jend=180): - - from numpy import arange, float64 - - if 'latitude' in self.dimensions(): return (self.dim('latitude'),self.dim('longitude'),) # already exists - - lons=-180+arange(360)*1.0+0.5 - lats=-90+arange(180)*1.0+0.5 - # - lats=lats[jstart:jend] - lons=lons[istart:iend] - # - self.automode() - dimlon=self.def_dim('longitude',lons.shape[0]) - dimlat=self.def_dim('latitude',lats.shape[0]) - - savedict=self.StandardVar(varname='latitude') - savedict['values']=lats.tolist() - savedict['actual_range']=(float(lats[0]),float(lats[-1])) - savedict['dims']=(dimlat,) - self.AddData(savedict) - - savedict=self.StandardVar(varname='longitude') - savedict['values']=lons.tolist() - savedict['actual_range']=(float(lons[0]),float(lons[-1])) - savedict['dims']=(dimlon,) - self.AddData(savedict) - - return (dimlat,dimlon,) - - def AddDateDim(self): - - self.automode() - if 'date' in self.dimensions(): return (self.dim('date'),) - return (self.def_dim('date',CDF.NC.UNLIMITED),) - - def AddDateDimFormat(self): - - self.automode() - if 'yyyymmddhhmmss' in self.dimensions(): return (self.dim('yyyymmddhhmmss'),) # already exists - return (self.def_dim('yyyymmddhhmmss',6),) - - def has_date(self,dd): - - if self.inq_unlimlen() > 0: - if dd in self.GetVariable('date').tolist(): - return True - else: - return False - else: - return False - - def GetVariable(self,varname): - """ get variable from ncf file""" - return array(self.var(varname).get()) - - def StandardVar(self,varname): - """ return properties of standard variables """ - import standardvariables - - if varname in standardvariables.standard_variables.keys(): - return standardvariables.standard_variables[varname] - else: - return standardvariables.standard_variables['unknown'] - - def AddData(self,datadict,nsets=1,silent=True): - """ add fields to file, at end of unlimited dimension""" - - existing_vars=self.variables() - - try: - next = datadict['count'] - except: - next=0 - - - if existing_vars.has_key(datadict['name']): - var = self.var(datadict['name']) - var[next:next+nsets]=datadict['values'] - else: - if not silent: print 'Creating new dataset: '+datadict['name'] - - if datadict.has_key('dtype'): - if datadict['dtype'] == 'int': - var = self.def_var(datadict['name'],CDF.NC.INT,datadict['dims']) - elif datadict['dtype'] == 'char': - var = self.def_var(datadict['name'],CDF.NC.CHAR,datadict['dims']) - elif datadict['dtype'] == 'double': - var = self.def_var(datadict['name'],CDF.NC.DOUBLE,datadict['dims']) - else: - var = self.def_var(datadict['name'],CDF.NC.FLOAT,datadict['dims']) - else: - var = self.def_var(datadict['name'],CDF.NC.FLOAT,datadict['dims']) - for k,v in datadict.iteritems(): - if k not in ['name','dims','values','_FillValue','count']: - setattr(var,k,v) - if var.isrecord(): - var[next:next+nsets]=datadict['values'] - else: - var[:]=datadict['values'] - -def GetVariable(file,varname): - """ get variable from HDF file""" - return array(file.select(varname).get()) - -def CreateDirs(rundat,dirname): - - dirname=os.path.join(rundat.outputdir,dirname) - if not os.path.exists(dirname): - print "Creating new output directory "+dirname - os.makedirs(dirname) - else: - print 'Writing files to directory: %s'%(dirname,) - return dirname - - -if __name__ == '__main__': - - try: - os.remove('test.nc') - except: - pass - ncf=CT_CDF('test.nc','create') - dimgrid=ncf.AddLatLonDim() - dimdate=ncf.AddDateDim() - dimidate=ncf.AddDateDimFormat() - diff --git a/da/ct/obs.py b/da/ct/obs.py deleted file mode 100755 index 97cda98887c586cb76684568b8778c944501bc74..0000000000000000000000000000000000000000 --- a/da/ct/obs.py +++ /dev/null @@ -1,249 +0,0 @@ -#!/usr/bin/env python -# obs.py - -""" -Author : peters - -Revision History: -File created on 28 Jul 2010. - -""" -import os -import sys -import logging -import datetime -sys.path.append(os.getcwd()) - -identifier = 'CarbonTracker CO2 mixing ratios' -version = '0.0' - -from da.baseclasses.obs import Observation - -################### Begin Class CtObservations ################### - -class CtObservations(Observation): - """ an object that holds data + methods and attributes needed to manipulate mixing ratio values """ - - def getid(self): - return identifier - - def getversion(self): - return version - - def Initialize(self,DaCycle): - - self.startdate = DaCycle['time.sample.start'] - self.enddate = DaCycle['time.sample.end'] - DaSystem = DaCycle.DaSystem - - sfname = DaSystem['obs.input.fname'] - filename = os.path.join(DaSystem['obs.input.dir'],sfname) - - if not os.path.exists(filename): - msg = 'Could not find the required observation input file (%s) ' % filename ; logging.error(msg) - raise IOError,msg - else: - self.ObsFilename = filename - - self.Data = MixingRatioList([]) - - return None - - def AddObs(self): - """ Returns a MixingRatioList holding individual MixingRatioSample objects for all obs in a file - - The CarbonTracker mixing ratio files are provided as one long list of obs for all possible dates. So we can - either: - - (1) read all, and the subselect the data we will use in the rest of this cycle - (2) Use nco to make a subset of the data - - For now, we will stick with option (1) - - """ - import Nio - import datetime as dtm - from string import strip, join - from numpy import array, logical_and - - ncf = Nio.open_file(self.ObsFilename) - idates = ncf.variables['date_components'].get_value() - 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) - - ids = ncf.variables['id'].get_value().take(subselect,axis=0) - sites = ncf.variables['site'].get_value().take(subselect,axis=0) - sites = [s.tostring().lower() for s in sites] - sites = map(strip,sites) - lats = ncf.variables['lat'].get_value().take(subselect,axis=0) - lons = ncf.variables['lon'].get_value().take(subselect,axis=0) - alts = ncf.variables['alt'].get_value().take(subselect,axis=0) - obs = ncf.variables['obs'].get_value().take(subselect,axis=0) - species = ncf.variables['species'].get_value().take(subselect,axis=0) - date = ncf.variables['date'].get_value().take(subselect,axis=0) - window = ncf.variables['sampling_strategy'].get_value().take(subselect,axis=0) - flags = ncf.variables['NOAA_QC_flags'].get_value().take(subselect,axis=0) - flags = [s.tostring().lower() for s in flags] - flags = map(strip,flags) - dummy = ncf.close() - - msg = "Successfully read data from obs file (%s)"%self.ObsFilename ; logging.debug(msg) - - - for n in range(len(dates)): - - self.Data.append( MixingRatioSample(dates[n],sites[n],obs[n],0.0,0.0,0.0,0.0,flags[n],alts[n],lats[n],lons[n],ids[n],0.0) ) - - msg = "Added %d observations to the Data list" % len(dates) ; logging.debug(msg) - - def AddSimulations(self,filename,silent=True): - """ Adds model simulated values to the mixing ratio objects """ - - import Nio - - ncf = Nio.open_file(filename,format='hdf') - ids = ncf.variables['id'].get_value() - simulated = ncf.variables['sampled_mole_frac'].get_value()*1e6 - dummy = ncf.close() - msg = "Successfully read data from model sample file (%s)"%filename ; logging.info(msg) - - obs_ids = self.Data.getvalues('id') - - obs_ids = obs_ids.tolist() - ids = map(int,ids) - - missing_samples = [] - - for id,val in zip(ids,simulated): - - if id in obs_ids: - - index = obs_ids.index(id) - dummy = self.Data[index].simulated = val[0] - - else: - - missing_samples.append(id) - - if not silent and missing_samples != []: - msg = 'Model samples were found that did not match any ID in the observation list. Skipping them...' ; logging.warning(msg) - #msg = '%s'%missing_samples ; logging.warning(msg) - - msg = "Added %d simulated values to the Data list" % (len(ids)-len(missing_samples) ) ; logging.debug(msg) - - def WriteSampleInfo(self, DaCycle): - """ - Write the information needed by the observation operator to a file. Return the filename that was written for later use - - For CarbonTracker, the input files are already in a format suitable for TM5 to read. All we need to do return the name of the file - that served as input for the Observation object. - - """ - import shutil - - obsinputfile = os.path.join(DaCycle['dir.input'],'observations.nc') - - dummy = shutil.copy(self.ObsFilename,obsinputfile) - msg = " [copying source] .... %s " % self.ObsFilename ; logging.debug(msg) - msg = " [to destination] .... %s " % obsinputfile ; logging.debug(msg) - msg = "Sample input file for obs operator now in place (%s)"%obsinputfile ; logging.info(msg) - - return obsinputfile - - - def AddModelDataMismatch(self, DaCycle): - """ - Get the model-data mismatch values for this cycle. - - """ - -################### End Class CtObservations ################### - - - -################### Begin Class MixingRatioSample ################### - -class MixingRatioSample(object): - """ - Holds the data that defines a Mixing Ratio 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,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',sdev=0.0): - - self.code = code.strip() # Site code - self.xdate = xdate # Date of obs - self.obs = obs # Value observed - self.simulated = simulated # Value simulated by model - self.resid = resid # Mixing ratio residuals - self.hphr = hphr # Mixing ratio prior uncertainty from fluxes and (HPH) and model data mismatch (R) - self.mdm = mdm # Model data mismatch - self.flag = flag # Flag - self.height = height # Sample height - self.lat = lat # Sample lat - self.lon = lon # Sample lon - self.id = evn # Event number - 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 - - def __str__(self): - day=self.xdate.strftime('%Y-%m-%d %H:%M:%S') - return ' '.join(map(str,[self.code,day,self.obs,self.flag,self.id])) - -################### End Class MixingRatioSample ################### - -################### Begin Class MixingRatioList ################### - -class MixingRatioList(list): - """ This is a special type of list that holds MixingRatioSample objects. It has methods to extract data from such a list easily """ - from numpy import array, ndarray - - def getvalues(self,name,constructor=array): - from numpy import ndarray - result = constructor([getattr(o,name) for o in self]) - if isinstance(result,ndarray): - return result.squeeze() - else: - return result - def unflagged(self): - return MixingRatioSample([o for o in self if o.flag == 0]) - def flagged(self): - return MixingRatioSample([o for o in self if o.flag != 0]) - def selectsite(self,site='mlo'): - l=[o for o in self if o.code == site] - return MixingRatioSample(l) - def selecthours(self,hours=range(1,24)): - l=[o for o in self if o.xdate.hour in hours] - return MixingRatioSample(l) - -################### End Class MixingRatioList ################### - -if __name__ == "__main__": - - from da.tools.initexit import StartLogger - from da.tools.pipeline import JobStart - from datetime import datetime - import logging - import sys,os - - sys.path.append(os.getcwd()) - - StartLogger() - - obs = CtObservations() - - DaCycle = JobStart(['-v'],{'rc':'da.rc'}) - DaCycle['time.sample.start'] = datetime(2000,1,1) - DaCycle['time.sample.end'] = datetime(2000,1,2) - - obs.Initialize(DaCycle) - obs.Validate() - obs.AddObs() - print(obs.Data.getvalues('obs')) - diff --git a/da/ct/optimizer.py b/da/ct/optimizer.py deleted file mode 100755 index b1746e5d63bafb52291371d5a5640f0ede518fcb..0000000000000000000000000000000000000000 --- a/da/ct/optimizer.py +++ /dev/null @@ -1,99 +0,0 @@ -#!/usr/bin/env python -# optimizer.py - -""" -Author : peters - -Revision History: -File created on 28 Jul 2010. - -""" - -import os -import sys -sys.path.append(os.getcwd()) -import logging -import datetime -from da.baseclasses.optimizer import Optimizer - - -identifier = 'Ensemble Square Root Filter' -version = '0.0' - -################### Begin Class CtOptimizer ################### - -class CtOptimizer(Optimizer): - """ - This creates an instance of a CarbonTracker optimization object. The base class it derives from is the optimizer object. - Additionally, this CtOptimizer implements a special localization option following the CT2007 method. - - All other methods are inherited from the base class Optimizer. - """ - def getid(self): - return identifier - - def getversion(self): - return version - - - def SetLocalization(self,type='None'): - """ determine which localization to use """ - - if type == 'CT2007': - self.localization = True - self.localizetype = 'CT2007' - else: - self.localization = False - self.localizetype = 'None' - - msg = "Current localization option is set to %s"%self.localizetype ; logging.info(msg) - - def Localize(self,n): - """ localize the Kalman Gain matrix """ - import numpy as np - - if not self.localization: return - if self.localizetype == 'CT2007': - - tvalue=1.97591 - if np.sqrt(self.R[n,n]) >= 1.5: - for r in range(self.nlag*self.nparams): - corr=np.corrcoef(self.HX_prime[n,:],self.X_prime[r,:].squeeze())[0,1] - prob=corr/np.sqrt((1.0-corr**2)/(self.nmembers-2)) - if abs(prob) < tvalue: - self.KG[r,n]=0.0 - - -################### End Class CtOptimizer ################### - - - - - -if __name__ == "__main__": - - import os - import sys - from da.tools.initexit import StartLogger - from da.tools.pipeline import JobStart - - sys.path.append(os.getcwd()) - - opts = ['-v'] - args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} - - StartLogger() - - DaCycle = JobStart(opts,args) - DaCycle.Initialize() - - opt = CtOptimizer() - - nobs = 100 - dims = ( int(DaCycle['time.nlag']), - int(DaCycle['forecast.nmembers']), - int(DaCycle.DaSystem['nparameters']), - nobs, ) - - opt.Initialize(dims) - opt.SetLocalization(type='CT2007') diff --git a/da/ct/standardvariables.py b/da/ct/standardvariables.py deleted file mode 100755 index 9eac8f64499b0e5e93c6f65c5de4e2f30681d877..0000000000000000000000000000000000000000 --- a/da/ct/standardvariables.py +++ /dev/null @@ -1,188 +0,0 @@ -standard_variables = { 'bio_flux_prior' : {'name' : 'bio_flux_prior',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, terrestrial vegetation, not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'bio_flux_opt' : {'name' : 'bio_flux_opt',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, terrestrial biosphere , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_prior' : {'name' : 'ocn_flux_prior',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, open ocean , not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_opt' : {'name' : 'ocn_flux_opt',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, open ocean , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'fossil_flux_imp' : {'name' : 'fossil_flux_imp',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, fossil fuel burning , imposed ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'fire_flux_imp' : {'name' : 'fire_flux_imp',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, biomass burning , imposed ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'bio_flux_prior_cov' : {'name' : 'bio_flux_prior_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, terrestrial vegetation , not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'bio_flux_opt_cov' : {'name' : 'bio_flux_opt_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, terrestrial vegetation , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_prior_cov' : {'name' : 'ocn_flux_prior_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, open ocean , not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_opt_cov' : {'name' : 'ocn_flux_opt_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, open ocean , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'decimal_date' : {'name' : 'decimal_date',\ - 'units' : 'years' ,\ - 'long_name' : 'dates and times', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'date', \ - 'dims' : (), \ - 'dtype' : 'double', \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'date' : {'name' : 'date',\ - 'units' : 'days since 2000-01-01 00:00:00 UTC' ,\ - 'long_name' : 'UTC dates and times', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'date', \ - 'dims' : (), \ - 'dtype' : 'double', \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'idate' : {'name' : 'idate',\ - 'units' : 'yyyy MM dd hh mm ss ' ,\ - 'long_name' : 'integer components of date and time', \ - 'standard_name' : 'calendar_components', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'dims' : (), \ - 'dtype' : 'int', \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'latitude' : {'name' : 'latitude',\ - 'units' : 'degrees_north ' ,\ - 'long_name' : 'latitude', \ - 'standard_name' : 'latitude', \ - 'comment' : 'center of interval',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'longitude' : {'name' : 'longitude',\ - 'units' : 'degrees_east ' ,\ - 'long_name' : 'longitude', \ - 'standard_name' : 'longitude', \ - 'comment' : 'center of interval',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'height' : {'name' : 'height',\ - 'units' : 'masl ' ,\ - 'long_name' : 'height_above_ground_level', \ - 'standard_name' : 'height_above_ground_level', \ - 'comment' : 'value is meters above sea level',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'co2' : {'name' : 'co2',\ - 'units' : 'micromol mol-1 ' ,\ - 'long_name' : 'mole_fraction_of_carbon_dioxide_in_air', \ - 'standard_name' : 'mole_fraction_of_carbon_dioxide_in_air', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'meanstate' : {'name' : 'statevectormean',\ - 'units' : 'unitless' ,\ - 'long_name' : 'mean_value_of_state_vector', \ - 'standard_name' : 'mean_value_of_state_vector', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ensemblestate': {'name' : 'statevectorensemble',\ - 'units' : 'unitless' ,\ - 'long_name' : 'ensemble_value_of_state_vector', \ - 'standard_name' : 'ensemble_value_of_state_vector', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'unknown' : {'name' : '',\ - 'units' : '' ,\ - 'long_name' : '', \ - 'standard_name' : '', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - } - - - - diff --git a/da/ct/statevector.py b/da/ct/statevector.py deleted file mode 100755 index 3ef5f858f0bd73e02bce17eebb3f1044f2249fd6..0000000000000000000000000000000000000000 --- a/da/ct/statevector.py +++ /dev/null @@ -1,167 +0,0 @@ -#!/usr/bin/env python -# ct_statevector_tools.py - -""" -Author : peters - -Revision History: -File created on 28 Jul 2010. - -""" - -import os -import sys -sys.path.append(os.getcwd()) - -import logging -import datetime -from da.baseclasses.statevector import EnsembleMember, StateVector -import numpy as np - -identifier = 'CarbonTracker Statevector ' -version = '0.0' - -################### Begin Class CtEnsembleMember ################### - -class CTEnsembleMember(EnsembleMember): - - def AddCustomFields(self,filehandle,DaCycle): - """ Placeholder to add more custom fields to the output of the ensemble member""" - from da.tools.io import CT_CDF, std_savedict - from da.ct.tools import StateToGrid - import numpy as np - - - data = StateToGrid(self.ParameterValues,DaCycle.DaSystem.regionmap) - - dimgrid = filehandle.AddLatLonDim() - - savedict = std_savedict.copy() - savedict['name'] = "parametermap" - savedict['long_name'] = "parametermap_for_member_%d"%self.membernumber - savedict['units'] = "unitless" - savedict['dims'] = dimgrid - savedict['values'] = data.tolist() - savedict['comment'] = 'These are gridded parameter values to use for member %d'%self.membernumber - dummy = filehandle.AddData(savedict) - - return None - -################### Begin Class CtStateVector ################### - -class CtStateVector(StateVector): - """ This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """ - - def getid(self): - return identifier - - def getversion(self): - return version - - def GetCovariance(self,DaSystem): - """ Make a new ensemble from specified matrices, the attribute lag refers to the position in the state vector. - Note that lag=1 means an index of 0 in python, hence the notation lag-1 in the indexing below. - The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag] - """ - - import Nio - - # Get the needed matrices from the specified covariance files - - file_ocn_cov = DaSystem['ocn.covariance'] - file_bio_cov = DaSystem['bio.covariance'] - - for file in [file_ocn_cov,file_bio_cov]: - - if not os.path.exists(file): - - msg = "Cannot find the specified file %s" % file ; logging.error(msg) - raise IOError,msg - else: - - msg = "Using covariance file: %s" % file ; logging.info(msg) - - f_ocn = Nio.open_file(file_ocn_cov,'r') - f_bio = Nio.open_file(file_bio_cov,'r') - - cov_ocn = f_ocn.variables['CORMAT'].get_value() - cov_bio = f_bio.variables['qprior'].get_value() - - dummy = f_ocn.close() - dummy = f_bio.close() - - dummy = logging.debug("Succesfully closed files after retrieving prior covariance matrices") - - # Once we have the matrices, we can start to make the full covariance matrix, and then decompose it - - fullcov = np.zeros((self.nparams,self.nparams),float) - - nocn = cov_ocn.shape[0] - nbio = cov_bio.shape[0] - - fullcov[0:nbio,0:nbio] = cov_bio - fullcov[nbio:nbio+nocn,nbio:nbio+nocn] = 0.16 *np.dot(cov_ocn,cov_ocn) - fullcov[nocn+nbio,nocn+nbio] = 1.e-10 - - # Visually inspect full covariance matrix if verbose - - try: - plt.imshow(fullcov) - plt.colorbar() - plt.savefig('fullcovariancematrix.png') - plt.close('all') - except: - pass - - return fullcov - - def GetNewMember(self,memberno): - """ Return an CTensemblemember object """ - - return CTEnsembleMember(memberno) - - -################### End Class CtStateVector ################### - - -if __name__ == "__main__": - - import os - import sys - from da.tools.initexit import StartLogger - from da.tools.pipeline import JobStart - import numpy as np - - sys.path.append(os.getcwd()) - - opts = ['-v'] - args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} - - StartLogger() - - DaCycle = JobStart(opts,args) - - DaCycle.Initialize() - - StateVector = CtStateVector() - - dims = ( int(DaCycle['time.nlag']), - int(DaCycle['forecast.nmembers']), - int(DaCycle.DaSystem['nparameters']), - ) - - StateVector.Initialize(dims) - - for n in range(dims[0]): - cov = StateVector.GetCovariance(DaCycle.DaSystem) - dummy = StateVector.MakeNewEnsemble(n+1,cov) - - StateVector.Propagate() - - savedir = DaCycle['dir.output'] - filename = os.path.join(savedir,'savestate.nc') - - dummy = StateVector.WriteToFile(DaCycle) - - StateVector.ReadFromFile(filename) - diff --git a/da/ct/tools.py b/da/ct/tools.py deleted file mode 100755 index 20cb223b8c4aa8d2930d9c070bbce845bcad5c1f..0000000000000000000000000000000000000000 --- a/da/ct/tools.py +++ /dev/null @@ -1,103 +0,0 @@ -#!/usr/bin/env python -# ct_tools.py - -""" -Author : peters - -Revision History: -File created on 12 Feb 2009. - -""" -import os -import sys -import logging -import datetime - -identifier = 'CarbonTracker CO2' - -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 - -if __name__ == "__main__": - - sys.path.append('../../') - - import os - import sys - from da.tools.general import StartLogger - from da.tools.initexit import CycleControl - import numpy as np - import da.tools.rc as rc - from pylab import * - - opts = ['-v'] - args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} - - StartLogger() - DaCycle = CycleControl(opts,args) - - DaCycle.Initialize() - print DaCycle - - a=arange(240)+100 - - b = StateToGrid(DaCycle.DaSystem,a) - - figure() - imshow(b) - colorbar() - print b.max() - - c = StateToGrid(DaCycle.DaSystem,b,reverse=True,avg=True) - - figure() - print c.max() - - plot(a,label='original') - plot(c,label='reconstructed') - legend(loc=0) - - - show() - - - - - diff --git a/da/doc/Makefile b/da/doc/Makefile deleted file mode 100644 index 4c501ee2bdb51292924be5024abbe4df4d2f5151..0000000000000000000000000000000000000000 --- a/da/doc/Makefile +++ /dev/null @@ -1,130 +0,0 @@ -# Makefile for Sphinx documentation -# - -# You can set these variables from the command line. -SPHINXOPTS = -SPHINXBUILD = sphinx-build-2.6 -PAPER = -BUILDDIR = build - -# Internal variables. -PAPEROPT_a4 = -D latex_paper_size=a4 -PAPEROPT_letter = -D latex_paper_size=letter -ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source - -.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest - -help: - @echo "Please use \`make <target>' where <target> is one of" - @echo " html to make standalone HTML files" - @echo " dirhtml to make HTML files named index.html in directories" - @echo " singlehtml to make a single large HTML file" - @echo " pickle to make pickle files" - @echo " json to make JSON files" - @echo " htmlhelp to make HTML files and a HTML help project" - @echo " qthelp to make HTML files and a qthelp project" - @echo " devhelp to make HTML files and a Devhelp project" - @echo " epub to make an epub" - @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter" - @echo " latexpdf to make LaTeX files and run them through pdflatex" - @echo " text to make text files" - @echo " man to make manual pages" - @echo " changes to make an overview of all changed/added/deprecated items" - @echo " linkcheck to check all external links for integrity" - @echo " doctest to run all doctests embedded in the documentation (if enabled)" - -clean: - -rm -rf $(BUILDDIR)/* - -html: - $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html - @echo - @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." - -dirhtml: - $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml - @echo - @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml." - -singlehtml: - $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml - @echo - @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml." - -pickle: - $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle - @echo - @echo "Build finished; now you can process the pickle files." - -json: - $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json - @echo - @echo "Build finished; now you can process the JSON files." - -htmlhelp: - $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp - @echo - @echo "Build finished; now you can run HTML Help Workshop with the" \ - ".hhp project file in $(BUILDDIR)/htmlhelp." - -qthelp: - $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp - @echo - @echo "Build finished; now you can run "qcollectiongenerator" with the" \ - ".qhcp project file in $(BUILDDIR)/qthelp, like this:" - @echo "# qcollectiongenerator $(BUILDDIR)/qthelp/CarbonTrackerDataAssimilationShell.qhcp" - @echo "To view the help file:" - @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/CarbonTrackerDataAssimilationShell.qhc" - -devhelp: - $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp - @echo - @echo "Build finished." - @echo "To view the help file:" - @echo "# mkdir -p $$HOME/.local/share/devhelp/CarbonTrackerDataAssimilationShell" - @echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/CarbonTrackerDataAssimilationShell" - @echo "# devhelp" - -epub: - $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub - @echo - @echo "Build finished. The epub file is in $(BUILDDIR)/epub." - -latex: - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex - @echo - @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." - @echo "Run \`make' in that directory to run these through (pdf)latex" \ - "(use \`make latexpdf' here to do that automatically)." - -latexpdf: - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex - @echo "Running LaTeX files through pdflatex..." - make -C $(BUILDDIR)/latex all-pdf - @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." - -text: - $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text - @echo - @echo "Build finished. The text files are in $(BUILDDIR)/text." - -man: - $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man - @echo - @echo "Build finished. The manual pages are in $(BUILDDIR)/man." - -changes: - $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes - @echo - @echo "The overview file is in $(BUILDDIR)/changes." - -linkcheck: - $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck - @echo - @echo "Link check complete; look for any errors in the above output " \ - "or in $(BUILDDIR)/linkcheck/output.txt." - -doctest: - $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest - @echo "Testing of doctests in the sources finished, look at the " \ - "results in $(BUILDDIR)/doctest/output.txt." diff --git a/da/doc/source/_static/banner1.png b/da/doc/source/_static/banner1.png deleted file mode 100644 index d505cf3edc998de70a54fac87cf735639554a004..0000000000000000000000000000000000000000 Binary files a/da/doc/source/_static/banner1.png and /dev/null differ diff --git a/da/doc/source/_templates/layout.html b/da/doc/source/_templates/layout.html deleted file mode 100644 index 15667a5d4dc7b877e22db86c24e1f43f12db75cb..0000000000000000000000000000000000000000 --- a/da/doc/source/_templates/layout.html +++ /dev/null @@ -1,7 +0,0 @@ -{% extends "!layout.html" %} - -{% block header %} -<div style="background-color: white; text-align: left; padding: 10px 10px 15px 15px"> -<a href="#"><img src="_static/banner1.png" border="0" alt="CarbonTracker"/></a> -</div> -{% endblock %} diff --git a/da/doc/source/conf.py b/da/doc/source/conf.py deleted file mode 100644 index 097333b4526d77dc09898f55fc7f7626bfbf6490..0000000000000000000000000000000000000000 --- a/da/doc/source/conf.py +++ /dev/null @@ -1,217 +0,0 @@ -# -*- coding: utf-8 -*- -# -# CarbonTracker Data Assimilation Shell documentation build configuration file, created by -# sphinx-quickstart on Sun Sep 26 13:39:23 2010. -# -# This file is execfile()d with the current directory set to its containing dir. -# -# Note that not all possible configuration values are present in this -# autogenerated file. -# -# All configuration values have a default; values that are commented out -# serve to show the default. - -import sys, os - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -sys.path.insert(0, os.path.abspath(os.path.join(os.getcwd(),'../../../'))) -print sys.path - -# -- General configuration ----------------------------------------------------- - -# If your documentation needs a minimal Sphinx version, state it here. -#needs_sphinx = '1.0' - -# Add any Sphinx extension module names here, as strings. They can be extensions -# coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.todo'] - -# Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] - -# The suffix of source filenames. -source_suffix = '.rst' - -# The encoding of source files. -#source_encoding = 'utf-8-sig' - -# The master toctree document. -master_doc = 'index' - -# General information about the project. -project = u'CarbonTracker Data Assimilation Shell' -copyright = u'2010, Wouter Peters' - -# The version info for the project you're documenting, acts as replacement for -# |version| and |release|, also used in various other places throughout the -# built documents. -# -# The short X.Y version. -version = '0.1' -# The full version, including alpha/beta/rc tags. -release = '0.1' - -# The language for content autogenerated by Sphinx. Refer to documentation -# for a list of supported languages. -#language = None - -# There are two options for replacing |today|: either, you set today to some -# non-false value, then it is used: -#today = '' -# Else, today_fmt is used as the format for a strftime call. -#today_fmt = '%B %d, %Y' - -# List of patterns, relative to source directory, that match files and -# directories to ignore when looking for source files. -exclude_patterns = [] - -# The reST default role (used for this markup: `text`) to use for all documents. -#default_role = None - -# If true, '()' will be appended to :func: etc. cross-reference text. -#add_function_parentheses = True - -# If true, the current module name will be prepended to all description -# unit titles (such as .. function::). -#add_module_names = True - -# If true, sectionauthor and moduleauthor directives will be shown in the -# output. They are ignored by default. -#show_authors = False - -# The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' - -# A list of ignored prefixes for module index sorting. -#modindex_common_prefix = [] - - -# -- Options for HTML output --------------------------------------------------- - -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -html_theme = 'sphinxdoc' - -# Theme options are theme-specific and customize the look and feel of a theme -# further. For a list of options available for each theme, see the -# documentation. -#html_theme_options = {} - -# Add any paths that contain custom themes here, relative to this directory. -#html_theme_path = [] - -# The name for this set of Sphinx documents. If None, it defaults to -# "<project> v<release> documentation". -html_title = 'CarbonTracker DAS' - -# A shorter title for the navigation bar. Default is the same as html_title. -#html_short_title = None - -# The name of an image file (relative to this directory) to place at the top -# of the sidebar. -#html_logo = './images/carbontracker.png' - -# The name of an image file (within the static path) to use as favicon of the -# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 -# pixels large. -#html_favicon = None - -# Add any paths that contain custom static files (such as style sheets) here, -# relative to this directory. They are copied after the builtin static files, -# so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] - -# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, -# using the given strftime format. -#html_last_updated_fmt = '%b %d, %Y' - -# If true, SmartyPants will be used to convert quotes and dashes to -# typographically correct entities. -#html_use_smartypants = True - -# Custom sidebar templates, maps document names to template names. -#html_sidebars = {} - -# Additional templates that should be rendered to pages, maps page names to -# template names. -#html_additional_pages = {} - -# If false, no module index is generated. -#html_domain_indices = True - -# If false, no index is generated. -#html_use_index = True - -# If true, the index is split into individual pages for each letter. -#html_split_index = False - -# If true, links to the reST sources are added to the pages. -#html_show_sourcelink = True - -# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. -#html_show_sphinx = True - -# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. -#html_show_copyright = True - -# If true, an OpenSearch description file will be output, and all pages will -# contain a <link> tag referring to it. The value of this option must be the -# base URL from which the finished HTML is served. -#html_use_opensearch = '' - -# This is the file name suffix for HTML files (e.g. ".xhtml"). -#html_file_suffix = None - -# Output file base name for HTML help builder. -htmlhelp_basename = 'CarbonTrackerDataAssimilationShelldoc' - - -# -- Options for LaTeX output -------------------------------------------------- - -# The paper size ('letter' or 'a4'). -#latex_paper_size = 'letter' - -# The font size ('10pt', '11pt' or '12pt'). -#latex_font_size = '10pt' - -# Grouping the document tree into LaTeX files. List of tuples -# (source start file, target name, title, author, documentclass [howto/manual]). -latex_documents = [ - ('index', 'CarbonTrackerDataAssimilationShell.tex', u'CarbonTracker Data Assimilation Shell Documentation', - u'Wouter Peters', 'manual'), -] - -# The name of an image file (relative to this directory) to place at the top of -# the title page. -#latex_logo = None - -# For "manual" documents, if this is true, then toplevel headings are parts, -# not chapters. -#latex_use_parts = False - -# If true, show page references after internal links. -#latex_show_pagerefs = False - -# If true, show URL addresses after external links. -#latex_show_urls = False - -# Additional stuff for the LaTeX preamble. -#latex_preamble = '' - -# Documents to append as an appendix to all manuals. -#latex_appendices = [] - -# If false, no module index is generated. -#latex_domain_indices = True - - -# -- Options for manual page output -------------------------------------------- - -# One entry per manual page. List of tuples -# (source start file, name, description, authors, manual section). -man_pages = [ - ('index', 'carbontrackerdataassimilationshell', u'CarbonTracker Data Assimilation Shell Documentation', - [u'Wouter Peters'], 1) -] diff --git a/da/doc/source/contents.rst b/da/doc/source/contents.rst deleted file mode 100644 index 1d8f599a7c2773359f53a2c593a64504bf99ed4d..0000000000000000000000000000000000000000 --- a/da/doc/source/contents.rst +++ /dev/null @@ -1,13 +0,0 @@ -.. _contents: - -Contents -======== - -.. toctree:: - :maxdepth: 2 - - System Requirements <systemrequirements> - Installing CTDAS <installing> - CTDAS Overview <overview> - CTDAS Tutorial <tutorial> - CTDAS Documentation <documentation> diff --git a/da/doc/source/docdesign.rtf b/da/doc/source/docdesign.rtf deleted file mode 100644 index d83f512dcda592e05d4bc302e91e55ad31df27a5..0000000000000000000000000000000000000000 --- a/da/doc/source/docdesign.rtf +++ /dev/null @@ -1,47 +0,0 @@ -{\rtf1\ansi\ansicpg1252\cocoartf1038\cocoasubrtf320 -{\fonttbl\f0\fswiss\fcharset0 Helvetica;} -{\colortbl;\red255\green255\blue255;} -\paperw11900\paperh16840\margl1440\margr1440\vieww9000\viewh8400\viewkind0 -\pard\tx566\tx1133\tx1700\tx2267\tx2834\tx3401\tx3968\tx4535\tx5102\tx5669\tx6236\tx6803\ql\qnatural\pardirnatural - -\f0\fs24 \cf0 Design of CT Documentation\ -\ -\ -Welcome\ -System Requirements\ -Installation\ -\ -Contents\ -\ -Overview of CTDAS\ - Design philosophy\ - Extensions\ -\ -Tutorial\ - Getting Started\ - Creating rc-files\ - Creating a run script\ - Setting up your platform\ - Running an inversion\ - Starting a new job\ - Restart from a crash\ - Extend an existing simulation\ - Modifying the system\ - Defining your own classes\ - Creating a new pipeline\ -\ -Reference/Documentation\ - Cycle Control\ - Platforms\ - DA System\ - Observations\ - State Vectors\ - Optimization\ - Observation Operator\ - Pipeline\ - Baseclasses\ - I/O\ -\ -Modules, Classes, Functions\ -\ - } \ No newline at end of file diff --git a/da/doc/source/documentation.rst b/da/doc/source/documentation.rst deleted file mode 100644 index 7d69a17624624cb2521b5596b895beeda10a87db..0000000000000000000000000000000000000000 --- a/da/doc/source/documentation.rst +++ /dev/null @@ -1,5 +0,0 @@ -.. _documentation: - -Documentation -====================== - diff --git a/da/doc/source/gettingstarted.rst b/da/doc/source/gettingstarted.rst deleted file mode 100644 index 18d11eac1bf5eb77fd914f086237a89f0e694c33..0000000000000000000000000000000000000000 --- a/da/doc/source/gettingstarted.rst +++ /dev/null @@ -1,7 +0,0 @@ -.. _gettingstarted: - -Getting Started with CTDAS -========================== - -.. automodule:: da.ct.obs - :members: CtObservations diff --git a/da/doc/source/images/carbontracker.png b/da/doc/source/images/carbontracker.png deleted file mode 100644 index e6da1e1af8e9f84166f8c5518014f351e66efd31..0000000000000000000000000000000000000000 Binary files a/da/doc/source/images/carbontracker.png and /dev/null differ diff --git a/da/doc/source/index.rst b/da/doc/source/index.rst deleted file mode 100644 index 05e8ca132f8f0c7e2f13a57a36741aaf472e33f2..0000000000000000000000000000000000000000 --- a/da/doc/source/index.rst +++ /dev/null @@ -1,30 +0,0 @@ -.. CarbonTracker Data Assimilation Shell documentation master file, created by - sphinx-quickstart on Sun Sep 26 13:39:23 2010. - You can adapt this file completely to your liking, but it should at least - contain the root `toctree` directive. - -Welcome to the CTDAS documentation! -=================================== - -This is the starting page for the documentation of the CarbonTracker Data Assimilation Shell (CTDAS). - -.. sidebar:: What is CTDAS? - - CTDAS is short for "CarbonTracker Data Assimilation Shell". - This is the implementation of an extendible data assimilation framework for CarbonTracker, - developed by NOAA ESRL & Wageningen University, in close cooperation with many partners around the world. - -The aim of the CTDAS system is to facilitate the use of CarbonTracker, and to -foster its development by its many international partners. At its most basic, the CTDAS is a simple -control system for CarbonTracker that deals with the running, optimization, analysis, and time stepping of the system. -For advanced users, the CTDAS provides an easy way to extend or modify the system by introducing new configurations -for the state vector, for the observations, or even for the transport model or optimization method. - -The CarbonTracker development team welcomes suggestions and feedback on the CTDAS system, and its documentation. - -Follow the link below to start navigating, or use the "next topic" button in the right-hand-side menu. - -.. toctree:: - :maxdepth: 1 - - List of contents <contents> diff --git a/da/doc/source/overview.rst b/da/doc/source/overview.rst deleted file mode 100644 index 9e3e0395b5fc6fe2867bf572115c6a2dd71ee3cf..0000000000000000000000000000000000000000 --- a/da/doc/source/overview.rst +++ /dev/null @@ -1,48 +0,0 @@ -.. _overview: - -.. index:: overview, philosophy, development centers - -Overview -======== - -This page gives a brief overview of CTDAS, with the intention to provide background information for (potential) new users. Afer reading the information below we recommend contacting one of the CTDAS development centers (Wageningen University and NOAA ESRL), or proceed to the :doc:`tutorial`. - -What is CTDAS? --------------- - -CTDAS is short for "CarbonTracker Data Assimilation Shell". -This is the implementation of an extendible data assimilation framework for CarbonTracker, -developed by NOAA ESRL & Wageningen University, in close cooperation with many partners around the world. - -The aim of the CTDAS system is to facilitate the use of CarbonTracker, and to -foster its development by its many international partners. - - -How can you use CTDAS? ----------------------- -At its most basic, the CTDAS is a simple -control system for CarbonTracker that deals with the running, optimization, analysis, and time stepping of the system. -For advanced users, the CTDAS provides an easy way to extend or modify the system by introducing new configurations -for the state vector, for the observations, or even for the transport model or optimization method. - -CTDAS Design philosophy ------------------------ - -CTDAS is implemented in python. Each component of the data assimilation system is written as a separate python class, which are -combined in a pipeline. Separate pipelines exist for an inverse simulation, or a forward run. Four classes compose the core of -the system: Observations, StateVector, ObservationOperator, and Optimizer. These are controlled by a CycleControl object, which -itself holds information from class PlatForm and class DaSystem. For each of these seven components, a "baseclass" exists that -describes the basic layout, and mandatory methods that each class needs to have. A specific implementation of a baseclass -"inherits" this basic behavior, and then extends or overwrites the methods. - -As a typical example: the baseclass class:`ObservationOperator` is a nearly empty object, that contains methods to Initialize() the ObservationOperator, to Validate() the input for it, to Run() the operator, and to SaveData() that is needed for a next cycle. The class TM5ObservationOperator is derived from this baseclass, and has a new implementation of each of these methods. When they are called in the pipeline, a TM5 model run is prepared and executed by the methods of the class, returning a set of samples for the optimizer to use in the -minimum least squares method. - -This design philosophy makes the pipeline completely independent of the specific implementation. It simply calls methods from the clasees it receives in a given order, and ensures that the flow of time and data works. The pipeline does not know (or need to know) whether it is running a gridded inversion with satellite data and the TM5 model without zoom, or whether it is running a methane inversion with TM3. The implementation is thus in the hands of the user. - -Extending CTDAS ---------------- -To extend CTDAS, one needs to write new classes that inherit from one of the baseclasses, and have specific functionality implemented under the methods called from the pipeline. - - - diff --git a/da/doc/source/systemrequirements.rst b/da/doc/source/systemrequirements.rst deleted file mode 100644 index 8abd0d4b0685a7d800c681ceff87b5b92028bc49..0000000000000000000000000000000000000000 --- a/da/doc/source/systemrequirements.rst +++ /dev/null @@ -1,53 +0,0 @@ -.. _systemrequirements: - -System Requirements -====================== - -The CarbonTracker DAS is programmed in python and uses it built-in functionality for many of its tasks. -Users of CTDAS are required to have a python installation on their system, with a small set of -open source add-ons. The minimum requirements are listed below, categorized by functionality. - -Getting the CTDAS code ----------------------- - -* A working ``subversion`` (SVN) installation is needed, to check your system type:: - - $ svn --version - - if the system returns an error, or a version < 1.5.0, please obtain svn from <http://subversion.tigris.org/> - - -Running CTDAS -------------- - -* A ``python2.4`` or later installation. - -.. note:: ``python3`` is not supported, one needs a python2.x version. To check your python version type:: - - $ python --version - -* The python module ``numpy``, not included with a standard installation. You can obtain numpy from <http://numpy.scipy.org> - -* The python module ``Nio`` created by NCAR/UCAR. This package is freely available as a binary download from UCAR's Earth System Grid at <http://www.earthsystemgrid.org> - -.. note:: Registration is required at NCAR/UCAR, their software is free for academic users. - - -Performing inversions ---------------------- - -* An important component of the CTDAS system is the observation operator, usually in the form of an atmospheric - transport model. This model is a stand-alone piece of code that can subsample your state vector, and return a - set of observations. **You therefore require an atmospheric transport model that can run indepently on your - platform.** - -.. note:: If you do not currently have access to such a transport model, or have no resources to run such a model - for long periods of time and for many ensemble configurations, you have probably stumbled onto this page not - understanding exactly what CTDAS is. We refer you to the :ref:`overview` for a more elaborate description. - -* An actual CarbonTracker inversion requires a large collection of input datasets, including - * Observations to be assimilated - * a-priori surface fluxes for your trace components - * a-priori uncertainty distributions for your unknowns - * ... and a few more items - diff --git a/da/doc/source/tutorial.rst b/da/doc/source/tutorial.rst deleted file mode 100644 index 1528abcfd4c1c982032c3a87febfcc7d32eda9da..0000000000000000000000000000000000000000 --- a/da/doc/source/tutorial.rst +++ /dev/null @@ -1,5 +0,0 @@ -.. _tutorial: - -Tutorial -====================== - diff --git a/da/platform/__init__.py b/da/platform/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/da/platform/jet.py b/da/platform/jet.py deleted file mode 100755 index 5194f1c4cf34c6a208584b575a11f7aaae1bfaa8..0000000000000000000000000000000000000000 --- a/da/platform/jet.py +++ /dev/null @@ -1,94 +0,0 @@ -#!/usr/bin/env python -# jet.py - -""" -Author : peters - -Revision History: -File created on 06 Sep 2010. - -""" - -import sys -import os -import logging -import subprocess - -from da.baseclasses.platform import PlatForm - -std_joboptions={'jobname':'test','jobaccount':'co2','jobnodes':'nserial 1','jobshell':'/bin/sh','depends':'','jobtime':'00:30:00','joblog':os.getcwd()} - -class JetPlatForm(PlatForm): - def __init__(self): - self.Identifier = 'NOAA jet' # the identifier gives the platform name - self.Version = '1.0' # the platform version used - - msg1 = '%s platform object initialized'%self.Identifier ; logging.debug(msg1) - msg2 = '%s version: %s'%(self.Identifier,self.Version) ; logging.debug(msg2) - - - def GetJobTemplate(self,joboptions={},block=False): - """ Return the job template for a given computing system, and fill it with options from the dictionary provided as argument""" - - template = """#$ -N jobname \n"""+ \ - """#$ -A jobaccount \n"""+ \ - """#$ -pe jobnodes \n"""+ \ - """#$ -l h_rt=jobtime \n"""+ \ - """#$ -S jobshell \n"""+ \ - """#$ -o joblog \n"""+ \ - """#$ -cwd\n"""+ \ - """#$ -r n\n"""+ \ - """#$ -j y\n""" - - if 'depends' in joboptions: - template += """#$ -hold_jid depends \n""" - - if block: - template += """#$ -sync y\n""" - - # First replace from passed dictionary - for k,v in joboptions.iteritems(): - while k in template: - template = template.replace(k,v) - - # Fill remaining values with std_options - for k,v in std_joboptions.iteritems(): - while k in template: - template = template.replace(k,v) - - return template - - def GetMyID(self): - try: - return os.environ['JOB_ID'] - except: - return os.getpid() - - def SubmitJob(self,jobfile): - """ This method submits a jobfile to the queue, and returns the queue ID """ - - cmd = ['qsub',jobfile] - output = subprocess.Popen(cmd, stdout=subprocess.PIPE).communicate()[0] ; logging.info(output) - retcode = output.split()[-1] - - msg = 'Job returned: %s'%retcode ; logging.debug(msg) - - return retcode - - def KillJob(self,jobid): - """ This method kills a running job """ - - output = subprocess.Popen(['qdel',jobid], stdout=subprocess.PIPE).communicate()[0] ; logging.info(output) - - return output - - def StatJob(self,jobid): - """ This method gets the status of a running job """ - import subprocess - - #output = subprocess.Popen(['sgestat'], stdout=subprocess.PIPE).communicate()[0] ; logging.info(output) - - return '' - -if __name__ == "__main__": - pass diff --git a/da/platform/maunaloa.py b/da/platform/maunaloa.py deleted file mode 100755 index a30dca84479dc3c3f483cac81df04d5288e06582..0000000000000000000000000000000000000000 --- a/da/platform/maunaloa.py +++ /dev/null @@ -1,29 +0,0 @@ -#!/usr/bin/env python -# maunaloa.py - -""" -Author : peters - -Revision History: -File created on 06 Sep 2010. - -""" - -import sys -import os -import logging -import subprocess - -from da.baseclasses.platform import PlatForm - -class MaunaloaPlatForm(PlatForm): - def __init__(self): - self.Identifier = 'WUR maunaloa' # the identifier gives the platform name - self.Version = '1.0' # the platform version used - - msg1 = '%s platform object initialized'%self.Identifier ; logging.debug(msg1) - msg2 = '%s version: %s'%(self.Identifier,self.Version) ; logging.debug(msg2) - - -if __name__ == "__main__": - pass diff --git a/da/test/test_optimizer.py b/da/test/test_optimizer.py deleted file mode 100755 index e86a7d6eddfbcd55367ce4b0bd2b540848e52c87..0000000000000000000000000000000000000000 --- a/da/test/test_optimizer.py +++ /dev/null @@ -1,243 +0,0 @@ -#!/usr/bin/env python -# test_optimizer.py - -""" -Author : peters - -Revision History: -File created on 04 Aug 2010. - -""" - -def SerialPyAgainstSerialFortran(): - """ Test the solution of the serial algorithm against the CT cy2 fortran generated one """ - - # get data from the savestate.hdf file from the first cycle of CarbonTracker 2009 release - - print "WARNING: The optimization algorithm has changed from the CT2009 release because of a bug" - print "WARNING: in the fortran code. Hence, the two solutions calculated are no longer the same." - print "WARNING: To change the python algorithm so that it corresponds to the fortran, change the" - print "WARNING: loop from m=n+1,nfcast to m=1,nfcast" - - savefile = '/data/CO2/peters/carbontracker/raw/ct09rc0i/20000108/savestate.hdf' - print savefile - - f = Nio.open_file(savefile,'r') - obs = f.variables['co2_obs_fcast'].get_value() - sel_obs = obs.shape[0] - - dims = ( int(DaCycle.da_settings['time.nlag']), - int(DaCycle.da_settings['forecast.nmembers']), - int(DaCycle.DaSystem.da_settings['nparameters']), - sel_obs, ) - - nlag,nmembers,nparams, nobs = dims - - optserial = CtOptimizer(dims) - opt = optserial - - opt.SetLocalization('CT2007') - - obs = f.variables['co2_obs_fcast'].get_value()[0:nobs] - opt.obs = obs - sim = f.variables['co2_sim_fcast'].get_value()[0:nobs] - opt.Hx = sim - error = f.variables['error_sim_fcast'].get_value()[0:nobs] - flags = f.variables['flag_sim_fcast'].get_value()[0:nobs] - opt.flags = flags - simana = f.variables['co2_sim_ana'].get_value()[0:nobs] - - for n in range(nobs): opt.R[n,n] = np.double(error[n]**2) - - xac=[] - adX=[] - for lag in range(nlag): - xpc = f.variables['xpc_%02d'%(lag+1)].get_value() - opt.x[lag*nparams:(lag+1)*nparams] = xpc - X = f.variables['pdX_%02d'%(lag+1)].get_value() - opt.X_prime[lag*nparams:(lag+1)*nparams,:] = np.transpose(X) - HX = f.variables['dF'][:,0:sel_obs] - opt.HX_prime[:,:] = np.transpose(HX) - - # Also create arrays of the analysis of the fortran code for later comparison - - xac.extend ( f.variables['xac_%02d'%(lag+1)].get_value()) - adX.append (f.variables['adX_%02d'%(lag+1)].get_value() ) - - xac=np.array(xac) - X_prime=np.array(adX).swapaxes(1,2).reshape((opt.nparams*opt.nlag,opt.nmembers)) - - opt.SerialMinimumLeastSquares() - - print "Maximum differences and correlation of 2 state vectors:" - print np.abs(xac-opt.x).max(),np.corrcoef(xac,opt.x)[0,1] - - plt.figure(1) - plt.plot(opt.x,label='SerialPy') - plt.plot(xac,label='SerialFortran') - plt.grid(True) - plt.legend(loc=0) - plt.title('Analysis of state vector') - - print "Maximum differences of 2 state vector deviations:" - print np.abs(X_prime-opt.X_prime).max() - - plt.figure(2) - plt.plot(opt.X_prime.flatten(),label='SerialPy') - plt.plot(X_prime.flatten(),label='SerialFortran') - plt.grid(True) - plt.legend(loc=0) - plt.title('Analysis of state vector deviations') - - print "Maximum differences and correlation of 2 simulated obs vectors:" - print np.abs(simana-opt.Hx).max(),np.corrcoef(simana,opt.Hx)[0,1] - - plt.figure(3) - plt.plot(opt.Hx,label='SerialPy') - plt.plot(simana,label='SerialFortran') - plt.grid(True) - plt.legend(loc=0) - plt.title('Analysis of CO2 mixing ratios') - plt.show() - - f.close() - -def SerialvsBulk(): - """ A test of the two algorithms currently implemented: serial vs bulk solution """ - - # get data from the savestate.hdf file from the first cycle of CarbonTracker 2009 release - - savefile = '/data/CO2/peters/carbontracker/raw/ct09rc0i/20000108/savestate.hdf' - print savefile - - f = Nio.open_file(savefile,'r') - obs = f.variables['co2_obs_fcast'].get_value() - - nobs = 77 - - dims = ( int(DaCycle.da_settings['time.nlag']), - int(DaCycle.da_settings['forecast.nmembers']), - int(DaCycle.DaSystem.da_settings['nparameters']), - nobs, ) - - nlag,nmembers,nparams, nobs = dims - - optbulk = CtOptimizer(dims) - optserial = CtOptimizer(dims) - - for o,opt in enumerate([optbulk,optserial]): - - opt.SetLocalization('CT2007') - - obs = f.variables['co2_obs_fcast'].get_value()[0:nobs] - opt.obs = obs - sim = f.variables['co2_sim_fcast'].get_value()[0:nobs] - opt.Hx = sim - error = f.variables['error_sim_fcast'].get_value()[0:nobs] - flags = f.variables['flag_sim_fcast'].get_value()[0:nobs] - opt.flags = flags - - for n in range(nobs): - opt.R[n,n] = np.double(error[n]**2) - - xac=[] - for lag in range(nlag): - xpc = f.variables['xpc_%02d'%(lag+1)].get_value() - opt.x[lag*nparams:(lag+1)*nparams] = xpc - X = f.variables['pdX_%02d'%(lag+1)].get_value() - opt.X_prime[lag*nparams:(lag+1)*nparams,:] = np.transpose(X) - HX = f.variables['dF'][:,0:nobs] - opt.HX_prime[:,:] = np.transpose(HX) - - if o == 0: - opt.BulkMinimumLeastSquares() - x1=opt.x - xp1=opt.X_prime - hx1=opt.Hx - hxp1=opt.HX_prime - hphr1=opt.HPHR - k1=opt.KG - if o == 1: - opt.SerialMinimumLeastSquares() - x2=opt.x - xp2=opt.X_prime - hx2=opt.Hx - hxp2=opt.HX_prime - hphr2=opt.HPHR - k2=opt.KG - - plt.figure() - - print "Maximum differences and correlation of 2 state vectors:" - print np.abs(x2-x1).max(),np.corrcoef(x2,x1)[0,1] - - plt.figure(1) - plt.plot(x1,label='Serial') - plt.plot(x2,label='Bulk') - plt.grid(True) - plt.legend(loc=0) - plt.title('Analysis of state vector') - - print "Maximum differences of 2 state vector deviations:" - print np.abs(xp2-xp1).max() - - plt.figure(2) - plt.plot(xp1.flatten(),label='Serial') - plt.plot(xp2.flatten(),label='Bulk') - plt.grid(True) - plt.legend(loc=0) - plt.title('Analysis of state vector deviations') - - print "Maximum differences and correlation of 2 simulated obs vectors:" - print np.abs(hx2-hx1).max(),np.corrcoef(hx2,hx1)[0,1] - - plt.figure(3) - plt.plot(hx1,label='Serial') - plt.plot(hx2,label='Bulk') - plt.title('Analysis of CO2 mixing ratios') - plt.grid(True) - plt.legend(loc=0) - - plt.show() - - f.close() - - - -if __name__ == "__main__": - - import sys - - sys.path.append('../../') - - import os - from da.tools.general import StartLogger - from da.tools.initexit import CycleControl - from da.ct.statevector import CtStateVector, CtMember, PrepareState - from da.ct.obs import CtObservations, MixingRatioSample - from da.ct.optimizer import CtOptimizer - import numpy as np - import da.tools.rc as rc - import datetime - import Nio - import matplotlib.pyplot as plt - - opts = ['-v'] - args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} - - StartLogger() - DaCycle = CycleControl(opts,args) - - dummy = CtMember(0) - DaCycle.Initialize() - StateVector = CtStateVector(DaCycle) - samples = CtObservations(DaCycle.DaSystem,datetime.datetime(2000,1,1)) - - # This is a test of the CarbonTracker release ct2009 fortran solution and the one implemented here - - #SerialPyAgainstSerialFortran() - #sys.exit(0) - - SerialvsBulk() - sys.exit(0) - diff --git a/da/tm5/__init__.py b/da/tm5/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/da/tm5/observationoperator.py b/da/tm5/observationoperator.py deleted file mode 100755 index 6fbf10d78e2fb4c82a9adcf901df23a513e8a43d..0000000000000000000000000000000000000000 --- a/da/tm5/observationoperator.py +++ /dev/null @@ -1,460 +0,0 @@ -#!/usr/bin/env python -# tm5_tools.py - -""" -Author : peters - -Revision History: -File created on 09 Feb 2009. -Major modifications to go to a class-based approach, July 2010. - -This module holds specific functions needed to use the TM5 model within the data assimilation shell. It uses the information -from the DA system in combination with the generic tm5.rc files. - -The TM5 model is now controlled by a python subprocess. This subprocess consists of an MPI wrapper (written in C) that spawns -a large number ( N= nmembers) of TM5 model instances under mpirun, and waits for them all to finish. - -The design of the system assumes that the tm5.x (executable) was pre-compiled with the normal TM5 tools, and is residing in a -directory specified by the ${RUNDIR} of a tm5 rc-file. This tm5 rc-file name is taken from the data assimilation rc-file. Thus, -this python shell does *not* compile the TM5 model for you! - -""" - -import os, sys -import subprocess -import logging -import shutil - -sys.path.append(os.getcwd()) - -identifier = 'TM5' -version = 'release 3.0' -mpi_shell_filename = 'tm5_mpi_wrapper' -mpi_shell_location = 'da/bin/' - -needed_rc_items = [ - 'rundir', - 'outputdir', - 'time.start', - 'time.break.nday' - ] - -################### Begin Class TM5 ################### - -from da.baseclasses.observationoperator import ObservationOperator - -class TM5ObservationOperator(ObservationOperator): - """ This class holds methods and variables that are needed to run the TM5 model. It is initiated with as only argument a TM5 rc-file - location. This rc-file will be used to figure out the settings for the run. - - *** This method of running TM5 assumes that a pre-compiled tm5.exe is present, and it will be run from time.start to time.final *** - - These settings can be modified later. To run a model version, simply compile the model using an existing TM5 rc-file, then - open python, and type: - - []> tm=TM5('/Users/peters/Modeling/TM5/tutorial.rc') - []> tm.WriteRc() - []> tm.WriteRunRc() - []> tm.Run() - - To use this class inside a data assimilation cycle, a stand-alone method "Initialize()" is included which modifies the TM5 - settings according to an external dictionary of values to overwrite, and then runs the TM5 model. - - - """ - - def __init__(self): - """ The instance of an TMObservationOperator is application dependent """ - self.Identifier = self.getid() # the identifier gives the model name - self.Version = self.getversion() # the model version used - - msg = 'Observation Operator object initialized: %s (%s)'%(self.Identifier, self.Version,) ; logging.info(msg) - - def getid(self): - return identifier - - def getversion(self): - return version - - def Initialize(self, DaCycle): - """ - Prepare a forward model TM5 run, this consists of: - - - reading the TM5 rc-file, - - validating it, - - modifying the values, - - Creating a tm5_runtime.rc file - - Removing the existing tm5.ok file if present - - """ - from da.tools.general import CreateLinks - - RcFileName = DaCycle['forecast.model.rc'] - dummy = self.LoadRc(RcFileName) # load the specified rc-file - dummy = self.ValidateRc() # validate the contents - -# Create a link from TM5 to the rundirectory of the das system - - sourcedir = self.tm_settings['rundir'] - targetdir = os.path.join(DaCycle['dir.exec'],'tm5') - dummy = CreateLinks(sourcedir,targetdir) - - DaCycle['dir.exec.tm5'] = targetdir - -# Write a modified TM5 model rc-file in which run/break times are defined by our da system - - NewItems = { - 'time.start' : DaCycle['time.sample.start'] , - 'time.final' : DaCycle['time.sample.end'] , - 'rundir' : DaCycle['dir.exec.tm5'] , - 'outputdir' : DaCycle['dir.output'] , - 'savedir' : DaCycle['dir.save'] , - 'das.input.dir' : DaCycle['dir.input'] - } - - if DaCycle['time.restart']: # If this is a restart from a previous cycle, the TM5 model should do a restart - NewItems['istart'] = 3 - if DaCycle['time.sample.window'] != 0: # If this is a restart from a previous time step wihtin the filter lag, the TM5 model should do a restart - NewItems['istart'] = 3 - - # If neither one is true, simply take the istart value from the tm5.rc file that was read - - self.ModifyRC(NewItems) - - self.WriteRc() - - self.WriteRunRc() - -# Copy the pre-compiled MPI wrapper to the execution directory - - targetdir = os.path.join(DaCycle['dir.exec.tm5']) - - if not os.path.exists(os.path.join(mpi_shell_location,mpi_shell_filename) ): - msg = "Cannot find the mpi_shell wrapper needed for completion (%s) in (%s)"% (mpi_shell_filename,mpi_shell_location) ; logging.error(msg) - msg = "Please see the %s/readme_wrapper.txt file for instructions to compile it"%mpi_shell_location ; logging.error(msg) - raise IOError - - shutil.copy(os.path.join(mpi_shell_location,mpi_shell_filename) ,os.path.join(targetdir,mpi_shell_filename) ) - - return 0 - - def LoadRc(self,RcFileName): - """ - This method loads a TM5 rc-file with settings for this simulation - """ - import da.tools.rc as rc - - self.tm_settings = rc.read(RcFileName) - self.RcFileName = RcFileName - self.Tm5RcLoaded = True - - msg = 'TM5 rc-file loaded successfully' ; logging.debug(msg) - - return True - - def ValidateRc(self): - """ - Validate the contents of the tm_settings dictionary and add extra values. The required items for the TM5 rc-file - are specified in the tm5_tools module, as dictionary variable "needed_rc_items". - """ - from da.tools.general import ToDatetime - - for k,v in self.tm_settings.iteritems(): - if v == 'True' : self.tm_settings[k] = True - if v == 'False': self.tm_settings[k] = False - if 'date' in k : self.tm_settings[k] = ToDatetime(v) - if 'time.start' in k : - self.tm_settings[k] = ToDatetime(v,fmt='TM5') - if 'time.final' in k : - self.tm_settings[k] = ToDatetime(v,fmt='TM5') - - for key in needed_rc_items: - - if not self.tm_settings.has_key(key): - status,msg = ( False,'Missing a required value in rc-file : %s' % key) - logging.error(msg) - raise IOError,msg - - status,msg = ( True,'rc-file has been validated succesfully' ) ; logging.debug(msg) - - - - - def ModifyRC(self,NewValues): - """ - Modify parts of the tm5 settings, for instance to give control of file locations to the DA shell - instead of to the tm5.rc script. - - Note that we replace these values in all {key,value} pairs of the tm5.rc file! - - """ - - for k,v in NewValues.iteritems(): - if self.tm_settings.has_key(k): - # keep previous value - v_orig= self.tm_settings[k] - #replace with new - self.tm_settings[k] = v - #replace all instances of old with new, but only if it concerns a name of a path!!! - if os.path.exists(str(v)): - for k_old,v_old in self.tm_settings.iteritems(): - if not isinstance(v_old,str): - continue - if str(v_orig) in str(v_old): - v_new = str(v_old).replace(str(v_orig),str(v)) - self.tm_settings[k_old] = v_new - - logging.debug('Replaced tm5 rc-item %s ' % k) - - else: - self.tm_settings[k] = v - logging.debug('Added new tm5 rc-item %s ' % k) - - - def WriteRc(self): - """ - Write the rc-file settings to a tm5.rc file in the rundir - """ - import da.tools.rc as rc - - tm5rcfilename = os.path.join(self.tm_settings['rundir'],'tm5.rc') - - dummy = rc.write(tm5rcfilename,self.tm_settings) - - msg = "Modified tm5.rc written (%s)" % tm5rcfilename ;logging.debug(msg) - - def WriteRunRc(self): - """ - Create the tm5-runtime.rc file which is read by initexit.F90 to control time loop and restart from save files - """ - import da.tools.rc as rc - - tm5rcfilename = os.path.join(self.tm_settings['rundir'],'tm5_runtime.rc') - - dummy = rc.write(tm5rcfilename,self.tm_settings) - - rc_runtm5={} - - rc_runtm5['year1'] = self.tm_settings['time.start'].year - rc_runtm5['month1'] = self.tm_settings['time.start'].month - rc_runtm5['day1'] = self.tm_settings['time.start'].day - rc_runtm5['hour1'] = self.tm_settings['time.start'].hour - rc_runtm5['minu1'] = 0 - rc_runtm5['sec1'] = 0 - - rc_runtm5['year2'] = self.tm_settings['time.final'].year - rc_runtm5['month2'] = self.tm_settings['time.final'].month - rc_runtm5['day2'] = self.tm_settings['time.final'].day - rc_runtm5['hour2'] = self.tm_settings['time.final'].hour - rc_runtm5['minu2'] = 0 - rc_runtm5['sec2'] = 0 - - rc_runtm5['istart'] = self.tm_settings['istart'] - rc_runtm5['savedir'] = self.tm_settings['savedir'] - rc_runtm5['outputdir'] = self.tm_settings['outputdir'] - - rc.write(tm5rcfilename,rc_runtm5) - - msg = "Modified tm5_runtime.rc written (%s)" % tm5rcfilename ;logging.debug(msg) - - - def ValidateInput(self,DaCycle): - """ - Make sure that parameter files are written to the TM5 inputdir, and that observation lists are present - """ - - datadir = self.tm_settings['das.input.dir'] - if not os.path.exists(datadir): - msg = "The specified input directory for the TM5 model to read from does not exist (%s), exiting..."%datadir ; logging.error(msg) - raise IOError,msg - - datafiles = os.listdir(datadir) - - obsfile = DaCycle['ObsOperator.inputfile'] - - if not os.path.exists(obsfile): - msg = "The specified obs input file for the TM5 model to read from does not exist (%s), exiting..."%obsfile ; logging.error(msg) - raise IOError,msg - - - for n in range(int(DaCycle['forecast.nmembers'])): - paramfile = 'parameters.%03d.nc'%n - if paramfile not in datafiles: - msg = "The specified parameter input file for the TM5 model to read from does not exist (%s), exiting..."%paramfile ; logging.error(msg) - raise IOError,msg - - # Next, make sure there is an actual model version compiled and ready to execute - - targetdir = os.path.join(self.tm_settings['rundir']) - self.Tm5Executable = os.path.join(targetdir,'tm5.x') - if not os.path.exists(self.Tm5Executable): - msg = "Required TM5 executable was not found %s"%self.Tm5Executable ; logging.error(msg) - msg = "Please compile the model with the specified rc-file and the regular TM5 scripts first" ; logging.error(msg) - raise IOError - - return 0 - - def Run(self,DaCycle): - """ - Start the TM5 executable. A new log file is started for the TM5 model IO, and then a subprocess is - spawned with the tm5_mpi_wrapper and the tm5.x executable. The exit code of the model is caught and - only if successfull on all processors will execution of the shell continue. - - """ - - cwd = os.getcwd() - - # From here on, several options should be implemented. - - # - # (1) Where an mpi process is forked to do a TM5 instance with N tracers, each an ensemble member - # - # (2) Where N processes are spawned, each being one TM5 instance representing one member - # - # (3) Where N/m processes are spawned, each being a TM5 instance that handles m ensemble members - # - # In principle, it is best to make these processes produce scripts that can be executed stand-alone, or - # be submitted to a queue. - # - - # Open logfile and spawn model, wait for finish and return code - - # Code for Option (1) - - code = self.TM5_under_mpirun(DaCycle) - - if code == 0: - logging.info('Finished model executable succesfully (%s)'%code) - self.Status = 'Success' - else: - logging.error('Error in model executable return code: %s ' % code) - self.Status = 'Failed' - raise OSError - - # Return to working directory - - os.chdir(cwd) - - return code - - def TM5_under_mpirun(self, DaCycle): - """ Method handles the case where a shell runs an MPI process that forks into N TM5 model instances """ - from string import join - import datetime - - DaPlatForm = DaCycle.DaPlatForm - - targetdir = os.path.join(self.tm_settings['rundir']) - - # Go to executable directory and start the subprocess, using a new logfile - - dummy = os.chdir(targetdir) - msg = 'Changing directory to %s ' % targetdir ;logging.debug(msg) - - # Remove the tm5.ok file from a previous run, placed back only if a successful TM5 run is executed - - okfile = 'tm5.ok' - if os.path.exists(okfile): - dummy = os.remove(okfile) - - nprocesses = DaCycle['forecast.nmembers'] - jobparams = {'jobname':'tm5', - 'jobnodes':'ncomp %d'%int(nprocesses), - 'jobtime':'00:30:00', - 'joblog':os.path.join(DaCycle['dir.jobs']) - } - - template = DaPlatForm.GetJobTemplate(jobparams,block=True) - template += 'cd %s\n'%targetdir - template += 'mpirun -np %d %s ./tm5.x\n'%(int(nprocesses),mpi_shell_filename,) - - jobfile = DaPlatForm.WriteJob(DaCycle,template,'tm5') - - msg = 'Submitting job at %s'%datetime.datetime.now() ; logging.debug(msg) - code = DaPlatForm.SubmitJob(jobfile) - msg = 'Resuming job at %s'%datetime.datetime.now() ; logging.debug(msg) - - if not os.path.exists(okfile): - code = -1 - else: - code = 0 - - return code - - - def TM5_With_N_tracers(self, nprocesses): - """ Method handles the case where one TM5 model instance with N tracers does the sampling of all ensemble members""" - from string import join - - modellogfile = open(self.ModelLogFilename,'w') - #cmd = ['openmpirun','-np', '10', mpi_shell_filename,'./tm5.x'] - cmd = ['mpirun','-np',str(nprocesses),'./tm5.x'] - msg = 'Starting model executable as subprocess (%s)'%join(cmd,' ') ; logging.info(msg) - - #code = subprocess.call(cmd,stdout=modellogfile,stderr=modellogfile) - - modellogfile.close() - - def SaveData(self): - """ Write the TM5 recovery data for the next cycle """ - - sourcedir = os.path.join(self.tm_settings['outputdir']) - targetdir = os.path.join(self.tm_settings['savedir']) - filter = ['save_%s'%self.tm_settings['time.final'].strftime('%Y%m%d')] - - msg = "Performing a full backup of TM5 save data" ; logging.debug(msg) - msg = " from directory: %s " % sourcedir ; logging.debug(msg) - msg = " to directory: %s " % targetdir ; logging.debug(msg) - msg = " with filter: %s " % filter ; logging.debug(msg) - - for file in os.listdir(sourcedir): - - file = os.path.join(sourcedir,file) - - if os.path.isdir(file): # skip dirs - skip = True - elif filter == []: # copy all - skip= False - else: # check filter - skip = True # default skip - for f in filter: - if f in file: - skip = False # unless in filter - break - - if skip: - msg = " [skip] .... %s " % file ; logging.debug(msg) - continue - - msg = " [copy] .... %s " % file ; logging.debug(msg) - dummy = shutil.copy(file,file.replace(sourcedir,targetdir) ) - -################### End Class TM5 ################### - - -if __name__ == "__main__": - - from da.tools.initexit import StartLogger - from da.tools.pipeline import JobStart - import datetime as dtm - import sys,os - - sys.path.append(os.getcwd()) - - StartLogger() - - DaCycle = JobStart([],{'rc':'da.rc'}) - DaCycle.Initialize() - DaCycle['time.sample.start'] = dtm.datetime(2000,1,1) - DaCycle['time.sample.end'] = dtm.datetime(2000,1,2) - DaCycle['time.sample.window'] = 0 - - tm=TM5ObservationOperator() - tm.Initialize(DaCycle) - tm.Run(DaCycle) - tm.SaveData() - - - - - diff --git a/da/tools/__init__.py b/da/tools/__init__.py deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000 diff --git a/da/tools/general.py b/da/tools/general.py deleted file mode 100755 index 93b06a093f9b9bd19ec0c95a6425b849699de3cb..0000000000000000000000000000000000000000 --- a/da/tools/general.py +++ /dev/null @@ -1,133 +0,0 @@ -#!/usr/bin/env python -# tools_da.py - -""" -Author : peters - -Revision History: -File created on 03 Oct 2008. - -Temporary module to hold classes and methods that are in development - -""" - -import logging -import os -import sys -import shutil -import datetime -import da.tools.rc as rc - -def ValidateRC(rcfile,needed_items): - """ validate the contents of an rc-file given a dictionary of required keys """ - - for k,v in rcfile.iteritems(): - if v == 'True' : rcfile[k] = True - if v == 'False': rcfile[k] = False - if 'date' in k : rcfile[k] = datetime.datetime.strptime(v,'%Y-%m-%d %H:%M:%S') - - for key in needed_items: - - if not rcfile.has_key(key): - status,msg = ( False,'Missing a required value in rc-file : %s' % key) - logging.error(msg) - raise IOError,msg - - status,msg = ( True,'rc-file has been validated succesfully' ) ; logging.debug(msg) - - -def CreateDirs(dirname,forceclean=False): - """ Create a directory and report success, only if non-existent """ - - if forceclean: - try: - shutil.rmtree(dirname) - except: - pass - - if not os.path.exists(dirname): - os.makedirs(dirname) - msg='Creating new directory %s' % dirname - logging.info(msg) - else: - msg='Using existing directory %s' % dirname - logging.debug(msg) - return None - -def CreateLinks(sourcedir,targetdir): - """ Create a symbolic link from, source to target and report success, note - that exisiting links are first removed and then recreated """ - - if os.path.exists(targetdir): - os.unlink(targetdir) - msg='Unlinking existing directory %s' % targetdir - logging.debug(msg) - try: - os.symlink(sourcedir,targetdir) - msg='Created new link from %s to %s' % (sourcedir, targetdir) - logging.debug(msg) - except OSError,msg: - msg='Failed to create link from %s to %s' % (sourcedir, targetdir) - logging.error(msg) - raise OSError - - return None - - - -def AdvanceTime(time_in,interval): - """ Advance time_in by a specified interval""" - - time_out=time_in - - if interval == 'month': # if monthly, this run will go to the first day of the next month - if time_in.month != 12: - time_out = datetime.datetime(time_in.year,time_in.month+1,1,time_in.hour,0,0) - else: - time_out = datetime.datetime(time_in.year+1,1,1,time_in.hour,0,0) # end of year provision - elif interval == 'week': - time_out = time_in + datetime.timedelta(days=7) - else: # assume that the interval specified is the number of days to run forward before resubmitting - time_out = time_in + datetime.timedelta(days=float(interval)) - - return time_out - - - -def ToDatetime(datestring,fmt=None): - """ convert a date string to a datetime object """ - import datetime - - if fmt == 'TM5': - datestring = '%04s-%02s-%02s %02s:%02s:00'%(datestring[0:4],datestring[4:6],datestring[6:8],datestring[8:10],datestring[10:12]) - - try: - return datetime.datetime.strptime(datestring,'%Y-%m-%d %H:%M:%S') - except: - date,time = datestring.split(' ') - year,month,day = map(int,date.split('-')) - hour,minute,second = map(int,time.split(':')) - return datetime.datetime(year,month,day,hour,minute,second) - -def StoreData(pickleobject,filename): - """ pickle object into a specified file name """ - import cPickle - - f = open(filename,'wb') - dummy = cPickle.dump(pickleobject,f,-1) - f.close() - - return None - -def RestoreData(filename): - """ unpickle object into a specified file name """ - import cPickle - - f = open(filename,'rb') - pickleobject = cPickle.load(f) - f.close() - - return pickleobject - -if __name__ == "__main__": - pass diff --git a/da/tools/initexit.py b/da/tools/initexit.py deleted file mode 100755 index 15126e76d18826d139596ab9c6acf051696c3a0b..0000000000000000000000000000000000000000 --- a/da/tools/initexit.py +++ /dev/null @@ -1,630 +0,0 @@ -#!/usr/bin/env python -# da_initexit.py - -""" -Author : peters - -Revision History: -File created on 13 May 2009. - -""" - -needed_da_items=[ - 'time.start', - 'time.finish', - 'time.nlag', - 'time.cycle', - 'dir.da_run', - 'forecast.model', - 'forecast.nmembers', - 'forecast.model.rc', - 'da.system'] - -# only needed in an earlier implemented where each substep was a separate job -# validprocesses = ['start','done','samplestate','advance','invert'] - -import logging -import os -import sys -import shutil -import datetime -import da.tools.rc as rc - - - -class CycleControl(dict): - """ - This object controls a data assimilation cycle. It is created using information from a da rc-file with settings, and - also using the command line options passed to the job control script that runs in the foreground, or through the queue. - This object is passed around throughout the code as it also contains information on directory structures, requested time - windows, and other resources. - """ - - def __init__(self,opts,args): - """ - Initialization occurs from passed options and arguments to the job control script. - """ - - self.LoadRc(args['rc']) - self.ValidateRC() - self.opts = opts - - # Add some useful variables to the rc-file dictionary - - self['jobrcfilename'] = self.RcFileName - self['dir.da_submit'] = os.getcwd() - self['da.crash.recover'] = '-r' in opts - self['verbose'] = '-v' in opts - self.DaSystem = None # to be filled later - - - def __str__(self): - """ - String representation of a CycleControl object - """ - - msg = "===============================================================" ; print msg - msg = "DA Cycle rc-file is %s" % self.RcFileName ; print msg - msg = "DA Cycle run directory is %s" % self['dir.da_run'] ; print msg - msg = "DA Cycle inverse system is %s" % self['da.system'] ; print msg - msg = "DA Cycle forecast model is %s" % self['forecast.model'] ; print msg - msg = "===============================================================" ; print msg - - return "" - - - def LoadRc(self,RcFileName): - """ - This method loads a DA Cycle rc-file with settings for this simulation - """ - - rcdata = rc.read(RcFileName) - for k,v in rcdata.iteritems(): - self[k] = v - self.RcFileName = RcFileName - self.DaRcLoaded = True - - msg = 'DA Cycle rc-file (%s) loaded successfully'%self.RcFileName ; logging.info(msg) - - return True - - - def ValidateRC(self): - """ - Validate the contents of the rc-file given a dictionary of required keys - """ - from da.tools.general import ToDatetime - - for k,v in self.iteritems(): - if v == 'True' : self[k] = True - if v == 'False': self[k] = False - if 'date' in k : self[k] = ToDatetime(v) - if 'time.start' in k : - self[k] = ToDatetime(v) - if 'time.end' in k : - self[k] = ToDatetime(v) - if 'time.finish' in k : - self[k] = ToDatetime(v) - - for key in needed_da_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 Cycle settings have been validated succesfully' ) ; logging.debug(msg) - - return None - - def ParseTimes(self): - """ - Parse time related parameters into datetime objects for later use - """ - from da.tools.general import AdvanceTime - - startdate = self['time.start'] - finaldate = self['time.finish'] - - if finaldate <= startdate: - msg = 'The start date (%s) is not greater than the end date (%s), please revise'%(startdate.strftime('%Y%m%d'),finaldate.strftime('%Y%m%d')) - logging.error(msg) - raise ValueError - # - cyclelength = self['time.cycle'] # get time step - -# Determine end date - - if cyclelength == 'infinite': - enddate = finaldate - else: - enddate = AdvanceTime(startdate,cyclelength) - - # - if enddate > finaldate: # do not run beyon finaldate - enddate = finaldate - - self['time.start'] = startdate - self['time.end'] = enddate - self['time.finish'] = finaldate - self['cyclelength'] = cyclelength - - msg = "===============================================================" ; logging.info(msg) - msg = "DA Cycle start date is %s" % startdate.strftime('%Y-%m-%d %H:%M') ; logging.info(msg) - msg = "DA Cycle end date is %s" % enddate.strftime('%Y-%m-%d %H:%M') ; logging.info(msg) - msg = "DA Cycle final date is %s" % finaldate.strftime('%Y-%m-%d %H:%M') ; logging.info(msg) - msg = "DA Cycle cycle length is %s" % cyclelength ; logging.info(msg) - msg = "DA Cycle restart is %s" % str(self['time.restart']) ; logging.info(msg) - msg = "===============================================================" ; logging.info(msg) - - return None - - def RandomSeed(self,action='read'): - """ Get the randomseed and save it, or read the random seed and set it""" - import cPickle - import numpy as np - - filename = os.path.join(self['dir.save'],'randomseed.pickle') - - if action == 'write': - f = open(filename,'wb') - seed = np.random.get_state() - dummy = cPickle.dump (seed,f,-1) - dummy = f.close() - - msg = "Saved the random seed generator values to file" - - if action == 'read': - f = open(filename,'rb') - seed = cPickle.load(f) - dummy = np.random.set_state(seed) - dummy = f.close() - - msg = "Retrieved the random seed generator values from file" - - logging.info(msg) - - return None - - def Initialize(self): - """ Determine how to proceed with this cycle: - - (a) recover from crash : get the save data from the one-ago folder and replace the da_runtime variables with those from the save dir - (b) restart cycle : use latest da_runtime variables from the exec dir - (c) fresh start : set up the required file structure for this simulation - """ - -# -# case 1: A recover from a previous crash, this is signaled by flag "-r" -# - if self['da.crash.recover']: - - msg = "Recovering simulation from data in: %s" % self['dir.da_run'] ; logging.info(msg) - - dummy = self.SetupFileStructure() - - dummy = self.RecoverRun() - - dummy = self.RandomSeed('read') -# -# case 2: A continuation, this is signaled by rc-item time.restart = True -# - elif self['time.restart']: - - msg = "Restarting filter from previous step" ; logging.info(msg) - - dummy = self.SetupFileStructure() - - dummy = self.RandomSeed('read') -# -# case 3: A fresh start, this is signaled by rc-item time.restart = False -# - elif not self['time.restart']: - msg = "First time step in filter sequence" ; logging.info(msg) - - dummy = self.SetupFileStructure() - - # expand jobrcfilename to include exec dir from now on. - - self['jobrcfilename'] = os.path.join(self['dir.exec'],self['jobrcfilename']) - - - self.ParseTimes() - self.WriteRC(self['jobrcfilename']) - - return None - - - def SetupFileStructure(self): - """ - Create file structure needed for data assimilation system. - In principle this looks like: - - ${da_rundir} - ${da_rundir}/input - ${da_rundir}/output - ${da_rundir}/exec - ${da_rundir}/diagnostics - ${da_rundir}/analysis - ${da_rundir}/jobs - ${da_rundir}/save - ${da_rundir}/save/one-ago - ${da_rundir}/save/two-ago - - Note that the exec dir will actually be a simlink to the directory where - the SIBCASA or model executable lives. This directory is passed through - the da.rc file. The observation input files will be placed in the exec dir, - and the resulting simulated values will be retrieved from there as well. - """ - from da.tools.general import CreateDirs - -# Create the run directory for this DA job, including I/O structure - - filtertime = self['time.start'].strftime('%Y%m%d') - - self['dir.exec'] = os.path.join(self['dir.da_run'],'exec') - self['dir.input'] = os.path.join(self['dir.da_run'],'input') - self['dir.output'] = os.path.join(self['dir.da_run'],'output',filtertime) - self['dir.diagnostics'] = os.path.join(self['dir.da_run'],'diagnostics') - self['dir.analysis'] = os.path.join(self['dir.da_run'],'analysis') - self['dir.jobs'] = os.path.join(self['dir.da_run'],'jobs') - self['dir.save'] = os.path.join(self['dir.da_run'],'save') - - CreateDirs(self['dir.da_run']) - CreateDirs(os.path.join(self['dir.exec'])) - CreateDirs(os.path.join(self['dir.input'])) - CreateDirs(os.path.join(self['dir.output'])) - CreateDirs(os.path.join(self['dir.diagnostics'])) - CreateDirs(os.path.join(self['dir.analysis'])) - CreateDirs(os.path.join(self['dir.jobs'])) - CreateDirs(os.path.join(self['dir.save'])) - CreateDirs(os.path.join(self['dir.save'],'one-ago')) - CreateDirs(os.path.join(self['dir.save'],'two-ago')) - - msg = 'Succesfully created the file structure for the assimilation job' ; logging.info(msg) - - - def RecoverRun(self): - """Prepare a recovery from a crashed run. This consists of: - - - copying all data from the save/one-ago folder, - - replacing all rc-file items with those from da_runtime.rc - - resetting the seed of the random number generator to the value it had before the crash - - """ - - # Move all data from the save/one-ago directory to the save/current directory - - dummy = self.MoveSaveData(io_option='restore',save_option='full',filter=[]) - - # Replace rc-items with those from the crashed run last rc-file - - file_rc_rec = os.path.join(self['dir.save'],'da_runtime.rc') - rc_rec = rc.read(file_rc_rec) - - for k,v in rc_rec.iteritems(): - self[k] = v - - self.ValidateRC() - - msg = "Replaced rc-items.... " ; logging.debug(msg) - msg = "Next cycle start date is %s" % self['time.start'] ; logging.debug(msg) - - return None - - def Finalize(self): - """ - Finalize the da cycle, this means writing the save data and rc-files for the next run - """ - - dummy = self.RandomSeed('write') - dummy = self.WriteNewRCfile() - dummy = self.CollectSaveData() - dummy = self.MoveSaveData(io_option='store',save_option='full') - dummy = self.SubmitNextCycle() - - def CollectSaveData(self): - """ Collect files needed for the restart of this cycle in case of a crash, or for the continuation of the next cycle. - All files needded are written to the save/ directory. - Currently, the following files are included - - -) The randomseed.pickle file, also already in the save directory - -) The da_runtime.rc file - -) The savestate.yyyymmdd.nc file - -) ... - - Note that we assume that the transport model restart files are directly written to the save/ directory and do not need to be collected - - """ - - filtertime = self['time.start'].strftime('%Y%m%d') - - file1 = os.path.join(self['dir.exec'],'da_runtime.rc') - file2 = os.path.join(self['dir.output'],'savestate.nc') - - FilesToSave = [file1,file2] - - targetdir = os.path.join(self['dir.save']) - - msg = "Collecting the required restart/save data" ; logging.info(msg) - msg = " to directory: %s " % targetdir ; logging.debug(msg) - - for file in FilesToSave: - - msg = " [copy] .... %s " % file ; logging.debug(msg) - dummy = shutil.copy(file,file.replace(os.path.split(file)[0],targetdir) ) - - - def MoveSaveData(self, io_option='restore', save_option='partial',filter=[]): - """ Store or restore model state to/from a temporary directory. Two save_option options are available: - - (1) save_option = 'partial' : Save the background concentration fields only, to use for a later advance of the background - (2) save_option = 'full' : Save all restart data, to use for a next cycle - - While also two IO options are available: - - (1) io_option = restore : Get data from a directory - (2) io_option = store : Save data to a directory - - In case of a 'store' command the target folder is re-created so that the contents are empty to begin with. - - There is currently room for the specific sub-models to also add save data to the mix. This is done through a - forecast.model dependend import statement, followed by a call to function: - - model.MoveSaveData(**args) - - """ - from da.tools.general import CreateDirs - - if io_option not in ['store','restore']: - raise ValueError,'Invalid option specified for io_option (%s)' % io_option - if save_option not in ['partial','full']: - raise ValueError,'Invalid option specified for save_option (%s)' % save_option - - - savedir = os.path.join(self['dir.save']) - oneagodir = os.path.join(self['dir.save'],'one-ago') - tempdir = os.path.join(self['dir.save'],'tmp') - - - if save_option == 'partial': # save background data to tmp dir - - if io_option == 'store': - - targetdir = tempdir - sourcedir = savedir - - elif io_option == 'restore': - - sourcedir = tempdir - targetdir = savedir - - elif save_option == 'full': # move all save data to one-ago dir - - if io_option == 'store': - - targetdir = oneagodir - sourcedir = savedir - - elif io_option == 'restore': - - sourcedir = oneagodir - targetdir = savedir - -# If "store" is requested, recreate target dir, cleaning the contents - - if io_option == 'store': - CreateDirs(os.path.join(targetdir),forceclean=True) - - msg = "Performing a %s %s of data" % (save_option,io_option) ; logging.debug(msg) - msg = " from directory: %s " % sourcedir ; logging.debug(msg) - msg = " to directory: %s " % targetdir ; logging.debug(msg) - msg = " with filter : %s " % filter ; logging.debug(msg) - - - for file in os.listdir(sourcedir): - - file = os.path.join(sourcedir,file) - - if os.path.isdir(file): # skip dirs - skip = True - elif filter == []: # copy all - skip= False - else: # check filter - skip = True # default skip - for f in filter: - if f in file: - skip = False # unless in filter - break - - if skip: - msg = " [skip] .... %s " % file ; logging.debug(msg) - continue - - - if io_option == 'store' and save_option == 'full': - #msg = " [move] .... %s " % file ; logging.debug(msg) - #dummy = shutil.move(file,file.replace(sourcedir,targetdir) ) - msg = " [copy] .... %s " % file ; logging.debug(msg) - dummy = shutil.copy(file,file.replace(sourcedir,targetdir) ) - - else: - msg = " [copy] .... %s " % file ; logging.debug(msg) - dummy = shutil.copy(file,file.replace(sourcedir,targetdir) ) - -# - def WriteNewRCfile(self): - """ Write the rc-file for the next DA cycle. Note that the start time for the next cycle is the end time of this one, while - the end time for the next cycle is the current end time + one cycle length. The resulting rc-file is written to: - - (1) to the dir.exec so that it can be used when resubmitting the next cycle - - """ - from da.tools.general import AdvanceTime - - # These first two lines advance the filter time for the next cycle - - self['time.start'] = self['time.end'] - self['time.end'] = AdvanceTime(self['time.end'],self['cyclelength']) - - # The rest is info needed for a system restart - - self['time.restart'] = True - #self['time.start'] = self['time.start'].strftime('%Y-%m-%d %H:%M:%S') - #self['time.finish'] = self['time.finish'].strftime('%Y-%m-%d %H:%M:%S') - - fname = os.path.join(self['dir.exec'],'da_runtime.rc') - self['da.restart.fname'] = fname - dummy = rc.write(fname,self) - - msg = 'Wrote new da_runtime.rc (%s)'%fname ; logging.debug(msg) - - - def WriteRC(self,fname): - """ Write RC file after each process to reflect updated info """ - - dummy = rc.write(fname,self) - msg = 'Wrote expanded rc-file (%s)'%(fname) ; logging.debug(msg) - return None - - def SubmitNextCycle(self): - """ - Submit next job in cycle, if needed. The end date is held against the final date. A subprocess is called that - submits the next job, without waiting for its result so this cycle can finish and exit. - """ - import subprocess - import os - - - DaPlatForm = self.DaPlatForm - - if self['time.start'] < self['time.finish']: - - jobparams = {'jobname':'das','jobaccount':'co2','jobnodes':'nserial 1','jobtime':'00:30:00','jobshell':'/bin/sh'} - template = DaPlatForm.GetJobTemplate(jobparams) - template += 'cd %s\n'%os.getcwd() - template += './das.py rc=%s'% self['da.restart.fname'] - - cd = self['time.end'].strftime('%Y%m%d') - jobfile = DaPlatForm.WriteJob(self,template,'cycle.%s.das'%cd) - jobid = DaPlatForm.SubmitJob(jobfile) - else: - logging.info('Final date reached, no new cycle started') - - return None - - def SubmitSubStep(self,stepname): - """ submit the next substep of a DA cycle""" - import subprocess - import os - from string import join - - - DaPlatForm = self.DaPlatForm - - jobparams = {'jobname':'das.%s'%stepname} - template = DaPlatForm.GetJobTemplate(jobparams) - template += 'cd %s\n'%os.getcwd() - template += './das.py rc=%s process=%s %s' % (self['jobrcfilename'],stepname,join(self.opts,''),) - jobfile = DaPlatForm.WriteJob(self,template,stepname) - jobid = DaPlatForm.SubmitJob(jobfile) - - return None - - -def StartLogger(): - """ start the logging of messages to screen""" - -# start the logging basic configuration by setting up a log file - - logging.basicConfig(level = logging.INFO, - format = ' [%(levelname)-7s] (%(asctime)s) py-%(module)-20s : %(message)s', - datefmt = '%Y-%m-%d %H:%M:%S') - -def ParseOptions(): - """ Function parses options from the command line and returns the arguments as a dictionary""" - import getopt - import sys - -# Parse keywords, the only option accepted so far is the "-h" flag for help - - opts=[] - args=[] - try: - opts, args = getopt.gnu_getopt(sys.argv[1:], "-hrv") - except getopt.GetoptError, msg: - logging.error('%s'%msg) - sys.exit(2) - - for options in opts: - options=options[0].lower() - if options == '-h': - print "" - print helptext - sys.exit(2) - if options == '-r': - logging.info('-r flag specified on command line: recovering from crash') - if options == '-v': - logging.info('-v flag specified on command line: extra verbose output') - dummy = logging.root.setLevel(logging.DEBUG) - -# Parse arguments and return as dictionary - - arguments={} - for item in args: - #item=item.lower() - -# Catch arguments that are passed not in "key=value" format - - if '=' in item: - key, arg = item.split('=') - else: - logging.error('%s'%'Argument passed without description (%s)' % item) - raise getopt.GetoptError,arg - - arguments[key]=arg - - if opts: opts=opts[0] - - return opts, arguments - -def ValidateOptsArgs(opts,args): - """ Validate the options and arguments passed from the command line before starting the cycle """ - - if not args.has_key("rc"): - msg = "There is no rc-file specified on the command line. Please use rc=yourfile.rc" ; logging.error(msg) - raise IOError,msg - elif not os.path.exists(args['rc']): - msg = "The specified rc-file (%s) does not exist " % args['rc'] ; logging.error(msg) - raise IOError,msg - - # WP not needed anymore - #if not args.has_key('process'): - # msg = "There is no process specified on the command line, assuming process=Start" ; logging.info(msg) - # args['process'] = 'start' - #if args['process'].lower() not in validprocesses: - # msg = "The specified process (%s) is not valid"%args['process'] ; logging.error(msg) - # raise IOError,msg - - return opts,args - -def CleanUpCycle(DaCycle): - """ - Move files, clean up rundir after a cycle has finished - """ - - # move log file to rundir/jobs - jobdir = os.path.join(DaCycle['dir.da_run'],"jobs") - logfile = 'das.o%s'%DaCycle.DaPlatForm.GetMyID() - - if os.path.exists(logfile): - joblogfile = os.path.join(jobdir,'cycle.%s.log'%DaCycle['time.start'].strftime('%Y%m%d')) - dummy = shutil.move(logfile,joblogfile) - msg = "....Moved %s to %s"%(logfile,joblogfile) ; logging.debug(msg) - - msg = "The complete log file is now at: %s"%(joblogfile) ; logging.info(msg) - - -if __name__ == "__main__": - pass - diff --git a/da/tools/io.py b/da/tools/io.py deleted file mode 100755 index 6cce56bec86e449c5ea5f3c073b460b045459caa..0000000000000000000000000000000000000000 --- a/da/tools/io.py +++ /dev/null @@ -1,390 +0,0 @@ -#!/usr/bin/env python -# io.py - -""" -Author : peters - -Revision History: -File created on 15 Oct 2008. -File modified for CT data assimilation system in July 2010, Wouter Peters - -""" -import datetime as dt -from numpy import array -import os - -disclaimer = "This data belongs to the CarbonTracker project" -email = "wouter.peters@wur.nl" -url = "http://carbontracker.wur.nl" -institution = "Wageningen University and Research Center" -source = "CarbonTracker release 2.0" -conventions = "CF-1.1" -historytext = 'Created on '+dt.datetime.now().strftime('%B %d, %Y')+' by %s'%os.environ['USER'] - -std_savedict={'name':'unknown','values':None,'dims':None,'units':'','long_name':'','_FillValue':float(-999.),'comment':''} - -class CT_CDF(object): - """ - - The CT_CDF object is a wrapper around the NIO functionality offered by NCAR's Nio library. It can read+write data in many formats, including - NetCDF, HDF, and GRIB. The CT_CDF class is able to make an instance of an Nio file object, and automatically add (meta)data to it in - a format that complies with CF-1.1 conventions. - - The user can specify a filename when creating an instance of the CT_CDF object, as well as all arguments used by the Nio.open_file function. - An object is returned that holds the Nio file handle in self.file, as well as some global attributes we want in each file coming out of - this data assimilation system. - - The methods of the CT_CDF object are: - - AddCTHeader() : Adds a standard header to the file with information about date, user, system, contact person, and a disclaimer. - AddParamsDim() : Add a dimension called 'nparams' - AddLatLonDim() : Add two dimension called 'latitude' and 'longitude' - AddMembersDim() : Add a dimension called 'nmembers' - AddObsDim() : Add a dimension called 'nobs' - AddLagDim() : Add a dimension called 'nlag', this one can be unlimited - AddDateDim() : Add a dimension called 'date', this one can be unlimited - AddDateDimFormat() : Add a dimension called 'datedimformat', this is a 6-integer format as in yyyy,mm,dd,HH,MM,SS - GetVariable() : return a variable from the dataset by name - StandardVar() : return a dictionary with precooked CF-1.1 information for a dataset to be written - inquire_unlimited() : get the length of the unlimited dimension - AddData() : Add data from a dictionary to the file - close() : close the file - - """ - def __init__(self,filename,*arguments): - import Nio - import os - - if os.path.exists(filename): - os.remove(filename) - self.file = Nio.open_file(filename,*arguments) - self.AddCTHeader() - - - def AddCTHeader(self): - - # - setattr(self.file,'Institution',institution) - setattr(self.file,'Contact',email) - setattr(self.file,'URL',url) - setattr(self.file,'Source',source) - setattr(self.file,'Convention',conventions) - setattr(self.file,'Disclaimer',disclaimer) - setattr(self,'History',historytext) - - def AddParamsDim(self,nparams): - - dimparams=self.file.create_dimension('nparameters',nparams) - - return ('nparameters',) - - def AddMembersDim(self,nmembers): - - dimmembers=self.file.create_dimension('nmembers',nmembers) - - return ('nmembers',) - - def AddLagDim(self,nlag=0,unlimited=True): - - if unlimited: - dimlag =self.file.create_dimension('nlag',None) - else: - dimlag=self.file.create_dimension('nlag',nlag) - - dimlag = ('nlag',) - - # Also create the variable of the same name so it can be queried for length - var = self.file.create_variable('nlag','i',dimlag) - - return ('nlag',) - - def AddObsDim(self,nobs): - - dimobs=self.file.create_dimension('nobs',nobs) - - return ('nobs',) - - def AddLatLonDim(self,istart=0,iend=360,jstart=0,jend=180): - """ - Add dimensions + data for latitude and longitude values of a rectangular 1x1 grid, the - scope of the arrays can be limited by using the istart,...,jend integers - - """ - - from numpy import arange, float64 - - if 'latitude' in self.file.dimensions.keys(): return ('latitude','longitude',) # already exists - - lons=-180+arange(360)*1.0+0.5 - lats=-90+arange(180)*1.0+0.5 - # - lats=lats[jstart:jend] - lons=lons[istart:iend] - # - dimlon=self.file.create_dimension('longitude',lons.shape[0]) - dimlon=('longitude',) - dimlat=self.file.create_dimension('latitude',lats.shape[0]) - dimlat=('latitude',) - - savedict=self.StandardVar(varname='latitude') - savedict['values']=lats.tolist() - savedict['actual_range']=(float(lats[0]),float(lats[-1])) - savedict['dims']=dimlat - self.AddData(savedict) - - savedict=self.StandardVar(varname='longitude') - savedict['values']=lons.tolist() - savedict['actual_range']=(float(lons[0]),float(lons[-1])) - savedict['dims']=dimlon - self.AddData(savedict) - - return dimlat+dimlon - - def AddDateDim(self,ndate=0,unlimited=True): - """ Add a date dimension, give it the unlimited length if requested """ - - if unlimited: - dimdate =self.file.create_dimension('date',None) - else: - dimdate=self.file.create_dimension('date',ndate) - - dimdate = ('date',) - - # Also create the variable of the same name so it can be queried for length - var = self.file.create_variable('date','d',dimdate) - - return dimdate - - def AddDateDimFormat(self): - """ Add a dimension representing a date format as yyyy/mm/dd/HH/MM/SS""" - - if 'yyyymmddhhmmss' in self.file.dimensions.keys(): - pass - else: - dummy = self.file.create_dimension('yyyymmddhhmmss',6) # already exists - dimdateformat = ('yyyymmddhhmmss',) - return dimdateformat - - def has_date(self,dd): - """ Check if a passed date (dd) is already present in the dates included in a file """ - - if self.file.variables['date'].shape[0] > 0: - if dd in self.file.variables['date'][:].tolist(): - return True - else: - return False - else: - return False - - def GetVariable(self,varname): - """ get variable from ncf file""" - return self.file.variables[varname][:] - - def StandardVar(self,varname): - """ return properties of standard variables """ - import standardvariables - - if varname in standardvariables.standard_variables.keys(): - return standardvariables.standard_variables[varname] - else: - var = std_savedict.copy() - var['name'] = varname - return var - - def inquire_unlimited(self): - """ Get the index and name of the unlimited dimension""" - - try: - index = self.file.dimensions.values().index(None) - unlimname = self.file.dimensions.keys()[index] - unlimlen = self.file.variables[unlimname].shape[0] - except: - unlimlen = -1 - unlimname = 'None' - - return (unlimname,unlimlen,) - - def AddData(self,datadict,silent=True): - """ - Add data to the Nio file object. This is achieved by passing a data dictionary to this routine that holds all metadata plus - the real data values. For example: - - savedict=std_savedict.copy() - savedict['name'] = 'testarray' - savedict['dims'] = ('date','latitude',) - savedict['count'] = unlimlen+n - savedict['values'] = np.arange(180)+n*1.5 - ncf.AddData(savedict,silent=False) - - this makes a copy of a standard variable dictionary first, and then populates it with (meta)data. The conventions are that - - ['dims'] = a tuple with names of the dimensions - ['count'] = an integer number that represents the first dimension along which to write (usually the unlimited dim) - ['values'] = a list or array that holds the data values to be written. Note that one cannot pass any other object. - - The counter for the unlimited dimension is needed because multiple datasets might be added with unlimited dimensions and - we don't want the unlimited counter to increase everytime. To get the value of the counter (needed when appending data for instance) - use - - unlimname,unlim_count = ncf.inquire_unlimited() - - where an integer value is returned for the length, and -1 indicates there is no unlimited dimension - - Note that there are special pre-cooked data dictionaries available in the module standardvariables, with many attributes - already specified. For instance: - - savedict=ncf.StandardVar(varname='date') - savedict['values']= [100.0+n] - savedict['count'] = unlimlen+n - savedict['dims']=('date',) - ncf.AddData(savedict,silent=False) - - has long_names, comments, and data ranges pre-specified. - - """ - import numpy as np - - # Make sure the passed dictionary has proper values for the important entries - - if None == datadict['values']: - raise IOError,"The passed data dictionary does not contain values, please provide data to write" - if None == datadict['dims']: - raise IOError,"The passed data dictionary does not contain valid dimensions, please provide names of dimensions for the dataset" - - if not isinstance(datadict['values'],(list,np.ndarray,)): - raise IOError,"Please pass a list or array to the AddData() method, not a %s object"%datadict['values'].__class__ - - # First, try to get the requested dataset number to place next in file. Note that this is an attribute of the - # savedict and has to be specified by the user. If not specified, the 'next' dataset is number 0 - - try: - next = datadict['count'] - except: - next=0 - - # Check if the requested variable already exists in the file, if so, get a variable handle directly - - existing_vars=self.file.variables.keys() - - if datadict['name'] in existing_vars: - - var = self.file.variables[datadict['name']] - - # If the variable name is new, create a new dataset - - else: - if not silent: print 'Creating new dataset: '+datadict['name'] - - if datadict.has_key('dtype'): # datatype is part of the data dictionary passed - - if datadict['dtype'] == 'int': - var = self.file.create_variable(datadict['name'],'i',datadict['dims']) - elif datadict['dtype'] == 'char': - var = self.file.create_variable(datadict['name'],'s1',datadict['dims']) - elif datadict['dtype'] == 'double': - var = self.file.create_variable(datadict['name'],'d',datadict['dims']) - else: - var = self.file.create_variable(datadict['name'],'f',datadict['dims']) - - else: - var = self.file.create_variable(datadict['name'],'d',datadict['dims']) # default is 'double' - - # All other keys in the datadict are assumed to be variable attributes, except for a few reserved ones - - for k,v in datadict.iteritems(): - if k not in ['name','dims','values','_FillValue','count']: - setattr(var,k,v) - - # Now that the variable instance is returned, start to add the data. First figure out if we are dealing with an unlimited dimension - - unlimname,unlimlen = self.inquire_unlimited() - has_unlimited_dim = (unlimname != 'None') - - # If so, check if the first dimension of the passed dataset corresponds to the unlimited dimension - - if has_unlimited_dim and var.dimensions[0] == unlimname: - - # Add the data to the unlimited dim, note the special occassion for passing only one value (a record) - - if len(datadict['values']) == 1: - var[next] = datadict['values'][0] - else: - var[next,:]=datadict['values'] - - # Add the data , note the special occassion for passing only one value (a record) - - else: - if len(datadict['values']) == 1: - var.assign_value(datadict['values'][0]) - else: - var[:]=datadict['values'] - - def close(self): - """ close the file object """ - - return self.file.close() - -def CreateDirs(rundat,dirname): - - dirname=os.path.join(rundat.outputdir,dirname) - if not os.path.exists(dirname): - print "Creating new output directory "+dirname - os.makedirs(dirname) - else: - print 'Writing files to directory: %s'%(dirname,) - return dirname - - -if __name__ == '__main__': - - import sys - import os - import numpy as np - - sys.path.append('../../') - - import da.tools.standardvariables - - try: - os.remove('test.nc') - except: - pass - - ncf=CT_CDF('test.nc','w') - dimgrid=ncf.AddLatLonDim() - dimdate=ncf.AddDateDim(unlimited=True) - #dimlag=ncf.AddLagDim(unlimited=True) - #dimidate=ncf.AddDateDimFormat() - - unlimname,unlimlen = ncf.inquire_unlimited() - - if unlimlen < 0: unlimlen=0 - - for n in range(100): - - savedict=ncf.StandardVar(varname='testarray') - savedict['dims'] = ('date','latitude',) - savedict['count'] = unlimlen+n - savedict['values'] = np.arange(180)+n*1.5 - ncf.AddData(savedict,silent=False) - - savedict=ncf.StandardVar(varname='testarray2') - savedict['dims'] = ('date','latitude',) - savedict['count'] = unlimlen+n - savedict['values'] = np.arange(180)-n*1.5 - ncf.AddData(savedict,silent=False) - - savedict=ncf.StandardVar(varname='date') - savedict['values']= [100.0+n] - savedict['count'] = unlimlen+n - savedict['dims']=dimdate - ncf.AddData(savedict,silent=False) - - #savedict=ncf.StandardVar(varname='testfail') - #savedict['values']=range(5) - #ncf.AddData(savedict,silent=False) - - print ncf.inquire_unlimited() - - ncf.file.close() - diff --git a/da/tools/pipeline.py b/da/tools/pipeline.py deleted file mode 100755 index dc09027b425f0a43c33e6299336ff59a6927395c..0000000000000000000000000000000000000000 --- a/da/tools/pipeline.py +++ /dev/null @@ -1,334 +0,0 @@ -#!/usr/bin/env python -# pipeline.py - -""" -Author : peters - -Revision History: -File created on 26 Aug 2010. - -""" -import logging -import os -import sys -import shutil -import datetime -import da.tools.rc as rc - -header = '\n\n *************************************** ' -footer = ' *************************************** \n ' - -def Main(PlatForm, DaSystem, Samples,StateVector,ObsOperator,Optimizer): - """ The main point of entry for the pipeline """ - -# Append current working dir to path - - sys.path.append(os.getcwd()) - -# Import methods and classes contained in this package - - from da.tools.initexit import CleanUpCycle - from da.tools.initexit import ValidateOptsArgs - from da.tools.initexit import ParseOptions - from da.tools.general import StoreData,RestoreData - - -# Parse options from the command line - - opts, args = ParseOptions() - -# Validate Options and arguments passed - - opts,args = ValidateOptsArgs(opts,args) - -# Create the Cycle Control object for this job - - DaCycle = JobStart(opts,args, DaSystem, PlatForm) - - msg = header+"Initializing Da Cycle"+footer ; logging.info(msg) - - dummy = DaCycle.Initialize() - - msg = header+"starting JobInput"+footer ; logging.info(msg) - - StateVector = JobInput(DaCycle, StateVector) - - msg = header+"starting SampleState"+footer ; logging.info(msg) - - dummy = SampleState(DaCycle,Samples,StateVector,ObsOperator) - - msg = header+"starting Invert"+footer ; logging.info(msg) - - dummy = Invert(DaCycle, StateVector, Optimizer) - - msg = header+"starting Advance"+footer ; logging.info(msg) - - dummy = Advance(DaCycle, Samples, StateVector, ObsOperator) - - msg = header+"starting SaveAndSubmit"+footer ; logging.info(msg) - - dummy = SaveAndSubmit(DaCycle, StateVector) - - msg = "Cycle finished...exiting" ; logging.info(msg) - - dummy = CleanUpCycle(DaCycle) - - sys.exit(0) - - -#################################################################################################### - -def JobStart(opts,args, DaSystem, DaPlatForm): - """ Set up the job specific directory structure and create an expanded rc-file """ - from da.tools.initexit import CycleControl - - # First create a CycleControl object that handles all the details of the da cycle - - DaCycle = CycleControl(opts,args) - - # Next fill the DaSystem object that contains all the details related to the DA system - dummy = DaSystem.Validate() - - dummy = DaSystem.Initialize() - - DaCycle.DaSystem = DaSystem - - DaCycle.DaPlatForm = DaPlatForm - - return DaCycle - - -def JobInput(DaCycle, StateVector): - """ Set up the input data for the forward model: obs and parameters/fluxes""" - - print DaCycle.DaSystem - - dims = ( int(DaCycle['time.nlag']), - int(DaCycle['forecast.nmembers']), - int(DaCycle.DaSystem['nparameters']), - ) - - dummy = StateVector.Initialize(dims) - - nlag = dims[0] - - # We now have an empty StateVector object that we need to populate with data. If this is a continuation from a previous cycle, we can read - # the previous StateVector values from a NetCDF file in the save/ directory. If this is the first cycle, we need to populate the StateVector - # with new values for each week. After we have constructed the StateVector, it will be propagated by one cycle length so it is ready to be used - # in the current cycle - - if not DaCycle['time.restart']: - - # Fill each week from n=1 to n=nlag with a new ensemble - - for n in range(0,nlag): - cov = StateVector.GetCovariance(DaCycle.DaSystem) - dummy = StateVector.MakeNewEnsemble(n+1,cov) - - else: - - # Read the StateVector data from file - - savedir = DaCycle['dir.save'] - filtertime = DaCycle['time.start'].strftime('%Y%m%d') - filename = os.path.join(savedir,'savestate.nc') - - StateVector.ReadFromFile(filename) - - # Now propagate the ensemble by one cycle to prepare for the current cycle - - dummy = StateVector.Propagate() - - return StateVector - -def SampleState(DaCycle,Samples,StateVector, ObservationOperator): - """ Sample the filter state for the inversion """ - - - # Before a forecast step, save all the data to a save/tmp directory so we can later recover it before the propagation step. - # This is especially important for: - # (i) The transport model restart data which holds the background mixing ratios. This is needed to run the model one cycle forward - # (ii) The random numbers (or the seed for the random number generator) so we can recreate the ensembles if needed - - dummy = DaCycle.MoveSaveData(io_option='store',save_option='partial',filter=[]) - msg = "All restart data have been copied to the save/tmp directory for future use" ; logging.debug(msg) - - nlag = int(DaCycle['time.nlag']) - msg = "Sampling model will be run over %d cycles"% nlag ; logging.info(msg) - - for lag in range(nlag): - - msg = header+".....Ensemble Kalman Filter at lag %d"% (int(lag)+1,) ; logging.info(msg) - - SampleOneCycle(DaCycle,Samples,StateVector, ObservationOperator,lag) - - # Optionally, post-processing of the model output can be added that deals for instance with - # sub-sampling of time series, vertical averaging, etc. - - return None - -def SampleOneCycle(DaCycle,Samples,StateVector, ObservationOperator,lag): - """ Perform all actions needed to sample one cycle """ - from da.tools.general import AdvanceTime - - # First set up the information for time start and time end of this sample - - cyclelen = DaCycle['cyclelength'] - - if lag == 0: - startdate = DaCycle['time.start'] - else: - startdate = DaCycle['time.sample.start'] - - - enddate = AdvanceTime(startdate,cyclelen) - - DaCycle['time.sample.start'] = startdate - DaCycle['time.sample.end'] = enddate - DaCycle['time.sample.window'] = lag - - msg = "New simulation interval set : " ; logging.info(msg) - msg = " start date : %s " % startdate.strftime('%F %H:%M') ; logging.info(msg) - msg = " end date : %s " % enddate.strftime('%F %H:%M') ; logging.info(msg) - - - # Implement something that writes the ensemble member parameter info to file, or manipulates them further into the - # type of info needed in your transport model - - dummy = WriteEnsembleData(DaCycle, StateVector, lag) - - dummy = Samples.Initialize(DaCycle) - - dummy = Samples.Validate() - - dummy = Samples.AddObs() - - filename = Samples.WriteSampleInfo(DaCycle) - - DaCycle['ObsOperator.inputfile'] = filename - - # Run the observation operator - - dummy = RunForecastModel(DaCycle,ObservationOperator) - - # Add model-data mismatch to all samples, this *might* use output from the ensemble in the future?? - - dummy = Samples.AddModelDataMismatch(DaCycle) - - msg = "Added Model Data Mismatch to all samples " ; logging.debug(msg) - - # Read forecast model samples that were written to NetCDF files by each member, also add the obs to the statevector - - silent = False - for Member in StateVector.EnsembleMembers[lag]: - filename = os.path.join(DaCycle['dir.output'],'samples.%03d.nc'%Member.membernumber) - dummy = Samples.AddSimulations(filename,silent) - Member.ModelSample = Samples - silent=True - - msg = "Added samples from the observation operator to each member of the StateVector" ; logging.debug(msg) - - # Advance the sample time - - startdate = AdvanceTime( startdate, cyclelen) - DaCycle['time.sample.start'] = startdate - -def Invert(DaCycle, StateVector, Optimizer ): - """ Perform the inverse calculation """ - - nobs = len(StateVector.EnsembleMembers[-1][0].ModelSample.Data.getvalues('obs')) - - dims = ( int(DaCycle['time.nlag']), - int(DaCycle['forecast.nmembers']), - int(DaCycle.DaSystem['nparameters']), - nobs, ) - - dummy = Optimizer.Initialize(dims) - - dummy = Optimizer.StateToMatrix(StateVector ) - - dummy = Optimizer.SetLocalization('None') - - if not DaCycle.DaSystem.has_key('opt.algorithm'): - msg = "There was no minimum least squares algorithm specified in the DA System rc file (key : opt.algorithm)" ; logging.info(msg) - msg = "...using serial algorithm as default..." ; logging.info(msg) - dummy = Optimizer.SerialMinimumLeastSquares() - elif DaCycle.DaSystem['opt.algorithm'] == 'serial': - msg = "Using the serial minimum least squares algorithm to solve ENKF equations" ; logging.info(msg) - dummy = Optimizer.SerialMinimumLeastSquares() - elif DaCycle.DaSystem['opt.algorithm'] == 'bulk': - msg = "Using the bulk minimum least squares algorithm to solve ENKF equations" ; logging.info(msg) - dummy = Optimizer.BulkMinimumLeastSquares() - - - dummy = Optimizer.MatrixToState(StateVector ) - - StateVector.isOptimized = True - - return None - -def Advance( DaCycle, Samples, StateVector, ObservationOperator): - """ Advance the filter state to the next step """ - - # This is the advance of the modeled CO2 state. Optionally, routines can be added to advance the state vector (mean+covariance) - - # Then, restore model state from the start of the filter - - dummy = DaCycle.MoveSaveData(io_option='restore',save_option='partial',filter=[]) - msg = "All restart data have been recovered from the save/tmp directory, this "+ \ - "resets the filter to %s "%DaCycle['time.start'] ; logging.debug(msg) - - msg = "Sampling model will be run over 1 cycle" ; logging.info(msg) - - SampleOneCycle(DaCycle,Samples,StateVector, ObservationOperator,0) - - return None - -def SaveAndSubmit( DaCycle, StateVector): - """ Save the model state and submit the next job """ - - savedir = DaCycle['dir.output'] - filename = os.path.join(savedir,'savestate.nc') - - dummy = StateVector.WriteToFile(filename) - dummy = DaCycle.Finalize() - - return None - -def RunForecastModel(DaCycle,ObsOperator): - """Prepare and execute a forecast step using an external Fortran model. Note that the flavor of model - used is determined in the very first line where the import statement of module "model" depends on a - setting in your da.rc file. After that choice, the module written specifically for a particular - application takes over. There are a few required variables and methods in such a module: - - version [variable, string] : defines the version number of the module, e.g. "1.0" - identifier [variable, string] : identifies the model used as a string variable, e.g. "TM5" - PrepareExe (method) : A method that preps the model simulation. It can for instance create an input - parameter file for the model (tm5.rc), or execute some scripts needed before the - actual execution starts (get_meteo). After this step, a call to the model - executable should start a succesfull simulation - Run (method) : Start the executable. How this is done depends on your model and might involve - submitting a string of jobs to the queue, or simply spawning a subprocess, or ... - """ - - status = ObsOperator.Initialize(DaCycle) - - status = ObsOperator.ValidateInput(DaCycle) - - status = ObsOperator.Run(DaCycle) - - dummy = ObsOperator.SaveData() - - return status - -def WriteEnsembleData(DaCycle, StateVector, lag): - """ Write the data needed by each model ensemble member to disk""" - members = StateVector.EnsembleMembers[lag] - for mem in members: - - mem.WriteToFile(DaCycle) - - msg = "Ensemble member parameter files were written to disk for lag = %d" % lag ; logging.info(msg) - return None - - diff --git a/da/tools/rc.py b/da/tools/rc.py deleted file mode 100755 index 02a1047d915e2a81052e7a963eca5b31bb19d50b..0000000000000000000000000000000000000000 --- a/da/tools/rc.py +++ /dev/null @@ -1,1147 +0,0 @@ -#! /usr/bin/env python -# rc.py - - -# ------------------------------------------------ -# help -# ------------------------------------------------ - -""" -Deal with model settings in `rc` format. - -RCFILES - - A rcfile is a text file with key/value pairs seperated by a ':', e.g. - - my.flag : T - my.answer : 42 - - The following functionality is supported: - - * Empty lines are ignored. - - * Comment lines are introduced by a '!' as first character. - - * Long values could be continued at the next line after a '\' as last character. - - * Include the key/value pairs from another file: - - #include an/other.rc - - * Substitute environment variables when available: - - tm.dir : ${HOME}/TM5/cy3 - - * Substitute values of other keys in a line: - - build.dir : ${tm.dir}/build - grid : glb300x200 - input.${grid}.path : /data/input/${grid} - - Substitions are allowed in both key names as well as values. - The substitutions are performed in a loop until nothing - has to be substituted anymore, or some substitutions could - not be applied at al; for the later an error is raised. - Values to be substituted could therefore be set before and - after they are used. - - Note that if a key has the same name as an environment variable, - the new value will be assigned to the key instead of the value - retrieved from the environment: - - HOME : /some/other/dir/ - - * Substitude some specials: - - ${pid} # evaluates to the current process id; - # useful for names of log files etc - ${script} # evaluates to the base name of the calling script, - # thus without .py etc - - * Instead of variables of the form '${..}' other patterns could be - specified with the optional 'marks' tupple (see below). - - * Old-style '#eval' lines are still supported: - - #eval RUNDIR = /path/to/mydir - tmdir : ${RUNDIR}/TM5/cy3 - - In this example, the value of RUNDIR will be evaluated and substituted - in all {key,value} pairs. This feature is obsolete and a warning will - be issued. The proper way to use this is with {key,value} pairs too: - - run.dir : /path/to/mydir - tmdir : ${run.dir}/TM5/cy3 - - * Comment starting with '!' is stripped from the values. - To have a value including exclamation marks, use '\!' but do - not expect that the rest of the value is scanned for comment too: - - my.value : -999 ! just an integer value - my.message : This value has 64 characters \! Count if you don't believe it ... - - * If you trust yourself you might try to use conditional expressions: - - #if ${my.number} == 1 - message : Welcome - #else - message : Whatever ... - #endif - - The conditions should be valid python expressions that evaluate to a boolean; - value substitutions are performed before evaluation. Examples: - - ${my.runmode} == 4 - "${my.tracer}" == "CH4" - - Keep it simple! Very complicated and nested if-statements might not be - resolved correctly, and are in any case not easy to understand for other users! - - In the example above, an exception could be raised by the special error expression; - everything behind the '#error' mark is displayed as an error message: - - #error No settings provided for value : ${my.value} - - -USAGE AS SCRIPT - - Called in script form, the following syntaxis is supported: - - rc.py [options] <rcfile> <key> - rc.py -h|--help - - The <rcfile> is read and the value defined by <key> is printed - to the standard output. - - Use the --help option for more documentation. - - -USAGE AS PYTHON MODULE - - Import the module with: - - import rc - - Initialiase by reading all settings in a rcfile, - supporting the functionality described in the 'RCFILES' section. - - rcf = RcFile( 'settings.rc' ) - - The initialisation accepts some optional arguments. - Set the silent flag to True to ignore warnings. - - rcf = RcFile( 'settings.rc', silent=False ) - - Use the optional 'marks' tupple to define that variables to be expanded - are marked other than '${..}' but rather '<mark1>..<mark2>' : - - rcf = RcFile( 'settings.rc', marks=('${','}') ) - - Test to see if a key is defined: - - if rcf.has_key('my.flag') : - print 'value of my flag is : ', rcf['my.flag'] - - Extract a list with all keys: - - rcf.keys() - - A 'get' function is provided to extract values: - - * by default, the 'get' function returns the value as a str type: - - s = rcf.get('my.value') - - * a second argument is the name of the python type to which - the value is converted to: - - i = rcf.get('my.flag','int') - - * if the return value should be a 'bool', the result is - True for values : 'True' , 'T', 'yes', or '1' , - and False for value : 'False', 'F', 'no' , or '0' ; - for other values an error is raised; - - * return a default value if the key is not found: - - rcf.get( 'my.flag', default=False ) - - * print a debug message to the logging system for each extracted key: - - rcf.get( 'my.flag', verbose=True ) - - Add a new value, comment is optional: - - rcf.add( 'my.iter', 2, comment='iteration number for restart' ) - - Assign a new value to an existing key: - - rcf.replace( 'my.flag', True ) - - Scan a character line for all occurances of ${<key>} and subsitute for - the rc value assigned to <key> : - - line = rcf.substitute( line ) - - Write the dictionary (with all variables expanded and included files included) - to new file: - - rcf.write('newfile.rc') - - -HISTORY - - 2008? Andy Jacobson, NOAA - Translation to python of original shell script 'go_readrc' . - 2009-06 Wouter Peters, WUR - Support substitution of previously defined variables. - 2009-06 Arjo Segers, TNO - Support include files. - 2009-09 Arjo Segers, TNO - Re-coded into class. - Implemented substitution loop. - 2009-11 Arjo Segers, JRC - Added main program to run this file as a shell script. - Added replace and substitute routines. - 2010-03 Arjo Segers, JRC - Support simple if-statements. - Support comment in values. - -""" - - -# ------------------------------------------------ -# classes -# ------------------------------------------------ - - -class RcFile(object) : - - """ - Class to store settings read from a rcfile. - """ - - def __init__( self, filename, silent=False, marks=('${','}') ) : - - """ - - Usage: - - rcf = RcFile( 'settings.rc' [,silent=False] [marks=('${','}')] ) - - Read an rc-file and expand all the keys and values into a dictionary. - Do not shout messages if silent is set to True. - The 2-item tupple (mark1,mark2) could be used to re-define the default - key pattern '${..}' into something else: - <mark1>...<mark2> - - """ - - # external: - import re - import os - import sys - import logging - - # info ... - logging.debug( 'reading rcfile %s ...' % filename ) - - # check ... - if not os.path.exists(filename) : - msg = 'rcfile not found : %s' % filename ; logging.error(msg) - raise IOError, msg - #endif - - # store file name: - self.filename = filename - # store rc-file root directory: - self.rootdir = os.path.split(filename)[0] - - # storage for processed rcfile: - self.outfile = [] - - # storage for key/value pairs: - self.values = {} - - # open the specified rc-file: - f = open(filename,'r') - # store all lines in a list: - inpfile = f.readlines() - # close: - f.close() - - # flags: - warned_for_eval = False - - # pass counter: - ipass = 1 - - # loop until all substitutions and inclusions are done: - while True : - - # start again with empty output file: - self.outfile = [] - # init current line: - line = '' - # assume nothing has to be done after this loop: - something_done = False - something_to_be_done = False - unresolved = [] - - # stack for conditional evaluation; - # each element is a tuple with elements: - # resolved (boolean) true if the if-statement is evaluated - # flag (boolean) true if the lines below the statement - # are to be included - # anyflag (boolean) used to check if any of the 'if' or 'elif' conditions - # in this sequence evaluated to True - # line (char) description of the line for messages - ifstack = [] - - #print '' - #print '---[pass %i]-------------------------------------' % ipass - #for line in inpfile : print line.strip() - - # loop over lines in input file: - iline = -1 - for inpline in inpfile : - - # line counter: - iline = iline + 1 - - # remove end-of-line character: - inpline = inpline.strip() - - ## DEBUG: display current line ... - #print '%4i | %s' % (iline,inpline) - - # - # empty lines - # - - # skip empty lines: - if len(inpline) == 0 : - # add empty line to output: - self.outfile.append('\n') - # next will be a new line: - line = '' - # next input line: - continue - #endif - - # - # comment lines - # - - # skip comment: - if inpline.startswith('!') : - # add copy to output file: - self.outfile.append( '%s\n' % inpline ) - # next will be a new line: - line = '' - # next input line: - continue - #endif - - # - # continuation lines - # - - # current line has continuation mark '\' at the end ? - # then add this input line: - if line.endswith('\\') : - # remove continuation character, add input line: - line = line[:-1]+inpline - else : - # bright new line: - line = inpline - #endif - - # continuation mark ? then next line of input file: - if line.endswith('\\') : continue - - # - # line info - # - - # line number and text for messages: - line_info = '%6i | %s' % (iline+1,line) - - # - # conditional settings (1) - # - - # is this the begin of a new condition ? - mark = '#if' - if line.startswith(mark) : - # push temporary flag to stack, will be replaced after evaluation of condition: - ifstack.append( ( False, True, False, line_info ) ) - # debug ... - #print 'xxx1 ', ifstack - #endif - - mark = '#elif' - if line.startswith(mark) : - # check ... - if len(ifstack) == 0 : - logging.error( 'found orphin elif in rcfile on line :' ) - logging.error( ' %s' % line_info ) - raise Exception - #endif - # remove current top from stack: - resolved,flag,anyflag,msg = ifstack.pop() - # push temporary flag to stack, will be replaced after evaluation of condition: - ifstack.append( ( resolved, True, anyflag, line_info ) ) - # debug ... - #print 'xxx1 ', ifstack - #endif - - mark = '#else' - if line.startswith(mark) : - # check ... - if len(ifstack) == 0 : - logging.error( 'found lonely else in rcfile on line :' ) - logging.error( ' %s' % line_info ) - raise Exception - #endif - # remove current top from stack: - resolved,flag,anyflag,msg = ifstack.pop() - # get higher level settings: - if len(ifstack) > 0 : - resolved_prev,flag_prev,anyflag_prev,msg_prev = ifstack[-1] - else : - flag_prev = True - #endif - # should next lines be included ? - new_flag = (not flag) and (not anyflag) and flag_prev - # add else block with reversed flag, take into acount higher level - ifstack.append( ( resolved, new_flag, False, line_info ) ) - # debug ... - #print 'xxx1 ', ifstack - # copy to output: - self.outfile.append( '%s\n' % line ) - # next input line: - continue - #endif - - # is this the end of a condition ? - mark = '#endif' - if line.startswith(mark) : - # check ... - if len(ifstack) == 0 : - logging.error( 'found lonely endif in rcfile on line :' ) - logging.error( ' %s' % line_info ) - raise Exception - #endif - # remove top from stack: - top = ifstack.pop() - # copy to output: - self.outfile.append( '%s\n' % line ) - # next input line: - continue - #endif - - # within if-statements ? - if len(ifstack) > 0 : - # extract current top: - resolved,flag,anyflag,msg = ifstack[-1] - # already resolved ? then check if this line should be skipped: - if resolved and (not flag) : - # skip this line, but include commented version in output: - self.outfile.append( '!%s\n' % line ) - # read the next input line: - continue - #endif - #endif - - # - # handle '#eval' lines - # - - mark = '#eval' - if line.startswith(mark): - # info .. - if not warned_for_eval : - if not silent: logging.warning( 'the #eval statements in rc-files are deprecated, use {key:value} pairs instead' ) - warned_for_eval = True - #endif - # add commented copy to output: - self.outfile.append( '!evaluated>>> '+line ) - # remove leading mark: - line = line.lstrip(mark) - # multiple settings are seperated by ';' : - evals = line.split(';') - # insert in output file: - for k in range(len(evals)) : - # split in key and value: - key,value = evals[k].split('=') - # insert: - self.outfile.append( '%s : %s' % (key,value) ) - #endfor - # next input line: - continue - #endif - - # - # replace ${..} values - # - - # ensure that common marks are evaluated correctly: - start_mark = marks[0].replace('{','\{').replace('<','\<').replace('$','\$') - close_mark = marks[1].replace('}','\}').replace('>','\>') - - # set syntax of keywords to be matched, e.g. '${...}' : - pattern = start_mark+'[A-Za-z0-9_.]+'+close_mark - - # make a regular expression that matches all variables: - rc_varpat = re.compile( pattern ) - - # search all matching paterns: - pats = re.findall(rc_varpat,line) - # counter for unexpanded substitutions: - ntobedone = 0 - # loop over matches: - for pat in pats : - # remove enclosing characters: - key = pat.lstrip(start_mark).rstrip(close_mark) - # test some dictionaries for matching key: - if self.values.has_key(key) : - # get previously defined value: - val = self.values[key] - # substitute value: - line = line.replace(pat,val) - # set flag: - something_done = True - elif os.environ.has_key(key) : - # get value from environment: - val = os.environ[key] - # substitute value: - line = line.replace(pat,val) - # set flag: - something_done = True - elif key == 'pid' : - # special value: process id; convert to character: - val = '%i' % os.getpid() - # substitute value: - line = line.replace(pat,val) - # set flag: - something_done = True - elif key == 'script' : - # special value: base name of the calling script, without extension: - script,ext = os.path.splitext(os.path.basename(sys.argv[0])) - # substitute value: - line = line.replace(pat,script) - # set flag: - something_done = True - else : - # could not substitute yet; set flag: - ntobedone = ntobedone + 1 - # continue with next substitution: - continue - #endif - #endfor # matched patterns - # not all substituted ? - if ntobedone > 0 : - # add line to output: - self.outfile.append(line) - # a new pass is needed: - something_to_be_done = True - # store for info messages: - unresolved.append(line) - # next input line: - continue - #endif - - # - # handle '#include' lines - # - - mark = '#include' - if line.startswith(mark) : - # remove leading mark, what remains is file to be included: - inc_file = line.lstrip(mark).strip() - # check ... - if not os.path.exists(inc_file) : - inc_file = os.path.join(self.rootdir,inc_file) - logging.debug( 'Added rootdir to requested include: %s' % inc_file ) - - if not os.path.exists(inc_file) : - logging.error( 'include file not found : %s' % inc_file ) - msg = 'ERROR - include file not found : %s' % inc_file - raise IOError,msg - #endif - # read content: - inc_f = open( inc_file, 'r' ) - inc_rc = inc_f.readlines() - inc_f.close() - # add extra comment for output file: - self.outfile.append( '! >>> %s >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n' % inc_file ) - self.outfile.extend( inc_rc ) - self.outfile.append( '! <<< %s <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n' % inc_file ) - # set flag: - something_done = True - # a new pass is needed: - something_to_be_done = True - # next input line: - continue - #endif - - - # - # conditional settings (2) - # - - # evaluate conditional expressions: - mark1 = '#if' - mark2 = '#elif' - if line.startswith(mark1) or line.startswith(mark2) : - # remove leading mark, what remains is logical expression: - expression = line.lstrip(mark1).strip() - expression = line.lstrip(mark2).strip() - # common mistake is to add a ':' as in python; remove this: - if expression.endswith(':') : expression = expression.rstrip(':').strip() - # evaluate: - try : - flag = eval( expression ) - except : - logging.error( 'could not evaluate expression:' ) - logging.error( ' %s' % expression ) - logging.error( 'on line:' ) - logging.error( line_info ) - sys.exit(1) - #endtry - # remove temporary top added before during this pass: - tmp_statement,tmp_flag,tmp_anyflag,tmp_msg = ifstack.pop() - # extract current top if necessary: - if len(ifstack) > 0 : - dummy_statement,prev_flag,dummy_anyflag,dummy_msg = ifstack[-1] - else : - prev_flag = True - #endif - # should next lines be included ? - new_flag = prev_flag and tmp_flag and flag - # any if/elif evaluated to true in this sequence ? - new_anyflag = tmp_anyflag or new_flag - # add to stack, now resolved, take into accout current flag: - ifstack.append( ( True, new_flag, new_anyflag, line_info ) ) - # debug ... - #print 'xxx2 ', ifstack - # copy to output: - self.outfile.append( '%s\n' % line ) - # next input line: - continue - #endif - - # - # error message - # - - # special command to rais an exception: - mark = '#error' - if line.startswith(mark) : - # remove leading mark, what remains is error message: - msg = line.lstrip(mark).strip() - # display: - logging.error( msg ) - # add info: - logging.error( 'error message found on line:' ) - logging.error( line_info ) - # stop: - raise Exception - #endif - - # - # checks - # - - # common mistake ... - if line.startswith('#') : - logging.error( 'line in rcfile starts with "#" but is not an "#include" or other special line;' ) - logging.error( 'if it is supposed to be comment, please start with "!" ...' ) - logging.error( ' rcfile : %s' % filename ) - logging.error( ' line : %s' % line ) - raise IOError - #endif - - # check ... - if ':' not in line : - logging.error( 'key/value line should contain a ":"' ) - logging.error( ' rcfile : %s' % filename ) - logging.error( ' line : %s' % line ) - raise IOError - #endif - - # - # add to output - # - - # add line to output: - self.outfile.append( '%s\n' % line ) - - # - # add key/value pair - # - - # not if inside an unresolved if-statement ... - if len(ifstack) > 0 : - # get top values: - resolved,flag,anyflag,msg = ifstack[-1] - # not resolved yet ? then continue: - if not resolved : continue - #endif - - # split in key and value; - # value might contain ':' too, so at maximum 1 split: - key,val = line.split(':',1) - - # remove comment from value: - if '!' in val : - # not if '\!' is in the value ... - if not '\!' in val : val,comment = val.split('!') - # replace all slash-comments: - val = val.replace('\!','!') - #endif - - # remove spaces: - key = key.strip() - val = val.strip() - - # already defined ? - if self.values.has_key(key) : - # no problem if values are the same, but otherwise ... - if self.values[key] != val : - logging.error( 'key found twice in "%s" :' % filename ) - logging.error( ' %s : %s' % (key,str(self.values[key])) ) - logging.error( ' %s : %s' % (key,str(val)) ) - raise Exception - #endif - else : - # store new value: - self.values[key] = val - # set flag: - something_done = True - #endif - - # display key and value ... - #print ' --> %s : %s' % (key,val) - - #endfor # loop over lines in text - - ## info ... - #print '~~~ outfile ~~~~~~~~~~~~~~~~~~~~~~~' - #for line in self.outfile : print line.strip() - #print '~~~ key/values ~~~~~~~~~~~~~~~~~~~~' - #for k,v in self.iteritems() : - # print '%s : %s' % (k,v) - ##endfor - #print '-------------------------------------------------' - #print '' - - # check ... - if len(ifstack) > 0 : - logging.error( 'unterminated if-statement ; current stack:' ) - for resolved,flag,anyflag,msg in ifstack : logging.error( msg ) - logging.error( 'please fix the rcfile or debug this script ...' ) - raise Exception - #endif - - # check ... - if something_to_be_done : - # check for unterminated loop ... - if not something_done : - logging.error( 'could not resolve the following lines in rcfile:' ) - for uline in unresolved : logging.error( ' %s' % uline ) - logging.error( 'please fix the rcfile or debug this script ...' ) - raise Exception - #endif - else : - # finished ... - break - #endif - - # for safety ... - if ipass == 100 : - logging.error( 'resolving rc file has reached pass %i ; something wrong ?' % ipass ) - raise Exception - #endif - - # new pass: - ipass = ipass + 1 - # renew input: - inpfile = self.outfile - - #endwhile # something to be done - - #enddef # __init__ - - - # *** - - - def has_key( self, key ) : - - # from dictionairy: - return self.values.has_key(key) - - #enddef - - - # *** - - - def keys( self ) : - - # from dictionairy: - return self.values.keys() - - #enddef - - - # *** - - - def get( self, key, totype='', default=None, verbose=False ) : - - """ - rcf.get( 'my.value' [,default=None] ) - Return element 'key' from the dictionairy. - If the element is not present but a default is specified, than return - the default value. - If 'verbose' is set to True, then print debug messages to the logging - about which values is returned for the given key. - The option argument 'totype' defines the conversion to a Python type. - If 'totype' is set to 'bool', the return value is the - boolean True for values 'T', 'True', 'yes', and '1', - and False for 'F', 'False', 'no', or '0' ; - for other values, an exception will be raised. - """ - - # external: - import logging - - # element found ? - if self.values.has_key(key) : - # copy value: - value = self.values[key] - # convert ? - if totype == 'bool' : - # convert to boolean: - if value in ['T','True','yes','1'] : - value = True - elif value in ['F','False','no','0'] : - value = False - else : - logging.error( "value of key '%s' is not a boolean : %s" % (key,str(value)) ) - raise Exception - #endif - elif len(totype) > 0 : - # convert to other type ... - value = eval( '%s(%s)' % (totype,value) ) - #endif - # for debugging ... - if verbose : logging.debug( 'rc setting "%s" : "%s"' % (key,str(value)) ) - else : - # default value specified ? - if default != None : - # copy default: - value = default - # for debugging ... - if verbose : logging.debug( 'rc setting "%s" : "%s" (deault)' % (key,str(value)) ) - else : - # something wrong ... - logging.error( "key '%s' not found in '%s' and no default specified" % (key,self.filename) ) - raise Exception - #endif - #endif - - # ok - return value - - #enddef - - - # *** - - - def replace( self, key, val ) : - - """ - Replace a key by a new value. - """ - - # external: - import logging - - # search for a line '<key> : <val>' - # loop over lines in output file: - found = False - for iline in range(len(self.outfile)) : - # extract: - line = self.outfile[iline] - # skip lines that are no key:value pair for sure ... - if ':' not in line : continue - # split once at first ':' - k,v = line.split(':',1) - # match ? - if k.strip() == key : - # replace line in original file: - self.outfile[iline] = '%s : %s\n' % (k,str(val)) - # replace value: - self.values[key] = val - # set flag: - found = True - # found, thus no need to continue: - break - #endif - #endfor # lines - # not found ? - if not found : - logging.error( 'could not replace key : %s' % key ) - raise Exception - #endif - - # ok - return - - #enddef - - - # *** - - - def add( self, key, val, comment='' ) : - - """Add a new key/value pair.""" - - # add lines: - self.outfile.append( '\n' ) - if len(comment) > 0 : self.outfile.append( '! %s\n' % comment ) - self.outfile.append( '%s : %s\n' % (key,str(val)) ) - - # add to dictionairy: - self.values[key] = val - - # ok - return - - #enddef - - - # *** - - - def substitute( self, line, marks=('${','}') ) : - - """ - Return a line with all '${..}' parts replaced by the corresponding rcfile values. - The 2-item tupple (mark1,mark2) could be used to re-define the default - key pattern '${..}' into something else: - <mark1>...<mark2> - """ - - # external: - import re - - # ensure that common marks are evaluated correctly: - start_mark = marks[0].replace('{','\{').replace('<','\<').replace('$','\$') - close_mark = marks[1].replace('}','\}').replace('>','\>') - - # set syntax of keywords to be matched, e.g. '${...}' : - pattern = start_mark+'[A-Za-z0-9_.]+'+close_mark - - # make a regular expression that matches all variables: - rc_varpat = re.compile( pattern ) - - # search all matching paterns: - pats = re.findall(rc_varpat,line) - # loop over matches: - for pat in pats : - # remove enclosing characters: - key = pat.lstrip(start_mark).rstrip(close_mark) - # test dictionary for matching key: - if self.values.has_key(key) : - # get previously defined value: - val = self.values[key] - # substitute value: - line = line.replace(pat,val) - #endif - #endfor # matched patterns - - # ok - return line - - #enddef - - - # *** - - - def WriteFile( self, filename ) : - - """ write the dictionary to file""" - - # open file for writing: - f = open(filename,'w') - - ## loop over key/value pairs: - #for k,v in self.iteritems(): - # # add line; at least the specified number of characters - # # is used for the key: - # f.write( '%-20s:%s\n' % (k,v) ) - ##endfor - - # write processed input: - f.writelines( self.outfile ) - - # close file: - f.close() - - #endif - - -#endclass # RcFile - -def read(rcfilename,silent=False): - """ - This method reads an rc-file by making an instance of the RcFile class, and then returns the dictionary of values only. This - makes it backwards compatible with older implementations of the rc.py module - """ - - rcdict = RcFile(rcfilename,silent=silent) - - return rcdict.values - -def write( filename, rcdict): - """ - This method writes an rc-file dictionary. This is included to make this module backwards compatible with - older implementations of the rc.py module - """ - - # open file for writing: - f = open(filename,'w') - - # loop over key/value pairs: - for k,v in rcdict.items(): - # add line; at least the specified number of characters - # is used for the key: - f.write( '%-20s:%s\n' % (k,v) ) - #endfor - - # close file: - f.close() - - - -# ------------------------------------------------ -# test -# ------------------------------------------------ - - -if __name__ == '__main__': - - # external ... - import sys - import optparse - import logging - - # extract arguments from sys.argv array: - # 0 = name of calling script, 1: = actual arguments - args = sys.argv[1:] - - # set text for 'usage' help line: - usage = "\n %prog <rcfile> <key> [-b|--bool] [--default<=value>]\n %prog <rcfile> -w|--write\n %prog -h|--help\n %prog -d|--doc" - - # initialise the option parser: - parser = optparse.OptionParser(usage=usage) - - # define options: - parser.add_option( "-d", "--doc", - help="print documentation", - dest="doc", action="store_true", default=False ) - parser.add_option( "-b", "--bool", - help="return 'True' for values 'T', 'True', 'yes', or '1', and 'False' for 'F', 'False', 'no', or '0'", - dest="boolean", action="store_true", default=False ) - parser.add_option( "--default", - help="default value returned if key is not found", - dest="default", action="store" ) - parser.add_option( "-w", "--write", - help="write pre-processed rcfile", - dest="write", action="store_true", default=False ) - - # now parse the actual arguments: - opts,args = parser.parse_args( args=args ) - - # print documentation ? - if opts.doc : - print __doc__ - sys.exit(0) - #endif - - # recfile argument should be provided: - if len(args) < 1 : - parser.error("no name of rcfile provided\n") - #endif - # extract: - rcfile = args[0] - - # read rcfile in dictionary: - try : - rcf = RcFile(rcfile) - except : - logging.error( sys.exc_info()[1] ) - sys.exit(1) - #endtry - - # print pre-processed file ? - if opts.write : - for line in rcf.outfile : print line.strip() - sys.exit(0) - #endif - - # key argument should be provided: - if len(args) < 2 : - parser.error("no name of rckey provided\n") - #endif - # extract: - rckey = args[1] - - # key present ? - if rcf.has_key(rckey) : - - # print requested value: - if opts.boolean : - # extract value: - flag = rcf.get(rckey,'bool') - # print result: - if flag : - print 'True' - else : - print 'False' - #endif - else : - # extract value: - value = rcf.get(rckey) - # display: - print value - #endif - - else : - - # default value provided ? - if opts.default != None : - # display: - print opts.default - else : - print 'ERROR - key "%s" not found in rcfile "%s" and no default specified' % (rckey,rcfile) - sys.exit(1) - #endif - - #endif - -#endif - - -# ------------------------------------------------ -# end -# ------------------------------------------------ - diff --git a/da/tools/standardvariables.py b/da/tools/standardvariables.py deleted file mode 100755 index 9eac8f64499b0e5e93c6f65c5de4e2f30681d877..0000000000000000000000000000000000000000 --- a/da/tools/standardvariables.py +++ /dev/null @@ -1,188 +0,0 @@ -standard_variables = { 'bio_flux_prior' : {'name' : 'bio_flux_prior',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, terrestrial vegetation, not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'bio_flux_opt' : {'name' : 'bio_flux_opt',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, terrestrial biosphere , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_prior' : {'name' : 'ocn_flux_prior',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, open ocean , not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_opt' : {'name' : 'ocn_flux_opt',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, open ocean , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'fossil_flux_imp' : {'name' : 'fossil_flux_imp',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, fossil fuel burning , imposed ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'fire_flux_imp' : {'name' : 'fire_flux_imp',\ - 'units' : 'mol m-2 s-1' ,\ - 'long_name' : 'Surface flux of carbon dioxide, biomass burning , imposed ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'surface_carbon_dioxide_mole_flux', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'bio_flux_prior_cov' : {'name' : 'bio_flux_prior_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, terrestrial vegetation , not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'bio_flux_opt_cov' : {'name' : 'bio_flux_opt_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, terrestrial vegetation , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_prior_cov' : {'name' : 'ocn_flux_prior_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, open ocean , not optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ocn_flux_opt_cov' : {'name' : 'ocn_flux_opt_cov',\ - 'units' : 'mol2 region-2 s-2' ,\ - 'long_name' : 'Covariance of surface flux of carbon dioxide, open ocean , optimized ', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : '', \ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'decimal_date' : {'name' : 'decimal_date',\ - 'units' : 'years' ,\ - 'long_name' : 'dates and times', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'date', \ - 'dims' : (), \ - 'dtype' : 'double', \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'date' : {'name' : 'date',\ - 'units' : 'days since 2000-01-01 00:00:00 UTC' ,\ - 'long_name' : 'UTC dates and times', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'standard_name' : 'date', \ - 'dims' : (), \ - 'dtype' : 'double', \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'idate' : {'name' : 'idate',\ - 'units' : 'yyyy MM dd hh mm ss ' ,\ - 'long_name' : 'integer components of date and time', \ - 'standard_name' : 'calendar_components', \ - 'comment' : 'time-interval average, centered on times in the date axis', \ - 'dims' : (), \ - 'dtype' : 'int', \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'latitude' : {'name' : 'latitude',\ - 'units' : 'degrees_north ' ,\ - 'long_name' : 'latitude', \ - 'standard_name' : 'latitude', \ - 'comment' : 'center of interval',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'longitude' : {'name' : 'longitude',\ - 'units' : 'degrees_east ' ,\ - 'long_name' : 'longitude', \ - 'standard_name' : 'longitude', \ - 'comment' : 'center of interval',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'height' : {'name' : 'height',\ - 'units' : 'masl ' ,\ - 'long_name' : 'height_above_ground_level', \ - 'standard_name' : 'height_above_ground_level', \ - 'comment' : 'value is meters above sea level',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'co2' : {'name' : 'co2',\ - 'units' : 'micromol mol-1 ' ,\ - 'long_name' : 'mole_fraction_of_carbon_dioxide_in_air', \ - 'standard_name' : 'mole_fraction_of_carbon_dioxide_in_air', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'meanstate' : {'name' : 'statevectormean',\ - 'units' : 'unitless' ,\ - 'long_name' : 'mean_value_of_state_vector', \ - 'standard_name' : 'mean_value_of_state_vector', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'ensemblestate': {'name' : 'statevectorensemble',\ - 'units' : 'unitless' ,\ - 'long_name' : 'ensemble_value_of_state_vector', \ - 'standard_name' : 'ensemble_value_of_state_vector', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - 'unknown' : {'name' : '',\ - 'units' : '' ,\ - 'long_name' : '', \ - 'standard_name' : '', \ - 'comment' : '',\ - 'dims' : (), \ - 'values' : [], \ - 'count' : 0 \ - } , \ - } - - - - diff --git a/dajet.rc b/dajet.rc deleted file mode 100644 index ab30afb5c9ad431c8534e0df6b00d1c96b126cbc..0000000000000000000000000000000000000000 --- a/dajet.rc +++ /dev/null @@ -1,19 +0,0 @@ -! Info on the data assimilation cycle - -time.restart : False -time.start : 2000-01-01 00:00:00 -time.finish : 2000-01-15 00:00:00 -time.cycle : 7 -time.nlag : 2 -dir.da_run : ${HOME}/tmp/test_da - -! Info on the DA system used - -da.system : CarbonTracker -da.system.rc : carbontrackerjet.rc - -! Info on the forward model to be used - -forecast.model : TM5 -forecast.model.rc : ${HOME}/TM/TM5_new/ct_new.rc -forecast.nmembers : 4 diff --git a/das.py b/das.py deleted file mode 100755 index 3e72236dadfac65326b309f550121a6e600ae608..0000000000000000000000000000000000000000 --- a/das.py +++ /dev/null @@ -1,47 +0,0 @@ -#!/usr/bin/env python -# das.py - -#$ -N das -#$ -A co2 -#$ -pe nserial 1 -#$ -l h_rt=00:30:00 -#$ -S /bin/sh -#$ -cwd -#$ -j y -#$ -r n - -from da.tools.initexit import StartLogger -from da.tools.pipeline import Main - -dummy = StartLogger() - -########################################################################################### -### IMPORT THE APPLICATION SPECIFIC MODULES HERE, TO BE PASSED INTO THE MAIN PIPELINE!!! ## -########################################################################################### - -from da.platform.maunaloa import MaunaloaPlatForm -from da.ct.dasystem import CtDaSystem -from da.ct.statevector import CtStateVector -from da.ct.obs import CtObservations -from da.tm5.observationoperator import TM5ObservationOperator -from da.ct.optimizer import CtOptimizer - -PlatForm = MaunaloaPlatForm() -DaSystem = CtDaSystem('carbontrackerjet.rc') -StateVector = CtStateVector() -Samples = CtObservations() -ObsOperator = TM5ObservationOperator() -Optimizer = CtOptimizer() - -########################################################################################## -################### ENTER THE PIPELINE WITH THE OBJECTS PASSED BY THE USER ############### -########################################################################################## - - -print "\n ********************************************************************************************************" -print " *************************************** Entering Pipeline ******************************************" -print " ********************************************************************************************************\n" - -Main(PlatForm, DaSystem, Samples,StateVector,ObsOperator,Optimizer) - - diff --git a/das_with_jobs.py b/das_with_jobs.py deleted file mode 100755 index f9a110630e96e0a3141c1c97bb08e9ac33dd0380..0000000000000000000000000000000000000000 --- a/das_with_jobs.py +++ /dev/null @@ -1,129 +0,0 @@ -#!/usr/bin/env python -# das.py - -""" -Author : peters - -Revision History: -File created on 29 Sep 2009. - -""" -if __name__ == "__main__": - import sys - import os - import logging - import shutil - -# Append current working dir to path - - sys.path.append(os.getcwd()) - -# Import methods and classes contained in this package - - from da.tools.initexit import CleanUpCycle - from da.tools.initexit import ValidateOptsArgs - from da.tools.initexit import ParseOptions - from da.tools.initexit import StartLogger - from da.tools.pipeline import * - from da.tools.general import StoreData,RestoreData - -# Start a logger for all that happens from here on - - dummy = StartLogger() - -# Parse options from the command line - - opts, args = ParseOptions() - -# Validate Options and arguments passed - - opts,args = ValidateOptsArgs(opts,args) - - process = args['process'].lower() - -# Create the Cycle Control object for this job - - DaCycle = JobStart(opts,args) - - -# Start the requested subprocesses - - if process == 'start': - - msg = header+"Initializing Da Cycle"+footer ; logging.info(msg) - - dummy = DaCycle.Initialize() - - msg = header+"starting JobInput"+footer ; logging.info(msg) - - StateVector = JobInput(DaCycle) - - savefilename = os.path.join(DaCycle.da_settings['dir.save'],'statevector.pickle') - - dummy = StoreData(StateVector,savefilename) - - dummy = DaCycle.SubmitSubStep('samplestate') - - sys.exit(0) - - - savefilename = os.path.join(DaCycle.da_settings['dir.save'],'statevector.pickle') - - if process == 'samplestate': - msg = header+"starting SampleState"+footer ; logging.info(msg) - - StateVector = RestoreData(savefilename) - - dummy = SampleState(DaCycle,StateVector) - - dummy = StoreData(StateVector,savefilename) - - dummy = DaCycle.SubmitSubStep('invert') - - dummy = CleanUpCycle(DaCycle) - - sys.exit(0) - - if process == 'invert': - msg = header+"starting Invert"+footer ; logging.info(msg) - - StateVector = RestoreData(savefilename) - - dummy = Invert(DaCycle, StateVector) - - StateVector = StoreData(StateVector,savefilename) - - dummy = DaCycle.SubmitSubStep('advance') - - dummy = CleanUpCycle(DaCycle) - - sys.exit(0) - - - if process == 'advance': - msg = header+"starting Advance"+footer ; logging.info(msg) - - StateVector = RestoreData(savefilename) - - dummy = Advance(DaCycle, StateVector) - - StateVector = StoreData(StateVector,savefilename) - - dummy = DaCycle.SubmitSubStep('done') - - dummy = CleanUpCycle(DaCycle) - - sys.exit(0) - - if process == 'done': - msg = header+"starting SaveAndSubmit"+footer ; logging.info(msg) - - StateVector = RestoreData(savefilename) - - dummy = SaveAndSubmit(DaCycle, StateVector) - - msg = "Cycle finished...exiting" ; logging.info(msg) - - dummy = CleanUpCycle(DaCycle) - - sys.exit(0)