Commit 1f31f4e1 authored by brunner's avatar brunner
Browse files

domain as 9 regions, fractional plant types

parent cee00f6a
! CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
! Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
! updates of the code. See also: http://www.carbontracker.eu.
!
! This program is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software Foundation,
! version 3. This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this
! program. If not, see <http://www.gnu.org/licenses/>.
! author: Wouter Peters
!
! This is a blueprint for an rc-file used in CTDAS. Feel free to modify it, and please go to the main webpage for further documentation.
!
! Note that rc-files have the convention that commented lines start with an exclamation mark (!), while special lines start with a hashtag (#).
!
! When running the script start_ctdas.sh, this /.rc file will be copied to your run directory, and some items will be replaced for you.
! The result will be a nearly ready-to-go rc-file for your assimilation job. The entries and their meaning are explained by the comments below.
!
!
! HISTORY:
!
! Created on August 20th, 2013 by Wouter Peters
!
!
! The time for which to start and end the data assimilation experiment in format YYYY-MM-DD HH:MM:SS
! the following 3 lines are for initial start
time.start : 2013-04-01 00:00:00
time.finish : 2013-04-07 23:00:00
time.end : 2013-04-07 23:00:00
abs.time.start : 2013-04-01 00:00:00
! Whether to restart the CTDAS system from a previous cycle, or to start the sequence fresh. Valid entries are T/F/True/False/TRUE/FALSE
time.restart : F
da.restart.tstamp : 2013-04-01 00:00:00
! The length of a cycle is given in days, such that the integer 7 denotes the typically used weekly cycle. Valid entries are integers > 1
time.cycle : 7
! The number of cycles of lag to use for a smoother version of CTDAS. CarbonTracker CO2 typically uses 5 weeks of lag. Valid entries are integers > 0
time.nlag : 2
! The directory under which the code, input, and output will be stored. This is the base directory for a run. The word
! '/' will be replaced through the start_ctdas.sh script by a user-specified folder name. DO NOT REPLACE
run.name : 9reg
dir.da_run : /scratch/snx3000/parsenov/${run.name}
restartmap.dir : ${dir.da_run}/input
! The resources used to complete the data assimilation experiment. This depends on your computing platform.
! The number of cycles per job denotes how many cycles should be completed before starting a new process or job, this
! allows you to complete many cycles before resubmitting a job to the queue and having to wait again for resources.
! Valid entries are integers > 0
da.resources.ncycles_per_job : 1
! The ntasks specifies the number of threads to use for the MPI part of the code, if relevant. Note that the CTDAS code
! itself is not parallelized and the python code underlying CTDAS does not use multiple processors. The chosen observation
! operator though might use many processors, like TM5. Valid entries are integers > 0
da.resources.ntasks : 1
! This specifies the amount of wall-clock time to request for each job. Its value depends on your computing platform and might take
! any form appropriate for your system. Typically, HPC queueing systems allow you a certain number of hours of usage before
! your job is killed, and you are expected to finalize and submit a next job before that time. Valid entries are strings.
da.resources.ntime : 44:00:00
! The resource settings above will cause the creation of a job file in which 2 cycles will be run, and 30 threads
! are asked for a duration of 4 hours
!
! Info on the DA system used, this depends on your application of CTDAS and might refer to for instance CO2, or CH4 optimizations.
!
da.system : CarbonTracker
! The specific settings for your system are read from a separate rc-file, which points to the data directories, observations, etc
da.system.rc : da/rc/carbontracker_cosmo.rc
! This flag should probably be moved to the da.system.rc file. It denotes which type of filtering to use in the optimizer
da.system.localization : CT2007
! Info on the observation operator to be used, these keys help to identify the settings for the transport model in this case
da.obsoperator : cosmo
!
! The TM5 transport model is controlled by an rc-file as well. The value below refers to the configuration of the TM5 model to
! be used as observation operator in this experiment.
!
da.obsoperator.home : /store/empa/em05/parsenov/cosmo_processing_chain
da.bio.input : /store/empa/em05/parsenov/cosmo_input/vprm/processed
da.bg.input : /store/empa/em05/parsenov/cosmo_input/icbc/processed
da.obsoperator.rc : ${da.obsoperator.home}/tm5-ctdas-ei-zoom.rc
!forward.savestate.exceptsam : TRUE
!
! The number of ensemble members used in the experiment. Valid entries are integers > 2
!
da.optimizer.nmembers : 21
nparameters : 181
! Finally, info on the archive task, if any. Archive tasks are run after each cycle to ensure that the results of each cycle are
! preserved, even if you run on scratch space or a temporary disk. Since an experiment can take multiple weeks to complete, moving
! your results out of the way, or backing them up, is usually a good idea. Note that the tasks are commented and need to be uncommented
! to use this feature.
! The following key identifies that two archive tasks will be executed, one called 'alldata' and the other 'resultsonly'.
!task.rsync : alldata onlyresults
! The specifics for the first task.
! 1> Which source directories to back up. Valid entry is a list of folders separated by spaces
! 2> Which destination directory to use. Valid entries are a folder name, or server and folder name in rsync format as below
! 3> Which flags to add to the rsync command
! The settings below will result in an rsync command that looks like:
!
! rsync -auv -e ssh ${dir.da_run} you@yourserver.com:/yourfolder/
!
!task.rsync.alldata.sourcedirs : ${dir.da_run}
!task.rsync.alldata.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.alldata.flags g -auv -e ssh
! Repeated for rsync task 2, note that we only back up the analysis and output dirs here
!task.rsync.onlyresults.sourcedirs : ${dir.da_run}/analysis ${dir.da_run}/output
!task.rsync.onlyresults.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.onlyresults.flags : -auv -e ssh
......@@ -50,7 +50,7 @@ class Optimizer(object):
self.nlag = dims[0]
self.nmembers = dims[1]
# self.nparams = dims[2]
self.nparams = 23
self.nparams = 181
self.nobs = dims[3]
self.create_matrices()
......@@ -109,9 +109,29 @@ class Optimizer(object):
for n in range(self.nlag):
samples = statevector.obs_to_assimilate[n]
members = statevector.ensemble_members[n]
# params_2d = members[0].param_values[:-1].reshape(9,22)
# params = np.zeros(shape=(9,20))
# for r in range(0,9):
# params[r,:] = np.delete(np.delete(params_2d[r,:],21),10)
# params = np.hstack((params.flatten(),members[0].param_values[-1])) # put back BG
# self.x[n * self.nparams:(n + 1) * self.nparams] = params
self.x[n * self.nparams:(n + 1) * self.nparams] = members[0].param_values
# params = np.zeros(shape=(self.nmembers,181))
# params_3d = np.zeros(shape=(self.nmembers,9,20))
# for m,mm in enumerate(members):
# params_temp = mm.param_values[:-1].reshape(9,22)
# params_bg = mm.param_values[-1]
# for r in range(0,9):
# params_3d[m,r,:] = np.delete(np.delete(params_temp[r,:],21),10)
# params[m,:] = np.hstack((params_3d[m,:,:].flatten(),params_bg))
# params = params_3d.reshape(self.nmembers,181)
# self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([params[m,:] for m,mm in enumerate(members)]))
self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([m.param_values for m in members]))
#self.X_prime[n * self.nparams:(n + 1) * self.nparams, :] = np.transpose(np.array([np.delete(np.delete(m.param_values,21),10) for m in members]))
if samples != None:
self.rejection_threshold = samples.rejection_threshold
......@@ -153,7 +173,17 @@ class Optimizer(object):
for n in range(self.nlag):
members = statevector.ensemble_members[n]
for m, mem in enumerate(members):
members[m].param_values[:] = self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]
members[m].param_values[:] = self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]
# params = (self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams]) # 181
# params_bg = params[-1] # BG
# params = params[:-1] # remove BG for now
# params = params.reshape(9,20)
# params_holder = np.zeros((9,22))
# for r in range(0,8):
# params_holder[r,:] = np.insert(np.insert(params[r,:],20,0),10,0)
# members[m].param_values[:] = np.hstack((params_holder.flatten(),params_bg))
#members[m].param_values[:] = np.insert(np.insert(self.X_prime[n * self.nparams:(n + 1) * self.nparams, m] + self.x[n * self.nparams:(n + 1) * self.nparams],10,0),21,0)
logging.debug('Returning optimized data to the StateVector, setting "StateVector.isOptimized = True" ')
......
......@@ -57,19 +57,31 @@ class CO2StateVector(StateVector):
10. None
"""
# fullcov = np.array([ \
# (1.00, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.36, 1.00, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 1.00, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 0.36, 1.00, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 0.16, 0.16, 1.00, 0.36, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 0.16, 0.16, 0.36, 1.00, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 1.00, 0.16, 0.16, 0.16, 0.00), \
#(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 1.00, 0.16, 0.16, 0.00), \
# (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 1.00, 0.16, 0.00), \
# (0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 1.00, 0.00), \
# # (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.e-10) ])
# (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00) ])
fullcov = np.array([ \
(1.00, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
(0.36, 1.00, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
(0.16, 0.16, 1.00, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
(0.16, 0.16, 0.36, 1.00, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
(0.16, 0.16, 0.16, 0.16, 1.00, 0.36, 0.04, 0.04, 0.04, 0.01, 0.00), \
(0.16, 0.16, 0.16, 0.16, 0.36, 1.00, 0.04, 0.04, 0.04, 0.01, 0.00), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 1.00, 0.16, 0.16, 0.16, 0.00), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 1.00, 0.16, 0.16, 0.00), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 1.00, 0.16, 0.00), \
(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 1.00, 0.00), \
# (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.e-10) ])
(0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00) ])
(1.00, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.36, 1.00, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 1.00, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.36, 1.00, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.16, 0.16, 1.00, 0.36, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.16, 0.16, 0.36, 1.00, 0.04, 0.04, 0.04, 0.01), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 1.00, 0.16, 0.16, 0.16), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 1.00, 0.16, 0.16), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 1.00, 0.16), \
(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 1.00) ])
return fullcov
......
......@@ -51,13 +51,16 @@ class ObservationOperator(object):
def prepare_run(self):
""" Prepare the running of the actual forecast model, for example compile code """
# Define the name of the file that will contain the modeled output of each observation
# Define the name of the file that will contain the modeled output of each observation
self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
def run(self,lag,dacycle,statevector):
members = statevector.ensemble_members[lag]
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
#self.nparams = int(dacycle.dasystem['nparameters'])
self.nparams = int(self.dacycle['nparameters'])
absolute_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
absolute_start_time_ch = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y-%m-%d'))
starth = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.start']).days)*24
......@@ -65,7 +68,7 @@ class ObservationOperator(object):
f = io.CT_CDF(self.simulated_file, method='create')
logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
dimid = f.createDimension('obs_num', size=None)
dimid = ('obs_num',)
savedict = io.std_savedict.copy()
......@@ -86,13 +89,14 @@ class ObservationOperator(object):
savedict['units'] = "ppm"
savedict['dims'] = dimid + dimmember
savedict['comment'] = "Simulated model value created by COSMO"
f.add_data(savedict,nsets=0)
# Open file with x,y,z,t of model samples that need to be sampled
# Open file with x,y,z,t of model samples that need to be sampled
f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read')
# Get simulated values and ID
# Get simulated values and ID
ids = f_in.get_variable('obs_num')
obs = f_in.get_variable('observed')
......@@ -123,34 +127,47 @@ class ObservationOperator(object):
# for ncfile in ncfilelist:
# infile = os.path.join(ncfile + '.nc')
# UNCOMMENT FROM HERE
# co2_bg = np.empty(self.forecast_nmembers)
#
# logging.info('Multiplying emissions with parameters for lag %d' % (lag))
# for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start']+timedelta(hours=24*lag*int(dacycle['time.cycle'])), until=dacycle['time.start']+timedelta(hours=(lag+1)*24*int(dacycle['time.cycle']))):
# for ens in range(0,self.forecast_nmembers):
# dthh = dt.strftime('%H')
# co2_bg[ens] = members[ens].param_values[-1]
# ens = str(ens).zfill(3)
# cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_gpp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_resp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# # logging.debug('Background CO2 params are (%s)' % co2_bg)
# if dthh=='00':
# ct(dt.strftime('%Y%m%d'), co2_bg)
#
# cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
# cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))
#
# os.chdir(dacycle['da.obsoperator.home'])
#
self.lambda_file = os.path.join(self.outputdir, 'lambda.%s.nc' % self.dacycle['time.sample.stamp'])
ofile = Dataset(self.lambda_file, mode='w')
opar = ofile.createDimension('param', self.nparams)
omem = ofile.createDimension('member', self.forecast_nmembers)#len(members.nmembers))
l = ofile.createVariable('lambda', np.float32, ('member','param'),fill_value=-999.99)
co2 = np.empty(shape=(self.forecast_nmembers,self.nparams))
for m in range(0,20):
co2[m,:] = members[m].param_values
l[:] = co2
ofile.close()
# UN/COMMENT FROM HERE
co2_bg = np.empty(self.forecast_nmembers)
logging.info('Multiplying emissions with parameters for lag %d' % (lag))
for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start']+timedelta(hours=24*lag*int(dacycle['time.cycle'])), until=dacycle['time.start']+timedelta(hours=(lag+1)*24*int(dacycle['time.cycle']))):
for ens in range(0,self.forecast_nmembers):
dthh = dt.strftime('%H')
co2_bg[ens] = members[ens].param_values[-1]
ens = str(ens).zfill(3)
cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_gpp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_resp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# logging.debug('Background CO2 params are (%s)' % co2_bg)
if dthh=='00':
ct(dt.strftime('%Y%m%d'), co2_bg)
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))
os.chdir(dacycle['da.obsoperator.home'])
if os.path.exists(dacycle['dir.da_run']+'/'+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
if os.path.exists(dacycle['dir.da_run']+"/non_opt_"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/old_non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
else:
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
os.system('python run_chain.py 21ens '+absolute_start_time_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
os.system('python run_chain.py '+self.dacycle['run.name']+' '+absolute_start_time_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
logging.info('COSMO done!')
os.chdir(dacycle['dir.da_run'])
args = [
......@@ -160,19 +177,17 @@ class ObservationOperator(object):
with Pool(self.forecast_nmembers) as pool:
pool.starmap(self.extract_model_data, args)
# pool.close()
# pool.join()
for i in range(0,self.forecast_nmembers):
idx = str(i).zfill(3)
# cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/51ens/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
# cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/OK_DONT_TOUCH/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp']) # last run with non-frac
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6 # in ppm
ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)):
f.variables['obs_num'][j] = data[0]
f.variables['obs_num'][j] = data[0]
f.variables['flask'][j,:] = model_data[:,j]
f.close()
#### WARNING ACHTUNG PAZNJA POZOR VNEMANIE data[2] is model data mismatch (=1000) by default in tools/io4.py!!! pavle
......@@ -235,7 +250,7 @@ class ObservationOperator(object):
site_height.main(cosmo_out, str(ens), ss, time_stamp)
cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("797", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc") # this needs changing to 1009 (797 + 212)
cdo.intlevel("1009", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc")
cdo.intlevel("3580", input = cosmo_out+"CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = cosmo_out+"CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc")
......
......@@ -48,7 +48,9 @@ class CO2Optimizer(Optimizer):
self.localization = True
self.localizetype = 'CT2007'
#T-test values for two-tailed student's T-test using 95% confidence interval for some options of nmembers
if self.nmembers == 50:
if self.nmembers == 21:
self.tvalue = 2.08
elif self.nmembers == 50:
self.tvalue = 2.0086
elif self.nmembers == 100:
self.tvalue = 1.9840
......
......@@ -150,7 +150,8 @@ class StateVector(object):
self.nlag = int(dacycle['time.nlag'])
self.nmembers = int(dacycle['da.optimizer.nmembers'])
# self.nparams = int(dacycle.dasystem['nparameters'])
self.nparams = 23
self.nparams = 181
#self.nparams = 198
self.nobs = 0
self.obs_to_assimilate = () # empty containter to hold observations to assimilate later on
......@@ -168,13 +169,11 @@ class StateVector(object):
# that maps onto it. From this map, a dictionary is created that allows a reverse look-up so that we can map parameters to a grid.
mapfile = os.path.join(dacycle.dasystem['regionsfile'])
# ncf = io.ct_read(mapfile, 'read')
# self.gridmap = ncf.get_variable('regions')
# self.tcmap = ncf.get_variable('transcom_regions')
# ncf.close()
regfile = Dataset(mapfile,mode='r')
self.gridmap = np.squeeze(regfile.variables['regions'])
self.tcmap = np.squeeze(regfile.variables['transcom_regions'])
self.lat = np.squeeze(regfile.variables['latitude'])
self.lon = np.squeeze(regfile.variables['longitude'])
regfile.close()
logging.debug("A TransCom map on 1x1 degree was read from file %s" % dacycle.dasystem['regionsfile'])
......@@ -182,33 +181,34 @@ class StateVector(object):
# Create a dictionary for state <-> gridded map conversions
nparams = 23
#nparams = self.gridmap.max()
# nparams = 198
nparams = 181
self.griddict = {}
for r in range(1, 11):
# for r in range(1, int(nparams) + 1):
sel = np.nonzero(self.gridmap.flat == r)
if len(sel[0]) > 0:
self.griddict[r] = sel
self.griddict[r+11] = sel # pavle - expand dictionary for nparam values because of RESP
for pft in range(1,18):
for r in range(1, 11):
sel = np.nonzero(self.gridmap[pft,:,:].flat == r)
if len(sel[0]) > 0:
self.griddict[pft,r] = sel
self.griddict[pft,r+11] = sel # pavle - expand dictionary for nparam values because of RESP
# pavle: sel sorts out regions by PFT
logging.debug("A dictionary to map grids to states and vice versa was created")
# Create a matrix for state <-> TransCom conversions
self.tcmatrix = np.zeros((self.nparams, 23), 'float')
self.tcmatrix = np.zeros((self.nparams, 9), 'float')
for r in range(1, self.nparams + 1):
sel = np.nonzero(self.gridmap.flat == r)
if len(sel[0]) < 1:
continue
else:
n_tc = set(self.tcmap.flatten().take(sel[0]))
if len(n_tc) > 1:
logging.error("Parameter %d seems to map to multiple TransCom regions (%s), I do not know how to handle this" % (r, n_tc))
raise ValueError
self.tcmatrix[r - 1, int(n_tc.pop()) - 1] = 1.0
for pft in range(0,17):
for r in range(1, self.nparams + 1):
sel = np.nonzero(self.gridmap[pft,:,:].flat == r)
if len(sel[0]) < 1:
continue
else:
n_tc = set(self.tcmap.flatten().take(sel[0]))
if len(n_tc) > 1:
logging.error("Parameter %d seems to map to multiple TransCom regions (%s), I do not know how to handle this" % (r, n_tc))
raise ValueError
self.tcmatrix[r - 1, int(n_tc.pop()) - 1] = 1.0
logging.debug("A matrix to map states to TransCom regions and vice versa was created")
......@@ -292,96 +292,17 @@ class StateVector(object):
# Create members 1:nmembers and add to ensemble_members list
for member in range(1, self.nmembers):
#rands = np.random.uniform(low=-0.5, high=0.5, size=self.nparams-1)
rands = np.random.uniform(low=-1., high=1., size=self.nparams-1)
rands_bg = np.random.uniform(low=-0.05, high=0.05, size=1)
# if member == 1 and lag == 0:
# rands=np.array([-0.5402309, 0.3021888, 0.346095, 0.4942706, 0.6276228, 0.7450897, 0.5155229, -0.3006688, 0.8386819, -0.1197997, -1, -0.9025682, -1, 0.5018581, 0.923963, -0.2253573, 0.4716174, 0.4637589, 0.7411034, 0.6526616, -0.3188075, -1, 0.01435228])
# if member == 2 and lag == 0:
# rands=np.array([0.1876527, -0.4323478, 0.2562822, -0.7032343, 0.83417, 0.1378682, 0.3969375, 0.2101818, 0.5018877, 0.1390382, -1, 0.9267841, 0.1679351, 0.1019922, 0.4781876, -0.08182608, 0.9295725, 0.05694205, -0.119032, -0.8025782, -0.8610154, -1, -0.03559723])
# if member == 3 and lag == 0:
# rands=np.array([-0.8030154, 0.08172598, -0.7622297, -0.4516532, 0.4969012, 0.8023733, -0.08237687, -0.8089882, -0.8048904, -0.4082271, -1, -0.1105984, -0.4859152, 0.3398715, -0.09020101, 0.462075, 0.5386102, -0.865896, 0.5063729, -0.01304542, -0.2682804, -1, 0.04435889])
# if member == 4 and lag == 0:
# rands=np.array([-0.555479, -0.6089986, -0.664399, -0.03338498, -1, -1, 0.1820136, 0.1071932, -0.9634671, -1, -1, -0.8547595, 0.273797, -0.802806, 0.1928162, -0.5642408, 0.6471559, 0.2174408, -0.05968733, -0.004075899, 0.9075829, -1, -0.02335107])
# if member == 5 and lag == 0:
# rands=np.array([-0.674361, 0.3454686, -0.4433188, -0.8486685, -0.7747395, 0.203481, -1, -0.09257206, -0.3892125, -0.02775251, -1, -0.6137948, -0.8518382, -0.8779979, -0.6405048, -0.8992591, -0.6312428, 0.4434872, -0.05931631, -0.2671082, -0.7183358, -1, -0.02273439])
# if member == 6 and lag == 0:
# rands=np.array([-0.4753621, -1, 0.6610132, -0.3525974, -0.987463, -0.2421917, -0.8256455, 0.1535834, 0.6603956, 0.2776085, -1, -0.6808311, -0.6428144, -1, -0.518185, -1, -1, 0.4222396, 0.7343495, 0.6562976, -0.2627398, -1, 0.00420737])
#if member == 7 and lag == 0:
# rands=np.array([0.9235092, -0.08976277, -0.693751, 0.7479168, -0.5282985, -0.325265, 0.5994204, 0.1087474, -0.8495662, 0.5450585, -1, -0.4888758, -0.6259907, 0.5597795, 0.7957715, 0.04843958, 0.2360165, -0.7671149, 0.2514055, -0.872794, -0.5365307, -1, -0.02834922])
# if member == 8 and lag == 0:
# rands=np.array([-0.2851942, 0.2424788, 0.6036716, 0.1861704, 0.2724828, -0.1842902, 0.4256651, 0.9268349, -0.5386767, -0.282816, -1, 0.1387041, 0.6110855, -0.5958148, -0.9803032, -0.5043366, -0.6089351, 0.2659722, -0.7406561, 0.409389, -0.07663137, -1, -0.04040922])
# if member == 9 and lag == 0:
# rands=np.array([0.7698311, 0.2813426, 0.4236566, -0.5569826, 0.06671345, 0.1696645, -0.5119113, 0.5269542, 0.4109541, 0.5580103, -1, -0.4272517, -0.556591, 0.1682963, 0.02460092, -0.2566992, -0.9167355, -0.7835014, 0.7101904, -0.5695945, -0.1746321, -1, 0.04223957])
# if member == 10 and lag == 0:
# rands=np.array([-0.6941028, -0.9384637, -0.9250847, -1, -0.959329, -1, 0.04787683, -0.896835, -0.4911715, -0.01541466, -1, 0.9626108, 0.6731427, 0.5944214, 1.05593, 1.0805, -0.1627321, -0.154229, 0.6146008, 0.420991, 0.7912573, -1, -0.02473733])
# if member == 11 and lag == 0:
# rands=np.array([0.96286, 0.9415601, -0.7406601, 0.119507, 0.5506695, 0.5813184, 0.9465097, -0.658381, -0.5169885, 0.3601404, -1, -0.1221153, -0.911568, 0.4350488, 0.983259, 0.4603367, 0.4221282, -0.3391019, -0.192896, 0.3793332, 0.8009401, -1, 0.03354285])
# if member == 12 and lag == 0:
# rands=np.array([0.4491677, -0.5364176, 0.1777764, 0.3457868, 0.4886276, 0.6140389, -0.9525508, 0.7758317, 0.7040095, 0.8019212, -1, 0.2732197, 0.994552, -0.199925, 0.7970242, 0.5645052, -0.6018554, 0.01864898, 0.0417577, 0.3765279, -0.2793771, -1, -0.008472145])
# if member == 13 and lag == 0:
# rands=np.array([-0.8257746, 0.5533802, -0.7910783, 0.4864411, -0.5154006, -0.5978034, -0.202793, -0.9381504, 0.2701632, 0.3198348, -1, 0.3341525, 1.001176, -0.2380723, 0.1933595, -0.6892359, -0.4271372, 0.6341484, -0.5026917, -0.3053397, 0.750175, -1, -0.02964199])
#if member == 14 and lag == 0:
# rands=np.array([0.701148, -0.1619855, -0.063532, -0.4197367, -0.6209376, -1, -0.8872383, -0.4309417, -1, -1, -1, 0.1452325, -0.2525761, 0.1687515, -0.3141938, -0.1218985, -0.4126042, 0.1641715, 0.3630136, -0.5947977, 0.7974821, -1, -0.004756214])
# if member == 15 and lag == 0:
# rands=np.array([-0.2590746, -0.5643597, -0.2787321, -0.2195021, -0.7023828, 0.5482118, 0.6232124, 0.5892509, 0.6535864, -0.5035124, -1, -0.6891443, 0.3709087, -0.1523427, -0.9473059, 0.1686965, 0.04820674, 0.07875252, 0.5189905, 0.2113448, -0.708707, -1, -0.01257968])
# if member == 16 and lag == 0:
# rands=np.array([0.4056181, 0.04536337, -0.6866099, -0.6401204, -0.2879847, -0.7785423, 0.6233089, 0.2740796, -0.8484439, 0.9647028, -1, 0.789863, 1.163193, -0.2164443, 0.08588136, 0.5577631, -0.3278782, -0.1751566, -0.7633625, 0.7689492, -0.819483, -1, -0.03903025])
# if member == 17 and lag == 0:
# rands=np.array([0.2155399, 0.1660601, 0.4649725, 0.5862293, -0.4452171, -0.7389292, -0.07185336, 0.9012561, -0.2824533, -0.7226023, -1, 0.4194106, -0.2461459, 0.109305, 0.2449136, -0.09785557, 0.7447503, -0.1461458, 0.3223822, 0.4597907, 0.9065241, -1, -0.0007600483])
# if member == 18 and lag == 0:
# rands=np.array([0.04116038, 0.7824114, 0.4636176, 0.1557975, 0.4527292, -0.4133459, 0.994531, -0.3580682, 0.5178631, 0.2745949, -1, 0.5335826, 0.265704, 0.1912572, 0.536295, -0.7093778, 0.07351156, -0.2118817, -0.7678121, 0.43079, 0.512942, -1, 0.03622137])
# if member == 19 and lag == 0:
# rands=np.array([0.3184969, 0.2478539, 0.1378153, -0.6565253, 0.6140021, -0.2936001, 0.6872, 0.8872706, 0.4869822, -0.2930091, -1, 0.2853186, 0.9346947, 0.7716588, -0.1080804, 0.3278734, -0.2541114, -0.7900179, -0.259804, 0.6708161, -0.8029982, -1, 0.000798168])
# if member == 20 and lag == 0:
# rands=np.array([-0.7732743, -0.3873762, -0.8547782, -0.07859742, -1, -0.4397394, 0.2048925, -0.3029643, 0.7081758, -0.06017716, -1, -0.466458, -0.2800255, 0.1250267, 0.632107, -0.7519485, 0.4637258, 0.5653031, 0.4806009, -0.01624972, 0.05281375, -1, 0.01537185])
#if member == 1 and lag == 1:
# rands=np.array([0.4061777, -0.5013675, -0.2862317, 0.6116837, 0.05803703, -0.3041813, 0.4205766, -0.04346657, -0.4199566, -0.2088348, -1, 0.8685411, -0.1089062, 0.8366363, -0.278753, -0.4832861, -0.147602, 0.3227859, 0.2312384, 0.357722, -0.03119155, -1, -0.03759673])
# if member == 2 and lag == 1:
# rands=np.array([-0.3528758, 0.2302615, 0.8945007, -0.2294943, -0.7223361, 0.4254291, 0.2018437, 0.555488, -0.03991819, 0.2198012, -1, 0.918696, 0.3752462, 1.003467, -0.3229041, 0.6925161, 1.036227, 0.7551609, 1.008353, 0.6003665, -0.6694191, -1, 0.03411322])
# if member == 3 and lag == 1:
# rands=np.array([0.732138, 0.0511689, 0.5308715, -0.3620436, 0.1138427, 0.3123921, 0.7225392, 0.8311706, -0.7413127, 0.02754115, -1, 0.2934451, -0.7651439, -0.574313, -0.9932647, -0.5833654, -0.5632083, 0.2126374, -0.4117518, -0.6983544, -0.5544686, -1, 0.03195827])
# if member == 4 and lag == 1:
# rands=np.array([-0.5030815, -0.02779967, -0.9290833, 0.1386678, -0.510511, 0.297296, -0.6832802, -0.1534745, -0.7899511, 0.4862387, -1, -0.991854, -0.2682375, -0.7712117, -0.8905275, 0.3489431, 0.6761383, 0.1792146, 0.8675478, -0.7689523, -0.8875247, -1, 0.01273509])
# if member == 5 and lag == 1:
# rands=np.array([0.3466186, 0.6708195, -0.4234087, 0.791083, 0.5170874, -0.4746583, 0.9374641, 0.9230136, -0.5804763, -0.02489379, -1, 0.9292546, 0.3582832, -0.2736852, -0.4632118, 1.018919, 0.09223775, 0.3311332, 0.2550537, 0.7294372, 0.4026192, -1, 0.03795173])
# if member == 6 and lag == 1:
# rands=np.array([0.3490212, 0.2809639, 0.08980897, -0.6552433, 0.3108646, 0.9600278, -0.7445735, -0.929427, 0.6871993, -0.1122855, -1, 0.2775505, 0.7744503, 0.5353921, 0.8993617, 0.5068132, 0.6897388, -0.5868602, -0.5443915, 0.2293534, -0.880945, -1, -0.04769203])
# if member == 7 and lag == 1:
# rands=np.array([0.7321754, 0.06772549, 0.2179533, 0.2919183, 1.062581, 0.7554626, 0.5785612, 1.011077, 0.6109563, -0.350915, -1, -0.579572, -0.07610841, -0.3409401, 0.7316601, -0.904253, -1, -0.01014853, 0.5529011, 0.534495, -0.3398449, -1, -0.007357216])
#if member == 8 and lag == 1:
# rands=np.array([-0.4717625, -0.881941, -0.1655097, 0.2257761, 0.8264135, -0.7091568, 0.2949457, -0.296845, 0.794414, -0.743104, -1, 0.9034545, -0.2101449, 0.3021894, -0.4184734, -0.4638716, 0.6389658, 0.9005892, -0.5002754, 0.6826805, -0.3747995, -1, 0.03224342])
# if member == 9 and lag == 1:
# rands=np.array([-0.5808461, -0.8636444, 0.7052112, 0.7909082, -0.8140187, -0.02347049, 0.5638158, -0.7550521, 0.5693654, 0.07545193, -1, -0.2417251, 0.6826159, 0.9508564, -0.3956808, 0.8499226, 0.09175592, 0.5388746, 0.8177778, 0.4029195, 0.4111119, -1, -0.01329621])
# if member == 10 and lag == 1:
# rands=np.array([0.2166863, -0.437274, -0.7773396, 0.3759225, -0.7887743, -0.5783228, 0.2769741, -0.8951628, 0.5877202, -0.7473742, -1, 0.8568208, -0.3967594, -0.02381929, 0.911634, 0.1854728, 0.347165, -0.7240146, 0.6684847, -0.05459603, -0.2197302, -1, -0.04627957])
# if member == 11 and lag == 1:
# rands=np.array([-0.5553841, 0.1981317, -0.01706193, 0.3010842, 0.3364336, 0.6311802, -0.4822653, -1, -0.2353798, -0.4258329, -1, -0.3820901, -0.1350132, 0.1199797, -0.4477318, -0.03276956, 0.2006144, 0.7612231, -0.4382411, 0.3140147, -0.7996647, -1, -0.0351413])
# if member == 12 and lag == 1:
# rands=np.array([0.1795571, 0.7117739, 0.3948795, 0.01603693, 0.4342803, 0.3915263, -0.1187998, 0.3347789, -0.1160258, -0.6199519, -1, 0.9331298, 0.2128856, 0.1029978, 0.5622161, -0.03326744, 0.1113134, 0.718336, -0.7506785, 0.8058456, -0.165978, -1, 0.03335277])
# if member == 13 and lag == 1:
# rands=np.array([0.1133219, 0.1868301, 0.6434304, 0.06698409, 0.7134196, 0.6549824, -0.3687601, 0.2397162, 0.3933871, -0.8959517, -1, -0.9083285, 0.1290343, -0.8501528, -0.1640191, -0.04551715, 0.2021651, 0.1527653, 0.4235286, 0.7385602, -0.1721378, -1, 0.008510991])
# if member == 14 and lag == 1:
# rands=np.array([0.5128248, -0.025835, -0.3690266, 0.5315114, -0.3562845, -0.8088912, 0.6070702, 0.8863306, 0.127363, 0.4851948, -1, 0.2463859, 0.5403649, 0.1009353, -0.3769973, 0.8082039, 0.642301, -0.3246575, -0.8921521, 0.485802, 0.06472221, -1, 0.042198])
#if member == 15 and lag == 1:
# rands=np.array([-0.6925536, -0.1888134, 0.7955419, -0.1718084, -0.4925962, -0.63966, 0.712366, 0.7542216, 0.6485603, 0.9584286, -1, 0.2036391, -0.547542, -0.1605409, 0.5345823, -0.8189099, -0.4620013, 0.2116281, 0.5273049, -0.5836464, -0.7646005, -1, 0.03012722])
# if member == 16 and lag == 1:
# rands=np.array([-0.9031112, -0.9747161, 0.6262707, 0.1299344, -0.656873, -1, 0.7236382, 0.9888599, 0.599555, -0.3344736, -1, 0.1587437, 0.9665449, 0.7327783, 0.2025174, 1.045059, 0.7680195, 0.202516, 0.3561991, -0.08263882, -0.7926366, -1, -0.0340635])
# if member == 17 and lag == 1:
# rands=np.array([-0.4468101, -0.6611752, 0.02045458, -0.818266, -0.06696166, -0.2870016, -0.8372102, -1, 0.6159038, 0.1100874, -1, -0.2941307, -0.8388062, 0.7502818, -0.02831579, -0.4171109, 0.3452238, 0.3286322, -0.1746102, 0.741174, -0.1962416, -1, -0.03235182])
# if member == 18 and lag == 1:
# rands=np.array([-0.117545, 0.6358077, -0.8399339, -0.4401485, -0.4653139, -0.1412946, 0.368727, 0.5386308, 0.07112807, 0.2995664, -1, -0.3923609, -0.8697481, 0.283159, -0.1124206, -0.2029478, -1, -0.9128151, 0.2968773, -0.4699179, -0.1645328, -1, -0.04468708])
# if member == 19 and lag == 1:
# rands=np.array([-0.3586554, -0.8109571, -0.4370517, 0.1238612, -0.996357, -0.6696278, 0.5044641, 0.6546122, -0.0149131, 0.8158504, -1, 0.8035823, 0.565383, 0.2036764, -0.6208936, 0.4221191, 0.4371626, -0.6936532, -0.6731901, -1, 0.1557615, -1, -0.02435511])
# if member == 20 and lag == 1:
# rands=np.array([-0.6186742, 0.1736774, 0.3401965, 0.8706512, -0.6920347, -0.2222075, 0.1594849, 0.3773506, 1.039722, 0.3430154, -1, 0.8352696, 1.13694, -0.20392, 0.1833831, 1.04152, 0.1299644, 0.07491443, -0.6444398, 0.2261413, 0.01262438, -1, -0.03116111])
newmember = EnsembleMember(member)
# newmember.param_values = rands + newmean
newmember.param_values = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel()
# newmember.param_values[10] = 0.
# newmember.param_values[21] = 0.
rands = rands.reshape([9,20])
randsC = np.zeros(shape=(9,20))
for r in range(0,9):
randsC[r,:] = np.hstack((np.dot(C, rands[r,0:10]),np.dot(C, rands[r,10:20])))
newmember.param_values = (np.hstack((randsC.flatten(),rands_bg)) + newmean).ravel() #THIS ONE
#newmember.param_values = (np.hstack((np.dot(C, rands[0:11]),np.dot(C, rands[11:22]), rands_bg)) + newmean).ravel() #THIS ONE
newmember.param_values[newmember.param_values<0.] = 0.
newmember.param_values[newmember.param_values>2.] = 2.
self.ensemble_members[lag].append(newmember)
......@@ -508,83 +429,20 @@ class StateVector(object):
logging.info('Successfully read the State Vector from file (%s) ' % filename)
def write_members_to_file(self, lag, outdir,endswith='.nc'):
"""
:param: lag: Which lag step of the filter to write, must lie in range [1,...,nlag]
:param: outdir: Directory where to write files
:param: endswith: Optional label to add to the filename, default is simply .nc
:rtyp