Include Emissionmodel in inversion

.gitignore

@@ -2,11 +2,5 @@
 *.py
 *.rc
 *.jb
-joblog
-# Exclude the logfiles
-*.log
-# Exclude files created by the slurm system
-slurm*.out
-# Exclude compiled files
-da/**/__pycache__/
-da/**/*.pyc
+/da/*/__pycache__
+/da/*/*.pyc

da/ccffdas/observationoperator.py

@@ -40,15 +40,15 @@ from da.tools.general import create_dirs, to_datetime
 from da.ccffdas.emissionmodel import EmisModel 
 from da.baseclasses.observationoperator import ObservationOperator
 
-#try: # Import memoization. This speeds up the functions a lot.
-#    from memoization import cached
-#except: # The memoization package is not always included. Import the memoization from #functools
-#    from functools import lru_cache as cached
-#    cached = cached()
+try: # Import memoization. This speeds up the functions a lot.
+    from memoization import cached
+except: # The memoization package is not always included. Import the memoization from functools
+    from functools import lru_cache as cached
+    cached = cached()
 import xlrd
 
 
-from da.ccffdas import category_info
+from da.ffdas import category_info
 categories = category_info.categories
 # global constants, which will be used in the following classes
 identifier = 'WRF-STILT'
@@ -86,12 +86,10 @@ epsilon_14CO2_gpp = -36.
 alpha_14CO2_gpp   = 1 + epsilon_14CO2_gpp/1000.                                  # = 0.964
 ################### Begin Class STILT ###################
 
-model_settings = rc.read('/projects/0/ctdas/awoude/develop/exec/da/rc/stilt/stilt.rc')
+model_settings = rc.read('/projects/0/ctdas/RINGO/inversions/ffdas/exec/da/rc/ffdas/stilt.rc')
 recalc_factors = [1, 1000]
 spname = ['CO2', 'CO']
-
-
-def run_STILT(dacycle, footprint, datepoint, i_species, path, i_member):
+def run_STILT(footprint, datepoint, site, i_species, path, i_member):
     """This function reads in STILT footprints and returns hourly concentration data
     Input:
         footprint: np.array: the footprint of the station
@@ -102,37 +100,39 @@ def run_STILT(dacycle, footprint, datepoint, i_species, path, i_member):
     Returns:
         float: the concentration increase at the respective location due to the emissions from the respective species"""
     # get the date:
-    spatial_emissions = get_spatial_emissions(dacycle, i_member, i_species, path, datepoint)
+    temp_profiles = get_temporal_profiles(datepoint, site, path,  i_member)
+    spatial_emissions = get_spatial_emissions(i_member, i_species, path)
 
     # multiply footprint with emissions for each hour in the back trajectory. Indices: Time, Category, lat, lon
-    foot_flux = (footprint[None, :, :, :] * spatial_emissions[:, :, :, :]).sum()
+    foot_flux = (footprint[:, None, :, :] * spatial_emissions[None, :, :, :] * temp_profiles[:, :, None, None]).sum()
     concentration_increase = float(recalc_factors[i_species]) * foot_flux
 
     return concentration_increase
 
-#@cached
-#def get_temporal_profiles(datepoint, station_name, path, i_member):
-#    """ Function that reads in the temporal profiles for the current timestep
-#    Input:
-#        datepoint: datepoint: datetime for which the temporal profile should be found
-#        i_station: int: the index of the location for which the c14 concentration should be calculated
-#        i_member: int: the index of the member for which the simulation is run
-#    Returns:
-#        np.array (2): the time profiles for all categories. Indices: time, category"""
-#    #read temporal profiles for the times within the footprint
-#    time_indices = get_time_indices(datepoint)
-#    temporal_prior = path + 'prior_temporal_{0}_{1:03}.nc'.format(station_name, i_member)
-#    temporal_prior = io.ct_read(temporal_prior, 'read')
-#    temp_profiles = []
-#    for category, values in categories.items():
-#        temporal_var_name = values['temporal']
-#        temporal_variation = temporal_prior.get_variable(temporal_var_name)[time_indices]
-#        temp_profiles.append(temporal_variation)
-#    #temporal_prior.close()
-#    temp_profiles = np.array(temp_profiles)
-#    return temp_profiles.T # Transpose to be similar to spatial data in dimensions
-
-def get_spatial_emissions(dacycle, i_member, i_species, path, datepoint):
+@cached
+def get_temporal_profiles(datepoint, station_name, path, i_member):
+    """ Function that reads in the temporal profiles for the current timestep
+    Input:
+        datepoint: datepoint: datetime for which the temporal profile should be found
+        i_station: int: the index of the location for which the c14 concentration should be calculated
+        i_member: int: the index of the member for which the simulation is run
+    Returns:
+        np.array (2): the time profiles for all categories. Indices: time, category"""
+    #read temporal profiles for the times within the footprint
+    time_indices = get_time_indices(datepoint)
+    temporal_prior = path + 'prior_temporal_{0}_{1:03}.nc'.format(station_name, i_member)
+    temporal_prior = io.ct_read(temporal_prior, 'read')
+    temp_profiles = []
+    for category, values in categories.items():
+        temporal_var_name = values['temporal']
+        temporal_variation = temporal_prior.get_variable(temporal_var_name)[time_indices]
+        temp_profiles.append(temporal_variation)
+    #temporal_prior.close()
+    temp_profiles = np.array(temp_profiles)
+    return temp_profiles.T # Transpose to be similar to spatial data in dimensions
+
+@cached
+def get_spatial_emissions(i_member,  i_species, path):
     """ Function that gets the spatial emissions
     Input:
         i_member: int: the index of the member for which the simulation is run
@@ -140,15 +140,9 @@ def get_spatial_emissions(dacycle, i_member, i_species, path, datepoint):
     Returns:
         np.ndarray (3d): the spatial emissions per category, lat and lon"""
     #read spatially distributed emissions calculated with the emission model
-    backtimes = 24 + int(dacycle['time.cycle']) * 24 # extra hours in backtime
-    backtimes = dtm.timedelta(hours=backtimes)
-    
-    start_time = (dacycle['time.sample.end'] - backtimes)
-    indices = get_time_indices(datepoint, start_time)
     emisfile = path +'prior_spatial_{0:03d}.nc'.format(i_member) #format: species[SNAP,lon,lat]
     f = io.ct_read(emisfile,method='read')
     emis=f.get_variable(spname[i_species])
-    emis = emis[:, indices, :, :] # index the interesting times
     f.close()
     return emis
 
@@ -166,17 +160,16 @@ def get_time_index_nc(time=None, startdate=None):
     # Note that this is in hours, and thus assumes that the flux files are hourly as well
     timediff = time - startdate
     timediff_hours = int(timediff.total_seconds()/3600) # 1hour could also be softcoded from time.cycle
-    time_index = int(timediff_hours)
+    time_index = int(timediff_hours) + DO_RINGO
     return time_index
 
-#@cached(ttl=5)
-
-def get_time_indices(datepoint, startdate=None, backtimes=24):
+@cached(ttl=5)
+def get_time_indices(datepoint, startdate=None):
     """Function that gets the time indices in the flux files
     Because if the footprint is for 24 hours back, we need the fluxes 24 hours back"""
 
     time_index = get_time_index_nc(datepoint, startdate=startdate)
-    return slice(time_index - backtimes, time_index)
+    return slice(time_index - int(model_settings['num_backtimes']), time_index)
 
 def get_biosphere_concentration(foot, gpp_mem, ter_mem):
     """Function that calculates the atmospheric increase due to the exchange of fluxes over the footprint
@@ -235,7 +228,7 @@ class STILTObservationOperator(ObservationOperator):
         self.btime = int(dacycle.dasystem['run.backtime'])
         self.categories = categories
         self.inputdir = dacycle.dasystem['inputdir']
-        self.paramdict = rc.read(dacycle.dasystem['paramdict'])
+        self.paramdict = rc.read('/projects/0/ctdas/RINGO/inversions/Data/paramdict.rc')
 
         biosphere_fluxes = nc.Dataset(dacycle.dasystem['biosphere_fluxdir'])
         self.gpp = biosphere_fluxes['GPP'][:]#[time_indices]
@@ -356,14 +349,10 @@ class STILTObservationOperator(ObservationOperator):
     def run_dynamic_emis_model(self):
         """Function that runs the dynamic emission model and writes the spatial emissions and time profiles to nc files
         based on emismodel.py"""
-        logging.info('Calculating the emissions; entering Emismodel.py')
         emismodel = EmisModel()
         emismodel.setup(self.dacycle)
+        emismodel.get_emis(self.dacycle, self.allsamples, do_pseudo=0)
         self.emismodel = emismodel
-        backtimes = 24 + int(self.dacycle['time.cycle']) * 24 # extra hours in backtime
-        indices = get_time_indices(self.dacycle['time.sample.end'], backtimes=backtimes)
-        emismodel.get_emis(self.dacycle, self.allsamples, indices=indices, do_pseudo=0)
-        logging.info('Emissions calculated')
 
     def prepare_biosphere(self, datepoint):
         """Function that prepares the biosphere fluxes on a map
@@ -438,7 +427,7 @@ class STILTObservationOperator(ObservationOperator):
         # Note that this is in hours, and thus assumes that the flux files are hourly as well
         timediff = time - startdate
         timediff_hours = int(timediff.total_seconds()/3600) # 1hour could also be softcoded from time.cycle
-        time_index = int(timediff_hours)
+        time_index = int(timediff_hours) + DO_RINGO
         return time_index
     
     def get_time_indices(self, datepoint, startdate=None):
@@ -525,6 +514,7 @@ class STILTObservationOperator(ObservationOperator):
         Input:
             i_member: int: member number for which the biosphere concentration should be calculated
             site: str: location for which the concentration should be calculated
+            datepoint: datetime.datetime: the datepoint for which the concentration should be calculated
             total: bool, optional: Whether to returned summed values or gpp and ter seperately as tuple
         Returns:
             if total == True:
@@ -604,7 +594,7 @@ class STILTObservationOperator(ObservationOperator):
         fluxmap = data[fluxname][:]
         data.close()
         return fluxmap
-
+    
     def run_forecast_model(self,dacycle,do_pseudo,adv):
         self.prepare_run(do_pseudo,adv)
         self.run(dacycle,do_pseudo)
@@ -632,7 +622,7 @@ class STILTObservationOperator(ObservationOperator):
             # Check if it is a new date
             self.sample = sample
             datepoint = sample.xdate
-            indices = get_time_indices(datepoint)
+
 
             if datepoint != previously_done_datepoint: 
                 # If it is a new date, we need a new biosphere prior
@@ -670,6 +660,12 @@ class STILTObservationOperator(ObservationOperator):
         # add the calculated concentrations to the object
         self.calculated_concentrations = calculated_concentrations
 
+
+        # Clear the cache from the functions, so that a new cycle does 
+        # (with different parameters) does not find the same results.
+        get_spatial_emissions.cache_clear()
+        get_temporal_profiles.cache_clear()
+
     def calc_concentrations(self, sample):
         """Function that calculates the concentration for a sample and all members.
         Gets the relevant information from the sample and then calculates the concentration
@@ -702,7 +698,7 @@ class STILTObservationOperator(ObservationOperator):
         if not 'C14' in species.upper():
             pool = Pool(self.forecast_nmembers)
             # Create the function that calculates the conentration
-            func = partial(run_STILT, self.dacycle, self.foot, datepoint, i_species, self.inputdir)
+            func = partial(run_STILT, self.foot, datepoint, site, i_species, self.inputdir)
             # We need to run over all members
             memberlist = list(range(0, self.forecast_nmembers))
             

da/ccffdas/pipeline.py

@@ -32,7 +32,7 @@ def ensemble_smoother_pipeline(dacycle, platform, dasystem, samples, statevector
     
     invert(dacycle, statevector, optimizer, obsoperator)
     
-    #advance(dacycle, samples, statevector, obsoperator)
+    advance(dacycle, samples, statevector, obsoperator)
         
     save_and_submit(dacycle, statevector)    
     logging.info("Cycle finished...exiting pipeline")

da/ccffdas/stilt-ops_urbanall.rc

 !!! Info for the CarbonTracker data assimilation system
-basepath : /projects/0/ctdas/awoude/develop/
+basepath : /projects/0/ctdas/RINGO/inversions/
 name     :  
-strategy : CO2
-datadir         :  /projects/0/ctdas/RINGO/inversions/Data
-inputdir        : ${basepath}/input/
-outputdir       : ${basepath}/output/
-restartdir      : ${basepath}/restart/
+strategy : CO2C14
+datadir         : ${basepath}/Data
+inputdir        : ${basepath}/${name}${strategy}/input/
+outputdir       : ${basepath}/${name}${strategy}/output/
+restartdir      : ${basepath}/${name}${strategy}/restart/
 ! list of all observation sites
 obs.input.id   : obsfiles.csv
 ! number of observation sites included; number of species included and to be used in inversion
@@ -15,17 +15,16 @@ obs.dir        : obsfiles${strategy}
 do.co         :      0
 do.c14integrated:    0
 do.c14targeted:      0
-obs.spec.name  : CO2
+obs.spec.name  : CO2,CO
 ! number of emission categories defined in the emission model
 obs.cat.nr     : 14
 ! For Rdam obs
-obs.sites.rc        : ${datadir}/sites_weights2.rc
+obs.sites.rc        : ${datadir}/sites_weights.rc
 ! number of parameters
 ! In the covmatrix and statevector, the ff parameters are first, then the bio parameters!
-nffparameters     : 24
-nbioparameters    : 6
-nparameters       : 30
-paramdict         : /projects/0/ctdas/RINGO/inversions/Data/paramdict.rc
+nffparameters     : 16
+nbioparameters    : 4 
+nparameters       : 20
 ! set fixed seed for random number generator, or use 0 if you want to use any random seed
 random.seed     : 4385
 !file with prior estimate of scaling factors (statevector) and covariances

da/platform/cartesius.py

@@ -100,7 +100,7 @@ class CartesiusPlatform(Platform):
                    """#SBATCH -n jobnodes \n""" + \
                    """#SBATCH -t jobtime \n""" + \
                    """#SBATCH -o joblog \n""" + \
-                   """module load Python/3.6.6-intel-2018b\n""" + \
+                   """module load python\n""" + \
                    """module load nco\n""" + \
            """\n"""
 

da/rc/stilt/stilt.rc

@@ -11,15 +11,17 @@
 ! You should have received a copy of the GNU General Public License along with this 
 ! program. If not, see <http://www.gnu.org/licenses/>. 
      homedir : /projects/0/ctdas/RINGO
-biosphere_fluxdir      : /projects/0/ctdas/RINGO/EmissionInventories/True/RINGO_ORCHIDEE_GPP_TER_dC14_old.nc
-nuclear_timeprofiledir : /projects/0/ctdas/RINGO/EmissionInventories/True/Power_TempDistributionFrance.xlsx
-nuclear_fluxdir        : /projects/0/ctdas/RINGO/EmissionInventories/True/RINGO_NuclearEmissions_14CO2.nc
+      sibdir : /projects/0/ctdas/RINGO/EmissionInventories/ 
       bgfdir : /projects/0/ctdas/RINGO/EmissionInventories/
      footdir : /projects/0/ctdas/RINGO/STILT_Output/
-inputdir: /projects/0/ctdas/RINGO/inversions/ffdas/input
+    bounddir : /projects/0/ctdas/RINGO/EmissionInventories/
+    inputdir : /home/awoude/ctdas_test/ctdas_stilt/input/
+      outdir : /home/awoude/ctdas_test/ctdas_stilt/output/
         nsam : 100
   ndayscycle : 10
   sysstadate : 2016-01-15 00:00:00
   cycstadate : 2016-01-20 00:00:00
 files_startdate : 2016-01-01 00:00:00
-num_backtimes: 24
+tracerdir    : /home/awoude/ctdas_test/ctdas_stilt/exec/da/rc/stilt/tracers
+! Now list the tracers separated by a comma, without a space!
+tracers      : co2,222Rn

da/tools/initexit.py

@@ -157,7 +157,7 @@ class CycleControl(dict):
                 self[k] = False
             if 'date' in k : 
                 self[k] = to_datetime(v)
-            if k in ['time.start', 'time.end', 'time.finish', 'da.restart.tstamp', 'time.fxstart']:
+            if k in ['time.start', 'time.end', 'time.finish', 'da.restart.tstamp']:
                 self[k] = to_datetime(v)
         for key in needed_da_items:
             if key not in self:
@@ -333,7 +333,7 @@ class CycleControl(dict):
             strippedname = os.path.split(self['jobrcfilename'])[-1]
             self['jobrcfilename'] = os.path.join(self['dir.exec'], strippedname)
 #            shutil.copy(os.path.join(self.dasystem['regionsfile']),os.path.join(self['dir.exec'],'da','analysis','copied_regions.nc'))
-            #logging.info('Copied regions file to the analysis directory: %s'%os.path.join(self.dasystem['regionsfile'])) 
+            logging.info('Copied regions file to the analysis directory: %s'%os.path.join(self.dasystem['regionsfile'])) 
             if 'extendedregionsfile' in self.dasystem:
 #                shutil.copy(os.path.join(self.dasystem['extendedregionsfile']),os.path.join(self['dir.exec'],'da','analysis','copied_regions_extended.nc')) 
                 logging.info('Copied extended regions file to the analysis directory: %s'%os.path.join(self.dasystem['extendedregionsfile'])) 

da/tools/rc.py

@@ -262,7 +262,7 @@ class RcFile(object) :
     Class to store settings read from a rcfile.
     """   
 
-    def __init__(self, filename, silent=False, marks=('${', '}'), do_eval=True):
+    def __init__(self, filename, silent=False, marks=('${', '}')) :
 
         """ 
         
@@ -836,6 +836,7 @@ class RcFile(object) :
                     #endif
                 else :
                     if do_eval:
+                        print(val)
                         if isinstance(val, str):
                             if any(operator in val for operator in '+-*/'):
                                 try: val = eval(val)