Changeset ff99eae in flex_extract.git
- Timestamp:
- Jun 1, 2018, 8:34:59 PM (6 years ago)
- Branches:
- master, ctbto, dev
- Children:
- e1228f3
- Parents:
- ccab809
- Files:
-
- 1 added
- 7 edited
- 8 moved
Legend:
- Unmodified
- Added
- Removed
-
python/ControlFile.py
r812283d rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # T ODOAP4 # ToDo AP 5 5 # - write a test class 6 6 #************************************************************************ … … 37 37 # - __init__ 38 38 # - __str__ 39 # - tolist 39 # - to_list 40 # 41 # @Class Attributes: 42 # - start_date 43 # - end_date 44 # - accuracy 45 # - omega 46 # - cwc 47 # - omegadiff 48 # - etadiff 49 # - level 50 # - levelist 51 # - step 52 # - maxstep 53 # - prefix 54 # - makefile 55 # - basetime 56 # - date_chunk 57 # - grib2flexpart 58 # - exedir 59 # - flexpart_root_scripts 60 # - ecmwfdatadir 40 61 # 41 62 #******************************************************************************* … … 46 67 import os 47 68 import inspect 69 48 70 # software specific module from flex_extract 49 import Tools 71 from tools import get_list_as_string, my_error 50 72 51 73 # ------------------------------------------------------------------------------ 52 74 # CLASS 53 75 # ------------------------------------------------------------------------------ 54 class ControlFile :76 class ControlFile(object): 55 77 ''' 56 78 Class containing the information of the flex_extract CONTROL file. … … 114 136 for d in dd: 115 137 data.append(d) 116 pass117 138 if len(data) == 2: 118 139 if '$' in data[1]: … … 126 147 data[1] = data[1][:i] + var + data[1][k+1:] 127 148 else: 128 Tools.myerror(None, 129 'Could not find variable ' + 130 data[1][j+1:k] + 131 ' while reading ' + 132 filename) 149 my_error(None, 'Could not find variable ' + 150 data[1][j+1:k] + ' while reading ' + 151 filename) 133 152 setattr(self, data[0].lower() + '_expanded', data[1]) 134 153 else: … … 160 179 if not hasattr(self, 'levelist'): 161 180 if not hasattr(self, 'level'): 162 print ('Warning: neither levelist nor level \163 specified in CONTROL file' )181 print 'Warning: neither levelist nor level \ 182 specified in CONTROL file' 164 183 else: 165 184 self.levelist = '1/to/' + self.level … … 224 243 return ', '.join("%s: %s" % item for item in attrs.items()) 225 244 226 def to list(self):245 def to_list(self): 227 246 ''' 228 247 @Description: … … 255 274 pass 256 275 else: 257 if type(item[1]) is list:276 if isinstance(item[1], list): 258 277 stot = '' 259 278 for s in item[1]: … … 265 284 266 285 return sorted(l) 286 287 # def to_dict(self): 288 # ''' 289 290 # ''' 291 # parameters_dict = vars(self) 292 293 # # remove unneeded parameter 294 # parameters_dict.pop('_expanded', None) 295 # parameters_dict.pop('exedir', None) 296 # parameters_dict.pop('flexpart_root_scripts', None) 297 # parameters_dict.pop('ecmwfdatadir', None) 298 299 # parameters_dict_str = {} 300 # for key, value in parameters_dict.iteritems(): 301 # if isinstance(value, list): 302 # parameters_dict_str[str(key)] = get_list_as_string(value, ' ') 303 # else: 304 # parameters_dict_str[str(key)] = str(value) 305 306 # return parameters_dict_str -
python/EcFlexpart.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # T ODOAP4 # ToDo AP 5 5 # - specifiy file header documentation 6 6 # - add class description in header information … … 21 21 # - extended with class Control 22 22 # - removed functions mkdir_p, daterange, years_between, months_between 23 # - added functions darain, dapoly, to paramId, init128, normalexit,24 # my error, cleanup, install_args_and_control,23 # - added functions darain, dapoly, to_param_id, init128, normal_exit, 24 # my_error, clean_up, install_args_and_control, 25 25 # interpret_args_and_control, 26 26 # - removed function __del__ in class EIFLexpart … … 40 40 # - applied PEP8 style guide 41 41 # - added documentation 42 # - removed function getFlexpartTime in class E CFlexpart42 # - removed function getFlexpartTime in class EcFlexpart 43 43 # - outsourced class ControlFile 44 44 # - outsourced class MarsRetrieval 45 # - changed class name from EIFlexpart to E CFlexpart45 # - changed class name from EIFlexpart to EcFlexpart 46 46 # - applied minor code changes (style) 47 47 # - removed "dead code" , e.g. retrieval of Q since it is not needed … … 70 70 # - deacc_fluxes 71 71 # 72 # @Class Attributes: 73 # - dtime 74 # - basetime 75 # - server 76 # - marsclass 77 # - stream 78 # - resol 79 # - accuracy 80 # - number 81 # - expver 82 # - glevelist 83 # - area 84 # - grid 85 # - level 86 # - levelist 87 # - types 88 # - dates 89 # - area 90 # - gaussian 91 # - params 92 # - inputdir 93 # - outputfilelist 94 # 72 95 #******************************************************************************* 73 96 #pylint: disable=unsupported-assignment-operation 97 # this is disabled because its an error in pylint for this specific case 98 #pylint: disable=consider-using-enumerate 99 # this is not useful in this case 74 100 # ------------------------------------------------------------------------------ 75 101 # MODULES … … 78 104 import shutil 79 105 import os 80 import sys81 import inspect82 106 import glob 83 import datetime 84 from numpy import * 85 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter 86 ecapi = True 87 try: 88 import ecmwfapi 89 except ImportError: 90 ecapi = False 91 from gribapi import * 107 from datetime import datetime, timedelta 108 import numpy as np 109 from gribapi import grib_set, grib_index_select, grib_new_from_index, grib_get,\ 110 grib_write, grib_get_values, grib_set_values, grib_release,\ 111 grib_index_release, grib_index_get 92 112 93 113 # software specific classes and modules from flex_extract 94 from GribTools import GribTools 95 from Tools import init128, toparamId, silentremove, product 96 from ControlFile import ControlFile 97 from MARSretrieval import MARSretrieval 98 import Disagg 114 from Gribtools import Gribtools 115 from tools import init128, to_param_id, silent_remove, product, my_error 116 from MarsRetrieval import MarsRetrieval 117 import disaggregation 99 118 100 119 # ------------------------------------------------------------------------------ 101 120 # CLASS 102 121 # ------------------------------------------------------------------------------ 103 class E CFlexpart:122 class EcFlexpart(object): 104 123 ''' 105 124 Class to retrieve FLEXPART specific ECMWF data. … … 111 130 ''' 112 131 @Description: 113 Creates an object/instance of E CFlexpart with the132 Creates an object/instance of EcFlexpart with the 114 133 associated settings of its attributes for the retrieval. 115 134 116 135 @Input: 117 self: instance of E CFlexpart136 self: instance of EcFlexpart 118 137 The current object of the class. 119 138 … … 142 161 # different mars types for retrieving data for flexpart 143 162 self.types = dict() 144 try: 145 if c.maxstep > len(c.type): # Pure forecast mode 146 c.type = [c.type[1]] 147 c.step = ['{:0>3}'.format(int(c.step[0]))] 148 c.time = [c.time[0]] 149 for i in range(1, c.maxstep + 1): 150 c.type.append(c.type[0]) 151 c.step.append('{:0>3}'.format(i)) 152 c.time.append(c.time[0]) 153 except: 154 pass 163 164 if c.maxstep > len(c.type): # Pure forecast mode 165 c.type = [c.type[1]] 166 c.step = ['{:0>3}'.format(int(c.step[0]))] 167 c.time = [c.time[0]] 168 for i in range(1, c.maxstep + 1): 169 c.type.append(c.type[0]) 170 c.step.append('{:0>3}'.format(i)) 171 c.time.append(c.time[0]) 155 172 156 173 self.inputdir = c.inputdir … … 175 192 btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0] 176 193 177 if mod(i, int(c.dtime)) == 0 and \ 178 (c.maxstep > 24 or i in btlist): 194 if i % int(c.dtime) == 0 and c.maxstep > 24 or i in btlist: 179 195 180 196 if ty not in self.types.keys(): … … 182 198 183 199 if ti not in self.types[ty]['times']: 184 if len(self.types[ty]['times']) > 0:200 if self.types[ty]['times']: 185 201 self.types[ty]['times'] += '/' 186 202 self.types[ty]['times'] += ti 187 203 188 204 if st not in self.types[ty]['steps']: 189 if len(self.types[ty]['steps']) > 0:205 if self.types[ty]['steps']: 190 206 self.types[ty]['steps'] += '/' 191 207 self.types[ty]['steps'] += st 192 208 i += 1 193 209 194 # Different grids need different retrievals 195 # SH = Spherical Harmonics, GG = Gaussian Grid, 196 # OG = Output Grid, ML = MultiLevel, SL = SingleLevel 197 self.params = {'SH__ML': '', 'SH__SL': '', 198 'GG__ML': '', 'GG__SL': '', 199 'OG__ML': '', 'OG__SL': '', 200 'OG_OROLSM_SL': '', 'OG_acc_SL': ''} 210 201 211 self.marsclass = c.marsclass 202 212 self.stream = c.stream … … 205 215 self.accuracy = c.accuracy 206 216 self.level = c.level 207 try: 217 218 if c.levelist: 208 219 self.levelist = c.levelist 209 e xcept:220 else: 210 221 self.levelist = '1/to/' + c.level 211 222 … … 213 224 self.glevelist = '1/to/' + c.level 214 225 215 try:226 if c.gaussian: 216 227 self.gaussian = c.gaussian 217 e xcept:228 else: 218 229 self.gaussian = '' 219 230 220 try:231 if c.expver: 221 232 self.expver = c.expver 222 e xcept:233 else: 223 234 self.expver = '1' 224 235 225 try:236 if c.number: 226 237 self.number = c.number 227 e xcept:238 else: 228 239 self.number = '0' 229 240 … … 250 261 # 4) Download also data for WRF 251 262 263 264 # Different grids need different retrievals 265 # SH = Spherical Harmonics, GG = Gaussian Grid, 266 # OG = Output Grid, ML = MultiLevel, SL = SingleLevel 267 self.params = {'SH__ML': '', 'SH__SL': '', 268 'GG__ML': '', 'GG__SL': '', 269 'OG__ML': '', 'OG__SL': '', 270 'OG_OROLSM_SL': '', 'OG_acc_SL': ''} 271 252 272 if fluxes is False: 253 273 self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF'] … … 255 275 self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \ 256 276 'SFC', '1', self.grid] 257 if len(c.addpar) > 0:277 if c.addpar: 258 278 if c.addpar[0] == '/': 259 279 c.addpar = c.addpar[1:] … … 277 297 self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF'] 278 298 else: 279 print ('Warning: This is a very costly parameter combination, \280 use only for debugging!' )299 print 'Warning: This is a very costly parameter combination, \ 300 use only for debugging!' 281 301 self.params['GG__SL'] = ['Q', 'ML', '1', \ 282 302 '{}'.format((int(self.resol) + 1) / 2)] … … 287 307 self.params['OG__ML'][0] += '/W' 288 308 289 try: 290 # add cloud water content if necessary 291 if c.cwc == '1': 292 self.params['OG__ML'][0] += '/CLWC/CIWC' 293 except: 294 pass 295 296 try: 297 # add vorticity and geopotential height for WRF if necessary 298 if c.wrf == '1': 299 self.params['OG__ML'][0] += '/Z/VO' 300 if '/D' not in self.params['OG__ML'][0]: 301 self.params['OG__ML'][0] += '/D' 302 #wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ / 303 # stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper() 304 wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \ 305 139/170/183/236/39/40/41/42'.upper() 306 lwrt_sfc = wrf_sfc.split('/') 307 for par in lwrt_sfc: 308 if par not in self.params['OG__SL'][0]: 309 self.params['OG__SL'][0] += '/' + par 310 except: 311 pass 309 # add cloud water content if necessary 310 if c.cwc == '1': 311 self.params['OG__ML'][0] += '/CLWC/CIWC' 312 313 # add vorticity and geopotential height for WRF if necessary 314 if c.wrf == '1': 315 self.params['OG__ML'][0] += '/Z/VO' 316 if '/D' not in self.params['OG__ML'][0]: 317 self.params['OG__ML'][0] += '/D' 318 #wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ / 319 # stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper() 320 wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \ 321 139/170/183/236/39/40/41/42'.upper() 322 lwrt_sfc = wrf_sfc.split('/') 323 for par in lwrt_sfc: 324 if par not in self.params['OG__SL'][0]: 325 self.params['OG__SL'][0] += '/' + par 326 312 327 else: 313 328 self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \ … … 328 343 329 344 @Input: 330 self: instance of E CFlexpart345 self: instance of EcFlexpart 331 346 The current object of the class. 332 347 … … 352 367 353 368 self.inputdir = c.inputdir 354 area = asarray(self.area.split('/')).astype(float)355 grid = asarray(self.grid.split('/')).astype(float)369 area = np.asarray(self.area.split('/')).astype(float) 370 grid = np.asarray(self.grid.split('/')).astype(float) 356 371 357 372 if area[1] > area[3]: 358 373 area[1] -= 360 359 zyk = abs((area[3] - area[1] - 360.) + grid[1]) < 1.e-6360 374 maxl = int((area[3] - area[1]) / grid[1]) + 1 361 375 maxb = int((area[0] - area[2]) / grid[0]) + 1 … … 364 378 f.write('&NAMGEN\n') 365 379 f.write(',\n '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb), 366 'mlevel = ' + self.level, 367 'mlevelist = ' + '"' + self.levelist + '"', 368 'mnauf = ' + self.resol, 'metapar = ' + '77', 369 'rlo0 = ' + str(area[1]), 'rlo1 = ' + str(area[3]), 370 'rla0 = ' + str(area[2]), 'rla1 = ' + str(area[0]), 371 'momega = ' + c.omega, 'momegadiff = ' + c.omegadiff, 372 'mgauss = ' + c.gauss, 'msmooth = ' + c.smooth, 373 'meta = ' + c.eta, 'metadiff = ' + c.etadiff, 374 'mdpdeta = ' + c.dpdeta])) 380 'mlevel = ' + self.level, 381 'mlevelist = ' + '"' + self.levelist + '"', 382 'mnauf = ' + self.resol, 383 'metapar = ' + '77', 384 'rlo0 = ' + str(area[1]), 385 'rlo1 = ' + str(area[3]), 386 'rla0 = ' + str(area[2]), 387 'rla1 = ' + str(area[0]), 388 'momega = ' + c.omega, 389 'momegadiff = ' + c.omegadiff, 390 'mgauss = ' + c.gauss, 391 'msmooth = ' + c.smooth, 392 'meta = ' + c.eta, 393 'metadiff = ' + c.etadiff, 394 'mdpdeta = ' + c.dpdeta])) 375 395 376 396 f.write('\n/\n') … … 386 406 387 407 @Input: 388 self: instance of E CFlexpart408 self: instance of EcFlexpart 389 409 The current object of the class. 390 410 … … 449 469 # ------ on demand path -------------------------------------------------- 450 470 if self.basetime is None: 451 MR = MARSretrieval(self.server, 452 marsclass=self.marsclass, stream=mfstream, 453 type=mftype, levtype=pv[1], levelist=pv[2], 454 resol=self.resol, gaussian=gaussian, 455 accuracy=self.accuracy, grid=pv[3], 456 target=mftarget, area=area, date=mfdate, 457 time=mftime, number=self.number, step=mfstep, 458 expver=self.expver, param=pv[0]) 459 460 MR.displayInfo() 461 MR.dataRetrieve() 471 MR = MarsRetrieval(self.server, 472 marsclass=self.marsclass, 473 stream=mfstream, 474 type=mftype, 475 levtype=pv[1], 476 levelist=pv[2], 477 resol=self.resol, 478 gaussian=gaussian, 479 accuracy=self.accuracy, 480 grid=pv[3], 481 target=mftarget, 482 area=area, 483 date=mfdate, 484 time=mftime, 485 number=self.number, 486 step=mfstep, 487 expver=self.expver, 488 param=pv[0]) 489 490 MR.display_info() 491 MR.data_retrieve() 462 492 # ------ operational path ------------------------------------------------ 463 493 else: … … 475 505 tm1 = -1 476 506 477 maxtime = datetime.datetime.strptime( 478 mfdate.split('/')[-1] + mftime.split('/')[tm1], 479 '%Y%m%d%H') + datetime.timedelta( 480 hours=int(mfstep.split('/')[sm1])) 481 elimit = datetime.datetime.strptime( 482 mfdate.split('/')[-1] + 483 self.basetime, '%Y%m%d%H') 507 maxdate = datetime.strptime(mfdate.split('/')[-1] + 508 mftime.split('/')[tm1], 509 '%Y%m%d%H') 510 istep = int(mfstep.split('/')[sm1]) 511 maxtime = maxdate + timedelta(hours=istep) 512 513 elimit = datetime.strptime(mfdate.split('/')[-1] + 514 self.basetime, '%Y%m%d%H') 484 515 485 516 if self.basetime == '12': … … 492 523 # if maxtime-elimit<12h reduce step for last time 493 524 # A split of the MARS job into 2 is likely necessary. 494 maxtime = elimit - datetime.timedelta(hours=24) 495 mfdate = '/'.join(('/'.join(mfdate.split('/')[:-1]), 496 datetime.datetime.strftime( 497 maxtime, '%Y%m%d'))) 498 499 MR = MARSretrieval(self.server, 500 marsclass=self.marsclass, 501 stream=self.stream, type=mftype, 502 levtype=pv[1], levelist=pv[2], 503 resol=self.resol, gaussian=gaussian, 504 accuracy=self.accuracy, grid=pv[3], 505 target=mftarget, area=area, 506 date=mfdate, time=mftime, 507 number=self.number, step=mfstep, 508 expver=self.expver, param=pv[0]) 509 510 MR.displayInfo() 511 MR.dataRetrieve() 512 513 maxtime = elimit - datetime.timedelta(hours=12) 514 mfdate = datetime.datetime.strftime(maxtime, 515 '%Y%m%d') 525 maxtime = elimit - timedelta(hours=24) 526 mfdate = '/'.join(['/'.join(mfdate.split('/')[:-1]), 527 datetime.strftime(maxtime, 528 '%Y%m%d')]) 529 530 MR = MarsRetrieval(self.server, 531 marsclass=self.marsclass, 532 stream=self.stream, 533 type=mftype, 534 levtype=pv[1], 535 levelist=pv[2], 536 resol=self.resol, 537 gaussian=gaussian, 538 accuracy=self.accuracy, 539 grid=pv[3], 540 target=mftarget, 541 area=area, 542 date=mfdate, 543 time=mftime, 544 number=self.number, 545 step=mfstep, 546 expver=self.expver, 547 param=pv[0]) 548 549 MR.display_info() 550 MR.data_retrieve() 551 552 maxtime = elimit - timedelta(hours=12) 553 mfdate = datetime.strftime(maxtime, '%Y%m%d') 516 554 mftime = '00' 517 555 mftarget = self.inputdir + "/" + ftype + pk + \ … … 519 557 '.' + str(os.getpid()) + ".grb" 520 558 521 MR = MARSretrieval(self.server, 522 marsclass=self.marsclass, 523 stream=self.stream, type=mftype, 524 levtype=pv[1], levelist=pv[2], 525 resol=self.resol, gaussian=gaussian, 526 accuracy=self.accuracy, grid=pv[3], 527 target=mftarget, area=area, 528 date=mfdate, time=mftime, 529 number=self.number, step=mfstep, 530 expver=self.expver, param=pv[0]) 531 532 MR.displayInfo() 533 MR.dataRetrieve() 559 MR = MarsRetrieval(self.server, 560 marsclass=self.marsclass, 561 stream=self.stream, 562 type=mftype, 563 levtype=pv[1], 564 levelist=pv[2], 565 resol=self.resol, 566 gaussian=gaussian, 567 accuracy=self.accuracy, 568 grid=pv[3], 569 target=mftarget, 570 area=area, 571 date=mfdate, 572 time=mftime, 573 number=self.number, 574 step=mfstep, 575 expver=self.expver, 576 param=pv[0]) 577 578 MR.display_info() 579 MR.data_retrieve() 534 580 # -------------- non flux data ------------------------ 535 581 else: 536 MR = MARSretrieval(self.server, 537 marsclass=self.marsclass, 538 stream=self.stream, type=mftype, 539 levtype=pv[1], levelist=pv[2], 540 resol=self.resol, gaussian=gaussian, 541 accuracy=self.accuracy, grid=pv[3], 542 target=mftarget, area=area, 543 date=mfdate, time=mftime, 544 number=self.number, step=mfstep, 545 expver=self.expver, param=pv[0]) 546 547 MR.displayInfo() 548 MR.dataRetrieve() 582 MR = MarsRetrieval(self.server, 583 marsclass=self.marsclass, 584 stream=self.stream, 585 type=mftype, 586 levtype=pv[1], 587 levelist=pv[2], 588 resol=self.resol, 589 gaussian=gaussian, 590 accuracy=self.accuracy, 591 grid=pv[3], 592 target=mftarget, 593 area=area, 594 date=mfdate, 595 time=mftime, 596 number=self.number, 597 step=mfstep, 598 expver=self.expver, 599 param=pv[0]) 600 601 MR.display_info() 602 MR.data_retrieve() 549 603 else: # basetime == 0 ??? #AP 550 604 551 maxtime = elimit - datetime.timedelta(hours=24)552 mfdate = datetime. datetime.strftime(maxtime,'%Y%m%d')605 maxtime = elimit - timedelta(hours=24) 606 mfdate = datetime.strftime(maxtime, '%Y%m%d') 553 607 mftimesave = ''.join(mftime) 554 608 … … 556 610 times = mftime.split('/') 557 611 while ((int(times[0]) + 558 int(mfstep.split('/')[0]) <= 12) and559 (pk != 'OG_OROLSM__SL') and 'acc' not in pk):612 int(mfstep.split('/')[0]) <= 12) and 613 (pk != 'OG_OROLSM__SL') and 'acc' not in pk): 560 614 times = times[1:] 561 615 if len(times) > 1: … … 564 618 mftime = times[0] 565 619 566 MR = MARSretrieval(self.server, 567 marsclass=self.marsclass, 568 stream=self.stream, type=mftype, 569 levtype=pv[1], levelist=pv[2], 570 resol=self.resol, gaussian=gaussian, 571 accuracy=self.accuracy, grid=pv[3], 572 target=mftarget, area=area, 573 date=mfdate, time=mftime, 574 number=self.number, step=mfstep, 575 expver=self.expver, param=pv[0]) 576 577 MR.displayInfo() 578 MR.dataRetrieve() 620 MR = MarsRetrieval(self.server, 621 marsclass=self.marsclass, 622 stream=self.stream, 623 type=mftype, 624 levtype=pv[1], 625 levelist=pv[2], 626 resol=self.resol, 627 gaussian=gaussian, 628 accuracy=self.accuracy, 629 grid=pv[3], 630 target=mftarget, 631 area=area, 632 date=mfdate, 633 time=mftime, 634 number=self.number, 635 step=mfstep, 636 expver=self.expver, 637 param=pv[0]) 638 639 MR.display_info() 640 MR.data_retrieve() 579 641 580 642 if (int(mftimesave.split('/')[0]) == 0 and 581 int(mfstep.split('/')[0]) == 0 and 582 pk != 'OG_OROLSM__SL'): 583 mfdate = datetime.datetime.strftime(elimit,'%Y%m%d') 643 int(mfstep.split('/')[0]) == 0 and 644 pk != 'OG_OROLSM__SL'): 645 646 mfdate = datetime.strftime(elimit, '%Y%m%d') 584 647 mftime = '00' 585 648 mfstep = '000' … … 588 651 '.' + str(os.getpid()) + ".grb" 589 652 590 MR = MARSretrieval(self.server, 591 marsclass=self.marsclass, 592 stream=self.stream, type=mftype, 593 levtype=pv[1], levelist=pv[2], 594 resol=self.resol, gaussian=gaussian, 595 accuracy=self.accuracy, grid=pv[3], 596 target=mftarget, area=area, 597 date=mfdate, time=mftime, 598 number=self.number, step=mfstep, 599 expver=self.expver, param=pv[0]) 600 601 MR.displayInfo() 602 MR.dataRetrieve() 603 604 print("MARS retrieve done... ") 653 MR = MarsRetrieval(self.server, 654 marsclass=self.marsclass, 655 stream=self.stream, 656 type=mftype, 657 levtype=pv[1], 658 levelist=pv[2], 659 resol=self.resol, 660 gaussian=gaussian, 661 accuracy=self.accuracy, 662 grid=pv[3], 663 target=mftarget, 664 area=area, 665 date=mfdate, 666 time=mftime, 667 number=self.number, 668 step=mfstep, 669 expver=self.expver, 670 param=pv[0]) 671 672 MR.display_info() 673 MR.data_retrieve() 674 675 print "MARS retrieve done... " 605 676 606 677 return … … 621 692 622 693 @Input: 623 self: instance of E CFlexpart694 self: instance of EcFlexpart 624 695 The current object of the class. 625 696 … … 642 713 ''' 643 714 644 print ('\n\nPostprocessing:\n Format: {}\n'.format(c.format))715 print '\n\nPostprocessing:\n Format: {}\n'.format(c.format) 645 716 646 717 if c.ecapi is False: … … 652 723 c.destination = os.getenv('DESTINATION') 653 724 print('ectrans: {}\n gateway: {}\n destination: {}\n ' 654 655 656 print ('Output filelist: \n')657 print (self.outputfilelist)725 .format(c.ectrans, c.gateway, c.destination)) 726 727 print 'Output filelist: \n' 728 print self.outputfilelist 658 729 659 730 if c.format.lower() == 'grib2': 660 731 for ofile in self.outputfilelist: 661 732 p = subprocess.check_call(['grib_set', '-s', 'edition=2, \ 662 663 733 productDefinitionTemplateNumber=8', 734 ofile, ofile + '_2']) 664 735 p = subprocess.check_call(['mv', ofile + '_2', ofile]) 665 736 … … 667 738 for ofile in self.outputfilelist: 668 739 p = subprocess.check_call(['ectrans', '-overwrite', '-gateway', 669 670 740 c.gateway, '-remote', c.destination, 741 '-source', ofile]) 671 742 print('ectrans:', p) 672 743 … … 674 745 for ofile in self.outputfilelist: 675 746 p = subprocess.check_call(['ecp', '-o', ofile, 676 747 os.path.expandvars(c.ecfsdir)]) 677 748 678 749 if c.outputdir != c.inputdir: … … 692 763 if '.' in fname[-1]: 693 764 l = fname[-1].split('.') 694 timestamp = datetime. datetime.strptime(l[0][-6:] + l[1],695 696 timestamp += datetime.timedelta(hours=int(l[2]))697 cdate = datetime. datetime.strftime(timestamp, '%Y%m%d')698 chms = datetime. datetime.strftime(timestamp, '%H%M%S')765 timestamp = datetime.strptime(l[0][-6:] + l[1], 766 '%y%m%d%H') 767 timestamp += timedelta(hours=int(l[2])) 768 cdate = datetime.strftime(timestamp, '%Y%m%d') 769 chms = datetime.strftime(timestamp, '%H%M%S') 699 770 else: 700 771 cdate = '20' + fname[-1][-8:-2] … … 708 779 # generate pathnames file 709 780 pwd = os.path.abspath(c.outputdir) 710 with open(pwd + '/pathnames', 'w') as f:781 with open(pwd + '/pathnames', 'w') as f: 711 782 f.write(pwd + '/Options/\n') 712 783 f.write(pwd + '/\n') … … 720 791 721 792 # read template COMMAND file 722 with open(os.path.expandvars( 723 os.path.expanduser(c.flexpart_root_scripts)) + 724 '/../Options/COMMAND', 'r') as f: 793 with open(os.path.expandvars(os.path.expanduser( 794 c.flexpart_root_scripts)) + '/../Options/COMMAND', 'r') as f: 725 795 lflist = f.read().split('\n') 726 796 … … 747 817 # afterwards switch back to the working dir 748 818 os.chdir(c.outputdir) 749 p = subprocess.check_call([os.path.expandvars( 750 os.path.expanduser(c.flexpart_root_scripts)) + 751 '/../FLEXPART_PROGRAM/grib2flexpart', 752 'useAvailable', '.']) 819 p = subprocess.check_call([ 820 os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts)) 821 + '/../FLEXPART_PROGRAM/grib2flexpart', 'useAvailable', '.']) 753 822 os.chdir(pwd) 754 823 … … 771 840 772 841 @Input: 773 self: instance of E CFlexpart842 self: instance of EcFlexpart 774 843 The current object of the class. 775 844 776 inputfiles: instance of U IOFiles845 inputfiles: instance of UioFiles 777 846 Contains a list of files. 778 847 … … 796 865 table128 = init128(c.ecmwfdatadir + 797 866 '/grib_templates/ecmwf_grib1_table_128') 798 wrfpars = to paramId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\867 wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\ 799 868 stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4', 800 869 table128) … … 802 871 index_keys = ["date", "time", "step"] 803 872 indexfile = c.inputdir + "/date_time_stepRange.idx" 804 silent remove(indexfile)805 grib = Grib Tools(inputfiles.files)873 silent_remove(indexfile) 874 grib = Gribtools(inputfiles.files) 806 875 # creates new index file 807 876 iid = grib.index(index_keys=index_keys, index_file=indexfile) … … 811 880 for key in index_keys: 812 881 index_vals.append(grib_index_get(iid, key)) 813 print (index_vals[-1])882 print index_vals[-1] 814 883 # index_vals looks for example like: 815 884 # index_vals[0]: ('20171106', '20171107', '20171108') ; date … … 842 911 # remove old fort.* files and open new ones 843 912 for k, f in fdict.iteritems(): 844 silent remove(c.inputdir + "/fort." + k)913 silent_remove(c.inputdir + "/fort." + k) 845 914 fdict[k] = open(c.inputdir + '/fort.' + k, 'w') 846 915 847 916 cdate = str(grib_get(gid, 'date')) 848 917 time = grib_get(gid, 'time') 849 type = grib_get(gid, 'type')850 918 step = grib_get(gid, 'step') 851 919 # create correct timestamp from the three time informations 852 920 # date, time, step 853 timestamp = datetime.datetime.strptime( 854 cdate + '{:0>2}'.format(time/100), '%Y%m%d%H') 855 timestamp += datetime.timedelta(hours=int(step)) 856 857 cdateH = datetime.datetime.strftime(timestamp, '%Y%m%d%H') 858 chms = datetime.datetime.strftime(timestamp, '%H%M%S') 921 timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100), 922 '%Y%m%d%H') 923 timestamp += timedelta(hours=int(step)) 924 925 cdateH = datetime.strftime(timestamp, '%Y%m%d%H') 859 926 860 927 if c.basetime is not None: 861 slimit = datetime.datetime.strptime( 862 c.start_date + '00', '%Y%m%d%H') 928 slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H') 863 929 bt = '23' 864 930 if c.basetime == '00': 865 931 bt = '00' 866 slimit = datetime.datetime.strptime( 867 c.end_date + bt, '%Y%m%d%H') - \ 868 datetime.timedelta(hours=12-int(c.dtime)) 932 slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\ 933 - timedelta(hours=12-int(c.dtime)) 869 934 if c.basetime == '12': 870 935 bt = '12' 871 slimit = datetime.datetime.strptime( 872 c.end_date + bt, '%Y%m%d%H') - \ 873 datetime.timedelta(hours=12-int(c.dtime)) 874 875 elimit = datetime.datetime.strptime( 876 c.end_date + bt, '%Y%m%d%H') 936 slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\ 937 - timedelta(hours=12-int(c.dtime)) 938 939 elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H') 877 940 878 941 if timestamp < slimit or timestamp > elimit: … … 901 964 paramId = grib_get(gid, 'paramId') 902 965 gridtype = grib_get(gid, 'gridType') 903 datatype = grib_get(gid, 'dataType')904 966 levtype = grib_get(gid, 'typeOfLevel') 905 967 if paramId == 133 and gridtype == 'reduced_gg': … … 939 1001 savedfields.append(paramId) 940 1002 else: 941 print ('duplicate ' + str(paramId) + ' not written')1003 print 'duplicate ' + str(paramId) + ' not written' 942 1004 943 1005 try: … … 963 1025 os.chdir(c.inputdir) 964 1026 if os.stat('fort.21').st_size == 0 and int(c.eta) == 1: 965 print ('Parameter 77 (etadot) is missing, most likely it is \966 not available for this type or date/time\n' )967 print ('Check parameters CLASS, TYPE, STREAM, START_DATE\n')968 my error(c, 'fort.21 is empty while parameter eta is set \1027 print 'Parameter 77 (etadot) is missing, most likely it is \ 1028 not available for this type or date/time\n' 1029 print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n' 1030 my_error(c, 'fort.21 is empty while parameter eta is set \ 969 1031 to 1 in CONTROL file') 970 1032 971 1033 # create the corresponding output file fort.15 972 1034 # (generated by CONVERT2) + fort.16 (paramId 167 and 168) 973 p = subprocess.check_call([os.path.expandvars( 974 os.path.expanduser(c.exedir)) + '/CONVERT2'], shell=True) 1035 p = subprocess.check_call( 1036 [os.path.expandvars(os.path.expanduser(c.exedir)) + 1037 '/CONVERT2'], shell=True) 975 1038 os.chdir(pwd) 976 1039 … … 983 1046 suffix = cdateH[2:10] 984 1047 fnout += suffix 985 print ("outputfile = " + fnout)1048 print "outputfile = " + fnout 986 1049 self.outputfilelist.append(fnout) # needed for final processing 987 1050 … … 1005 1068 open(c.inputdir + '/fort.25', 'rb'), fout) 1006 1069 1007 try: 1008 if c.wrf == '1': 1009 fwrf.close() 1010 except: 1011 pass 1070 if c.wrf == '1': 1071 fwrf.close() 1012 1072 1013 1073 grib_index_release(iid) … … 1025 1085 1026 1086 @Input: 1027 self: instance of E CFlexpart1087 self: instance of EcFlexpart 1028 1088 The current object of the class. 1029 1089 1030 inputfiles: instance of U IOFiles1090 inputfiles: instance of UioFiles 1031 1091 Contains a list of files. 1032 1092 … … 1050 1110 table128 = init128(c.ecmwfdatadir + 1051 1111 '/grib_templates/ecmwf_grib1_table_128') 1052 pars = to paramId(self.params['OG_acc_SL'][0], table128)1112 pars = to_param_id(self.params['OG_acc_SL'][0], table128) 1053 1113 index_keys = ["date", "time", "step"] 1054 1114 indexfile = c.inputdir + "/date_time_stepRange.idx" 1055 silent remove(indexfile)1056 grib = Grib Tools(inputfiles.files)1115 silent_remove(indexfile) 1116 grib = Gribtools(inputfiles.files) 1057 1117 # creates new index file 1058 1118 iid = grib.index(index_keys=index_keys, index_file=indexfile) … … 1061 1121 index_vals = [] 1062 1122 for key in index_keys: 1063 key_vals = grib_index_get(iid, key)1064 print (key_vals)1123 key_vals = grib_index_get(iid, key) 1124 print key_vals 1065 1125 # have to sort the steps for disaggregation, 1066 1126 # therefore convert to int first … … 1083 1143 stepsdict[str(p)] = [] 1084 1144 1145 print 'maxstep: ', c.maxstep 1146 1085 1147 for prod in product(*index_vals): 1086 1148 # e.g. prod = ('20170505', '0', '12') … … 1092 1154 1093 1155 gid = grib_new_from_index(iid) 1094 # do convert2 program if gid at this time is not None,1095 # therefore save in hid1096 hid = gid1097 1156 if gid is not None: 1098 1157 cdate = grib_get(gid, 'date') 1099 1158 time = grib_get(gid, 'time') 1100 type = grib_get(gid, 'type')1101 1159 step = grib_get(gid, 'step') 1102 1160 # date+time+step-2*dtime 1103 1161 # (since interpolated value valid for step-2*dtime) 1104 sdate = datetime .datetime(year=cdate/10000,1105 month=mod(cdate,10000)/100,1106 day=mod(cdate,100),1107 1108 fdate = sdate + datetime.timedelta(1109 hours=step-2*int(c.dtime))1110 sdates = sdate + datetime.timedelta(hours=step)1162 sdate = datetime(year=cdate/10000, 1163 month=(cdate % 10000)/100, 1164 day=(cdate % 100), 1165 hour=time/100) 1166 fdate = sdate + timedelta(hours=step-2*int(c.dtime)) 1167 sdates = sdate + timedelta(hours=step) 1168 elimit = None 1111 1169 else: 1112 1170 break … … 1126 1184 else: 1127 1185 fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H') 1128 gnout = c.inputdir + '/flux' + (fdate+datetime.timedelta( 1129 hours = int(c.dtime))).strftime('%Y%m%d%H') 1186 gnout = c.inputdir + '/flux' + (fdate + 1187 timedelta(hours=int(c.dtime)) 1188 ).strftime('%Y%m%d%H') 1130 1189 hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H') 1131 1190 g = open(gnout, 'w') 1132 1191 h = open(hnout, 'w') 1133 1192 1134 print ("outputfile = " + fnout)1193 print "outputfile = " + fnout 1135 1194 f = open(fnout, 'w') 1136 1195 … … 1156 1215 fak = 3600. 1157 1216 1158 values = ( reshape(values, (nj, ni))).flatten() / fak1217 values = (np.reshape(values, (nj, ni))).flatten() / fak 1159 1218 vdp.append(values[:]) # save the accumulated values 1160 1219 if step <= int(c.dtime): … … 1164 1223 1165 1224 print(cparamId, atime, step, len(values), 1166 values[0], std(values))1225 values[0], np.std(values)) 1167 1226 # save the 1/3-hourly or specific values 1168 1227 # svdp.append(values[:]) … … 1172 1231 if len(svdp) > 3: 1173 1232 if cparamId == '142' or cparamId == '143': 1174 values = Disagg.darain(svdp)1233 values = disaggregation.darain(svdp) 1175 1234 else: 1176 values = Disagg.dapoly(svdp)1235 values = disaggregation.dapoly(svdp) 1177 1236 1178 1237 if not (step == c.maxstep and c.maxstep > 12 \ … … 1197 1256 1198 1257 if c.basetime is not None: 1199 elimit = datetime.datetime.strptime(c.end_date + 1200 c.basetime, 1201 '%Y%m%d%H') 1258 elimit = datetime.strptime(c.end_date + 1259 c.basetime, '%Y%m%d%H') 1202 1260 else: 1203 elimit = sdate + datetime.timedelta(2*int(c.dtime))1261 elimit = sdate + timedelta(2*int(c.dtime)) 1204 1262 1205 1263 # squeeze out information of last two steps contained 1206 1264 # in svdp 1207 1265 # if step+int(c.dtime) == c.maxstep and c.maxstep>12 1208 # or sdates+ datetime.timedelta(hours = int(c.dtime))1266 # or sdates+timedelta(hours = int(c.dtime)) 1209 1267 # >= elimit: 1210 1268 # Note that svdp[0] has not been popped in this case 1211 1269 1212 if (step == c.maxstep and c.maxstep > 121213 or sdates == elimit):1270 if step == c.maxstep and c.maxstep > 12 or \ 1271 sdates == elimit: 1214 1272 1215 1273 values = svdp[3] 1216 1274 grib_set_values(gid, values) 1217 1275 grib_set(gid, 'step', 0) 1218 truedatetime = fdate + datetime.timedelta(1219 hours=2*int(c.dtime))1276 truedatetime = fdate + timedelta(hours= 1277 2*int(c.dtime)) 1220 1278 grib_set(gid, 'time', truedatetime.hour * 100) 1221 1279 grib_set(gid, 'date', truedatetime.year * 10000 + … … 1226 1284 #values = (svdp[1]+svdp[2])/2. 1227 1285 if cparamId == '142' or cparamId == '143': 1228 values = Disagg.darain(list(reversed(svdp)))1286 values = disaggregation.darain(list(reversed(svdp))) 1229 1287 else: 1230 values = Disagg.dapoly(list(reversed(svdp))) 1231 1232 grib_set(gid, 'step',0) 1233 truedatetime = fdate + datetime.timedelta( 1234 hours=int(c.dtime)) 1288 values = disaggregation.dapoly(list(reversed(svdp))) 1289 1290 grib_set(gid, 'step', 0) 1291 truedatetime = fdate + timedelta(hours=int(c.dtime)) 1235 1292 grib_set(gid, 'time', truedatetime.hour * 100) 1236 1293 grib_set(gid, 'date', truedatetime.year * 10000 + … … 1250 1307 grib_index_release(iid) 1251 1308 1309 exit() 1252 1310 return -
python/GribTools.py
rccab809 rff99eae 1 1 #!/usr/bin/env python 2 2 # -*- coding: utf-8 -*- 3 #************************************************************************4 # TODO AP5 # - GribTools name möglicherweise etwas verwirrend.6 # - change self.filename in self.filenames!!!7 # - bis auf --init-- und index wird keine Funktion verwendet!?8 #************************************************************************9 3 #******************************************************************************* 10 4 # @Author: Anne Fouilloux (University of Oslo) … … 36 30 # @Class Content: 37 31 # - __init__ 38 # - get keys39 # - set keys32 # - get_keys 33 # - set_keys 40 34 # - copy 41 35 # - index 36 # 37 # @Class Attributes: 38 # - filenames 42 39 # 43 40 #******************************************************************************* … … 47 44 # ------------------------------------------------------------------------------ 48 45 import os 49 from gribapi import * 46 from gribapi import grib_new_from_file, grib_is_defined, grib_get, \ 47 grib_release, grib_set, grib_write, grib_index_read, \ 48 grib_index_new_from_file, grib_index_add_file, \ 49 grib_index_write 50 50 51 51 # ------------------------------------------------------------------------------ 52 52 # CLASS 53 53 # ------------------------------------------------------------------------------ 54 class Grib Tools:54 class Gribtools(object): 55 55 ''' 56 56 Class for GRIB utilities (new methods) based on GRIB API … … 62 62 ''' 63 63 @Description: 64 Initialise an object of Grib Tools and assign a list64 Initialise an object of Gribtools and assign a list 65 65 of filenames. 66 66 … … 73 73 ''' 74 74 75 self.filename = filenames75 self.filenames = filenames 76 76 77 77 return 78 78 79 79 80 def get keys(self, keynames, wherekeynames=[], wherekeyvalues=[]):80 def get_keys(self, keynames, wherekeynames=[], wherekeyvalues=[]): 81 81 ''' 82 82 @Description: … … 88 88 List of keynames. 89 89 90 wherekeynames: list of ???, optional90 wherekeynames: list of strings, optional 91 91 Default value is an empty list. 92 92 93 wherekeyvalues: list of ???, optional93 wherekeyvalues: list of strings, optional 94 94 Default value is an empty list. 95 95 … … 99 99 ''' 100 100 101 fileid = open(self.filename , 'r')101 fileid = open(self.filenames, 'r') 102 102 103 103 return_list = [] … … 136 136 137 137 138 def set keys(self, fromfile, keynames, keyvalues, wherekeynames=[],139 wherekeyvalues=[], strict=False, filemode='w'):138 def set_keys(self, fromfile, keynames, keyvalues, wherekeynames=[], 139 wherekeyvalues=[], strict=False, filemode='w'): 140 140 ''' 141 141 @Description: … … 150 150 Filename of the input file to read the grib messages from. 151 151 152 keynames: list of ???152 keynames: list of strings 153 153 List of keynames. Default is an empty list. 154 154 155 keyvalues: list of ???155 keyvalues: list of strings 156 156 List of keynames. Default is an empty list. 157 157 158 wherekeynames: list of ???, optional158 wherekeynames: list of strings, optional 159 159 Default value is an empty list. 160 160 161 wherekeyvalues: list of ???, optional161 wherekeyvalues: list of strings, optional 162 162 Default value is an empty list. 163 163 … … 174 174 175 175 ''' 176 fout = open(self.filename , filemode)176 fout = open(self.filenames, filemode) 177 177 fin = open(fromfile) 178 178 … … 196 196 i += 1 197 197 198 #AP is it secured that the order of keynames is equal to keyvalues?199 198 if select: 200 199 i = 0 … … 203 202 i += 1 204 203 205 #AP this is a redundant code section 206 # delete the if/else : 207 # 208 # grib_write(gid_in, fout) 209 # 210 if strict: 211 if select: 212 grib_write(gid_in, fout) 213 else: 214 grib_write(gid_in, fout) 215 #AP end 204 grib_write(gid_in, fout) 205 216 206 grib_release(gid_in) 217 207 … … 238 228 function. Default is True. 239 229 240 keynames: list of ???, optional230 keynames: list of strings, optional 241 231 List of keynames. Default is an empty list. 242 232 243 keyvalues: list of ???, optional233 keyvalues: list of strings, optional 244 234 List of keynames. Default is an empty list. 245 235 … … 252 242 253 243 fin = open(filename_in) 254 fout = open(self.filename , filemode)244 fout = open(self.filenames, filemode) 255 245 256 246 while 1: … … 307 297 Grib index id. 308 298 ''' 309 print ("... index will be done")310 self.iid = None311 312 if (os.path.exists(index_file)):313 self.iid = grib_index_read(index_file)314 print ("Use existing index file: %s " % (index_file))299 print "... index will be done" 300 iid = None 301 302 if os.path.exists(index_file): 303 iid = grib_index_read(index_file) 304 print "Use existing index file: %s " % (index_file) 315 305 else: 316 for file in self.filename:317 print ("Inputfile: %s " % (file))318 if self.iid is None:319 self.iid = grib_index_new_from_file(file, index_keys)306 for filename in self.filenames: 307 print "Inputfile: %s " % (filename) 308 if iid is None: 309 iid = grib_index_new_from_file(filename, index_keys) 320 310 else: 321 311 print 'in else zweig' 322 grib_index_add_file(self.iid, file) 323 324 if self.iid is not None: 325 grib_index_write(self.iid, index_file) 326 327 print('... index done') 328 329 return self.iid 330 331 332 333 334 312 grib_index_add_file(iid, filename) 313 314 if iid is not None: 315 grib_index_write(iid, index_file) 316 317 print '... index done' 318 319 return iid -
python/MarsRetrieval.py
r991df6a rff99eae 1 1 #!/usr/bin/env python 2 2 # -*- coding: utf-8 -*- 3 #************************************************************************4 # TODO AP5 # -6 #************************************************************************7 3 #******************************************************************************* 8 4 # @Author: Anne Fouilloux (University of Oslo) … … 13 9 # 14 10 # November 2015 - Leopold Haimberger (University of Vienna): 15 # - optimized display Info16 # - optimized data Retrieve and seperate between python and shell11 # - optimized display_info 12 # - optimized data_retrieve and seperate between python and shell 17 13 # script call 18 14 # … … 37 33 # @Class Content: 38 34 # - __init__ 39 # - displayInfo 40 # - dataRetrieve 35 # - display_info 36 # - data_retrieve 37 # 38 # @Class Attributes: 39 # - server 40 # - marsclass 41 # - dtype 42 # - levtype 43 # - levelist 44 # - repres 45 # - date 46 # - resol 47 # - stream 48 # - area 49 # - time 50 # - step 51 # - expver 52 # - number 53 # - accuracy 54 # - grid 55 # - gaussian 56 # - target 57 # - param 41 58 # 42 59 #******************************************************************************* … … 48 65 import os 49 66 50 ecapi = True51 try:52 import ecmwfapi53 except ImportError:54 ecapi = False55 56 67 # ------------------------------------------------------------------------------ 57 68 # CLASS 58 69 # ------------------------------------------------------------------------------ 59 class M ARSretrieval:70 class MarsRetrieval(object): 60 71 ''' 61 72 Class for submitting MARS retrievals. … … 68 79 ''' 69 80 70 def __init__(self, server, marsclass = "ei", type = "", levtype ="",71 levelist = "", repres = "", date = "", resol = "", stream ="",72 area = "", time = "", step = "", expver = "1", number ="",73 accuracy = "", grid = "", gaussian = "", target ="",74 param =""):81 def __init__(self, server, marsclass="ei", dtype="", levtype="", 82 levelist="", repres="", date="", resol="", stream="", 83 area="", time="", step="", expver="1", number="", 84 accuracy="", grid="", gaussian="", target="", 85 param=""): 75 86 ''' 76 87 @Description: 77 Initialises the instance of the M ARSretrieval class and88 Initialises the instance of the MarsRetrieval class and 78 89 defines and assigns a set of the necessary retrieval parameters 79 90 for the FLEXPART input data. … … 84 95 85 96 @Input: 86 self: instance of M ARSretrieval97 self: instance of MarsRetrieval 87 98 For description see class documentation. 88 99 … … 96 107 Default is the ERA-Interim dataset "ei". 97 108 98 type: string, optional109 dtype: string, optional 99 110 Determines the type of fields to be retrieved. 100 111 Selects between observations, images or fields. … … 275 286 self.server = server 276 287 self.marsclass = marsclass 277 self. type =type288 self.dtype = dtype 278 289 self.levtype = levtype 279 290 self.levelist = levelist … … 296 307 297 308 298 def display Info(self):309 def display_info(self): 299 310 ''' 300 311 @Description: … … 302 313 303 314 @Input: 304 self: instance of M ARSretrieval315 self: instance of MarsRetrieval 305 316 For description see class documentation. 306 317 … … 314 325 # with their corresponding values 315 326 for item in attrs.items(): 316 if item[0] in ('server'):327 if item[0] in 'server': 317 328 pass 318 329 else: 319 print (item[0] + ': ' + str(item[1]))330 print item[0] + ': ' + str(item[1]) 320 331 321 332 return 322 333 323 def data Retrieve(self):334 def data_retrieve(self): 324 335 ''' 325 336 @Description: … … 330 341 331 342 @Input: 332 self: instance of M ARSretrieval343 self: instance of MarsRetrieval 333 344 For description see class documentation. 334 345 … … 344 355 s = 'ret' 345 356 for k, v in attrs.iteritems(): 346 if k in ('server'):357 if k in 'server': 347 358 continue 348 359 if k == 'marsclass': … … 360 371 self.server.execute(s, target) 361 372 except: 362 print ('MARS Request failed, \363 have you already registered at apps.ecmwf.int?' )373 print 'MARS Request failed, \ 374 have you already registered at apps.ecmwf.int?' 364 375 raise IOError 365 376 if os.stat(target).st_size == 0: 366 print ('MARS Request returned no data - please check request')377 print 'MARS Request returned no data - please check request' 367 378 raise IOError 368 379 # MARS request via extra process in shell … … 373 384 stderr=subprocess.PIPE, bufsize=1) 374 385 pout = p.communicate(input=s)[0] 375 print (pout.decode())386 print pout.decode() 376 387 377 388 if 'Some errors reported' in pout.decode(): 378 print ('MARS Request failed - please check request')389 print 'MARS Request failed - please check request' 379 390 raise IOError 380 391 381 392 if os.stat(target).st_size == 0: 382 print ('MARS Request returned no data - please check request')393 print 'MARS Request returned no data - please check request' 383 394 raise IOError 384 395 -
python/UioFiles.py
r991df6a rff99eae 1 1 #!/usr/bin/env python 2 2 # -*- coding: utf-8 -*- 3 #************************************************************************4 # TODO AP5 # - checken welche regelmässigen methoden auf diese Files noch angewendet werden6 # und dann hier implementieren7 # cleanup hier rein8 #************************************************************************9 3 #******************************************************************************* 10 4 # @Author: Anne Fouilloux (University of Oslo) … … 15 9 # 16 10 # November 2015 - Leopold Haimberger (University of Vienna): 17 # - modified method list Files to work with glob instead of listdir18 # - added pattern search in method list Files11 # - modified method list_files to work with glob instead of listdir 12 # - added pattern search in method list_files 19 13 # 20 14 # February 2018 - Anne Philipp (University of Vienna): 21 15 # - applied PEP8 style guide 22 16 # - added documentation 23 # - optimisation of method list Files since it didn't work correctly17 # - optimisation of method list_files since it didn't work correctly 24 18 # for sub directories 25 # - additional speed up of method list Files19 # - additional speed up of method list_files 26 20 # - modified the class so that it is initiated with a pattern instead 27 21 # of suffixes. Gives more precision in selection of files. … … 40 34 # @Class Content: 41 35 # - __init__ 42 # - listFiles 43 # - deleteFiles 36 # - list_files 37 # - delete_files 38 # 39 # @Class Attributes: 40 # - pattern 41 # - files 44 42 # 45 43 #******************************************************************************* … … 49 47 # ------------------------------------------------------------------------------ 50 48 import os 51 import glob52 49 import fnmatch 53 import time54 50 55 51 # software specific module from flex_extract 56 import profiling57 from Tools import silentremove52 #import profiling 53 from tools import silent_remove 58 54 59 55 # ------------------------------------------------------------------------------ … … 61 57 # ------------------------------------------------------------------------------ 62 58 63 class U IOFiles:59 class UioFiles(object): 64 60 ''' 65 61 Class to manipulate files. At initialisation it has the attribute … … 76 72 77 73 @Input: 78 self: instance of U IOFiles74 self: instance of UioFiles 79 75 Description see class documentation. 80 76 … … 87 83 88 84 self.pattern = pattern 85 self.files = None 89 86 90 87 return 91 88 92 89 #@profiling.timefn 93 def list Files(self, path, callid=0):90 def list_files(self, path, callid=0): 94 91 ''' 95 92 @Description: … … 98 95 99 96 @Input: 100 self: instance of U IOFiles97 self: instance of UioFiles 101 98 Description see class documentation. 102 99 … … 132 129 if subdirs: 133 130 for subdir in subdirs: 134 self.list Files(os.path.join(path, subdir), callid=1)131 self.list_files(os.path.join(path, subdir), callid=1) 135 132 136 133 return 137 134 138 def delete Files(self):135 def delete_files(self): 139 136 ''' 140 137 @Description: … … 142 139 143 140 @Input: 144 self: instance of U IOFiles141 self: instance of UioFiles 145 142 Description see class documentation. 146 143 … … 149 146 ''' 150 147 151 for fin self.files:152 silent remove(f)148 for old_file in self.files: 149 silent_remove(old_file) 153 150 154 151 return -
python/disaggregation.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # T ODOAP4 # ToDo AP 5 5 # - check alist of size 4 ? 6 6 # - write a test, IMPORTANT … … 20 20 # - added structured documentation 21 21 # - outsourced the disaggregation functions dapoly and darain 22 # to a new module named Disagg22 # to a new module named disaggregation 23 23 # 24 24 # @License: … … 29 29 # 30 30 # @Module Description: 31 # Disaggregation of deaccumulated flux data from an ECMWF model FG field.31 # disaggregationregation of deaccumulated flux data from an ECMWF model FG field. 32 32 # Initially the flux data to be concerned are: 33 33 # - large-scale precipitation … … 70 70 using a cubic polynomial solution which conserves the integrals 71 71 of the fluxes within each timespan. 72 Disaggregation is done for 4 accumluated timespans which generates72 disaggregationregation is done for 4 accumluated timespans which generates 73 73 a new, disaggregated value which is output at the central point 74 74 of the 4 accumulation timespans. This new point is used for linear … … 111 111 field using a modified linear solution which conserves the integrals 112 112 of the fluxes within each timespan. 113 Disaggregation is done for 4 accumluated timespans which generates113 disaggregationregation is done for 4 accumluated timespans which generates 114 114 a new, disaggregated value which is output at the central point 115 115 of the 4 accumulation timespans. This new point is used for linear … … 144 144 145 145 return nfield 146 147 148 149 150 151 152 153 154 155 -
python/get_mars_data.py
r991df6a rff99eae 1 1 #!/usr/bin/env python 2 2 # -*- coding: utf-8 -*- 3 #************************************************************************4 # TODO AP5 # - add function docstrings!!!!6 #************************************************************************7 3 #******************************************************************************* 8 4 # @Author: Anne Fouilloux (University of Oslo) … … 13 9 # 14 10 # November 2015 - Leopold Haimberger (University of Vienna): 15 # - moved the getEIdata program into a function "get MARSdata"11 # - moved the getEIdata program into a function "get_mars_data" 16 12 # - moved the AgurmentParser into a seperate function 17 13 # - adatpted the function for the use in flex_extract 18 # - renamed file to get MARSdata14 # - renamed file to get_mars_data 19 15 # 20 16 # February 2018 - Anne Philipp (University of Vienna): … … 42 38 # @Program Content: 43 39 # - main 44 # - get MARSdata40 # - get_mars_data 45 41 # 46 42 #******************************************************************************* … … 54 50 import inspect 55 51 try: 56 ecapi =True52 ecapi = True 57 53 import ecmwfapi 58 54 except ImportError: 59 ecapi=False 55 ecapi = False 56 57 # software specific classes and modules from flex_extract 58 from tools import my_error, normal_exit, interpret_args_and_control 59 from EcFlexpart import EcFlexpart 60 from UioFiles import UioFiles 60 61 61 62 # add path to pythonpath so that python finds its buddies 62 localpythonpath= os.path.dirname(os.path.abspath(63 LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath( 63 64 inspect.getfile(inspect.currentframe()))) 64 if localpythonpath not in sys.path: 65 sys.path.append(localpythonpath) 66 67 # software specific classes and modules from flex_extract 68 from ControlFile import ControlFile 69 from Tools import myerror, normalexit, \ 70 interpret_args_and_control 71 from ECFlexpart import ECFlexpart 72 from UIOFiles import UIOFiles 65 if LOCAL_PYTHON_PATH not in sys.path: 66 sys.path.append(LOCAL_PYTHON_PATH) 73 67 74 68 # ------------------------------------------------------------------------------ … … 78 72 ''' 79 73 @Description: 80 If get MARSdata is called from command line, this function controls74 If get_mars_data is called from command line, this function controls 81 75 the program flow and calls the argumentparser function and 82 the get MARSdata function for retrieving EC data.76 the get_mars_data function for retrieving EC data. 83 77 84 78 @Input: … … 89 83 ''' 90 84 args, c = interpret_args_and_control() 91 get MARSdata(args,c)92 normal exit(c)85 get_mars_data(c) 86 normal_exit(c) 93 87 94 88 return 95 89 96 def get MARSdata(args,c):90 def get_mars_data(c): 97 91 ''' 98 92 @Description: … … 103 97 104 98 @Input: 105 args: instance of ArgumentParser106 Contains the commandline arguments from script/program call.107 108 99 c: instance of class ControlFile 109 100 Contains all the parameters of CONTROL file, which are e.g.: … … 126 117 os.makedirs(c.inputdir) 127 118 128 print ("Retrieving EC data!")129 print ("start date %s " % (c.start_date))130 print ("end date %s " % (c.end_date))119 print "Retrieving EC data!" 120 print "start date %s " % (c.start_date) 121 print "end date %s " % (c.end_date) 131 122 132 123 if ecapi: … … 160 151 # -------------- flux data ------------------------------------------------ 161 152 print 'removing old flux content of ' + c.inputdir 162 tobecleaned = U IOFiles('*_acc_*.' + str(os.getppid()) + '.*.grb')163 tobecleaned.list Files(c.inputdir)164 tobecleaned.delete Files()153 tobecleaned = UioFiles('*_acc_*.' + str(os.getppid()) + '.*.grb') 154 tobecleaned.list_files(c.inputdir) 155 tobecleaned.delete_files() 165 156 166 157 # if forecast for maximum one day (upto 23h) are to be retrieved, … … 172 163 while day < endp1: 173 164 # retrieve MARS data for the whole period 174 flexpart = E CFlexpart(c, fluxes=True)165 flexpart = EcFlexpart(c, fluxes=True) 175 166 tmpday = day + datechunk - datetime.timedelta(days=1) 176 167 if tmpday < endp1: … … 186 177 flexpart.retrieve(server, dates, c.inputdir) 187 178 except IOError: 188 my error(c, 'MARS request failed')179 my_error(c, 'MARS request failed') 189 180 190 181 day += datechunk … … 199 190 while day <= end: 200 191 # retrieve MARS data for the whole period 201 flexpart = E CFlexpart(c, fluxes=True)192 flexpart = EcFlexpart(c, fluxes=True) 202 193 tmpday = day + datechunk - datetime.timedelta(days=1) 203 194 if tmpday < end: … … 213 204 flexpart.retrieve(server, dates, c.inputdir) 214 205 except IOError: 215 my error(c, 'MARS request failed')206 my_error(c, 'MARS request failed') 216 207 217 208 day += datechunk … … 219 210 # -------------- non flux data -------------------------------------------- 220 211 print 'removing old non flux content of ' + c.inputdir 221 tobecleaned = U IOFiles('*__*.' + str(os.getppid()) + '.*.grb')222 tobecleaned.list Files(c.inputdir)223 tobecleaned.delete Files()212 tobecleaned = UioFiles('*__*.' + str(os.getppid()) + '.*.grb') 213 tobecleaned.list_files(c.inputdir) 214 tobecleaned.delete_files() 224 215 225 216 day = start 226 217 while day <= end: 227 228 flexpart = ECFlexpart(c, fluxes=False)229 230 231 232 233 234 235 236 237 238 239 240 241 242 myerror(c, 'MARS request failed')243 244 218 # retrieve all non flux MARS data for the whole period 219 flexpart = EcFlexpart(c, fluxes=False) 220 tmpday = day + datechunk - datetime.timedelta(days=1) 221 if tmpday < end: 222 dates = day.strftime("%Y%m%d") + "/to/" + \ 223 tmpday.strftime("%Y%m%d") 224 else: 225 dates = day.strftime("%Y%m%d") + "/to/" + \ 226 end.strftime("%Y%m%d") 227 228 print "retrieve " + dates + " in dir " + c.inputdir 229 230 try: 231 flexpart.retrieve(server, dates, c.inputdir) 232 except IOError: 233 my_error(c, 'MARS request failed') 234 235 day += datechunk 245 236 246 237 return … … 248 239 if __name__ == "__main__": 249 240 main() 250 -
python/install.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # TODO AP 5 # - localpythonpath should not be set in module load section! 4 # ToDo AP 6 5 # - create a class Installation and divide installation in 3 subdefs for 7 6 # ecgate, local and cca seperatly … … 18 17 # - applied PEP8 style guide 19 18 # - added documentation 19 # - moved install_args_and_control in here 20 20 # 21 21 # @License: … … 45 45 # MODULES 46 46 # ------------------------------------------------------------------------------ 47 import datetime48 47 import os 49 48 import sys … … 51 50 import subprocess 52 51 import inspect 53 from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter 54 55 # add path to pythonpath so that python finds its buddies 56 localpythonpath = os.path.dirname(os.path.abspath( 57 inspect.getfile(inspect.currentframe()))) 58 if localpythonpath not in sys.path: 59 sys.path.append(localpythonpath) 52 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter 60 53 61 54 # software specific classes and modules from flex_extract 62 55 from ControlFile import ControlFile 56 57 # add path to pythonpath so that python finds its buddies 58 LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath( 59 inspect.getfile(inspect.currentframe()))) 60 if LOCAL_PYTHON_PATH not in sys.path: 61 sys.path.append(LOCAL_PYTHON_PATH) 63 62 64 63 # ------------------------------------------------------------------------------ … … 78 77 ''' 79 78 80 os.chdir( localpythonpath)79 os.chdir(LOCAL_PYTHON_PATH) 81 80 args, c = install_args_and_control() 82 81 … … 84 83 install_via_gateway(c, args.install_target) 85 84 else: 86 print ('Please specify installation target (local|ecgate|cca)')87 print ('use -h or --help for help')85 print 'Please specify installation target (local|ecgate|cca)' 86 print 'use -h or --help for help' 88 87 89 88 sys.exit() … … 152 151 try: 153 152 c = ControlFile(args.controlfile) 154 except :155 print ('Could not read CONTROL file "' + args.controlfile + '"')156 print ('Either it does not exist or its syntax is wrong.')157 print ('Try "' + sys.argv[0].split('/')[-1] +158 ' -h" to print usage information' )153 except IOError: 154 print 'Could not read CONTROL file "' + args.controlfile + '"' 155 print 'Either it does not exist or its syntax is wrong.' 156 print 'Try "' + sys.argv[0].split('/')[-1] + \ 157 ' -h" to print usage information' 159 158 exit(1) 160 159 161 160 if args.install_target != 'local': 162 if (args.ecgid is None or args.ecuid is None or args.gateway is None163 or args.destination is None):164 print ('Please enter your ECMWF user id and group id as well as \161 if args.ecgid is None or args.ecuid is None or args.gateway is None \ 162 or args.destination is None: 163 print 'Please enter your ECMWF user id and group id as well as \ 165 164 the \nname of the local gateway and the ectrans \ 166 destination ' )167 print ('with command line options --ecuid --ecgid \168 --gateway --destination' )169 print ('Try "' + sys.argv[0].split('/')[-1] +170 ' -h" to print usage information' )171 print ('Please consult ecaccess documentation or ECMWF user support \172 for further details' )165 destination ' 166 print 'with command line options --ecuid --ecgid \ 167 --gateway --destination' 168 print 'Try "' + sys.argv[0].split('/')[-1] + \ 169 ' -h" to print usage information' 170 print 'Please consult ecaccess documentation or ECMWF user support \ 171 for further details' 173 172 sys.exit(1) 174 173 else: … … 178 177 c.destination = args.destination 179 178 180 try:179 if args.makefile: 181 180 c.makefile = args.makefile 182 except:183 pass184 181 185 182 if args.install_target == 'local': … … 240 237 c.flexpart_root_scripts 241 238 else: 242 data ='export FLEXPART_ROOT_SCRIPTS=$HOME'239 data = 'export FLEXPART_ROOT_SCRIPTS=$HOME' 243 240 if target.lower() != 'local': 244 241 if '--workdir' in data: … … 292 289 fo.close() 293 290 294 295 296 291 if target.lower() == 'local': 297 292 # compile CONVERT2 298 293 if c.flexpart_root_scripts is None or c.flexpart_root_scripts == '../': 299 print ('Warning: FLEXPART_ROOT_SCRIPTS has not been specified')300 print ('Only CONVERT2 will be compiled in ' + ecd + '/../src')294 print 'Warning: FLEXPART_ROOT_SCRIPTS has not been specified' 295 print 'Only CONVERT2 will be compiled in ' + ecd + '/../src' 301 296 else: 302 297 c.flexpart_root_scripts = os.path.expandvars(os.path.expanduser( 303 298 c.flexpart_root_scripts)) 304 299 if os.path.abspath(ecd) != os.path.abspath(c.flexpart_root_scripts): 305 300 os.chdir('/') … … 311 306 try: 312 307 os.makedirs(c.flexpart_root_scripts + '/ECMWFDATA7.1') 313 except:308 finally: 314 309 pass 315 310 os.chdir(c.flexpart_root_scripts + '/ECMWFDATA7.1') … … 329 324 p = subprocess.check_call(['rm'] + flist) 330 325 try: 331 print (('Using makefile: ' + makefile))326 print 'Using makefile: ' + makefile 332 327 p = subprocess.check_call(['make', '-f', makefile]) 333 p = subprocess.check_call(['ls', '-l','CONVERT2']) 334 except: 335 print('compile failed - please edit ' + makefile + 336 ' or try another Makefile in the src directory.') 337 print('most likely GRIB_API_INCLUDE_DIR, GRIB_API_LIB ' 338 'and EMOSLIB must be adapted.') 339 print('Available Makefiles:') 340 print(glob.glob('Makefile*')) 341 328 p = subprocess.check_call(['ls', '-l', 'CONVERT2']) 329 except subprocess.CalledProcessError as e: 330 print 'compile failed with the following error:' 331 print e.output 332 print 'please edit ' + makefile + \ 333 ' or try another Makefile in the src directory.' 334 print 'most likely GRIB_API_INCLUDE_DIR, GRIB_API_LIB \ 335 and EMOSLIB must be adapted.' 336 print 'Available Makefiles:' 337 print glob.glob('Makefile*') 342 338 elif target.lower() == 'ecgate': 343 339 os.chdir('/') … … 352 348 'ecgate:/home/ms/' + c.ecgid + '/' + 353 349 c.ecuid + '/ECMWFDATA7.1.tar']) 354 except: 355 print('ecaccess-file-put failed! Probably the eccert key has expired.') 350 except subprocess.CalledProcessError as e: 351 print 'ecaccess-file-put failed! \ 352 Probably the eccert key has expired.' 356 353 exit(1) 357 p = subprocess.check_call(['ecaccess-job-submit', 358 '-queueName', 359 target, 360 ecd + 'python/compilejob.ksh']) 361 print('compilejob.ksh has been submitted to ecgate for ' 362 'installation in ' + c.ec_flexpart_root_scripts + 363 '/ECMWFDATA7.1') 364 print('You should get an email with subject flexcompile within ' 365 'the next few minutes') 354 355 try: 356 p = subprocess.check_call(['ecaccess-job-submit', 357 '-queueName', 358 target, 359 ecd + 'python/compilejob.ksh']) 360 print 'compilejob.ksh has been submitted to ecgate for \ 361 installation in ' + c.ec_flexpart_root_scripts + \ 362 '/ECMWFDATA7.1' 363 print 'You should get an email with subject flexcompile within \ 364 the next few minutes' 365 except subprocess.CalledProcessError as e: 366 print 'ecaccess-job-submit failed!' 367 exit(1) 366 368 367 369 elif target.lower() == 'cca': … … 377 379 'cca:/home/ms/' + c.ecgid + '/' + 378 380 c.ecuid + '/ECMWFDATA7.1.tar']) 379 except :380 print ('ecaccess-file-put failed! '381 'Probably the eccert key has expired.')381 except subprocess.CalledProcessError as e: 382 print 'ecaccess-file-put failed! \ 383 Probably the eccert key has expired.' 382 384 exit(1) 383 385 384 p=subprocess.check_call(['ecaccess-job-submit', 385 '-queueName', 386 target, 387 ecd + 'python/compilejob.ksh']) 388 print('compilejob.ksh has been submitted to cca for installation in ' + 389 c.ec_flexpart_root_scripts + '/ECMWFDATA7.1') 390 print('You should get an email with subject flexcompile ' 391 'within the next few minutes') 386 try: 387 p = subprocess.check_call(['ecaccess-job-submit', 388 '-queueName', 389 target, 390 ecd + 'python/compilejob.ksh']) 391 print 'compilejob.ksh has been submitted to cca for installation in ' +\ 392 c.ec_flexpart_root_scripts + '/ECMWFDATA7.1' 393 print 'You should get an email with subject flexcompile \ 394 within the next few minutes' 395 except subprocess.CalledProcessError as e: 396 print 'ecaccess-job-submit failed!' 397 exit(1) 392 398 393 399 else: 394 print ('ERROR: unknown installation target ', target)395 print ('Valid targets: ecgate, cca, local')400 print 'ERROR: unknown installation target ', target 401 print 'Valid targets: ecgate, cca, local' 396 402 397 403 return -
python/plot_retrieved.py
rccab809 rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # T ODOAP4 # ToDo AP 5 5 # - documentation der Funktionen 6 6 # - docu der progam functionality … … 19 19 # - created function main and moved the two function calls for 20 20 # arguments and plotting into it 21 # - added function get Basics to extract the boundary conditions21 # - added function get_basics to extract the boundary conditions 22 22 # of the data fields from the first grib file it gets. 23 23 # … … 33 33 # @Program Content: 34 34 # - main 35 # - get Basics36 # - plot Retrieved37 # - plot TS38 # - plot Map39 # - get PlotArgs35 # - get_basics 36 # - plot_retrieved 37 # - plot_timeseries 38 # - plot_map 39 # - get_plot_args 40 40 # 41 41 #******************************************************************************* … … 49 49 import inspect 50 50 import sys 51 import glob52 51 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter 53 54 from matplotlib.pylab import * 55 import matplotlib.patches as mpatches 56 from mpl_toolkits.basemap import Basemap, addcyclic 57 import matplotlib.colors as mcolors 58 from matplotlib.font_manager import FontProperties 59 from matplotlib.patches import Polygon 60 import matplotlib.cm as cmx 61 import matplotlib.colors as colors 62 #from rasotools.utils import stats 52 import matplotlib 53 import matplotlib.pyplot as plt 54 from mpl_toolkits.basemap import Basemap 55 from eccodes import codes_grib_new_from_file, codes_get, codes_release, \ 56 codes_get_values 57 import numpy as np 58 59 # software specific classes and modules from flex_extract 60 from ControlFile import ControlFile 61 from UioFiles import UioFiles 62 63 # add path to pythonpath so that python finds its buddies 64 LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath( 65 inspect.getfile(inspect.currentframe()))) 66 if LOCAL_PYTHON_PATH not in sys.path: 67 sys.path.append(LOCAL_PYTHON_PATH) 63 68 64 69 font = {'family': 'monospace', 'size': 12} 65 70 matplotlib.rcParams['xtick.major.pad'] = '20' 66 67 71 matplotlib.rc('font', **font) 68 69 from eccodes import *70 import numpy as np71 72 # add path to pythonpath so that python finds its buddies73 localpythonpath = os.path.dirname(os.path.abspath(74 inspect.getfile(inspect.currentframe())))75 if localpythonpath not in sys.path:76 sys.path.append(localpythonpath)77 78 # software specific classes and modules from flex_extract79 from Tools import silentremove80 from ControlFile import ControlFile81 #from GribTools import GribTools82 from UIOFiles import UIOFiles83 84 72 # ------------------------------------------------------------------------------ 85 73 # FUNCTION … … 88 76 ''' 89 77 @Description: 90 If plot Retrieved is called from command line, this function controls78 If plot_retrieved is called from command line, this function controls 91 79 the program flow and calls the argumentparser function and 92 the plot Retrieved function for plotting the retrieved GRIB data.80 the plot_retrieved function for plotting the retrieved GRIB data. 93 81 94 82 @Input: … … 98 86 <nothing> 99 87 ''' 100 args, c = get PlotArgs()101 plot Retrieved(args,c)88 args, c = get_plot_args() 89 plot_retrieved(c) 102 90 103 91 return 104 92 105 def get Basics(ifile, verb=False):93 def get_basics(ifile, verb=False): 106 94 """ 107 95 @Description: … … 113 101 ifile: string 114 102 Contains the full absolute path to the ECMWF grib file. 103 115 104 verb (opt): bool 116 105 Is True if there should be extra output in verbose mode. … … 127 116 'iDirectionIncrementInDegrees' 128 117 """ 129 from eccodes import *130 118 131 119 data = {} … … 141 129 142 130 # information needed from grib message 143 keys = [ 144 'Ni', 131 keys = ['Ni', 145 132 'Nj', 146 133 'latitudeOfFirstGridPointInDegrees', … … 149 136 'longitudeOfLastGridPointInDegrees', 150 137 'jDirectionIncrementInDegrees', 151 'iDirectionIncrementInDegrees' 152 ] 153 154 if verb: print('\nInformations are: ')138 'iDirectionIncrementInDegrees'] 139 140 if verb: 141 print '\nInformations are: ' 155 142 for key in keys: 156 143 # Get the value of the key in a grib message. 157 data[key] = codes_get(gid,key) 158 if verb: print "%s = %s" % (key,data[key]) 159 if verb: print '\n' 144 data[key] = codes_get(gid, key) 145 if verb: 146 print "%s = %s" % (key, data[key]) 147 if verb: 148 print '\n' 160 149 161 150 # Free the memory for the message referred as gribid. … … 164 153 return data 165 154 166 def get FilesPerDate(files, datelist):155 def get_files_per_date(files, datelist): 167 156 ''' 168 157 @Description: … … 171 160 172 161 @Input: 173 files: instance of U IOFiles162 files: instance of UioFiles 174 163 For description see class documentation. 175 164 It contains the attribute "files" which is a list of pathes … … 185 174 186 175 filelist = [] 187 for file in files:188 f date = file[-8:]189 ddate = datetime.datetime.strptime(f date, '%y%m%d%H')176 for filename in files: 177 filedate = filename[-8:] 178 ddate = datetime.datetime.strptime(filedate, '%y%m%d%H') 190 179 if ddate in datelist: 191 filelist.append(file )180 filelist.append(filename) 192 181 193 182 return filelist 194 183 195 def plot Retrieved(args,c):184 def plot_retrieved(c): 196 185 ''' 197 186 @Description: 198 187 Reads GRIB data from a specified time period, a list of levels 199 188 and a specified list of parameter. 200 201 @Input:202 args: instance of ArgumentParser203 Contains the commandline arguments from script/program call.204 205 c: instance of class ControlFile206 Contains all necessary information of a CONTROL file. The parameters207 are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,208 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,209 RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,210 SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,211 ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR,212 OUTPUTDIR, FLEXPART_ROOT_SCRIPTS213 For more information about format and content of the parameter see214 documentation.215 216 @Return:217 <nothing>218 '''219 start = datetime.datetime.strptime(c.start_date, '%Y%m%d%H')220 end = datetime.datetime.strptime(c.end_date, '%Y%m%d%H')221 222 # create datelist between start and end date223 datelist = [start] # initialise datelist with first date224 run_date = start225 while run_date < end:226 run_date += datetime.timedelta(hours=int(c.dtime))227 datelist.append(run_date)228 229 print 'datelist: ', datelist230 231 c.paramIds = asarray(c.paramIds, dtype='int')232 c.levels = asarray(c.levels, dtype='int')233 c.area = asarray(c.area)234 235 # index_keys = ["date", "time", "step"]236 # indexfile = c.inputdir + "/date_time_stepRange.idx"237 # silentremove(indexfile)238 # grib = GribTools(inputfiles.files)239 # # creates new index file240 # iid = grib.index(index_keys=index_keys, index_file=indexfile)241 242 # # read values of index keys243 # index_vals = []244 # for key in index_keys:245 # index_vals.append(grib_index_get(iid, key))246 # print(index_vals[-1])247 # # index_vals looks for example like:248 # # index_vals[0]: ('20171106', '20171107', '20171108') ; date249 # # index_vals[1]: ('0', '1200', '1800', '600') ; time250 # # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange251 252 253 254 255 #index_keys = ["date", "time", "step", "paramId"]256 #indexfile = c.inputdir + "/date_time_stepRange.idx"257 #silentremove(indexfile)258 259 files = UIOFiles(c.prefix+'*')260 files.listFiles(c.inputdir)261 ifiles = getFilesPerDate(files.files, datelist)262 ifiles.sort()263 264 gdict = getBasics(ifiles[0], verb=False)265 266 fdict = dict()267 fmeta = dict()268 fstamp = dict()269 for p in c.paramIds:270 for l in c.levels:271 key = '{:0>3}_{:0>3}'.format(p, l)272 fdict[key] = []273 fmeta[key] = []274 fstamp[key] = []275 276 for file in ifiles:277 f = open(file)278 print( "Opening file for reading data --- %s" % file)279 fdate = datetime.datetime.strptime(file[-8:], "%y%m%d%H")280 281 # Load in memory a grib message from a file.282 gid = codes_grib_new_from_file(f)283 while(gid is not None):284 gtype = codes_get(gid, 'type')285 paramId = codes_get(gid, 'paramId')286 parameterName = codes_get(gid, 'parameterName')287 level = codes_get(gid, 'level')288 289 if paramId in c.paramIds and level in c.levels:290 key = '{:0>3}_{:0>3}'.format(paramId, level)291 print 'key: ', key292 if len(fstamp[key]) != 0 :293 for i in range(len(fstamp[key])):294 if fdate < fstamp[key][i]:295 fstamp[key].insert(i, fdate)296 fmeta[key].insert(i, [paramId, parameterName, gtype,297 fdate, level])298 fdict[key].insert(i, flipud(reshape(299 codes_get_values(gid),300 [gdict['Nj'], gdict['Ni']])))301 break302 elif fdate > fstamp[key][i] and i == len(fstamp[key])-1:303 fstamp[key].append(fdate)304 fmeta[key].append([paramId, parameterName, gtype,305 fdate, level])306 fdict[key].append(flipud(reshape(307 codes_get_values(gid),308 [gdict['Nj'], gdict['Ni']])))309 break310 elif fdate > fstamp[key][i] and i != len(fstamp[key])-1 \311 and fdate < fstamp[key][i+1]:312 fstamp[key].insert(i, fdate)313 fmeta[key].insert(i, [paramId, parameterName, gtype,314 fdate, level])315 fdict[key].insert(i, flipud(reshape(316 codes_get_values(gid),317 [gdict['Nj'], gdict['Ni']])))318 break319 else:320 pass321 else:322 fstamp[key].append(fdate)323 fmeta[key].append((paramId, parameterName, gtype,324 fdate, level))325 fdict[key].append(flipud(reshape(326 codes_get_values(gid), [gdict['Nj'], gdict['Ni']])))327 328 codes_release(gid)329 330 # Load in memory a grib message from a file.331 gid = codes_grib_new_from_file(f)332 333 f.close()334 #print 'fstamp: ', fstamp335 #exit()336 337 for k in fdict.keys():338 print 'fmeta: ', len(fmeta),fmeta339 fml = fmeta[k]340 fdl = fdict[k]341 print 'fm1: ', len(fml), fml342 #print 'fd1: ', fdl343 #print zip(fdl, fml)344 for fd, fm in zip(fdl, fml):345 print fm346 ftitle = fm[1] + ' {} '.format(fm[-1]) + \347 datetime.datetime.strftime(fm[3], '%Y%m%d%H') #+ ' ' + stats(fd)348 pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \349 datetime.datetime.strftime(fm[3], '%Y%m%d%H')350 plotMap(c, fd, fm, gdict, ftitle, pname, 'png')351 352 for k in fdict.keys():353 fml = fmeta[k]354 fdl = fdict[k]355 fsl = fstamp[k]356 if fdl:357 fm = fml[0]358 fd = fdl[0]359 ftitle = fm[1] + ' {} '.format(fm[-1]) + \360 datetime.datetime.strftime(fm[3], '%Y%m%d%H') #+ ' ' + stats(fd)361 pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \362 datetime.datetime.strftime(fm[3], '%Y%m%d%H')363 lat = -20364 lon = 20365 plotTS(c, fdl, fml, fsl, lat, lon,366 gdict, ftitle, pname, 'png')367 368 return369 370 def plotTS(c, flist, fmetalist, ftimestamps, lat, lon,371 gdict, ftitle, filename, fending, show=False):372 '''373 @Description:374 375 @Input:376 c:377 378 flist:379 The actual data values to be plotted from the grib messages.380 381 fmetalist: list of strings382 Contains some meta date for the data field to be plotted:383 parameter id, parameter Name, grid type, datetime, level384 385 ftimestamps: list of datetime386 Contains the time stamps.387 388 lat:389 390 lon:391 392 gdict:393 394 ftitle: string395 The title of the timeseries.396 397 filename: string398 The time series is stored in a file with this name.399 400 fending: string401 Contains the type of plot, e.g. pdf or png402 403 show: boolean404 Decides if the plot is shown after plotting or hidden.405 406 @Return:407 <nothing>408 '''409 print 'plotting timeseries'410 411 t1 = time.time()412 413 llx = gdict['longitudeOfFirstGridPointInDegrees']414 if llx > 180. :415 llx -= 360.416 lly = gdict['latitudeOfLastGridPointInDegrees']417 dxout = gdict['iDirectionIncrementInDegrees']418 dyout = gdict['jDirectionIncrementInDegrees']419 urx = gdict['longitudeOfLastGridPointInDegrees']420 ury = gdict['latitudeOfFirstGridPointInDegrees']421 numxgrid = gdict['Ni']422 numygrid = gdict['Nj']423 424 farr = asarray(flist)425 #(time, lat, lon)426 427 lonindex = linspace(llx, urx, numxgrid)428 latindex = linspace(lly, ury, numygrid)429 #print len(lonindex), len(latindex), farr.shape430 431 #latindex = (lat + 90) * 180 / (gdict['Nj'] - 1)432 #lonindex = (lon + 179) * 360 / gdict['Ni']433 #print latindex, lonindex434 435 436 ts = farr[:, 0, 0]#latindex[0], lonindex[0]]437 438 fig = plt.figure(figsize=(12, 6.7))439 440 plt.plot(ftimestamps, ts)441 plt.title(ftitle)442 443 plt.savefig(c.outputdir+'/'+filename+'_TS.'+fending, facecolor=fig.get_facecolor(), edgecolor='none',format=fending)444 print 'created ', c.outputdir + '/' + filename445 if show == True:446 plt.show()447 fig.clf()448 plt.close(fig)449 450 print time.time() - t1, 's'451 452 return453 454 def plotMap(c, flist, fmetalist, gdict, ftitle, filename, fending, show=False):455 '''456 @Description:457 189 458 190 @Input: … … 468 200 documentation. 469 201 470 flist471 472 fmetalist473 474 gdict: dict475 Contains basic informations of the ECMWF grib files, e.g.476 'Ni', 'Nj', 'latitudeOfFirstGridPointInDegrees',477 'longitudeOfFirstGridPointInDegrees',478 'latitudeOfLastGridPointInDegrees',479 'longitudeOfLastGridPointInDegrees',480 'jDirectionIncrementInDegrees',481 'iDirectionIncrementInDegrees'482 483 ftitle: string484 The titel of the plot.485 486 filename: string487 The plot is stored in a file with this name.488 489 fending: string490 Contains the type of plot, e.g. pdf or png491 492 show: boolean493 Decides if the plot is shown after plotting or hidden.494 495 202 @Return: 496 203 <nothing> 497 204 ''' 498 print 'plotting map' 499 500 t1 = time.time() 501 502 fig = plt.figure(figsize=(12, 6.7)) 503 #mbaxes = fig.add_axes([0.05, 0.15, 0.8, 0.7]) 504 505 llx = gdict['longitudeOfFirstGridPointInDegrees'] #- 360 506 if llx > 180. : 507 llx -= 360. 508 lly = gdict['latitudeOfLastGridPointInDegrees'] 509 dxout = gdict['iDirectionIncrementInDegrees'] 510 dyout = gdict['jDirectionIncrementInDegrees'] 511 urx = gdict['longitudeOfLastGridPointInDegrees'] 512 ury = gdict['latitudeOfFirstGridPointInDegrees'] 513 numxgrid = gdict['Ni'] 514 numygrid = gdict['Nj'] 515 516 m = Basemap(projection='cyl', llcrnrlon=llx, llcrnrlat=lly, 517 urcrnrlon=urx, urcrnrlat=ury,resolution='i') 518 519 lw = 0.5 520 m.drawmapboundary() 521 x = linspace(llx, urx, numxgrid) 522 y = linspace(lly, ury, numygrid) 523 524 xx, yy = m(*meshgrid(x, y)) 525 526 #s = m.contourf(xx, yy, flist) 527 528 s = plt.imshow(flist.T, 529 extent=(llx, urx, lly, ury), 530 alpha=1.0, 531 interpolation='nearest' 532 #vmin=vn, 533 #vmax=vx, 534 #cmap=my_cmap, 535 #levels=levels, 536 #cmap=my_cmap, 537 #norm=LogNorm(vn,vx) 538 ) 539 540 title(ftitle, y=1.08) 541 cb = m.colorbar(s, location="right", pad="10%") 542 #cb.set_label('Contribution per cell (ng m$^{-3}$)',size=14) 543 544 thickline = np.arange(lly,ury+1,10.) 545 thinline = np.arange(lly,ury+1,5.) 546 m.drawparallels(thickline,color='gray',dashes=[1,1],linewidth=0.5,labels=[1,1,1,1], xoffset=1.) # draw parallels 547 m.drawparallels(np.setdiff1d(thinline,thickline),color='lightgray',dashes=[1,1],linewidth=0.5,labels=[0,0,0,0]) # draw parallels 548 549 thickline = np.arange(llx,urx+1,10.) 550 thinline = np.arange(llx,urx+1,5.) 551 m.drawmeridians(thickline,color='gray',dashes=[1,1],linewidth=0.5,labels=[1,1,1,1],yoffset=1.) # draw meridians 552 m.drawmeridians(np.setdiff1d(thinline,thickline),color='lightgray',dashes=[1,1],linewidth=0.5,labels=[0,0,0,0]) # draw meridians 553 554 m.drawcoastlines() 555 m.drawcountries() 556 557 plt.savefig(c.outputdir+'/'+filename+'_MAP.'+fending, facecolor=fig.get_facecolor(), edgecolor='none',format=fending) 558 print 'created ', c.outputdir + '/' + filename 559 if show == True: 560 plt.show() 561 fig.clf() 562 plt.close(fig) 563 564 print time.time() - t1, 's' 205 start = datetime.datetime.strptime(c.start_date, '%Y%m%d%H') 206 end = datetime.datetime.strptime(c.end_date, '%Y%m%d%H') 207 208 # create datelist between start and end date 209 datelist = [start] # initialise datelist with first date 210 run_date = start 211 while run_date < end: 212 run_date += datetime.timedelta(hours=int(c.dtime)) 213 datelist.append(run_date) 214 215 print 'datelist: ', datelist 216 217 c.paramIds = np.asarray(c.paramIds, dtype='int') 218 c.levels = np.asarray(c.levels, dtype='int') 219 c.area = np.asarray(c.area) 220 221 files = UioFiles(c.prefix+'*') 222 files.list_files(c.inputdir) 223 ifiles = get_files_per_date(files.files, datelist) 224 ifiles.sort() 225 226 gdict = get_basics(ifiles[0], verb=False) 227 228 fdict = dict() 229 fmeta = dict() 230 fstamp = dict() 231 for p in c.paramIds: 232 for l in c.levels: 233 key = '{:0>3}_{:0>3}'.format(p, l) 234 fdict[key] = [] 235 fmeta[key] = [] 236 fstamp[key] = [] 237 238 for filename in ifiles: 239 f = open(filename) 240 print "Opening file for reading data --- %s" % filename 241 fdate = datetime.datetime.strptime(filename[-8:], "%y%m%d%H") 242 243 # Load in memory a grib message from a file. 244 gid = codes_grib_new_from_file(f) 245 while gid is not None: 246 gtype = codes_get(gid, 'type') 247 paramId = codes_get(gid, 'paramId') 248 parameterName = codes_get(gid, 'parameterName') 249 level = codes_get(gid, 'level') 250 251 if paramId in c.paramIds and level in c.levels: 252 key = '{:0>3}_{:0>3}'.format(paramId, level) 253 print 'key: ', key 254 if fstamp[key]: 255 for i in range(len(fstamp[key])): 256 if fdate < fstamp[key][i]: 257 fstamp[key].insert(i, fdate) 258 fmeta[key].insert(i, [paramId, parameterName, gtype, 259 fdate, level]) 260 fdict[key].insert(i, np.flipud(np.reshape( 261 codes_get_values(gid), 262 [gdict['Nj'], gdict['Ni']]))) 263 break 264 elif fdate > fstamp[key][i] and i == len(fstamp[key])-1: 265 fstamp[key].append(fdate) 266 fmeta[key].append([paramId, parameterName, gtype, 267 fdate, level]) 268 fdict[key].append(np.flipud(np.reshape( 269 codes_get_values(gid), 270 [gdict['Nj'], gdict['Ni']]))) 271 break 272 elif fdate > fstamp[key][i] and i != len(fstamp[key])-1 \ 273 and fdate < fstamp[key][i+1]: 274 fstamp[key].insert(i, fdate) 275 fmeta[key].insert(i, [paramId, parameterName, gtype, 276 fdate, level]) 277 fdict[key].insert(i, np.flipud(np.reshape( 278 codes_get_values(gid), 279 [gdict['Nj'], gdict['Ni']]))) 280 break 281 else: 282 pass 283 else: 284 fstamp[key].append(fdate) 285 fmeta[key].append((paramId, parameterName, gtype, 286 fdate, level)) 287 fdict[key].append(np.flipud(np.reshape( 288 codes_get_values(gid), [gdict['Nj'], gdict['Ni']]))) 289 290 codes_release(gid) 291 292 # Load in memory a grib message from a file. 293 gid = codes_grib_new_from_file(f) 294 295 f.close() 296 297 for k in fdict.iterkeys(): 298 print 'fmeta: ', len(fmeta), fmeta 299 fml = fmeta[k] 300 fdl = fdict[k] 301 print 'fm1: ', len(fml), fml 302 for fd, fm in zip(fdl, fml): 303 print fm 304 ftitle = fm[1] + ' {} '.format(fm[-1]) + \ 305 datetime.datetime.strftime(fm[3], '%Y%m%d%H') 306 pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \ 307 datetime.datetime.strftime(fm[3], '%Y%m%d%H') 308 plot_map(c, fd, fm, gdict, ftitle, pname, 'png') 309 310 for k in fdict.iterkeys(): 311 fml = fmeta[k] 312 fdl = fdict[k] 313 fsl = fstamp[k] 314 if fdl: 315 fm = fml[0] 316 fd = fdl[0] 317 ftitle = fm[1] + ' {} '.format(fm[-1]) + \ 318 datetime.datetime.strftime(fm[3], '%Y%m%d%H') 319 pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \ 320 datetime.datetime.strftime(fm[3], '%Y%m%d%H') 321 lat = -20. 322 lon = 20. 323 plot_timeseries(c, fdl, fml, fsl, lat, lon, gdict, 324 ftitle, pname, 'png') 565 325 566 326 return 567 327 568 def getPlotArgs(): 328 def plot_timeseries(c, flist, fmetalist, ftimestamps, lat, lon, 329 gdict, ftitle, filename, fending, show=False): 569 330 ''' 570 331 @Description: 571 Assigns the command line arguments and reads CONTROL file 572 content. Apply default values for non mentioned arguments. 332 Creates a timeseries plot for a given lat/lon position. 573 333 574 334 @Input: 575 <nothing>576 577 @Return:578 args: instance of ArgumentParser579 Contains the commandline arguments from script/program call.580 581 335 c: instance of class ControlFile 582 336 Contains all necessary information of a CONTROL file. The parameters … … 589 343 For more information about format and content of the parameter see 590 344 documentation. 345 346 flist: numpy array, 2d 347 The actual data values to be plotted from the grib messages. 348 349 fmetalist: list of strings 350 Contains some meta date for the data field to be plotted: 351 parameter id, parameter Name, grid type, datetime, level 352 353 ftimestamps: list of datetime 354 Contains the time stamps. 355 356 lat: float 357 The latitude for which the timeseries should be plotted. 358 359 lon: float 360 The longitude for which the timeseries should be plotted. 361 362 gdict: dict 363 Contains basic informations of the ECMWF grib files, e.g. 364 'Ni', 'Nj', 'latitudeOfFirstGridPointInDegrees', 365 'longitudeOfFirstGridPointInDegrees', 366 'latitudeOfLastGridPointInDegrees', 367 'longitudeOfLastGridPointInDegrees', 368 'jDirectionIncrementInDegrees', 369 'iDirectionIncrementInDegrees' 370 371 ftitle: string 372 The title of the timeseries. 373 374 filename: string 375 The time series is stored in a file with this name. 376 377 fending: string 378 Contains the type of plot, e.g. pdf or png 379 380 show: boolean 381 Decides if the plot is shown after plotting or hidden. 382 383 @Return: 384 <nothing> 385 ''' 386 print 'plotting timeseries' 387 388 t1 = time.time() 389 390 #llx = gdict['longitudeOfFirstGridPointInDegrees'] 391 #if llx > 180. : 392 # llx -= 360. 393 #lly = gdict['latitudeOfLastGridPointInDegrees'] 394 #dxout = gdict['iDirectionIncrementInDegrees'] 395 #dyout = gdict['jDirectionIncrementInDegrees'] 396 #urx = gdict['longitudeOfLastGridPointInDegrees'] 397 #ury = gdict['latitudeOfFirstGridPointInDegrees'] 398 #numxgrid = gdict['Ni'] 399 #numygrid = gdict['Nj'] 400 401 farr = np.asarray(flist) 402 #(time, lat, lon) 403 404 #lonindex = linspace(llx, urx, numxgrid) 405 #latindex = linspace(lly, ury, numygrid) 406 407 ts = farr[:, 0, 0] 408 409 fig = plt.figure(figsize=(12, 6.7)) 410 411 plt.plot(ftimestamps, ts) 412 plt.title(ftitle) 413 414 plt.savefig(c.outputdir + '/' + filename + '_TS.' + fending, 415 facecolor=fig.get_facecolor(), 416 edgecolor='none', 417 format=fending) 418 print 'created ', c.outputdir + '/' + filename 419 if show: 420 plt.show() 421 fig.clf() 422 plt.close(fig) 423 424 print time.time() - t1, 's' 425 426 return 427 428 def plot_map(c, flist, fmetalist, gdict, ftitle, filename, fending, show=False): 429 ''' 430 @Description: 431 Creates a basemap plot with imshow for a given data field. 432 433 @Input: 434 c: instance of class ControlFile 435 Contains all necessary information of a CONTROL file. The parameters 436 are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 437 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST, 438 RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, 439 SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, 440 ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR, 441 OUTPUTDIR, FLEXPART_ROOT_SCRIPTS 442 For more information about format and content of the parameter see 443 documentation. 444 445 flist: numpy array, 2d 446 The actual data values to be plotted from the grib messages. 447 448 fmetalist: list of strings 449 Contains some meta date for the data field to be plotted: 450 parameter id, parameter Name, grid type, datetime, level 451 452 gdict: dict 453 Contains basic informations of the ECMWF grib files, e.g. 454 'Ni', 'Nj', 'latitudeOfFirstGridPointInDegrees', 455 'longitudeOfFirstGridPointInDegrees', 456 'latitudeOfLastGridPointInDegrees', 457 'longitudeOfLastGridPointInDegrees', 458 'jDirectionIncrementInDegrees', 459 'iDirectionIncrementInDegrees' 460 461 ftitle: string 462 The titel of the plot. 463 464 filename: string 465 The plot is stored in a file with this name. 466 467 fending: string 468 Contains the type of plot, e.g. pdf or png 469 470 show: boolean 471 Decides if the plot is shown after plotting or hidden. 472 473 @Return: 474 <nothing> 475 ''' 476 print 'plotting map' 477 478 t1 = time.time() 479 480 fig = plt.figure(figsize=(12, 6.7)) 481 #mbaxes = fig.add_axes([0.05, 0.15, 0.8, 0.7]) 482 483 llx = gdict['longitudeOfFirstGridPointInDegrees'] #- 360 484 if llx > 180.: 485 llx -= 360. 486 lly = gdict['latitudeOfLastGridPointInDegrees'] 487 #dxout = gdict['iDirectionIncrementInDegrees'] 488 #dyout = gdict['jDirectionIncrementInDegrees'] 489 urx = gdict['longitudeOfLastGridPointInDegrees'] 490 ury = gdict['latitudeOfFirstGridPointInDegrees'] 491 #numxgrid = gdict['Ni'] 492 #numygrid = gdict['Nj'] 493 494 m = Basemap(projection='cyl', llcrnrlon=llx, llcrnrlat=lly, 495 urcrnrlon=urx, urcrnrlat=ury, resolution='i') 496 497 #lw = 0.5 498 m.drawmapboundary() 499 #x = linspace(llx, urx, numxgrid) 500 #y = linspace(lly, ury, numygrid) 501 502 #xx, yy = m(*meshgrid(x, y)) 503 504 #s = m.contourf(xx, yy, flist) 505 506 s = plt.imshow(flist.T, 507 extent=(llx, urx, lly, ury), 508 alpha=1.0, 509 interpolation='nearest' 510 #vmin=vn, 511 #vmax=vx, 512 #cmap=my_cmap, 513 #levels=levels, 514 #cmap=my_cmap, 515 #norm=LogNorm(vn,vx) 516 ) 517 518 plt.title(ftitle, y=1.08) 519 cb = m.colorbar(s, location="right", pad="10%") 520 cb.set_label('label', size=14) 521 522 thickline = np.arange(lly, ury+1, 10.) 523 thinline = np.arange(lly, ury+1, 5.) 524 m.drawparallels(thickline, 525 color='gray', 526 dashes=[1, 1], 527 linewidth=0.5, 528 labels=[1, 1, 1, 1], 529 xoffset=1.) 530 m.drawparallels(np.setdiff1d(thinline, thickline), 531 color='lightgray', 532 dashes=[1, 1], 533 linewidth=0.5, 534 labels=[0, 0, 0, 0]) 535 536 thickline = np.arange(llx, urx+1, 10.) 537 thinline = np.arange(llx, urx+1, 5.) 538 m.drawmeridians(thickline, 539 color='gray', 540 dashes=[1, 1], 541 linewidth=0.5, 542 labels=[1, 1, 1, 1], 543 yoffset=1.) 544 m.drawmeridians(np.setdiff1d(thinline, thickline), 545 color='lightgray', 546 dashes=[1, 1], 547 linewidth=0.5, 548 labels=[0, 0, 0, 0]) 549 550 m.drawcoastlines() 551 m.drawcountries() 552 553 plt.savefig(c.outputdir + '/' + filename + '_MAP.' + fending, 554 facecolor=fig.get_facecolor(), 555 edgecolor='none', 556 format=fending) 557 print 'created ', c.outputdir + '/' + filename 558 if show: 559 plt.show() 560 fig.clf() 561 plt.close(fig) 562 563 print time.time() - t1, 's' 564 565 return 566 567 def get_plot_args(): 568 ''' 569 @Description: 570 Assigns the command line arguments and reads CONTROL file 571 content. Apply default values for non mentioned arguments. 572 573 @Input: 574 <nothing> 575 576 @Return: 577 args: instance of ArgumentParser 578 Contains the commandline arguments from script/program call. 579 580 c: instance of class ControlFile 581 Contains all necessary information of a CONTROL file. The parameters 582 are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 583 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST, 584 RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, 585 SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, 586 ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR, 587 OUTPUTDIR, FLEXPART_ROOT_SCRIPTS 588 For more information about format and content of the parameter see 589 documentation. 591 590 ''' 592 591 parser = ArgumentParser(description='Plot retrieved GRIB data from ' + \ … … 597 596 parser.add_argument("--start_date", dest="start_date", 598 597 help="start date YYYYMMDD") 599 parser.add_argument( 600 598 parser.add_argument("--end_date", dest="end_date", 599 help="end_date YYYYMMDD") 601 600 602 601 parser.add_argument("--start_step", dest="start_step", 603 602 help="start step in hours") 604 parser.add_argument( 605 603 parser.add_argument("--end_step", dest="end_step", 604 help="end step in hours") 606 605 607 606 # some arguments that override the default in the CONTROL file … … 635 634 except IOError: 636 635 try: 637 c = ControlFile(localpythonpath + args.controlfile) 638 639 except: 636 c = ControlFile(LOCAL_PYTHON_PATH + args.controlfile) 637 except IOError: 640 638 print 'Could not read CONTROL file "' + args.controlfile + '"' 641 639 print 'Either it does not exist or its syntax is wrong.' … … 682 680 if __name__ == "__main__": 683 681 main() 684 -
python/prepare_flexpart.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # TODO AP 5 # - wieso cleanup in main wenn es in prepareflexpart bereits abgefragt wurde? 6 # doppelt gemoppelt? 4 # ToDo AP 7 5 # - wieso start=startm1 wenn basetime = 0 ? wenn die fluxes nicht mehr 8 6 # relevant sind? verstehe ich nicht … … 29 27 # - added documentation 30 28 # - minor changes in programming style for consistence 31 # - BUG: removed call of clean up-Function after call of29 # - BUG: removed call of clean_up-Function after call of 32 30 # prepareFlexpart in main since it is already called in 33 31 # prepareFlexpart at the end! 34 32 # - created function main and moved the two function calls for 35 # arguments and prepare FLEXPARTinto it33 # arguments and prepare_flexpart into it 36 34 # 37 35 # @License: … … 44 42 # This program prepares the final version of the grib files which are 45 43 # then used by FLEXPART. It converts the bunch of grib files extracted 46 # via get MARSdata by doing for example the necessary conversion to get44 # via get_mars_data by doing for example the necessary conversion to get 47 45 # consistent grids or the disaggregation of flux data. Finally, the 48 46 # program combines the data fields in files per available hour with the … … 52 50 # @Program Content: 53 51 # - main 54 # - prepare FLEXPART52 # - prepare_flexpart 55 53 # 56 54 #******************************************************************************* … … 59 57 # MODULES 60 58 # ------------------------------------------------------------------------------ 61 import shutil62 59 import datetime 63 #import time64 60 import os 65 61 import inspect 66 62 import sys 67 63 import socket 68 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter69 64 70 hostname = socket.gethostname() 71 ecapi = 'ecmwf' not in hostname 65 # software specific classes and modules from flex_extract 66 from UioFiles import UioFiles 67 from tools import interpret_args_and_control, clean_up 68 from EcFlexpart import EcFlexpart 69 70 ecapi = 'ecmwf' not in socket.gethostname() 72 71 try: 73 72 if ecapi: … … 77 76 78 77 # add path to pythonpath so that python finds its buddies 79 localpythonpath= os.path.dirname(os.path.abspath(78 LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath( 80 79 inspect.getfile(inspect.currentframe()))) 81 if localpythonpathnot in sys.path:82 sys.path.append( localpythonpath)80 if LOCAL_PYTHON_PATH not in sys.path: 81 sys.path.append(LOCAL_PYTHON_PATH) 83 82 84 # software specific classes and modules from flex_extract85 from UIOFiles import UIOFiles86 from Tools import interpret_args_and_control, cleanup87 from ECFlexpart import ECFlexpart88 83 # ------------------------------------------------------------------------------ 89 84 # FUNCTION … … 92 87 ''' 93 88 @Description: 94 If prepare FLEXPARTis called from command line, this function controls89 If prepare_flexpart is called from command line, this function controls 95 90 the program flow and calls the argumentparser function and 96 the prepare FLEXPARTfunction for preparation of GRIB data for FLEXPART.91 the prepare_flexpart function for preparation of GRIB data for FLEXPART. 97 92 98 93 @Input: … … 103 98 ''' 104 99 args, c = interpret_args_and_control() 105 prepare FLEXPART(args, c)100 prepare_flexpart(args, c) 106 101 107 102 return 108 103 109 def prepare FLEXPART(args, c):104 def prepare_flexpart(args, c): 110 105 ''' 111 106 @Description: 112 Lists all grib files retrieved from MARS with get MARSdata and107 Lists all grib files retrieved from MARS with get_mars_data and 113 108 uses prepares data for the use in FLEXPART. Specific data fields 114 109 are converted to a different grid and the flux data are going to be … … 157 152 # one day after the end date is needed 158 153 startm1 = start - datetime.timedelta(days=1) 159 endp1 = end + datetime.timedelta(days=1)154 # endp1 = end + datetime.timedelta(days=1) 160 155 161 156 # get all files with flux data to be deaccumulated 162 inputfiles = U IOFiles('*OG_acc_SL*.' + c.ppid + '.*')163 inputfiles.list Files(c.inputdir)157 inputfiles = UioFiles('*OG_acc_SL*.' + c.ppid + '.*') 158 inputfiles.list_files(c.inputdir) 164 159 165 160 # create output dir if necessary … … 168 163 169 164 # deaccumulate the flux data 170 flexpart = E CFlexpart(c, fluxes=True)165 flexpart = EcFlexpart(c, fluxes=True) 171 166 flexpart.write_namelist(c, 'fort.4') 172 167 flexpart.deacc_fluxes(inputfiles, c) 173 168 174 print ('Prepare ' + start.strftime("%Y%m%d") +175 "/to/" + end.strftime("%Y%m%d") )169 print 'Prepare ' + start.strftime("%Y%m%d") + \ 170 "/to/" + end.strftime("%Y%m%d") 176 171 177 172 # get a list of all files from the root inputdir 178 inputfiles = U IOFiles('????__??.*' + c.ppid + '.*')179 inputfiles.list Files(c.inputdir)173 inputfiles = UioFiles('????__??.*' + c.ppid + '.*') 174 inputfiles.list_files(c.inputdir) 180 175 181 176 # produce FLEXPART-ready GRIB files and … … 185 180 start = startm1 186 181 187 flexpart = E CFlexpart(c, fluxes=False)182 flexpart = EcFlexpart(c, fluxes=False) 188 183 flexpart.create(inputfiles, c) 189 184 flexpart.process_output(c) … … 192 187 # otherwise delete temporary files 193 188 if int(c.debug) != 0: 194 print ('Temporary files left intact')189 print 'Temporary files left intact' 195 190 else: 196 clean up(c)191 clean_up(c) 197 192 198 193 return -
python/profiling.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # T ODOAP4 # ToDo AP 5 5 # - check of license of book content 6 6 #************************************************************************ … … 66 66 result = fn(*args, **kwargs) 67 67 t2 = time.time() 68 print ("@timefn:" + fn.func_name + " took " + str(t2 - t1) + " seconds")68 print "@timefn:" + fn.func_name + " took " + str(t2 - t1) + " seconds" 69 69 70 70 return result -
python/submit.py
r991df6a rff99eae 1 1 #!/usr/bin/env python 2 2 # -*- coding: utf-8 -*- 3 #************************************************************************4 # TODO AP5 #6 # - Change History ist nicht angepasst ans File! Original geben lassen7 # - dead code ? what to do?8 # - seperate operational and reanlysis for clarification9 # - divide in two submits , ondemand und operational10 # -11 #************************************************************************12 13 3 #******************************************************************************* 14 4 # @Author: Anne Fouilloux (University of Oslo) … … 37 27 # program flow. 38 28 # If it is supposed to work locally then it works through the necessary 39 # functions get MARSdata and prepareFlexpart. Otherwise it prepares29 # functions get_mars_data and prepareFlexpart. Otherwise it prepares 40 30 # a shell job script which will do the necessary work on the 41 31 # ECMWF server and is submitted via ecaccess-job-submit. … … 52 42 import os 53 43 import sys 54 import glob55 44 import subprocess 56 45 import inspect 57 # add path to pythonpath so that python finds its buddies58 localpythonpath = os.path.dirname(os.path.abspath(59 inspect.getfile(inspect.currentframe())))60 if localpythonpath not in sys.path:61 sys.path.append(localpythonpath)62 46 63 47 # software specific classes and modules from flex_extract 64 from Tools import interpret_args_and_control, normalexit 65 from getMARSdata import getMARSdata 66 from prepareFLEXPART import prepareFLEXPART 48 from tools import interpret_args_and_control, normal_exit 49 from get_mars_data import get_mars_data 50 from prepare_flexpart import prepare_flexpart 51 52 # add path to pythonpath so that python finds its buddies 53 LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath( 54 inspect.getfile(inspect.currentframe()))) 55 if LOCAL_PYTHON_PATH not in sys.path: 56 sys.path.append(LOCAL_PYTHON_PATH) 67 57 68 58 # ------------------------------------------------------------------------------ … … 84 74 <nothing> 85 75 ''' 86 calledfromdir = os.getcwd() 76 77 called_from_dir = os.getcwd() 87 78 args, c = interpret_args_and_control() 79 88 80 # on local side 89 81 if args.queue is None: 90 82 if c.inputdir[0] != '/': 91 c.inputdir = os.path.join(called fromdir, c.inputdir)83 c.inputdir = os.path.join(called_from_dir, c.inputdir) 92 84 if c.outputdir[0] != '/': 93 c.outputdir = os.path.join(called fromdir, c.outputdir)94 get MARSdata(args, c)95 prepare FLEXPART(args, c)96 normal exit(c)85 c.outputdir = os.path.join(called_from_dir, c.outputdir) 86 get_mars_data(args, c) 87 prepare_flexpart(args, c) 88 normal_exit(c) 97 89 # on ECMWF server 98 90 else: … … 139 131 140 132 # put all parameters of ControlFile instance into a list 141 clist = c.to list()133 clist = c.to_list() # ondemand 142 134 colist = [] # operational 143 135 mt = 0 144 136 145 #AP wieso 2 for loops?146 137 for elem in clist: 147 138 if 'maxstep' in elem: … … 152 143 elem = 'start_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}' 153 144 if 'end_date' in elem: 154 #AP Fehler?! Muss end_date heissen 155 elem = 'start_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}' 145 elem = 'end_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}' 156 146 if 'base_time' in elem: 157 147 elem = 'base_time ' + '${MSJ_BASETIME}' … … 173 163 p = subprocess.check_call(['ecaccess-job-submit', '-queueName', 174 164 queue, 'job.ksh']) 175 except: 176 print('ecaccess-job-submit failed, probably eccert has expired') 165 except subprocess.CalledProcessError as e: 166 print 'ecaccess-job-submit failed!' 167 print 'Error Message: ' 168 print e.output 177 169 exit(1) 178 170 179 print ('You should get an email with subject flex.hostname.pid')171 print 'You should get an email with subject flex.hostname.pid' 180 172 181 173 return -
python/test_suite.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # TODO AP 5 # 4 # ToDo AP 6 5 # - provide more tests 7 6 # - provide more documentation 8 # -9 7 #************************************************************************ 10 8 … … 28 26 # @Program Functionality: 29 27 # This script triggers the ECMWFDATA test suite. Call with 30 # test suite.py [test group]28 # test_suite.py [test group] 31 29 # 32 30 # @Program Content: … … 46 44 # ------------------------------------------------------------------------------ 47 45 try: 48 taskfile = open('test suite.json')49 except :50 print 'could not open suite definition file test suite.json'46 taskfile = open('test_suite.json') 47 except IOError: 48 print 'could not open suite definition file test_suite.json' 51 49 exit() 52 50 … … 73 71 try: 74 72 tk, tv = g, tasks[g] 75 except:76 continue73 finally: 74 pass 77 75 garglist = [] 78 76 for ttk, ttv in tv.iteritems(): … … 80 78 if ttk != 'script': 81 79 garglist.append('--' + ttk) 82 if '$' == ttv[0]:80 if ttv[0] == '$': 83 81 garglist.append(os.path.expandvars(ttv)) 84 82 else: … … 89 87 for tttk, tttv in ttv.iteritems(): 90 88 if isinstance(tttv, basestring): 91 92 93 94 95 89 arglist.append('--' + tttk) 90 if '$' in tttv[0]: 91 arglist.append(os.path.expandvars(tttv)) 92 else: 93 arglist.append(tttv) 96 94 print 'Command: ', ' '.join([tv['script']] + garglist + arglist) 97 95 o = '../test/' + tk + '_' + ttk + '_' + '_'.join(ttv.keys()) … … 101 99 p = subprocess.check_call([tv['script']] + garglist + arglist, 102 100 stdout=f, stderr=f) 103 except :101 except subprocess.CalledProcessError as e: 104 102 f.write('\nFAILED\n') 105 103 print 'FAILED' … … 111 109 print str(jobcounter-jobfailed) + ' successful, ' + str(jobfailed) + ' failed' 112 110 print 'If tasks have been submitted via ECACCESS please check emails' 113 114 #print tasks -
python/tools.py
r991df6a rff99eae 2 2 # -*- coding: utf-8 -*- 3 3 #************************************************************************ 4 # T ODOAP5 # - check my error6 # - check normal exit7 # - check get ListAsString4 # ToDo AP 5 # - check my_error 6 # - check normal_exit 7 # - check get_list_as_string 8 8 # - seperate args and control interpretation 9 9 #************************************************************************ … … 15 15 # @Change History: 16 16 # October 2014 - Anne Fouilloux (University of Oslo) 17 # - created functions silent remove and product (taken from ECMWF)17 # - created functions silent_remove and product (taken from ECMWF) 18 18 # 19 19 # November 2015 - Leopold Haimberger (University of Vienna) 20 # - created functions: interpret_args_and_control, clean up21 # my error, normalexit, init128, toparamId20 # - created functions: interpret_args_and_control, clean_up 21 # my_error, normal_exit, init128, to_param_id 22 22 # 23 23 # April 2018 - Anne Philipp (University of Vienna): 24 24 # - applied PEP8 style guide 25 25 # - added documentation 26 # - moved all functions from file Flexpart Tools to this file Tools27 # - added function get ListAsString26 # - moved all functions from file Flexparttools to this file tools 27 # - added function get_list_as_string 28 28 # 29 29 # @License: … … 39 39 # @Module Content: 40 40 # - interpret_args_and_control 41 # - clean up42 # - my error43 # - normal exit41 # - clean_up 42 # - my_error 43 # - normal_exit 44 44 # - product 45 # - silent remove45 # - silent_remove 46 46 # - init128 47 # - to paramId48 # - get ListAsString47 # - to_param_id 48 # - get_list_as_string 49 49 # 50 50 #******************************************************************************* … … 57 57 import sys 58 58 import glob 59 import inspect 60 import subprocess 59 61 import traceback 60 from numpy import *61 from gribapi import *62 62 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter 63 import numpy as np 63 64 64 65 # software specific class from flex_extract … … 123 124 parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts", 124 125 help="FLEXPART root directory (to find grib2flexpart \ 125 and COMMAND file)\n \Normally ECMWFDATA resides in \126 and COMMAND file)\n Normally ECMWFDATA resides in \ 126 127 the scripts directory of the FLEXPART distribution") 127 128 128 # this is only used by prepare FLEXPART.py to rerun a postprocessing step129 # this is only used by prepare_flexpart.py to rerun a postprocessing step 129 130 parser.add_argument("--ppid", dest="ppid", 130 131 help="Specify parent process id for \ 131 rerun of prepare FLEXPART")132 rerun of prepare_flexpart") 132 133 133 134 # arguments for job submission to ECMWF, only needed by submit.py … … 152 153 except IOError: 153 154 try: 154 c = ControlFile(localpythonpath + args.controlfile) 155 except: 156 print('Could not read CONTROL file "' + args.controlfile + '"') 157 print('Either it does not exist or its syntax is wrong.') 158 print('Try "' + sys.argv[0].split('/')[-1] + 159 ' -h" to print usage information') 155 LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath( 156 inspect.getfile(inspect.currentframe()))) 157 c = ControlFile(LOCAL_PYTHON_PATH + args.controlfile) 158 except IOError: 159 print 'Could not read CONTROL file "' + args.controlfile + '"' 160 print 'Either it does not exist or its syntax is wrong.' 161 print 'Try "' + sys.argv[0].split('/')[-1] + \ 162 ' -h" to print usage information' 160 163 exit(1) 161 164 162 165 # check for having at least a starting date 163 166 if args.start_date is None and getattr(c, 'start_date') is None: 164 print ('start_date specified neither in command line nor \165 in CONTROL file ' + args.controlfile )166 print ('Try "' + sys.argv[0].split('/')[-1] +167 ' -h" to print usage information' )167 print 'start_date specified neither in command line nor \ 168 in CONTROL file ' + args.controlfile 169 print 'Try "' + sys.argv[0].split('/')[-1] + \ 170 ' -h" to print usage information' 168 171 exit(1) 169 172 … … 206 209 l = args.area.split('/') 207 210 if afloat: 208 for i in range(len(l)):209 l[i] = str(int(float(l[i]) * 1000))211 for i, item in enumerate(l): 212 item = str(int(float(item) * 1000)) 210 213 c.upper, c.left, c.lower, c.right = l 211 214 … … 218 221 if 'to' in args.step.lower(): 219 222 if 'by' in args.step.lower(): 220 ilist = arange(int(l[0]), int(l[2]) + 1, int(l[4]))223 ilist = np.arange(int(l[0]), int(l[2]) + 1, int(l[4])) 221 224 c.step = ['{:0>3}'.format(i) for i in ilist] 222 225 else: 223 my error(None, args.step + ':\n' +224 'please use "by" as well if "to" is used')226 my_error(None, args.step + ':\n' + 227 'please use "by" as well if "to" is used') 225 228 else: 226 229 c.step = l … … 239 242 240 243 241 def clean up(c):244 def clean_up(c): 242 245 ''' 243 246 @Description: … … 263 266 ''' 264 267 265 print ("cleanup")268 print "clean_up" 266 269 267 270 cleanlist = glob.glob(c.inputdir + "/*") 268 271 for cl in cleanlist: 269 272 if c.prefix not in cl: 270 silent remove(cl)273 silent_remove(cl) 271 274 if c.ecapi is False and (c.ectrans == '1' or c.ecstorage == '1'): 272 silent remove(cl)273 274 print ("Done")275 silent_remove(cl) 276 277 print "Done" 275 278 276 279 return 277 280 278 281 279 def my error(c, message='ERROR'):282 def my_error(c, message='ERROR'): 280 283 ''' 281 284 @Description: … … 303 306 <nothing> 304 307 ''' 305 # uncomment if user wants email notification directly from python 306 #try: 307 #target = c.mailfail 308 #except AttributeError: 309 #target = os.getenv('USER') 310 311 #if(type(target) is not list): 312 #target = [target] 313 314 print(message) 315 316 # uncomment if user wants email notification directly from python 317 #for t in target: 318 #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)], 319 # stdin = subprocess.PIPE, stdout = subprocess.PIPE, 320 # stderr = subprocess.PIPE, bufsize = 1) 321 #tr = '\n'.join(traceback.format_stack()) 322 #pout = p.communicate(input = message+'\n\n'+tr)[0] 323 #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode() 308 309 print message 310 311 # comment if user does not want email notification directly from python 312 try: 313 target = [] 314 target.extend(c.mailfail) 315 except AttributeError: 316 target = [] 317 target.extend(os.getenv('USER')) 318 319 for t in target: 320 p = subprocess.Popen(['mail', '-s ECMWFDATA v7.0 ERROR', 321 os.path.expandvars(t)], 322 stdin=subprocess.PIPE, 323 stdout=subprocess.PIPE, 324 stderr=subprocess.PIPE, 325 bufsize=1) 326 tr = '\n'.join(traceback.format_stack()) 327 pout = p.communicate(input=message + '\n\n' + tr)[0] 328 print 'Email sent to ' + os.path.expandvars(t) + ' ' + pout.decode() 324 329 325 330 exit(1) … … 328 333 329 334 330 def normal exit(c, message='Done!'):335 def normal_exit(c, message='Done!'): 331 336 ''' 332 337 @Description: … … 354 359 355 360 ''' 356 # Uncomment if user wants notification directly from python357 #try: 358 #target = c.mailops359 #if(type(target) is not list):360 #target = [target]361 #for t in target:362 #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit',363 # os.path.expandvars(t)],364 # stdin = subprocess.PIPE,365 # stdout =subprocess.PIPE,366 # stderr = subprocess.PIPE, bufsize = 1)367 #pout = p.communicate(input = message+'\n\n')[0]368 #print pout.decode()369 #except:370 #pass371 372 print(message)361 print message 362 363 # comment if user does not want notification directly from python 364 try: 365 target = [] 366 target.extend(c.mailops) 367 for t in target: 368 p = subprocess.Popen(['mail', '-s ECMWFDATA v7.0 normal exit', 369 os.path.expandvars(t)], 370 stdin=subprocess.PIPE, 371 stdout=subprocess.PIPE, 372 stderr=subprocess.PIPE, 373 bufsize=1) 374 pout = p.communicate(input=message+'\n\n')[0] 375 print 'Email sent to ' + os.path.expandvars(t) + ' ' + pout.decode() 376 finally: 377 pass 373 378 374 379 return … … 411 416 412 417 413 def silent remove(filename):418 def silent_remove(filename): 414 419 ''' 415 420 @Description: … … 435 440 436 441 437 def init128(f n):442 def init128(filepath): 438 443 ''' 439 444 @Description: … … 441 446 442 447 @Input: 443 f n: string448 filepath: string 444 449 Path to file of ECMWF grib table number 128. 445 450 … … 451 456 ''' 452 457 table128 = dict() 453 with open(f n) as f:458 with open(filepath) as f: 454 459 fdata = f.read().split('\n') 455 460 for data in fdata: … … 460 465 461 466 462 def to paramId(pars, table):467 def to_param_id(pars, table): 463 468 ''' 464 469 @Description: … … 486 491 ipar = [] 487 492 for par in cpar: 488 found = False489 493 for k, v in table.iteritems(): 490 494 if par == k or par == v: 491 495 ipar.append(int(k)) 492 found = True493 496 break 494 if found is False:495 print ('Warning: par ' + par + ' not found in table 128')497 else: 498 print 'Warning: par ' + par + ' not found in table 128' 496 499 497 500 return ipar 498 501 499 def get ListAsString(listObj):502 def get_list_as_string(list_obj, concatenate_sign=', '): 500 503 ''' 501 504 @Description: … … 503 506 504 507 @Input: 505 list Obj: list508 list_obj: list 506 509 A list with arbitrary content. 507 510 508 @Return: 509 strOfList: string 511 concatenate_sign: string, optional 512 A string which is used to concatenate the single 513 list elements. Default value is ", ". 514 515 @Return: 516 str_of_list: string 510 517 The content of the list as a single string. 511 518 ''' 512 519 513 str OfList = ", ".join( str(l) for l in listObj)514 515 return str OfList520 str_of_list = concatenate_sign.join(str(l) for l in list_obj) 521 522 return str_of_list
Note: See TracChangeset
for help on using the changeset viewer.