improved documentation

da/ffdas/observationoperator.py

@@ -4,10 +4,11 @@
 """
 Author : W. He
 Adapted by Ingrid Super, last revisions 14-8-2018
+Adapted by Auke van der Woude, 09-2019 -- ...
 
 Revision History:
-Basing on Wouter's codes, replace the TM5 model with the STILT model, April 2015.
-
+| -- Basing on Wouter's codes, replace the TM5 model with the STILT model, April 2015.
+| -- Expanding code to include C14 and a biosphere module 
 This module holds specific functions needed to use the STILT model within the data assimilation shell. It uses the information
 from the DA system in combination with the generic stilt.rc files.
 
@@ -242,19 +243,6 @@ class STILTObservationOperator(ObservationOperator):
         while dumdate<(self.timefinishkey): # self.timefinishkey : dt
             self.datelist.append(dumdate)
             dumdate=dumdate+dt
-        self.prepare_background(do_pseudo, self.timefinishkey)
-
-    def prepare_background(self, do_pseudo, datepoint):
-        background_concs = []
-        indices = self.get_time_indices(datepoint)
-        for species in self.spname:
-            filename = self.bgfile.strip('.nc') + '_' + species.lower() + '.nc'
-            bf = io.ct_read(filename, 'read')
-            background_conc = bf.get_variable('obs')[indices]
-            background_concs.append(background_conc)
-        background_concs = np.array(background_concs)
-        self.bg = background_concs
-        bf.close()
 
     ### !!!! Only cache if made dependent on datepoint/datetime !!!!
     def get_time_index_nc(self, time=None):
@@ -288,7 +276,12 @@ class STILTObservationOperator(ObservationOperator):
         """Function that gets the footprint for the current time and site.
         Returns a 3D np.array with dims (time, lat, lon)
         Input:
-            site: str of 3 letter sitename. Is converted to uppercase. Example: 'hei'"""
+            i_loc: number of location. Requires object to have a 'sitenames' property
+            datepoint: datetime of the footprint to be found.
+        Returns:
+            np.array (3d): Footprint with as indices time, lat, lon
+            """
+
         site = self.sitenames[i_loc].upper()
         path = self.model_settings['footdir'] + '/' + site.upper() + '24'
         fname = path + '/footprint_{0}_{1}x{2:02d}x{3:02d}x{4:02d}*.nc'.format(site.upper(),
@@ -305,11 +298,17 @@ class STILTObservationOperator(ObservationOperator):
     @cached 
     def get_background(self, i_species, i_loc, datepoint):
         """Function that finds the center of mass of the first footprint and the time corresponding to it.
+            and finds the concentration in the center of mass. This is used as the background.
         Input: 
-            fp: 3d np.array of the footprint
-        time is taken from the observationoperator self."""
+            i_species: int: the index of the species for which the background should be found
+            i_loc    : int: the index of the locatin for which the background should be found
+            datepoint: datetime.datetime: the datetime of the background concentration
+        Returns:
+            float: background concentration
+            """
         from scipy import ndimage
         
+        # First, get the footprint and find the first time of influence
         foot = self.get_foot(i_loc, datepoint)
         start_influence = -1 # In backtimes
         total_influence = 0
@@ -317,13 +316,15 @@ class STILTObservationOperator(ObservationOperator):
             start_influence += 1
             total_influence = fp[start_influence].sum()
         
+        # Get the center of mass of the first influence
         center_of_mass = ndimage.measurements.center_of_mass(fp[start_influence])
         center_of_mass = np.array(np.rint(center_of_mass), dtype=int)
 
         species_name = self.spname[i_species]
 
         index  = self.get_time_index_nc() - (int(self.model_settings['num_backtimes']) - start_influence)
-
+        
+        # Find the concentration at the center of mass at the correct time.
         background_file = self.model_settings['bgfdir'] + '/background.nc'
         background_map  = nc.Dataset(background_file)[species_name]
         background_conc = background_map[index, center_of_mass[0], center_of_mass[1]]
@@ -332,15 +333,32 @@ class STILTObservationOperator(ObservationOperator):
     
     @cached
     def get_background_orig(self, i_species, i_datepoint):
+        """Function that finds the background concentration, non-time dependent and hard-coded.
+        Input:
+            i_species: int: the index of the species for which the background should be found
+            i_loc    : int: the index of the locatin for which the background should be found
+            datepoint: datetime.datetime: the datetime of the background concentration
+        Returns:
+            float: background concentration
+            """
+
         backgrounds = 406.15, 127.2
         return backgrounds[i_species]
 
     @cached
     def get_c14_concentration(self, i_loc, i_datepoint,  gpp, ter, ff_flux, background):
-        """Function that gets the c14 concentration based on the emissions by the biosphere and
-        emissions by nuclear plants. Based Ch4 of Deniza's thesis: Modeling $\Delta^{14}CO_2 for Western Europe
-        DobsCO2obs = Dbg(CO2bg + CO2p + CO2r) + DffCO2ff + Dn14CO2n"""
-        
+        """Function that gets the c14 concentration based on the emissions by the biosphere, 
+            fossil fuels and nuclear power. The constants are defined at the top of this script.
+        Input: 
+            i_loc: int: the index of the location for which the c14 concentration should be calculated
+            i_datepoint: int: the index of the datepoint for which the c14 should be calculated
+            gpp: float: The gross primary production in umol/s
+            ter: float: The total ecosystem respiration in umol/s
+            ff_flux: float: The fossil fuel flux in umol/s
+            background: float: The background CO2 concentration
+        Returns: 
+            float: The C14 in ppm
+            float: Delta(14C) """
         datepoint = self.datelist[i_datepoint]
         # First, get the footprint
         site = self.sitenames[i_loc]
@@ -393,6 +411,12 @@ class STILTObservationOperator(ObservationOperator):
 
     @cached
     def get_nc_variable(self, file, fluxname):
+        """Helper function that gets the values from an nc file
+        Input: 
+            file: str: filename of the nc file of which the value should be taken
+            fluxname: str: name of the variable in the file
+        Returns:
+            np.ndarray: the values in the nc file"""
         data = nc.Dataset(file, 'r')
         fluxmap = data[fluxname][:]
         data.close()
@@ -402,15 +426,19 @@ class STILTObservationOperator(ObservationOperator):
     def get_biosphere_concentration(self,i_loc, datepoint, total=False):
         """Function that calculates the atmospheric increase due to the exchange of fluxes over the footprint
         Input:
-            site: str of 3 chars with the location. Example: 'hei'
-            tracer_info: dict with information on the tracer:
-                path to fluxes,
-                half-life,
-                recalculation factor"""
+            i_loc: int: index of the 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:
+                float: The total biosphere flux in umol/s
+            else: 
+                tuple of 2 floats: GPP and TER in umol/s """
+
         # First, get the footprint
         site = self.sitenames[i_loc]
         foot = self.get_foot(i_loc, datepoint)
-
+        # Get the indices 
         indices = self.get_time_indices(datepoint)
         # Get the tracer fluxes
         file = self.model_settings['biosphere_fluxdir']
@@ -424,8 +452,14 @@ class STILTObservationOperator(ObservationOperator):
         else: return gpp_increase, ter_increase
 
     @cached
-    def get_temporal_profiles(self, i_datepoint, i_station,  i_member, do_pseudo=0):
-        """ Function that reads in the temporal profiles for the current timestep"""
+    def get_temporal_profiles(self, i_datepoint, i_station,  i_member):
+        """ Function that reads in the temporal profiles for the current timestep
+        Input:
+            i_datepoint: int: the index of the datepoint for which the c14 should be calculated
+            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 = self.get_time_indices(self.datelist[i_datepoint])
         station_name = self.sitenames[i_station]
@@ -442,6 +476,12 @@ class STILTObservationOperator(ObservationOperator):
 
     @cached
     def get_spatial_emissions(self, i_member, i_species):
+        """ Function that gets the spatial emissions
+        Input:
+            i_member: int: the index of the member for which the simulation is run
+            i_species: int: the index of the species for which the simulation is run
+        Returns:
+            np.ndarray (3d): the spatial emissions per category, lat and lon"""
         #read spatially distributed emissions calculated with the emission model
         emisfile = os.path.join(self.model_settings['datadir'],'prior_spatial_%03d.nc'%i_member) #format: species[SNAP,lon,lat]
         f = io.ct_read(emisfile,method='read')
@@ -459,7 +499,7 @@ class STILTObservationOperator(ObservationOperator):
             i_species  : int: index of the species
             do_pseudo  : bool: whether to do pseudo-observations or not
         Returns:
-            None"""
+            float: the concentration increase at the respective location due to the emissions from the respective species"""
         # get the date:
         datepoint = self.datelist[i_datepoint]
 
@@ -493,9 +533,9 @@ class STILTObservationOperator(ObservationOperator):
             self.save_obs()
 
     def run(self,dacycle,do_pseudo):
-        """ This function calls the OPS and STILT functions for each locations, species, ensemble member and time step"""
+        """Function that loops over all members, locations, species and datepoints and calculates the concentration increase at that point (in time)
+        Adds the data as attribute to this object."""
 
-        #### do STILT
         totser=np.array(int(self.dacycle['da.optimizer.nmembers'])*[(self.nrloc*self.nrspc*len(self.datelist))*[0.]])
         for i_member in range(int(self.dacycle['da.optimizer.nmembers'])):
             conc_STILT=[]