[5763793] | 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" |
---|
[0e29ef4] | 355 | silentremove(indexfile) |
---|
[5763793] | 356 | grib=GribTools(inputfiles.files) |
---|
| 357 | iid=grib.index(index_keys=index_keys, index_file = indexfile) |
---|
| 358 | |
---|
[0e29ef4] | 359 | print 'index done...' |
---|
[5763793] | 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) |
---|
[0e29ef4] | 393 | # print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime |
---|
[5763793] | 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+'/') |
---|
[0e29ef4] | 425 | print "outputdir = " + outputdir+'/'+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime |
---|
[5763793] | 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 | |
---|
[0e29ef4] | 433 | def __del__(self): |
---|
[5763793] | 434 | print "clean" |
---|
[0e29ef4] | 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") |
---|
[5763793] | 444 | |
---|