Changeset 092aaf1 in flex_extract.git


Ignore:
Timestamp:
Dec 13, 2018, 11:53:49 PM (5 months ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
dev
Children:
5c0a578
Parents:
145fbe0
Message:

modularized the init function; implemented disaggregation with new interpolation method

File:
1 edited

Legend:

Unmodified
Added
Removed
  • source/python/classes/EcFlexpart.py

    re5884c9 r092aaf1  
    9090sys.path.append('../')
    9191import _config
    92 from GribTools import GribTools
     92from GribUtil import GribUtil
    9393from mods.tools import (init128, to_param_id, silent_remove, product,
    9494                        my_error, make_dir, get_informations, get_dimensions)
     
    125125
    126126        '''
    127         # set a counter for the number of mars requests generated
     127        # set a counter for the number of generated mars requests
    128128        self.mreq_count = 0
    129 
    130         # different mars types for retrieving data for flexpart
    131         self.types = dict()
    132 
    133         # Pure forecast mode
    134         if c.purefc:
    135             c.type = [c.type[0]]
    136             c.step = ['{:0>3}'.format(int(c.step[0]))]
    137             c.time = [c.time[0]]
    138             for i in range(1, c.maxstep + 1):
    139                 c.type.append(c.type[0])
    140                 c.step.append('{:0>3}'.format(i))
    141                 c.time.append(c.time[0])
    142129
    143130        self.inputdir = c.inputdir
     
    145132        self.basetime = c.basetime
    146133        self.dtime = c.dtime
    147         i = 0
    148         if fluxes and not c.purefc:
    149             self.types[str(c.acctype)] = {'times': str(c.acctime),
    150                                           'steps': '{}/to/{}/by/{}'.format(
    151                                               c.dtime, c.accmaxstep, c.dtime)}
    152         else:
    153             for ty, st, ti in zip(c.type, c.step, c.time):
    154                 btlist = range(24)
    155                 if c.basetime == '12':
    156                     btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
    157                 if c.basetime == '00':
    158                     btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
    159 
    160                 if ((ty.upper() == 'AN' and (int(c.time[i]) % int(c.dtime)) == 0) or
    161                     (ty.upper() != 'AN' and (int(c.step[i]) % int(c.dtime)) == 0 and
    162                      (int(c.step[i]) % int(c.dtime) == 0)) ) and \
    163                     (int(c.time[i]) in btlist or c.purefc):
    164 
    165                     if ty not in self.types.keys():
    166                         self.types[ty] = {'times': '', 'steps': ''}
    167 
    168                     if ti not in self.types[ty]['times']:
    169                         if self.types[ty]['times']:
    170                             self.types[ty]['times'] += '/'
    171                         self.types[ty]['times'] += ti
    172 
    173                     if st not in self.types[ty]['steps']:
    174                         if self.types[ty]['steps']:
    175                             self.types[ty]['steps'] += '/'
    176                         self.types[ty]['steps'] += st
    177                 i += 1
    178 
     134        self.acctype = c.acctype
     135        self.acctime = c.acctime
     136        self.accmaxstep = c.accmaxstep
    179137        self.marsclass = c.marsclass
    180138        self.stream = c.stream
     
    182140        self.resol = c.resol
    183141        self.accuracy = c.accuracy
     142        self.addpar = c.addpar
    184143        self.level = c.level
    185144        self.expver = c.expver
    186145        self.levelist = c.levelist
    187         # for gaussian grid retrieval
    188         self.glevelist = '1/to/' + c.level
     146        self.glevelist = '1/to/' + c.level # in case of gaussian grid
    189147        self.gaussian = c.gaussian
    190148        self.grid = c.grid
    191149        self.area = c.area
     150        self.purefc = c.purefc
    192151        self.outputfilelist = []
    193152
    194 
    195         # Now comes the nasty part that deals with the different
    196         # scenarios we have:
    197         # 1) Calculation of etadot on
    198         #    a) Gaussian grid
    199         #    b) Output grid
    200         #    c) Output grid using parameter 77 retrieved from MARS
    201         # 3) Calculation/Retrieval of omega
    202         # 4) Download also data for WRF
    203 
    204         # Different grids need different retrievals
    205         # SH = Spherical Harmonics, GG = Gaussian Grid,
    206         # OG = Output Grid, ML = MultiLevel, SL = SingleLevel
    207         self.params = {'SH__ML': '', 'SH__SL': '',
    208                        'GG__ML': '', 'GG__SL': '',
    209                        'OG__ML': '', 'OG__SL': '',
    210                        'OG_OROLSM_SL': '', 'OG_acc_SL': ''}
    211         # the self.params dictionary stores a list of
    212         # [param, levtype, levelist, grid] per key
    213 
    214         if fluxes is False:
    215             self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
    216             #                        "SD/MSL/TCC/10U/10V/2T/2D/129/172"
    217             self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \
    218                                      'SFC', '1', self.grid]
    219             if c.addpar:
    220                 if c.addpar[0] == '/':
    221                     c.addpar = c.addpar[1:]
    222                 self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
    223 
    224             if c.marsclass.upper() == 'EA' or c.marsclass.upper() == 'EP':
    225                 self.params['OG_OROLSM__SL'] = ["160/27/28/244",
    226                                                 'SFC', '1', self.grid]
    227             else:
    228                 self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
    229                                                 'SFC', '1', self.grid]
    230 
    231             self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
    232 
    233             #if c.gauss == '0' and c.eta == '1':
    234             if not c.gauss and c.eta:
    235                 # the simplest case
    236                 self.params['OG__ML'][0] += '/U/V/77'
    237             #elif c.gauss == '0' and c.eta == '0':
    238             elif not c.gauss and not c.eta:
    239             # this is not recommended (inaccurate)
    240                 self.params['OG__ML'][0] += '/U/V'
    241             #elif c.gauss == '1' and c.eta == '0':
    242             elif c.gauss and not c.eta:
    243                 # this is needed for data before 2008, or for reanalysis data
    244                 self.params['GG__SL'] = ['Q', 'ML', '1', \
    245                                          '{}'.format((int(self.resol) + 1) / 2)]
    246                 self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
    247             else:
    248                 print('Warning: This is a very costly parameter combination, \
    249                       use only for debugging!')
    250                 self.params['GG__SL'] = ['Q', 'ML', '1', \
    251                                          '{}'.format((int(self.resol) + 1) / 2)]
    252                 self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
    253                                          '{}'.format((int(self.resol) + 1) / 2)]
    254 
    255             if c.omega:
    256                 self.params['OG__ML'][0] += '/W'
    257 
    258             # add cloud water content if necessary
    259             if c.cwc:
    260                 self.params['OG__ML'][0] += '/CLWC/CIWC'
    261 
    262             # add vorticity and geopotential height for WRF if necessary
    263             if c.wrf:
    264                 self.params['OG__ML'][0] += '/Z/VO'
    265                 if '/D' not in self.params['OG__ML'][0]:
    266                     self.params['OG__ML'][0] += '/D'
    267                 #wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ /
    268                 #           stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper()
    269                 wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \
    270                            139/170/183/236/39/40/41/42'
    271                 lwrt_sfc = wrf_sfc.split('/')
    272                 for par in lwrt_sfc:
    273                     if par not in self.params['OG__SL'][0]:
    274                         self.params['OG__SL'][0] += '/' + par
    275 
     153        # Define the different types of field combinations (type, time, step)
     154        self.types = {}
     155        # Define the parameters and their level types, level list and
     156        # grid resolution for the retrieval job
     157        self.params = {}
     158
     159        if fluxes:
     160            self._create_params_fluxes()
    276161        else:
    277             self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
    278                                         'SFC', '1', self.grid]
    279 
     162            self._create_params(c.gauss, c.eta, c.omega, c.cwc, c.wrf)
     163
     164        if fluxes and not c.purefc:
     165            self._create_field_types_fluxes()
     166        else:
     167            self._create_field_types(c.type, c.time, c.step)
     168        return
     169
     170    def _create_field_types(self, ftype, ftime, fstep):
     171        '''Create the combination of field type, time and forecast step.
     172
     173        Parameters:
     174        -----------
     175        ftype : :obj:`list` of :obj:`string`
     176            List of field types.
     177
     178        ftime : :obj:`list` of :obj:`string`
     179            The time in hours of the field.
     180
     181        fstep : :obj:`string`
     182            Specifies the forecast time step from forecast base time.
     183            Valid values are hours (HH) from forecast base time.
     184
     185        Return
     186        ------
     187
     188        '''
     189        i = 0
     190        for ty, st, ti in zip(ftype, fstep, ftime):
     191            btlist = range(24)
     192            if self.basetime == '12':
     193                btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
     194            if self.basetime == '00':
     195                btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
     196
     197            # if ((ty.upper() == 'AN' and (int(c.time[i]) % int(c.dtime)) == 0) or
     198                # (ty.upper() != 'AN' and (int(c.step[i]) % int(c.dtime)) == 0 and
     199                 # (int(c.step[i]) % int(c.dtime) == 0)) ) and \
     200                 # (int(c.time[i]) in btlist or c.purefc):
     201
     202            if (int(ti) in btlist) or self.purefc:
     203
     204                if ((ty.upper() == 'AN' and (int(ti) % int(self.dtime)) == 0) or
     205                    (ty.upper() != 'AN' and (int(st) % int(self.dtime)) == 0)):
     206
     207                    if ty not in self.types.keys():
     208                        self.types[ty] = {'times': '', 'steps': ''}
     209
     210                    if ti not in self.types[ty]['times']:
     211                        if self.types[ty]['times']:
     212                            self.types[ty]['times'] += '/'
     213                        self.types[ty]['times'] += ti
     214
     215                    if st not in self.types[ty]['steps']:
     216                        if self.types[ty]['steps']:
     217                            self.types[ty]['steps'] += '/'
     218                        self.types[ty]['steps'] += st
     219            i += 1
     220        return
     221
     222    def _create_field_types_fluxes(self):
     223        '''Create the combination of field type, time and forecast step
     224        for the flux data.
     225
     226        Parameters:
     227        -----------
     228
     229        Return
     230        ------
     231
     232        '''
     233        self.types[str(self.acctype)] = {'times': str(self.acctime),
     234                                         'steps': '{}/to/{}/by/{}'.format(
     235                                             self.dtime,
     236                                             self.accmaxstep,
     237                                             self.dtime)}
     238        return
     239
     240    def _create_params(self, gauss, eta, omega, cwc, wrf):
     241        '''Define the specific parameter settings for retrievment.
     242
     243        The different parameters need specific grid types and level types
     244        for retrievement. We might get following combination of types
     245        (depending on selection and availability):
     246        (These are short cuts for the grib file names (leading sequence)
     247        SH__ML, OG__ML, GG__ML, SH__SL, OG__SL, GG__SL, OG_OROLSM_SL
     248        where:
     249            SH = Spherical Harmonics, GG = Gaussian Grid, OG = Output Grid,
     250            ML = Model Level, SL = Surface Level
     251
     252        For each of this combination there is a list of parameter names,
     253        the level type, the level list and the grid resolution.
     254
     255        There are different scenarios for data extraction from MARS:
     256        1) Retrieval of etadot
     257           eta=1, gauss=0, omega=0
     258        2) Calculation of etadot from divergence
     259           eta=0, gauss=1, omega=0
     260        3) Calculation of etadot from omega (for makes sense for debugging)
     261           eta=0, gauss=0, omega=1
     262        4) Retrieval and Calculation of etadot (only for debugging)
     263           eta=1, gauss=1, omega=0
     264        5) Download also specific model and surface level data for FLEXPART-WRF
     265
     266        Parameters:
     267        -----------
     268        gauss : :obj:`integer`
     269            Gaussian grid is retrieved.
     270
     271        eta : :obj:`integer`
     272            Etadot parameter will be directly retrieved.
     273
     274        omega : :obj:`integer`
     275            The omega paramterwill be retrieved.
     276
     277        cwc : :obj:`integer`
     278            The cloud liquid and ice water content will be retrieved.
     279
     280        wrf : :obj:`integer`
     281            Additional model level and surface level data will be retrieved for
     282            WRF/FLEXPART-WRF simulations.
     283
     284        Return
     285        ------
     286
     287        '''
     288        # SURFACE FIELDS
     289        #-----------------------------------------------------------------------
     290        self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
     291        self.params['OG__SL'] = ['SD/MSL/TCC/10U/10V/2T/2D/Z/LSM', \
     292                                 'SFC', '1', self.grid]
     293        if self.addpar:
     294            self.params['OG__SL'][0] += self.addpar
     295
     296        if self.marsclass.upper() == 'EA' or self.marsclass.upper() == 'EP':
     297            self.params['OG_OROLSM__SL'] = ["SDOR/CVL/CVH/FSR",
     298                                            'SFC', '1', self.grid]
     299        else:
     300            self.params['OG_OROLSM__SL'] = ["SDOR/CVL/CVH/SR", \
     301                                            'SFC', '1', self.grid]
     302
     303        # MODEL LEVEL FIELDS
     304        #-----------------------------------------------------------------------
     305        self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
     306
     307        if not gauss and eta:
     308            self.params['OG__ML'][0] += '/U/V/ETADOT'
     309        elif gauss and not eta:
     310            self.params['GG__SL'] = ['Q', 'ML', '1', \
     311                                     '{}'.format((int(self.resol) + 1) / 2)]
     312            self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
     313        elif not gauss and not eta:
     314            self.params['OG__ML'][0] += '/U/V'
     315        else:
     316            print('Warning: Collecting etadot and parameters for gaussian grid \
     317                            is a very costly parameter combination, \
     318                            use this combination only for debugging!')
     319            self.params['GG__SL'] = ['Q', 'ML', '1', \
     320                                     '{}'.format((int(self.resol) + 1) / 2)]
     321            self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \
     322                                     '{}'.format((int(self.resol) + 1) / 2)]
     323
     324        if omega:
     325            self.params['OG__ML'][0] += '/W'
     326
     327        if cwc:
     328            self.params['OG__ML'][0] += '/CLWC/CIWC'
     329
     330        # ADDITIONAL FIELDS FOR FLEXPART-WRF MODEL (IF QUESTIONED)
     331        #-----------------------------------------------------------------------
     332        if wrf:
     333            self.params['OG__ML'][0] += '/Z/VO'
     334            if '/D' not in self.params['OG__ML'][0]:
     335                self.params['OG__ML'][0] += '/D'
     336            wrf_sfc = ['SP','SKT','SST','CI','STL1','STL2', 'STL3','STL4',
     337                       'SWVL1','SWVL2','SWVL3','SWVL4']
     338            for par in wrf_sfc:
     339                if par not in self.params['OG__SL'][0]:
     340                    self.params['OG__SL'][0] += '/' + par
     341
     342        return
     343
     344
     345    def _create_params_fluxes(self):
     346        '''Define the parameter setting for flux data.
     347
     348        Flux data are accumulated fields in time and are stored on the
     349        surface level. The leading short cut name for the grib files is:
     350        "OG_acc_SL" with OG for Regular Output Grid, SL for Surface Level, and
     351        acc for Accumulated Grid.
     352        The params dictionary stores a list of parameter names, the level type,
     353        the level list and the grid resolution.
     354
     355        The flux data are: LSP/CP/SSHF/EWSS/NSSS/SSR
     356
     357        Parameters:
     358        -----------
     359
     360        Return
     361        ------
     362
     363        '''
     364        self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \
     365                                    'SFC', '1', self.grid]
    280366        return
    281367
     
    401487        indexfile = os.path.join(inputdir, _config.FILE_GRIB_INDEX)
    402488        silent_remove(indexfile)
    403         grib = GribTools(inputfiles.files)
     489        grib = GribUtil(inputfiles.files)
    404490        # creates new index file
    405491        iid = grib.index(index_keys=index_keys, index_file=indexfile)
     
    749835        # index_vals looks like e.g.:
    750836        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
    751         # index_vals[1]: ('0', '1200', '1800', '600') ; time
     837        # index_vals[1]: ('0', '600', '1200', '1800') ; time
    752838        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
    753839
    754840        if c.rrint:
    755             start_date = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
    756             end_date = datetime.strptime(c.end_date + '23', '%Y%m%d%H')
     841            if not c.purefc:
     842                start_date = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
     843                end_date = datetime.strptime(c.end_date + '23', '%Y%m%d%H')
     844            else:
     845                sdate_str = c.start_date + '{:0>2}'.format(index_vals[1][0])
     846                start_date = datetime.strptime(sdate_str, '%Y%m%d%H')
     847                edate_str = c.end_date + '{:0>2}'.format(index_vals[1][-1])
     848                end_date = datetime.strptime(edate_str, '%Y%m%d%H')
     849                end_date = end_date + timedelta(hours=c.maxstep)
     850
    757851            info = get_informations(os.path.join(c.inputdir,
    758852                                                 inputfiles.files[0]))
    759             dims = get_dimensions(c, info)
     853            dims = get_dimensions(info, c.purefc, c.dtime, index_vals,
     854                                  start_date, end_date)
    760855            # create numpy array
    761856            lsp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
     
    763858            it_lsp = 0
    764859            it_cp = 0
     860            date_list = []
     861            step_list = []
    765862
    766863        # initialize dictionaries to store values
     
    792889            # create correct timestamp from the three time informations
    793890            cdate = str(codes_get(gid, 'date'))
    794             ctime = '{:0>2}'.format(codes_get(gid, 'time')/100)
    795             cstep = '{:0>3}'.format(codes_get(gid, 'step'))
     891            time = codes_get(gid, 'time')/100 # integer
     892            step = codes_get(gid, 'step') # integer
     893            ctime = '{:0>2}'.format(time)
     894            cstep = '{:0>3}'.format(step)
     895
    796896            t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
    797             t_dt = t_date + timedelta(hours=int(cstep))
    798             t_m1dt = t_date + timedelta(hours=int(cstep)-int(c.dtime))
    799             t_m2dt = t_date + timedelta(hours=int(cstep)-2*int(c.dtime))
     897            t_dt = t_date + timedelta(hours=step)
     898            t_m1dt = t_date + timedelta(hours=step-int(c.dtime))
     899            t_m2dt = t_date + timedelta(hours=step-2*int(c.dtime))
    800900            t_enddate = None
    801901
     
    823923            g_handle = open(gnout, 'w')
    824924
    825             # read message for message and store relevant data fields
     925            # read message for message and store relevant data fields, where
    826926            # data keywords are stored in pars
    827927            while True:
    828928                if not gid:
    829929                    break
    830                 parId = codes_get(gid, 'paramId')
    831                 step = codes_get(gid, 'step')
    832                 time = codes_get(gid, 'time')
    833                 ni = codes_get(gid, 'Ni')
    834                 nj = codes_get(gid, 'Nj')
     930                parId = codes_get(gid, 'paramId') # integer
     931                step = codes_get(gid, 'step') # integer
     932                time = codes_get(gid, 'time') # integer
     933                ni = codes_get(gid, 'Ni') # integer
     934                nj = codes_get(gid, 'Nj') # integer
    835935                if parId not in orig_vals.keys():
    836936                    # parameter is not a flux, find next one
     
    861961                # store precipitation if new disaggregation method is selected
    862962                # only the exact days are needed
    863                 if start_date <= t_dt <= end_date:
    864                     print 'DATE: ', t_dt
    865                     if c.rrint and parId == 142:
    866                         print len(deac_vals[parId])
    867                         lsp_np[:,it_lsp] = deac_vals[parId][-1][:]
    868                         it_lsp += 1
    869                     elif c.rrint and parId == 143:
    870                         cp_np[:,it_cp] = deac_vals[parId][-1][:]
    871                         it_cp += 1
     963                if c.rrint:
     964                    if start_date <= t_dt <= end_date:
     965                        if not c.purefc:
     966                            if t_dt not in date_list:
     967                                date_list.append(t_dt)
     968                                step_list = [0]
     969                        else:
     970                            if t_date not in date_list:
     971                                date_list.append(t_date)
     972                            if step not in step_list:
     973                                step_list.append(step)
     974                        if c.rrint and parId == 142:
     975                            lsp_np[:,it_lsp] = deac_vals[parId][-1][:]
     976                            it_lsp += 1
     977                        elif c.rrint and parId == 143:
     978                            cp_np[:,it_cp] = deac_vals[parId][-1][:]
     979                            it_cp += 1
    872980
    873981                # information printout
     
    875983
    876984                # length of deac_vals[parId] corresponds to the
    877                 # number of time steps (max. 4)
     985                # number of time steps, max. 4 are needed for disaggegration
     986                # with the old and original method
    878987                # run over all grib messages and perform
    879                 # shifting in time for old disaggegration method
     988                # shifting in time
    880989                if len(deac_vals[parId]) >= 3:
    881990                    if len(deac_vals[parId]) > 3:
     
    9241033                            if step == c.maxstep and c.purefc or \
    9251034                               t_dt == t_enddate:
    926 
    927                                 values = deac_vals[parId][3]
    928                                 codes_set_values(gid, values)
    929                                 codes_set(gid, 'stepRange', 0)
    930                                 truedatetime = t_m2dt + timedelta(hours=
    931                                                                   2*int(c.dtime))
    932                                 codes_set(gid, 'time', truedatetime.hour * 100)
    933                                 codes_set(gid, 'date', truedatetime.year * 10000 +
    934                                           truedatetime.month * 100 +
    935                                           truedatetime.day)
    936                                 codes_write(gid, h_handle)
    937 
    938                                 #values = (deac_vals[parId][1]+deac_vals[parId][2])/2.
     1035                                # last step
     1036                                if c.purefc:
     1037                                    values = deac_vals[parId][3]
     1038                                    codes_set_values(gid, values)
     1039                                    codes_set(gid, 'stepRange', step)
     1040                                    #truedatetime = t_m2dt + timedelta(hours=2*int(c.dtime))
     1041                                    codes_write(gid, h_handle)
     1042                                else:
     1043                                    values = deac_vals[parId][3]
     1044                                    codes_set_values(gid, values)
     1045                                    codes_set(gid, 'stepRange', 0)
     1046                                    truedatetime = t_m2dt + timedelta(hours=2*int(c.dtime))
     1047                                    codes_set(gid, 'time', truedatetime.hour * 100)
     1048                                    codes_set(gid, 'date', int(truedatetime.strftime('%Y%m%d')))
     1049                                    codes_write(gid, h_handle)
     1050
    9391051                                if parId == 142 or parId == 143:
    9401052                                    values = disaggregation.darain(list(reversed(deac_vals[parId])))
     
    9421054                                    values = disaggregation.dapoly(list(reversed(deac_vals[parId])))
    9431055
    944                                 codes_set(gid, 'stepRange', 0)
    945                                 truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
    946                                 codes_set(gid, 'time', truedatetime.hour * 100)
    947                                 codes_set(gid, 'date', truedatetime.year * 10000 +
    948                                           truedatetime.month * 100 +
    949                                           truedatetime.day)
    950                                 codes_set_values(gid, values)
    951                                 codes_write(gid, g_handle)
     1056                                # step before last step
     1057                                if c.purefc:
     1058                                    codes_set(gid, 'stepRange', step-int(c.dtime))
     1059                                    #truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
     1060                                    codes_set_values(gid, values)
     1061                                    codes_write(gid, g_handle)
     1062                                else:
     1063                                    codes_set(gid, 'stepRange', 0)
     1064                                    truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
     1065                                    codes_set(gid, 'time', truedatetime.hour * 100)
     1066                                    codes_set(gid, 'date', int(truedatetime.strftime('%Y%m%d')))
     1067                                    codes_set_values(gid, values)
     1068                                    codes_write(gid, g_handle)
    9521069
    9531070                codes_release(gid)
     
    9611078        codes_index_release(iid)
    9621079
    963         # disaggregation.IA3(lsp_np)
    964         # pro gitterpunkt die zeitserie disagg, oder geht es auch mit feld
    965         # wenn zurück dann über zeitschritte itterieren und rausschreiben.
    966         # original in die flux
    967         # subgrid wohin? extrafiles die beide oder jeweils einen subgrid beinhalten?
    968         # kann time hh und mm fassen?
    969 
    970         print lsp_np.shape
    971         print lsp_np
    972 
     1080        if c.rrint:
     1081            self._create_rr_grib_dummy(inputfiles.files[0], c.inputdir)
     1082
     1083            self._prep_new_rrint(ni, nj, lsp_np.shape[1], lsp_np,
     1084                                 cp_np, date_list, step_list, c)
    9731085
    9741086        return
    9751087
     1088    def _prep_new_rrint(self, ni, nj, nt, lsp_np, cp_np, date_list, step_list, c):
     1089        '''Calculates and writes out the disaggregated precipitation fields.
     1090
     1091        Disaggregation is done in time and original times are written to the
     1092        flux files, while the additional subgrid times are written to
     1093        extra files.
     1094
     1095        Parameters
     1096        ----------
     1097        ni : :obj:`integer`
     1098            Amount of zonal grid points.
     1099
     1100        nj : :obj:`integer`
     1101            Amount of meridional grid points.
     1102
     1103        nt : :obj:`integer`
     1104            Number of time steps.
     1105
     1106        lsp_np : :obj:`numpy array` of :obj:`float`
     1107            The large scale precipitation fields for each time step.
     1108            Shape (ni * nj, nt).
     1109
     1110        cp_np : :obj:`numpy array` of :obj:`float`
     1111            The convective precipitation fields for each time step.
     1112            Shape (ni * nj, nt).
     1113
     1114        date_list : :obj:`list` of :obj:`datetime`
     1115            The list of dates for which the disaggregation is to be done.
     1116
     1117        step_list : :obj:`list` of :obj:`integer`
     1118            The list of steps for a single forecast time.
     1119            Only necessary for pure forecasts.
     1120
     1121        c : :obj:`ControlFile`
     1122            Contains all the parameters of CONTROL file and
     1123            command line.
     1124
     1125        Return
     1126        ------
     1127
     1128        '''
     1129        print('... disaggregation or precipitation with new method.')
     1130        lsp_new_np = np.zeros((ni * nj, nt * 3), dtype=np.float64)
     1131        cp_new_np = np.zeros((ni * nj, nt * 3), dtype=np.float64)
     1132
     1133        # do the disaggregation, but neglect the last value of the
     1134        # time series. This one corresponds for example to 24 hour,
     1135        # which we don't need.
     1136        for ix in range(ni*nj):
     1137            lsp_new_np[ix,:] = disaggregation.IA3(lsp_np[ix,:])[:-1]
     1138            cp_new_np[ix,:] = disaggregation.IA3(cp_np[ix,:])[:-1]
     1139
     1140        # write to grib files (full/orig times to flux file and inbetween
     1141        # times into seperate end files)
     1142        print('... write disaggregated precipitation to files.')
     1143        it = 0
     1144        for date in date_list:
     1145            for step in step_list:
     1146                tmpfile = os.path.join(c.inputdir, 'rr_grib_dummy.grb')
     1147
     1148                if c.purefc:
     1149                    fluxfilename = 'flux' + date.strftime('%Y%m%d.%H') + \
     1150                                   '.{:0>3}'.format(step)
     1151                    filename1 = c.prefix + date.strftime('%y%m%d.%H') + \
     1152                                '.{:0>3}'.format(step) + '_1'
     1153                    filename2 = c.prefix + date.strftime('%y%m%d.%H') + \
     1154                                '.{:0>3}'.format(step) + '_2'
     1155                else:
     1156                    fluxfilename = 'flux' + date.strftime('%Y%m%d%H')
     1157                    filename1 = c.prefix + date.strftime('%y%m%d%H') + '_1'
     1158                    filename2 = c.prefix + date.strftime('%y%m%d%H') + '_2'
     1159
     1160                # collect for final processing
     1161                self.outputfilelist.append(os.path.basename(fluxfilename))
     1162                self.outputfilelist.append(os.path.basename(filename1))
     1163                self.outputfilelist.append(os.path.basename(filename2))
     1164
     1165                # write original time step to flux file as usual
     1166                fluxfile = GribUtil(os.path.join(c.inputdir, fluxfilename))
     1167                fluxfile.set_keys(tmpfile, filemode='a', strict=True,
     1168                                  wherekeynames=['paramId'], wherekeyvalues=[142],
     1169                                  keynames=['date','time','stepRange','values'],
     1170                                  keyvalues=[int(date.strftime('%Y%m%d')),
     1171                                             date.hour*100, step, lsp_new_np[:,it]],
     1172                                 )
     1173                fluxfile.set_keys(tmpfile, filemode='a', strict=True,
     1174                                  wherekeynames=['paramId'], wherekeyvalues=[143],
     1175                                  keynames=['date','time','stepRange','values'],
     1176                                  keyvalues=[int(date.strftime('%Y%m%d')),
     1177                                             date.hour*100, step, cp_new_np[:,it]]
     1178                                 )
     1179
     1180                # write first subgrid time step
     1181                endfile1 = GribUtil(os.path.join(c.inputdir, filename1))
     1182                endfile1.set_keys(tmpfile, filemode='w', strict=True,
     1183                                  wherekeynames=['paramId'], wherekeyvalues=[142],
     1184                                  keynames=['date','time','stepRange','values'],
     1185                                  keyvalues=[int(date.strftime('%Y%m%d')),
     1186                                             date.hour*100, step, lsp_new_np[:,it+1]]
     1187                                  )
     1188                endfile1.set_keys(tmpfile, filemode='a', strict=True,
     1189                                  wherekeynames=['paramId'], wherekeyvalues=[143],
     1190                                  keynames=['date','time','stepRange','values'],
     1191                                  keyvalues=[int(date.strftime('%Y%m%d')),
     1192                                             date.hour*100, step, cp_new_np[:,it+1]]
     1193                                 )
     1194
     1195                # write second subgrid time step
     1196                endfile2 = GribUtil(os.path.join(c.inputdir, filename2))
     1197                endfile2.set_keys(tmpfile, filemode='w', strict=True,
     1198                                  wherekeynames=['paramId'], wherekeyvalues=[142],
     1199                                  keynames=['date','time','stepRange','values'],
     1200                                  keyvalues=[int(date.strftime('%Y%m%d')),
     1201                                             date.hour*100, step, lsp_new_np[:,it+2]]
     1202                                 )
     1203                endfile2.set_keys(tmpfile, filemode='a', strict=True,
     1204                                  wherekeynames=['paramId'], wherekeyvalues=[143],
     1205                                  keynames=['date','time','stepRange','values'],
     1206                                  keyvalues=[int(date.strftime('%Y%m%d')),
     1207                                             date.hour*100, step, cp_new_np[:,it+2]]
     1208                                 )
     1209                it = it + 3 # jump to next original time step
     1210        return
     1211
     1212    def _create_rr_grib_dummy(self, ifile, inputdir):
     1213        '''Creates a grib file with a dummy message for the two precipitation
     1214        types lsp and cp each.
     1215
     1216        Parameters
     1217        ----------
     1218        ifile : :obj:`string`
     1219            Filename of the input file to read the grib messages from.
     1220
     1221        inputdir : :obj:`string`, optional
     1222            Path to the directory where the retrieved data is stored.
     1223
     1224        Return
     1225        ------
     1226
     1227        '''
     1228
     1229        gribfile = GribUtil(os.path.join(inputdir,'rr_grib_dummy.grb'))
     1230
     1231        gribfile.copy_dummy_msg(ifile, keynames=['paramId'],
     1232                      keyvalues=[142], filemode='w')
     1233
     1234        gribfile.copy_dummy_msg(ifile, keynames=['paramId'],
     1235                      keyvalues=[143], filemode='a')
     1236
     1237        return
    9761238
    9771239    def create(self, inputfiles, c):
     
    10291291        # index_vals looks like e.g.:
    10301292        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
    1031         # index_vals[1]: ('0', '1200', '1800', '600') ; time
     1293        # index_vals[1]: ('0', '600', '1200', '1800') ; time
    10321294        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
    10331295
     
    11171379                        codes_set(gid, 'paramId', 201031)
    11181380                        codes_write(gid, fdict['22'])
     1381                        scwc = None
    11191382                elif c.wrf and paramId in [129, 138, 155] and \
    11201383                      levtype == 'hybrid': # Z, VO, D
     
    13161579        # read template COMMAND file
    13171580        with open(os.path.expandvars(os.path.expanduser(
    1318             c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f:
     1581            c.flexpartdir)) + '/../Options/COMMAND', 'r') as f:
    13191582            lflist = f.read().split('\n')
    13201583
     
    13421605        os.chdir(c.outputdir)
    13431606        p = subprocess.check_call([
    1344             os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))
     1607            os.path.expandvars(os.path.expanduser(c.flexpartdir))
    13451608            + '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.'])
    13461609        os.chdir(pwd)
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG