Changeset ff99eae in flex_extract.git for python/EcFlexpart.py
- Timestamp:
- Jun 1, 2018, 8:34:59 PM (6 years ago)
- Branches:
- master, ctbto, dev
- Children:
- e1228f3
- Parents:
- ccab809
- File:
-
- 1 moved
-
python/EcFlexpart.py (moved) (moved from python/ECFlexpart.py) (57 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
