Changeset c274d9a in flex_extract.git
 Timestamp:
 Dec 7, 2018, 3:30:13 PM (15 months ago)
 Branches:
 dev
 Children:
 4d3b052
 Parents:
 d2febd4
 Location:
 source/python/classes
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

source/python/classes/ControlFile.py
r3f36e42 rc274d9a 150 150 self.public = 0 151 151 self.ecapi = None 152 self.rrint = 0 152 153 153 154 self.logicals = ['gauss', 'omega', 'omegadiff', 'eta', 'etadiff', 154 155 'dpdeta', 'cwc', 'wrf', 'grib2flexpart', 'ecstorage', 155 'ectrans', 'debug', 'request', 'public' ]156 'ectrans', 'debug', 'request', 'public', 'rrint'] 156 157 157 158 self.__read_controlfile__() … … 478 479 self.accmaxstep='12' 479 480 480 481 481 self.grid = check_grid(self.grid) 482 482 
source/python/classes/EcFlexpart.py
r3f36e42 rc274d9a 92 92 from GribTools import GribTools 93 93 from mods.tools import (init128, to_param_id, silent_remove, product, 94 my_error, make_dir )94 my_error, make_dir, get_informations, get_dimensions) 95 95 from MarsRetrieval import MarsRetrieval 96 96 import mods.disaggregation as disaggregation … … 724 724  725 725 inputfiles : :obj:`UioFiles` 726 Contains a list of files.726 Contains the list of files that contain flux data. 727 727 728 728 c : :obj:`ControlFile` … … 752 752 # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange 753 753 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 = {} 757 769 for p in pars: 758 valsdict[str(p)] = []759 svalsdict[str(p)] = []770 orig_vals[p] = [] 771 deac_vals[p] = [] 760 772 761 773 # "product" genereates each possible combination between the … … 813 825 # read message for message and store relevant data fields 814 826 # data keywords are stored in pars 815 while 1:827 while True: 816 828 if not gid: 817 829 break 818 cparamId = str(codes_get(gid, 'paramId'))830 parId = codes_get(gid, 'paramId') 819 831 step = codes_get(gid, 'step') 820 832 time = codes_get(gid, 'time') 821 833 ni = codes_get(gid, 'Ni') 822 834 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 deaccumulation needed 854 deac_vals[parId].append(values[:] / int(c.dtime)) 855 else: 856 # do deaccumulation 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) 830 893 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] 856 898 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)): 862 902 codes_set_values(gid, values) 863 903 … … 872 912 873 913 if c.basetime: 874 t_enddate = datetime.strptime(c.end_date + 875 c.basetime, 914 t_enddate = datetime.strptime(c.end_date + c.basetime, 876 915 '%Y%m%d%H') 877 916 else: 878 917 t_enddate = t_date + timedelta(2*int(c.dtime)) 879 918 880 # squeeze out information of last two steps contained881 # in svdp882 # if step+int(c.dtime) == c.maxstep and c.maxstep>12883 # or t_dt+timedelta(hours = int(c.dtime))884 # >= t_enddate:885 # Note that svdp[0] has not been popped in this case886 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) 915 954 916 955 codes_release(gid) … … 923 962 924 963 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 925 975 926 976 return
Note: See TracChangeset
for help on using the changeset viewer.