Changeset c274d9a in flex_extract.git


Ignore:
Timestamp:
Dec 7, 2018, 3:30:13 PM (3 months ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
dev
Children:
4d3b052
Parents:
d2febd4
Message:

added rrint param and updated deacc_fluxes function to store rain fields

Location:
source/python/classes
Files:
2 edited

Legend:

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

    r3f36e42 rc274d9a  
    150150        self.public = 0
    151151        self.ecapi = None
     152        self.rrint = 0
    152153
    153154        self.logicals = ['gauss', 'omega', 'omegadiff', 'eta', 'etadiff',
    154155                         'dpdeta', 'cwc', 'wrf', 'grib2flexpart', 'ecstorage',
    155                          'ectrans', 'debug', 'request', 'public']
     156                         'ectrans', 'debug', 'request', 'public', 'rrint']
    156157
    157158        self.__read_controlfile__()
     
    478479            self.accmaxstep='12'
    479480
    480 
    481481        self.grid = check_grid(self.grid)
    482482
  • source/python/classes/EcFlexpart.py

    r3f36e42 rc274d9a  
    9292from GribTools import GribTools
    9393from mods.tools import (init128, to_param_id, silent_remove, product,
    94                         my_error, make_dir)
     94                        my_error, make_dir, get_informations, get_dimensions)
    9595from MarsRetrieval import MarsRetrieval
    9696import mods.disaggregation as disaggregation
     
    724724        ----------
    725725        inputfiles : :obj:`UioFiles`
    726             Contains a list of files.
     726            Contains the list of files that contain flux data.
    727727
    728728        c : :obj:`ControlFile`
     
    752752        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
    753753
    754         # initialize dictionaries
    755         valsdict = {}
    756         svalsdict = {}
     754        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')
     757            info = get_informations(os.path.join(c.inputdir,
     758                                                 inputfiles.files[0]))
     759            dims = get_dimensions(c, info)
     760            # create numpy array
     761            lsp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
     762            cp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
     763            it_lsp = 0
     764            it_cp = 0
     765
     766        # initialize dictionaries to store values
     767        orig_vals = {}
     768        deac_vals = {}
    757769        for p in pars:
    758             valsdict[str(p)] = []
    759             svalsdict[str(p)] = []
     770            orig_vals[p] = []
     771            deac_vals[p] = []
    760772
    761773        # "product" genereates each possible combination between the
     
    813825            # read message for message and store relevant data fields
    814826            # data keywords are stored in pars
    815             while 1:
     827            while True:
    816828                if not gid:
    817829                    break
    818                 cparamId = str(codes_get(gid, 'paramId'))
     830                parId = codes_get(gid, 'paramId')
    819831                step = codes_get(gid, 'step')
    820832                time = codes_get(gid, 'time')
    821833                ni = codes_get(gid, 'Ni')
    822834                nj = codes_get(gid, 'Nj')
    823                 if cparamId in valsdict.keys():
    824                     values = codes_get_values(gid)
    825                     vdp = valsdict[cparamId]
    826                     svdp = svalsdict[cparamId]
    827 
    828                     if cparamId == '142' or cparamId == '143':
    829                         fak = 1. / 1000.
     835                if parId not in orig_vals.keys():
     836                    # parameter is not a flux, find next one
     837                    continue
     838
     839                # define conversion factor
     840                if parId == 142 or parId == 143:
     841                    fak = 1. / 1000.
     842                else:
     843                    fak = 3600.
     844
     845                # get parameter values and reshape
     846                values = codes_get_values(gid)
     847                values = (np.reshape(values, (nj, ni))).flatten() / fak
     848
     849                # save the original and accumulated values
     850                orig_vals[parId].append(values[:])
     851
     852                if c.marsclass.upper() == 'EA' or step <= int(c.dtime):
     853                    # no de-accumulation needed
     854                    deac_vals[parId].append(values[:] / int(c.dtime))
     855                else:
     856                    # do de-accumulation
     857                    deac_vals[parId].append(
     858                        (orig_vals[parId][-1] - orig_vals[parId][-2]) /
     859                         int(c.dtime))
     860
     861                # store precipitation if new disaggregation method is selected
     862                # 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
     872
     873                # information printout
     874                print(parId, time, step, len(values), values[0], np.std(values))
     875
     876                # length of deac_vals[parId] corresponds to the
     877                # number of time steps (max. 4)
     878                # run over all grib messages and perform
     879                # shifting in time for old disaggegration method
     880                if len(deac_vals[parId]) >= 3:
     881                    if len(deac_vals[parId]) > 3:
     882                        if not c.rrint and (parId == 142 or parId == 143):
     883                            values = disaggregation.darain(deac_vals[parId])
     884                        else:
     885                            values = disaggregation.dapoly(deac_vals[parId])
     886
     887                        if not (step == c.maxstep and c.maxstep > 12 \
     888                                or t_dt == t_enddate):
     889                            # remove first time step in list to shift
     890                            # time line
     891                            orig_vals[parId].pop(0)
     892                            deac_vals[parId].pop(0)
    830893                    else:
    831                         fak = 3600.
    832 
    833                     values = (np.reshape(values, (nj, ni))).flatten() / fak
    834                     vdp.append(values[:])  # save the accumulated values
    835                     if c.marsclass.upper() == 'EA' or \
    836                        step <= int(c.dtime):
    837                         svdp.append(values[:] / int(c.dtime))
    838                     else:  # deaccumulate values
    839                         svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
    840 
    841                     print(cparamId, time, step, len(values),
    842                           values[0], np.std(values))
    843 
    844                     # len(svdp) corresponds to the time
    845                     if len(svdp) >= 3:
    846                         if len(svdp) > 3:
    847                             if cparamId == '142' or cparamId == '143':
    848                                 values = disaggregation.darain(svdp)
    849                             else:
    850                                 values = disaggregation.dapoly(svdp)
    851 
    852                             if not (step == c.maxstep and c.maxstep > 12 \
    853                                     or t_dt == t_enddate):
    854                                 vdp.pop(0)
    855                                 svdp.pop(0)
     894                        # if the third time step is read (per parId),
     895                        # write out the first one as a boundary value
     896                        if c.maxstep > 12:
     897                            values = deac_vals[parId][1]
    856898                        else:
    857                             if c.maxstep > 12:
    858                                 values = svdp[1]
    859                             else:
    860                                 values = svdp[0]
    861 
     899                            values = deac_vals[parId][0]
     900
     901                    if not (c.rrint and (parId == 142 or parId == 143)):
    862902                        codes_set_values(gid, values)
    863903
     
    872912
    873913                        if c.basetime:
    874                             t_enddate = datetime.strptime(c.end_date +
    875                                                           c.basetime,
     914                            t_enddate = datetime.strptime(c.end_date + c.basetime,
    876915                                                          '%Y%m%d%H')
    877916                        else:
    878917                            t_enddate = t_date + timedelta(2*int(c.dtime))
    879918
    880                         # squeeze out information of last two steps contained
    881                         # in svdp
    882                         # if step+int(c.dtime) == c.maxstep and c.maxstep>12
    883                         # or t_dt+timedelta(hours = int(c.dtime))
    884                         # >= t_enddate:
    885                         # Note that svdp[0] has not been popped in this case
    886 
    887                         if step == c.maxstep and c.maxstep > 12 or \
    888                            t_dt == t_enddate:
    889 
    890                             values = svdp[3]
    891                             codes_set_values(gid, values)
    892                             codes_set(gid, 'stepRange', 0)
    893                             truedatetime = t_m2dt + timedelta(hours=
    894                                                               2*int(c.dtime))
    895                             codes_set(gid, 'time', truedatetime.hour * 100)
    896                             codes_set(gid, 'date', truedatetime.year * 10000 +
    897                                       truedatetime.month * 100 +
    898                                       truedatetime.day)
    899                             codes_write(gid, h_handle)
    900 
    901                             #values = (svdp[1]+svdp[2])/2.
    902                             if cparamId == '142' or cparamId == '143':
    903                                 values = disaggregation.darain(list(reversed(svdp)))
    904                             else:
    905                                 values = disaggregation.dapoly(list(reversed(svdp)))
    906 
    907                             codes_set(gid, 'stepRange', 0)
    908                             truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
    909                             codes_set(gid, 'time', truedatetime.hour * 100)
    910                             codes_set(gid, 'date', truedatetime.year * 10000 +
    911                                       truedatetime.month * 100 +
    912                                       truedatetime.day)
    913                             codes_set_values(gid, values)
    914                             codes_write(gid, g_handle)
     919                            # squeeze out information of last two steps contained
     920                            # in deac_vals[parId]
     921                            # if step+int(c.dtime) == c.maxstep and c.maxstep>12
     922                            # or t_dt+timedelta(hours = int(c.dtime))
     923                            # >= t_enddate:
     924                            # Note that deac_vals[parId][0] has not been popped in this case
     925
     926                            if step == c.maxstep and c.maxstep > 12 or \
     927                               t_dt == t_enddate:
     928
     929                                values = deac_vals[parId][3]
     930                                codes_set_values(gid, values)
     931                                codes_set(gid, 'stepRange', 0)
     932                                truedatetime = t_m2dt + timedelta(hours=
     933                                                                  2*int(c.dtime))
     934                                codes_set(gid, 'time', truedatetime.hour * 100)
     935                                codes_set(gid, 'date', truedatetime.year * 10000 +
     936                                          truedatetime.month * 100 +
     937                                          truedatetime.day)
     938                                codes_write(gid, h_handle)
     939
     940                                #values = (deac_vals[parId][1]+deac_vals[parId][2])/2.
     941                                if parId == 142 or parId == 143:
     942                                    values = disaggregation.darain(list(reversed(deac_vals[parId])))
     943                                else:
     944                                    values = disaggregation.dapoly(list(reversed(deac_vals[parId])))
     945
     946                                codes_set(gid, 'stepRange', 0)
     947                                truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
     948                                codes_set(gid, 'time', truedatetime.hour * 100)
     949                                codes_set(gid, 'date', truedatetime.year * 10000 +
     950                                          truedatetime.month * 100 +
     951                                          truedatetime.day)
     952                                codes_set_values(gid, values)
     953                                codes_write(gid, g_handle)
    915954
    916955                codes_release(gid)
     
    923962
    924963        codes_index_release(iid)
     964
     965        # disaggregation.IA3(lsp_np)
     966        # pro gitterpunkt die zeitserie disagg, oder geht es auch mit feld
     967        # wenn zurück dann über zeitschritte itterieren und rausschreiben.
     968        # original in die flux
     969        # subgrid wohin? extrafiles die beide oder jeweils einen subgrid beinhalten?
     970        # kann time hh und mm fassen?
     971
     972        print lsp_np.shape
     973        print lsp_np
     974
    925975
    926976        return
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG