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 | # |
---|
11 | import subprocess |
---|
12 | import shutil |
---|
13 | import os, errno |
---|
14 | import datetime |
---|
15 | from dateutil.relativedelta import relativedelta |
---|
16 | import re |
---|
17 | |
---|
18 | from string import atoi |
---|
19 | from numpy import * |
---|
20 | from ecmwfapi import ECMWFDataServer |
---|
21 | from gribapi import * |
---|
22 | from GribTools import GribTools |
---|
23 | |
---|
24 | def 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 | ############################################################### |
---|
38 | def 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 | ############################################################### |
---|
48 | def 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 | ############################################################## |
---|
58 | def 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 | ############################################################## |
---|
69 | def 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 | ############################################################## |
---|
82 | def 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 | ############################################################## |
---|
96 | class 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 | ############################################################## |
---|
208 | class 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 | |
---|