Changeset 45b99e6 in flex_extract.git

Timestamp:

Feb 15, 2019, 3:48:31 PM (5 years ago)

Author:

Anne Philipp <anne.philipp@…>

Branches:

master, ctbto, dev

Children:

b9f1214

Parents:

ae2756e

Message:

added possibility to extract ensemble members with ELDA stream (and ENFO)

Location:

source/python

Files:

• classes/EcFlexpart.py (modified) (10 diffs)
• mods/get_mars_data.py (modified) (2 diffs)
• mods/prepare_flexpart.py (modified) (1 diff)

source/python/classes/EcFlexpart.py

@@ -67 +67 @@
                      codes_get_values, codes_set_values, codes_set,
                      codes_write, codes_release, codes_new_from_index,
-                     codes_index_release, codes_index_get)
+                     codes_index_release, codes_index_get, codes_get_array,
+                     codes_set_array, codes_grib_new_from_file)
 
 # software specific classes and modules from flex_extract
@@ -77 +78 @@
                         execute_subprocess)
 from MarsRetrieval import MarsRetrieval
+from UioFiles import UioFiles
 import mods.disaggregation as disaggregation
 
@@ -663 +665 @@
 
         for ftype in self.types:
-            # fk contains field types such as
+            # ftype contains field types such as
             #     [AN, FC, PF, CV]
-            # fv contains all of the items of the belonging key
-            #     [times, steps]
             for pk, pv in self.params.iteritems():
                 # pk contains one of these keys of params
@@ -707 +707 @@
                     retr_param_dict['area'] = ""
                     retr_param_dict['gaussian'] = 'reduced'
+                if ftype.upper() == 'FC' and \
+                    'acc' not in retr_param_dict['target']:
+                    if (int(retr_param_dict['time'][0]) +
+                        int(retr_param_dict['step'][0])) > 23:
+                        dates = retr_param_dict['date'].split('/')
+                        sdate = datetime.strptime(dates[0], '%Y%m%d%H')
+                        sdate = sdate - timedelta(days=1)
+                        retr_param_dict['date'] = '/'.join(
+                            [sdate.strftime("%Y%m%d")] +
+                            retr_param_dict['date'][1:])
+
+                        print('CHANGED FC start date to ' +
+                              sdate.strftime("%Y%m%d") +
+                              ' to accomodate TIME=' +
+                              retr_param_dict['time'][0] +
+                              ', STEP=' +
+                              retr_param_dict['time'][0])
 
     # ------  on demand path  --------------------------------------------------
@@ -760 +777 @@
                             times = retr_param_dict['time'].split('/')
                             steps = retr_param_dict['step'].split('/')
-                            print 'times', times, int(times[0]), times[1:]
-                            print 'steps', steps, int(steps[0])
+
                             while int(times[0]) + int(steps[0]) <= 12:
-                                print 'HELLO'
                                 times = times[1:]
-                                print 'in while 1 ', times
-
                                 if len(times) > 1:
                                     retr_param_dict['time'] = '/'.join(times)
                                 else:
                                     retr_param_dict['time'] = times[0]
-
-                                print 'in while 2 ', times
-                                print retr_param_dict['time']
 
                         if (pk != 'OG_OROLSM__SL' and
@@ -916 +926 @@
         # get the values of the keys which are used for distinct access
         # of grib messages via product
-        index_keys = ["date", "time", "step"]
+        if '/' in self.number:
+            # more than one ensemble member is selected
+            index_keys = ["number", "date", "time", "step"]
+        else:
+            index_keys = ["date", "time", "step"]
         iid, index_vals = self._mk_index_values(c.inputdir,
                                                 inputfiles,
@@ -991 +1005 @@
                 t_enddate = t_date + timedelta(2*int(c.dtime))
 
+            # if necessary, add ensemble member number to filename suffix
+            # otherwise, add empty string
+            if 'number' in index_keys:
+                index_number = index_keys.index('number')
+                if len(index_vals[index_number]) > 1:
+                    numbersuffix = '.N{:0>3}'.format(int(prod[index_number]))
+            else:
+                numbersuffix = ''
+
             if c.purefc:
                 fnout = os.path.join(c.inputdir, 'flux' +
                                      t_date.strftime('%Y%m%d.%H') +
-                                     '.{:0>3}'.format(step-2*int(c.dtime)))
+                                     '.{:0>3}'.format(step-2*int(c.dtime)) +
+                                     numbersuffix)
                 gnout = os.path.join(c.inputdir, 'flux' +
                                      t_date.strftime('%Y%m%d.%H') +
-                                     '.{:0>3}'.format(step-int(c.dtime)))
+                                     '.{:0>3}'.format(step-int(c.dtime)) +
+                                     numbersuffix)
                 hnout = os.path.join(c.inputdir, 'flux' +
                                      t_date.strftime('%Y%m%d.%H') +
-                                     '.{:0>3}'.format(step))
+                                     '.{:0>3}'.format(step) +
+                                     numbersuffix)
             else:
                 fnout = os.path.join(c.inputdir, 'flux' +
-                                     t_m2dt.strftime('%Y%m%d%H'))
+                                     t_m2dt.strftime('%Y%m%d%H') + numbersuffix)
                 gnout = os.path.join(c.inputdir, 'flux' +
-                                     t_m1dt.strftime('%Y%m%d%H'))
+                                     t_m1dt.strftime('%Y%m%d%H') + numbersuffix)
                 hnout = os.path.join(c.inputdir, 'flux' +
-                                     t_dt.strftime('%Y%m%d%H'))
+                                     t_dt.strftime('%Y%m%d%H') + numbersuffix)
 
             print("outputfile = " + fnout)
@@ -1370 +1396 @@
         # get the values of the keys which are used for distinct access
         # of grib messages via product
-        index_keys = ["date", "time", "step"]
+        if '/' in self.number:
+            # more than one ensemble member is selected
+            index_keys = ["number", "date", "time", "step"]
+        else:
+            index_keys = ["date", "time", "step"]
         iid, index_vals = self._mk_index_values(c.inputdir,
                                                 inputfiles,
@@ -1524 +1554 @@
             else:
                 suffix = cdate_hour[2:10]
+
+            # if necessary, add ensemble member number to filename suffix
+            if 'number' in index_keys:
+                index_number = index_keys.index('number')
+                if len(index_vals[index_number]) > 1:
+                    suffix = suffix + '.N{:0>3}'.format(int(prod[index_number]))
+
             fnout = os.path.join(c.inputdir, c.prefix + suffix)
             print("outputfile = " + fnout)
@@ -1553 +1590 @@
 
         codes_index_release(iid)
+
+        return
+
+
+    def calc_extra_elda(self, path, prefix):
+        ''' Calculates extra ensemble members for ELDA - Stream.
+
+        Parameters
+        ----------
+        path : str
+            Path to the output files.
+
+        prefix : str
+            The prefix of the output filenames as defined in Control file.
+
+        Return
+        ------
+
+        '''
+        # max number
+        maxnum = int(self.number.split('/')[-1])
+
+        # get a list of all prepared output files with control forecast (CF)
+        CF_filelist = UioFiles(path, prefix + '*.N000')
+
+        for cffile in CF_filelist.files:
+            with open(cffile, 'rb') as f:
+                cfvalues=[]
+                while True:
+                    fid = codes_grib_new_from_file(f)
+                    if fid is None:
+                        break
+                    cfvalues.append(codes_get_array(fid, 'values'))
+                    codes_release(fid)
+
+            filename = cffile.split('N000')[0]
+            for i in range(1, maxnum+1): # max number nehmen
+
+                # read an ensemble member
+                g = open(filename + 'N{:0>3}'.format(i), 'rb')
+                # create file for newly calculated ensemble member
+                h = open(filename + 'N{:0>3}'.format(i+maxnum), 'wb')
+                # number of message in grib file
+                j = 0
+                while True:
+                    gid = codes_grib_new_from_file(g)
+                    if gid is None:
+                        break
+                    values = codes_get_array(gid, 'values')
+                    codes_set_array(gid, 'values',
+                                    values-2*(values-cfvalues[j]))
+                    codes_set(gid, 'number', i+maxnum)
+                    codes_write(gid, h)
+                    codes_release(gid)
+                    j += 1
+
+                g.close()
+                h.close()
+                print('wrote ' + filename + 'N{:0>3}'.format(i+maxnum))
+                self.outputfilelist.append(
+                    os.path.basename(filename + 'N{:0>3}'.format(i+maxnum)))
 
         return

source/python/mods/get_mars_data.py

@@ -209 +209 @@
 
 
+def check_dates_for_nonflux_fc_times(types, times):
+    '''
+    '''
+    for ty, ti in zip(types,times):
+        if ty.upper() == 'FC' and int(ti) == 18:
+            return True
+    return False
+
+
 def mk_dates(c, fluxes):
     '''Prepares start and end date depending on flux or non flux data.
@@ -258 +267 @@
         start = start - timedelta(days=1)
         end = end + timedelta(days=1)
+
+    # if we have non-flux forecast data starting at 18 UTC
+    # we need to start retrieving data one day in advance
+    if not fluxes and check_dates_for_nonflux_fc_times(c.type, c.time):
+        start = start - timedelta(days=1)
 
     return start, end, chunk

source/python/mods/prepare_flexpart.py

@@ -177 +177 @@
     flexpart = EcFlexpart(c, fluxes=False)
     flexpart.create(inputfiles, c)
+    if c.stream.lower() == 'elda':
+        flexpart.calc_extra_elda(os.path.join(c.inputdir, c.prefix))
     flexpart.process_output(c)