diff --git a/source/python/classes/EcFlexpart.py b/source/python/classes/EcFlexpart.py index 6576c892b69525b8fb9a5c587272316c2c2fc5a2..707ea80210618f58a30b779a93f2a0482b23bb59 100644 --- a/source/python/classes/EcFlexpart.py +++ b/source/python/classes/EcFlexpart.py @@ -81,14 +81,10 @@ import subprocess from datetime import datetime, timedelta import numpy as np -from gribapi import (grib_set, grib_index_select, grib_new_from_index, grib_get, - grib_write, grib_get_values, grib_set_values, grib_release, - grib_index_release, grib_index_get) - -# from eccodes import (codes_index_select, codes_new_from_index, codes_get, - # codes_get_values, codes_set_values, codes_set, - # codes_write, codes_release, codes_new_from_index, - # codes_index_release, codes_index_get) +from eccodes import (codes_index_select, codes_new_from_index, codes_get, + codes_get_values, codes_set_values, codes_set, + codes_write, codes_release, codes_new_from_index, + codes_index_release, codes_index_get) # software specific classes and modules from flex_extract sys.path.append('../') @@ -402,7 +398,7 @@ class EcFlexpart(object): Return ------ - iid : :obj:`grib_index` + iid : :obj:`codes_index` This is a grib specific index structure to access messages in a file. @@ -426,9 +422,7 @@ class EcFlexpart(object): # read the values of index keys index_vals = [] for key in index_keys: - #index_vals.append(grib_index_get(iid, key)) - #print(index_vals[-1]) - key_vals = grib_index_get(iid, key) + key_vals = codes_index_get(iid, key) print(key_vals) # have to sort the steps for disaggregation, # therefore convert to int first @@ -770,10 +764,10 @@ class EcFlexpart(object): print('current product: ', prod) for i in range(len(index_keys)): - grib_index_select(iid, index_keys[i], prod[i]) + codes_index_select(iid, index_keys[i], prod[i]) # get first id from current product - gid = grib_new_from_index(iid) + gid = codes_new_from_index(iid) # if there is no data for this specific time combination / product # skip the rest of the for loop and start with next timestep/product @@ -781,9 +775,9 @@ class EcFlexpart(object): continue # create correct timestamp from the three time informations - cdate = str(grib_get(gid, 'date')) - ctime = '{:0>2}'.format(grib_get(gid, 'time')/100) - cstep = '{:0>3}'.format(grib_get(gid, 'step')) + cdate = str(codes_get(gid, 'date')) + ctime = '{:0>2}'.format(codes_get(gid, 'time')/100) + cstep = '{:0>3}'.format(codes_get(gid, 'step')) t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H') t_dt = t_date + timedelta(hours=int(cstep)) t_m1dt = t_date + timedelta(hours=int(cstep)-int(c.dtime)) @@ -809,19 +803,22 @@ class EcFlexpart(object): t_dt.strftime('%Y%m%d%H')) print("outputfile = " + fnout) + f_handle = open(fnout, 'w') + h_handle = open(hnout, 'w') + g_handle = open(gnout, 'w') # read message for message and store relevant data fields # data keywords are stored in pars while 1: if not gid: break - cparamId = str(grib_get(gid, 'paramId')) - step = grib_get(gid, 'step') - time = grib_get(gid, 'time') - ni = grib_get(gid, 'Ni') - nj = grib_get(gid, 'Nj') + cparamId = str(codes_get(gid, 'paramId')) + step = codes_get(gid, 'step') + time = codes_get(gid, 'time') + ni = codes_get(gid, 'Ni') + nj = codes_get(gid, 'Nj') if cparamId in valsdict.keys(): - values = grib_get_values(gid) + values = codes_get_values(gid) vdp = valsdict[cparamId] svdp = svalsdict[cparamId] @@ -841,7 +838,7 @@ class EcFlexpart(object): print(cparamId, time, step, len(values), values[0], np.std(values)) - # len(svdp) correspond to the time + # len(svdp) corresponds to the time if len(svdp) >= 3: if len(svdp) > 3: if cparamId == '142' or cparamId == '143': @@ -859,17 +856,16 @@ class EcFlexpart(object): else: values = svdp[0] - grib_set_values(gid, values) + codes_set_values(gid, values) if c.maxstep > 12: - grib_set(gid, 'step', max(0, step-2*int(c.dtime))) + codes_set(gid, 'step', max(0, step-2*int(c.dtime))) else: - grib_set(gid, 'step', 0) - grib_set(gid, 'time', t_m2dt.hour*100) - grib_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d'))) + codes_set(gid, 'step', 0) + codes_set(gid, 'time', t_m2dt.hour*100) + codes_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d'))) - with open(fnout, 'w') as f_handle: - grib_write(gid, f_handle) + codes_write(gid, f_handle) if c.basetime: t_enddate = datetime.strptime(c.end_date + @@ -889,16 +885,15 @@ class EcFlexpart(object): t_dt == t_enddate: values = svdp[3] - grib_set_values(gid, values) - grib_set(gid, 'step', 0) + codes_set_values(gid, values) + codes_set(gid, 'step', 0) truedatetime = t_m2dt + timedelta(hours= 2*int(c.dtime)) - grib_set(gid, 'time', truedatetime.hour * 100) - grib_set(gid, 'date', truedatetime.year * 10000 + + codes_set(gid, 'time', truedatetime.hour * 100) + codes_set(gid, 'date', truedatetime.year * 10000 + truedatetime.month * 100 + truedatetime.day) - with open(hnout, 'w') as h_handle: - grib_write(gid, h_handle) + codes_write(gid, h_handle) #values = (svdp[1]+svdp[2])/2. if cparamId == '142' or cparamId == '143': @@ -906,21 +901,24 @@ class EcFlexpart(object): else: values = disaggregation.dapoly(list(reversed(svdp))) - grib_set(gid, 'step', 0) + codes_set(gid, 'step', 0) truedatetime = t_m2dt + timedelta(hours=int(c.dtime)) - grib_set(gid, 'time', truedatetime.hour * 100) - grib_set(gid, 'date', truedatetime.year * 10000 + - truedatetime.month * 100 + - truedatetime.day) - grib_set_values(gid, values) - with open(gnout, 'w') as g_handle: - grib_write(gid, g_handle) + codes_set(gid, 'time', truedatetime.hour * 100) + codes_set(gid, 'date', truedatetime.year * 10000 + + truedatetime.month * 100 + + truedatetime.day) + codes_set_values(gid, values) + codes_write(gid, g_handle) + + codes_release(gid) - grib_release(gid) + gid = codes_new_from_index(iid) - gid = grib_new_from_index(iid) + f_handle.close() + g_handle.close() + h_handle.close() - grib_index_release(iid) + codes_index_release(iid) return @@ -991,10 +989,10 @@ class EcFlexpart(object): print('current product: ', prod) for i in range(len(index_keys)): - grib_index_select(iid, index_keys[i], prod[i]) + codes_index_select(iid, index_keys[i], prod[i]) # get first id from current product - gid = grib_new_from_index(iid) + gid = codes_new_from_index(iid) # if there is no data for this specific time combination / product # skip the rest of the for loop and start with next timestep/product @@ -1009,9 +1007,9 @@ class EcFlexpart(object): fdict[k] = open(fortfile, 'w') # create correct timestamp from the three time informations - cdate = str(grib_get(gid, 'date')) - ctime = '{:0>2}'.format(grib_get(gid, 'time')/100) - cstep = '{:0>3}'.format(grib_get(gid, 'step')) + cdate = str(codes_get(gid, 'date')) + ctime = '{:0>2}'.format(codes_get(gid, 'time')/100) + cstep = '{:0>3}'.format(codes_get(gid, 'step')) timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H') timestamp += timedelta(hours=int(cstep)) cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H') @@ -1039,34 +1037,34 @@ class EcFlexpart(object): while 1: if not gid: break - paramId = grib_get(gid, 'paramId') - gridtype = grib_get(gid, 'gridType') - levtype = grib_get(gid, 'typeOfLevel') + paramId = codes_get(gid, 'paramId') + gridtype = codes_get(gid, 'gridType') + levtype = codes_get(gid, 'typeOfLevel') if paramId == 77: # ETADOT - grib_write(gid, fdict['21']) + codes_write(gid, fdict['21']) elif paramId == 130: # T - grib_write(gid, fdict['11']) + codes_write(gid, fdict['11']) elif paramId == 131 or paramId == 132: # U, V wind component - grib_write(gid, fdict['10']) + codes_write(gid, fdict['10']) elif paramId == 133 and gridtype != 'reduced_gg': # Q - grib_write(gid, fdict['17']) + codes_write(gid, fdict['17']) elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian - grib_write(gid, fdict['18']) + codes_write(gid, fdict['18']) elif paramId == 135: # W - grib_write(gid, fdict['19']) + codes_write(gid, fdict['19']) elif paramId == 152: # LNSP - grib_write(gid, fdict['12']) + codes_write(gid, fdict['12']) elif paramId == 155 and gridtype == 'sh': # D - grib_write(gid, fdict['13']) + codes_write(gid, fdict['13']) elif paramId == 246 or paramId == 247: # CLWC, CIWC # sum cloud liquid water and ice if scwc is None: - scwc = grib_get_values(gid) + scwc = codes_get_values(gid) else: - scwc += grib_get_values(gid) - grib_set_values(gid, scwc) - grib_set(gid, 'paramId', 201031) - grib_write(gid, fdict['22']) + scwc += codes_get_values(gid) + codes_set_values(gid, scwc) + codes_set(gid, 'paramId', 201031) + codes_write(gid, fdict['22']) elif c.wrf and paramId in [129, 138, 155] and \ levtype == 'hybrid': # Z, VO, D # do not do anything right now @@ -1076,7 +1074,7 @@ class EcFlexpart(object): if paramId not in savedfields: # SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR # and all ADDPAR parameter - grib_write(gid, fdict['16']) + codes_write(gid, fdict['16']) savedfields.append(paramId) else: print('duplicate ' + str(paramId) + ' not written') @@ -1086,15 +1084,15 @@ class EcFlexpart(object): # model layer if levtype == 'hybrid' and \ paramId in [129, 130, 131, 132, 133, 138, 155]: - grib_write(gid, fwrf) + codes_write(gid, fwrf) # sfc layer elif paramId in wrfpars: - grib_write(gid, fwrf) + codes_write(gid, fwrf) except AttributeError: pass - grib_release(gid) - gid = grib_new_from_index(iid) + codes_release(gid) + gid = codes_new_from_index(iid) for f in fdict.values(): f.close() @@ -1148,7 +1146,7 @@ class EcFlexpart(object): if c.wrf: fwrf.close() - grib_index_release(iid) + codes_index_release(iid) return @@ -1189,7 +1187,7 @@ class EcFlexpart(object): ofile = os.path.join(self.inputdir, ofile) if c.format.lower() == 'grib2': - p = subprocess.check_call(['grib_set', '-s', 'edition=2,', + p = subprocess.check_call(['codes_set', '-s', 'edition=2,', 'productDefinitionTemplateNumber=8', ofile, ofile + '_2']) p = subprocess.check_call(['mv', ofile + '_2', ofile]) diff --git a/source/python/classes/GribTools.py b/source/python/classes/GribTools.py index b2157032e2cf91fa73cf80c232b7a57bdc33edef..b85625259045fc0aa03401a19e0b664a4c7b65e9 100644 --- a/source/python/classes/GribTools.py +++ b/source/python/classes/GribTools.py @@ -43,16 +43,11 @@ # MODULES # ------------------------------------------------------------------------------ import os -from gribapi import grib_new_from_file, grib_is_defined, grib_get, \ - grib_release, grib_set, grib_write, grib_index_read, \ - grib_index_new_from_file, grib_index_add_file, \ - grib_index_write - -# from eccodes import (codes_grib_new_from_file, codes_is_defined, codes_get, - # codes_release, codes_set, codes_write, codes_index_read, - # codes_index_new_from_file, codes_index_add_file, - # codes_index_write) +from eccodes import (codes_grib_new_from_file, codes_is_defined, codes_get, + codes_release, codes_set, codes_write, codes_index_read, + codes_index_new_from_file, codes_index_add_file, + codes_index_write) # ------------------------------------------------------------------------------ # CLASS @@ -108,7 +103,7 @@ class GribTools(object): return_list = [] while 1: - gid = grib_new_from_file(fileid) + gid = codes_new_from_file(fileid) if gid is None: break @@ -120,20 +115,20 @@ class GribTools(object): select = True i = 0 for wherekey in wherekeynames: - if not grib_is_defined(gid, wherekey): + if not codes_is_defined(gid, wherekey): raise Exception("where key was not defined") select = (select and (str(wherekeyvalues[i]) == - str(grib_get(gid, wherekey)))) + str(codes_get(gid, wherekey)))) i += 1 if select: llist = [] for key in keynames: - llist.extend([str(grib_get(gid, key))]) + llist.extend([str(codes_get(gid, key))]) return_list.append(llist) - grib_release(gid) + codes_release(gid) fileid.close() @@ -181,7 +176,7 @@ class GribTools(object): fin = open(fromfile) while 1: - gid = grib_new_from_file(fin) + gid = codes_new_from_file(fin) if gid is None: break @@ -192,22 +187,22 @@ class GribTools(object): select = True i = 0 for wherekey in wherekeynames: - if not grib_is_defined(gid, wherekey): + if not codes_is_defined(gid, wherekey): raise Exception("where Key was not defined") select = (select and (str(wherekeyvalues[i]) == - str(grib_get(gid, wherekey)))) + str(codes_get(gid, wherekey)))) i += 1 if select: i = 0 for key in keynames: - grib_set(gid, key, keyvalues[i]) + codes_set(gid, key, keyvalues[i]) i += 1 - grib_write(gid, fout) + codes_write(gid, fout) - grib_release(gid) + codes_release(gid) fin.close() fout.close() @@ -249,7 +244,7 @@ class GribTools(object): fout = open(self.filenames, filemode) while 1: - gid = grib_new_from_file(fin) + gid = codes_new_from_file(fin) if gid is None: break @@ -260,21 +255,21 @@ class GribTools(object): select = True i = 0 for key in keynames: - if not grib_is_defined(gid, key): + if not codes_is_defined(gid, key): raise Exception("Key was not defined") if selectWhere: select = (select and (str(keyvalues[i]) == - str(grib_get(gid, key)))) + str(codes_get(gid, key)))) else: select = (select and (str(keyvalues[i]) != - str(grib_get(gid, key)))) + str(codes_get(gid, key)))) i += 1 if select: - grib_write(gid, fout) + codes_write(gid, fout) - grib_release(gid) + codes_release(gid) fin.close() fout.close() @@ -305,18 +300,18 @@ class GribTools(object): iid = None if os.path.exists(index_file): - iid = grib_index_read(index_file) + iid = codes_index_read(index_file) print("Use existing index file: %s " % (index_file)) else: for filename in self.filenames: print("Inputfile: %s " % (filename)) if iid is None: - iid = grib_index_new_from_file(filename, index_keys) + iid = codes_index_new_from_file(filename, index_keys) else: - grib_index_add_file(iid, filename) + codes_index_add_file(iid, filename) if iid is not None: - grib_index_write(iid, index_file) + codes_index_write(iid, index_file) print('... index done')