source: flexpart.git/preproc/python/FlexpartTools.py @ 0e29ef4

10.4.1_peseiFPv9.3.1FPv9.3.1b_testingFPv9.3.2GFS_025bugfixes+enhancementsdevfp9.3.1-20161214-nc4grib2nc4_repairrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 0e29ef4 was 0e29ef4, checked in by Anne Fouilloux <annefou@…>, 9 years ago

bug fix when looping over months/dates

  • Property mode set to 100644
File size: 17.0 KB
Line 
1#
2# (C) Copyright 2014 UIO.
3#
4# This software is licensed under the terms of the Apache Licence Version 2.0
5# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6#
7#
8# Creation: October 2014 - Anne Fouilloux - University of Oslo
9#
10#
11import subprocess
12import shutil
13import os, errno
14import datetime
15from dateutil.relativedelta import relativedelta
16import re
17
18from string import atoi
19from numpy  import *
20from ecmwfapi import ECMWFDataServer
21from gribapi import *
22from GribTools import GribTools
23
24def product(*args, **kwds):
25    # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
26    # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
27    pools = map(tuple, args) * kwds.get('repeat', 1)
28    result = [[]]
29    for pool in pools:
30        result = [x+[y] for x in result for y in pool]
31    for prod in result:
32        yield tuple(prod)
33
34###############################################################
35# utility to remove a file if it exists
36# it does not fail if the file does not exist
37###############################################################
38def silentremove(filename):
39    try:
40        os.remove(filename)
41    except OSError as e: # this would be "except OSError, e:" before Python 2.6
42        if e.errno != errno.ENOENT: # errno.ENOENT = no such file or directory
43            raise # re-raise exception if a different error occured
44
45###############################################################
46# Utility to create directories
47###############################################################
48def mkdir_p(file):
49    path = os.path.dirname(file)
50    try:
51        os.stat(path)
52    except:
53        os.makedirs(path)
54
55##############################################################
56# function to iterate over dates
57##############################################################
58def daterange( start_date, end_date ):
59    if start_date <= end_date:
60        for n in range( ( end_date - start_date ).days + 1 ):
61            yield start_date + datetime.timedelta( n )
62    else:
63        for n in range( ( start_date - end_date ).days + 1 ):
64            yield start_date - datetime.timedelta( n )
65
66##############################################################
67# function to iterate over years
68##############################################################
69def years_between(start_date,end_date):
70    years = []
71    cursor = start_date
72
73    while cursor <= end_date:
74        if cursor.year not in years:
75            years.append(cursor.year)
76        cursor += relativedelta(years=1)
77
78    return years
79##############################################################
80# function to iterate over months
81##############################################################
82def months_between(start_date,end_date):
83    months = []
84    cursor = start_date
85
86    while cursor <= end_date:
87        if cursor.month not in months:
88            months.append(cursor.month)
89        cursor += datetime.timedelta(days=1)
90
91    return months
92
93##############################################################
94# MARSretrieval class
95##############################################################
96class MARSretrieval:
97    'class for MARS retrievals'
98
99    def __init__(self,server, dataset="interim",marsclass="ei",type="",levtype="",levelist="", 
100                 repres="", date="",resol="",stream="",area="",time="",step="",expver="1",
101                 number="",accuracy="", grid="", gaussian="", target="",param=""):
102        self.dataset=dataset
103        self.marsclass=marsclass
104        self.type=type
105        self.levtype=levtype
106        self.levelist=levelist
107        self.repres=repres
108        self.date=date
109        self.resol=resol
110        self.stream=stream
111        self.area=area
112        self.time=time
113        self.step=step
114        self.expver=expver
115        self.target=target
116        self.param=param
117        self.number=number
118        self.accuracy=accuracy
119        self.grid=grid
120        self.gaussian=gaussian
121        self.server=server
122
123    def displayInfo(self):
124        if self.dataset:
125            print "dataset: ", self.dataset
126        if self.marsclass:
127            print "class: ", self.marsclass
128        if self.type:
129             print "type: ", self.type
130        if self.levtype:
131            print "levtype: ", self.levtype
132        if self.levelist:
133            print "levelist: ", self.levelist
134        if self.repres:
135            print "repres: ", self.repres
136        if self.date:
137            print "date: ", self.date
138        if self.resol:
139            print "resol: ", self.resol
140        if self.stream:
141            print "stream: ", self.stream
142        if self.area:
143            print "area: ", self.area
144        if self.time:
145            print "time: ", self.time
146        if self.step:
147            print "step: ", self.step
148        if self.expver:
149            print "expver: ", self.expver
150        if self.target:
151            print "target: ", self.target
152        if self.param:
153            print "param: ", self.param
154        if self.number:
155            print "number: ", self.number
156        if self.accuracy:
157            print "accuracy: ", self.accuracy
158        if self.grid:
159            print "grid: ", self.grid
160        if self.gaussian:
161            print "gaussian: ", self.gaussian
162
163    def dataRetrieve(self):
164        dicolist = {}
165        if self.dataset:
166                dicolist['dataset']="%s"%(self.dataset)
167        if self.marsclass:
168                dicolist['class']="%s"%(self.marsclass)
169        if self.date:
170                dicolist['date']="%s"%(self.date)
171        if self.time:
172                dicolist['time']="%s"%(self.time)
173        if self.expver:
174                dicolist['expver']="%s"%(self.expver)
175        if self.param:
176                dicolist['param']="%s"%(self.param)
177        if self.type:
178                dicolist['type']="%s"%(self.type)
179        if self.levtype:
180                dicolist['levtype']="%s"%(self.levtype)
181        if self.levelist:
182                dicolist['levelist']="%s"%(self.levelist)
183        if self.repres:
184                dicolist['repres']="%s"%(self.repres)
185        if self.resol:
186                dicolist['resol']="%s"%(self.resol)
187        if self.stream:
188                dicolist['stream']="%s"%(self.stream)
189        if self.area:
190                dicolist['area']="%s"%(self.area)
191        if self.step:
192                dicolist['step']="%s"%(self.step)
193        if self.target:
194                dicolist['target']="%s"%(self.target)
195        if self.number:
196                dicolist['number']="%s"%(self.number)
197        if self.accuracy:
198                dicolist['accuracy']="%s"%(self.accuracy)
199        if self.grid:
200                dicolist['grid']="%s"%(self.grid)
201        if self.gaussian:
202                dicolist['gaussian']="%s"%(self.gaussian)
203
204        self.server.retrieve(dicolist)
205
206
207##############################################################
208class EIFlexpart:
209    'class to retrieve Era Interim data for running FLEXPART'
210##############################################################
211
212    def __init__(self):
213# different mars types for retrieving reanalysis data for flexpart
214       self.types=["an","fc"]
215
216       self.mars={}
217# set type (an/fc) and times and steps for EI
218#         (real) time    [type     time  step ] as in mars  (step for forecast only)     
219       self.mars["00"] = ["an",    "00"]
220       self.mars["01"] = ["fc",    "00", "01"]
221       self.mars["02"] = ["fc",    "00", "02"]
222       self.mars["03"] = ["fc",    "00", "03"]
223       self.mars["04"] = ["fc",    "00", "04"]
224       self.mars["05"] = ["fc",    "00", "05"]
225       self.mars["06"] = ["an",    "06"]
226       self.mars["07"] = ["fc",    "00", "07"]
227       self.mars["08"] = ["fc",    "00", "08"]
228       self.mars["09"] = ["fc",    "00", "09"]
229       self.mars["10"] = ["fc",    "00", "10"]
230       self.mars["11"] = ["fc",    "00", "11"]
231       self.mars["12"] = ["an",    "12"]
232       self.mars["13"] = ["fc",    "12", "01"]
233       self.mars["14"] = ["fc",    "12", "02"]
234       self.mars["15"] = ["fc",    "12", "03"]
235       self.mars["16"] = ["fc",    "12", "04"]
236       self.mars["17"] = ["fc",    "12", "05"]
237       self.mars["18"] = ["an",    "18"]
238       self.mars["19"] = ["fc",    "12", "07"]
239       self.mars["20"] = ["fc",    "12", "08"]
240       self.mars["21"] = ["fc",    "12", "09"]
241       self.mars["22"] = ["fc",    "12", "10"]
242       self.mars["23"] = ["fc",    "12", "11"]
243       
244    def retrieve(self, server, dates,times, area, levels, outputdir):
245        self.dates=dates
246        self.times=times
247        self.server=server
248        self.area=area
249        self.levels=levels
250        if outputdir=="":
251          self.outputdir='.'
252        else:
253          self.outputdir=outputdir
254        # surface data
255        dataset="interim"
256        marsclass="ei"
257        expver="0001"
258        levels="1/to/"+ self.levels
259# Retrieve Q not for using Q but as a template for a reduced gaussian grid one date and time is enough
260# Take analysis at 00
261        qdate=self.dates
262        idx=qdate.find("/")
263        if (idx >0):
264          qdate=self.dates[:idx]
265       
266        Q= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type="an", levtype="ML", levelist="1",
267                         gaussian="reduced",grid="80", resol="159",accuracy="24",target=self.outputdir+"/"+"Q.grb",
268                         date=qdate, time="00",expver=expver, param="133.128")
269        Q.displayInfo()
270        Q.dataRetrieve()
271
272        for type in self.types:
273            stime=set()
274            sstep=set()
275            for time in self.times.split('/'):
276# not so nice... to review later!
277                if self.mars[time][0] == type:
278                  stime.add(self.mars[time][1])
279                  if (len(self.mars[time]) > 2):
280                     sstep.add(self.mars[time][2])
281         
282            rtime="/".join(stime)
283            if len(sstep) > 0:
284              rstep="/".join(sstep)
285            else:
286              rstep=""
287            print type
288            print rtime
289            if rstep!="":
290               print rstep
291            print "MARS retrieve start... "
292            paramUVD="131.128/132.128/155.128"
293            paramT="130.128"
294            paramLNSP="152.128"
295            param2D2T="167.128/168.128"
296            if rstep!="":
297                UVD= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="ML", levelist=levels,
298                         resol="159", accuracy="24",grid="OFF",target=self.outputdir+"/"+type+"UVD.grb",
299                         area=self.area, date=self.dates, time=rtime,step=rstep, expver=expver, param=paramUVD)
300                T= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="ML", levelist=levels,
301                         resol="159", accuracy="24",grid="1.0/1.0",target=self.outputdir+"/"+type+"T.grb",
302                         area=self.area, date=self.dates, time=rtime,step=rstep, expver=expver, param=paramT)
303                LNSP= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="ML", levelist="1",
304                         accuracy="24",target=self.outputdir+"/"+type + "LNSP.grb",
305                         area=self.area, date=self.dates, time=rtime,step=rstep, expver=expver, param=paramLNSP)
306                DT= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="SFC", levelist="OFF",
307                         resol="159",accuracy="24",grid="1.0/1.0",target=self.outputdir+"/"+type + "2D2T.grb",
308                         area=self.area, date=self.dates, time=rtime,step=rstep, expver=expver, param=param2D2T)
309                print "do not retrieve Q for forecast"
310            else:
311                UVD= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="ML", levelist=levels,
312                         resol="159", accuracy="24",grid="OFF",target=self.outputdir+"/"+type + "UVD.grb",
313                         area=self.area, date=self.dates, time=rtime,expver=expver, param=paramUVD)
314                T= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="ML", levelist=levels,
315                         resol="159", accuracy="24",grid="1.0/1.0",target=self.outputdir+"/"+type + "T.grb",
316                         area=self.area, date=self.dates, time=rtime,expver=expver, param=paramT)
317                LNSP= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="ML", levelist="1",
318                         accuracy="24",target=self.outputdir+"/"+type + "LNSP.grb",
319                         area=self.area, date=self.dates, time=rtime,expver=expver, param=paramLNSP)
320                DT= MARSretrieval(self.server, dataset=dataset, marsclass=marsclass, stream="oper", type=type, levtype="SFC", levelist="OFF",
321                         resol="159",accuracy="24",grid="1.0/1.0",target=self.outputdir+"/"+type + "2D2T.grb",
322                         area=self.area, date=self.dates, time=rtime,expver=expver, param=param2D2T)
323
324            UVD.displayInfo()
325            UVD.dataRetrieve()
326            T.dataRetrieve()
327            LNSP.dataRetrieve()
328            DT.dataRetrieve()
329        print "MARS retrieve done... "
330       
331    def getFlexpartTime(self, type,step, time):
332        if int(step) < 10:
333            cstep = '0' +  str(step)
334        else:
335            cstep = str(step)
336        if int(time/100) < 10:
337            ctime = '0' +  str(int(time/100))
338        else:
339            ctime = str(int(time/100))
340
341        ctype = str(type)
342        if ctype == 'fc':
343            myinfo = [ctype,ctime, cstep]
344        else:
345            myinfo = [ctype,ctime]
346        cflextime = None
347        for t, marsinfo in self.mars.items():
348            if myinfo == marsinfo:
349                cflextime=t
350        return cflextime
351
352    def create(self, inputfiles, outputdir):
353        index_keys=["date","time","stepRange"]
354        indexfile="date_time_stepRange.idx"
355        silentremove(indexfile)
356        grib=GribTools(inputfiles.files)
357        iid=grib.index(index_keys=index_keys, index_file = indexfile)
358
359        print 'index done...' 
360        silentremove("fort.10")
361        silentremove("fort.11")
362        silentremove("fort.12")
363        silentremove("fort.13")
364        silentremove("fort.18")
365        index_vals = []
366        for key in index_keys:
367            key_vals = grib_index_get(iid,key)
368
369            index_vals.append(key_vals)
370
371
372        for prod in product(*index_vals):
373            for i in range(len(index_keys)):
374                grib_index_select(iid,index_keys[i],prod[i])
375
376            f10 = open('fort.10','w')
377            f11 = open('fort.11','w')
378            f12 = open('fort.12','w')
379            f13 = open('fort.13','w')
380            f16 = open('fort.16','w')
381            gid = grib_new_from_index(iid)
382            hid = gid
383            cflextime = None
384            if gid is not None: 
385                cdate = str(grib_get(gid, 'date'))
386                cyear = cdate[:4]
387                cmonth = cdate[4:6]
388                cday = cdate[6:8]
389                time = grib_get(gid, 'time')
390                type = grib_get(gid, 'type')
391                step = grib_get(gid, 'stepRange')
392                cflextime = self.getFlexpartTime(type,step, time)
393#                print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
394            while 1: 
395                if gid is None: break
396                paramId = grib_get(gid, 'paramId')
397                if paramId == 133:
398# Relative humidity (Q.grb) is used as a template only so we need the first we "meet"
399                    fout=open('fort.18','w')
400                    grib_write(gid,fout)
401                    fout.close()
402                if paramId == 131 or paramId == 132:
403                    grib_write(gid, f10)
404                if paramId == 130:
405                    grib_write(gid,f11)
406                if paramId == 152:
407                    grib_write(gid,f12)
408                if paramId == 155:
409                    grib_write(gid,f13)
410                if paramId == 167 or paramId == 168:
411                    grib_write(gid, f16)                   
412                grib_release(gid)
413                gid = grib_new_from_index(iid)
414# call for CONVERT2
415           
416            f10.close()
417            f11.close()
418            f12.close()
419            f13.close()
420            f16.close()
421            if hid is not None:
422                p=subprocess.check_call(['CONVERT2'])
423# create the corresponding output file  fort.15 (generated by CONVERT2) + fort.16 (paramId 167 and paramId 168) 
424                mkdir_p(outputdir+'/'+cyear+'/'+cmonth+'/')
425                print "outputdir = " + outputdir+'/'+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime
426                fout = open(outputdir+'/'+cyear+'/'+cmonth+'/EI'+cyear[2:4]+cmonth+cday+cflextime,'wb')   
427                shutil.copyfileobj(open('fort.15','rb'), fout)     
428                shutil.copyfileobj(open('fort.16','rb'), fout)
429                fout.close()
430
431        grib_index_release(iid)
432
433    def __del__(self):
434        print "clean"
435        silentremove("fort.10")
436        silentremove("fort.11")
437        silentremove("fort.12")
438        silentremove("fort.13")
439        silentremove("fort.15")
440        silentremove("fort.16")
441        silentremove("fort.18")
442        silentremove("VERTICAL.EC")
443        silentremove("date_time_stepRange.idx")
444
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG