Changeset ff99eae in flex_extract.git for python/plot_retrieved.py


Ignore:
Timestamp:
Jun 1, 2018, 8:34:59 PM (6 years ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
master, ctbto, dev
Children:
e1228f3
Parents:
ccab809
Message:

completed application of pep8 style guide and pylint investigations. added documentation almost everywhere

File:
1 edited

Legend:

Unmodified
Added
Removed

  • python/plot_retrieved.py

    rccab809 rff99eae  
    22# -*- coding: utf-8 -*-
    33#************************************************************************
    4 # TODO AP
     4# ToDo AP
    55# - documentation der Funktionen
    66# - docu der progam functionality
    … …  
    1919#        - created function main and moved the two function calls for
    2020#          arguments and plotting into it
    21 #        - added function getBasics to extract the boundary conditions
     21#        - added function get_basics to extract the boundary conditions
    2222#          of the data fields from the first grib file it gets.
    2323#
    … …  
    3333# @Program Content:
    3434#    - main
    35 #    - getBasics
    36 #    - plotRetrieved
    37 #    - plotTS
    38 #    - plotMap
    39 #    - getPlotArgs
     35#    - get_basics
     36#    - plot_retrieved
     37#    - plot_timeseries
     38#    - plot_map
     39#    - get_plot_args
    4040#
    4141#*******************************************************************************
    … …  
    4949import inspect
    5050import sys
    51 import glob
    5251from 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
     52import matplotlib
     53import matplotlib.pyplot as plt
     54from mpl_toolkits.basemap import Basemap
     55from eccodes import codes_grib_new_from_file, codes_get, codes_release, \
     56                    codes_get_values
     57import numpy as np
     58
     59# software specific classes and modules from flex_extract
     60from ControlFile import ControlFile
     61from UioFiles import UioFiles
     62
     63# add path to pythonpath so that python finds its buddies
     64LOCAL_PYTHON_PATH = os.path.dirname(os.path.abspath(
     65    inspect.getfile(inspect.currentframe())))
     66if LOCAL_PYTHON_PATH not in sys.path:
     67    sys.path.append(LOCAL_PYTHON_PATH)
    6368
    6469font = {'family': 'monospace', 'size': 12}
    6570matplotlib.rcParams['xtick.major.pad'] = '20'
    66 
    6771matplotlib.rc('font', **font)
    68 
    69 from eccodes import *
    70 import numpy as np
    71 
    72 # add path to pythonpath so that python finds its buddies
    73 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_extract
    79 from Tools import silentremove
    80 from ControlFile import ControlFile
    81 #from GribTools import GribTools
    82 from UIOFiles import UIOFiles
    83 
    8472# ------------------------------------------------------------------------------
    8573# FUNCTION
    … …  
    8876    '''
    8977    @Description:
    90         If plotRetrieved is called from command line, this function controls
     78        If plot_retrieved is called from command line, this function controls
    9179        the program flow and calls the argumentparser function and
    92         the plotRetrieved function for plotting the retrieved GRIB data.
     80        the plot_retrieved function for plotting the retrieved GRIB data.
    9381
    9482    @Input:
    … …  
    9886        <nothing>
    9987    '''
    100     args, c = getPlotArgs()
    101     plotRetrieved(args, c)
     88    args, c = get_plot_args()
     89    plot_retrieved(c)
    10290
    10391    return
    10492
    105 def getBasics(ifile, verb=False):
     93def get_basics(ifile, verb=False):
    10694    """
    10795    @Description:
    … …  
    113101        ifile: string
    114102            Contains the full absolute path to the ECMWF grib file.
     103
    115104        verb (opt): bool
    116105            Is True if there should be extra output in verbose mode.
    … …  
    127116            'iDirectionIncrementInDegrees'
    128117    """
    129     from eccodes import *
    130118
    131119    data = {}
    … …  
    141129
    142130        # information needed from grib message
    143         keys = [
    144                 'Ni',
     131        keys = ['Ni',
    145132                'Nj',
    146133                'latitudeOfFirstGridPointInDegrees',
    … …  
    149136                'longitudeOfLastGridPointInDegrees',
    150137                'jDirectionIncrementInDegrees',
    151                 'iDirectionIncrementInDegrees'
    152                ]
    153 
    154         if verb: print('\nInformations are: ')
     138                'iDirectionIncrementInDegrees']
     139
     140        if verb:
     141            print '\nInformations are: '
    155142        for key in keys:
    156143            # 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'
    160149
    161150        # Free the memory for the message referred as gribid.
    … …  
    164153    return data
    165154
    166 def getFilesPerDate(files, datelist):
     155def get_files_per_date(files, datelist):
    167156    '''
    168157    @Description:
    … …  
    171160
    172161    @Input:
    173         files: instance of UIOFiles
     162        files: instance of UioFiles
    174163            For description see class documentation.
    175164            It contains the attribute "files" which is a list of pathes
    … …  
    185174
    186175    filelist = []
    187     for file in files:
    188         fdate = file[-8:]
    189         ddate = datetime.datetime.strptime(fdate, '%y%m%d%H')
     176    for filename in files:
     177        filedate = filename[-8:]
     178        ddate = datetime.datetime.strptime(filedate, '%y%m%d%H')
    190179        if ddate in datelist:
    191             filelist.append(file)
     180            filelist.append(filename)
    192181
    193182    return filelist
    194183
    195 def plotRetrieved(args, c):
     184def plot_retrieved(c):
    196185    '''
    197186    @Description:
    198187        Reads GRIB data from a specified time period, a list of levels
    199188        and a specified list of parameter.
    200 
    201     @Input:
    202         args: instance of ArgumentParser
    203             Contains the commandline arguments from script/program call.
    204 
    205         c: instance of class ControlFile
    206             Contains all necessary information of a CONTROL file. The parameters
    207             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_SCRIPTS
    213             For more information about format and content of the parameter see
    214             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 date
    223     datelist = [start] # initialise datelist with first date
    224     run_date = start
    225     while run_date < end:
    226         run_date += datetime.timedelta(hours=int(c.dtime))
    227         datelist.append(run_date)
    228 
    229     print 'datelist: ', datelist
    230 
    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 file
    240     # iid = grib.index(index_keys=index_keys, index_file=indexfile)
    241 
    242     # # read values of index keys
    243     # 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') ; date
    249         # # index_vals[1]: ('0', '1200', '1800', '600') ; time
    250         # # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
    251 
    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: ', key
    292                 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                             break
    302                         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                             break
    310                         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                             break
    319                         else:
    320                             pass
    321                 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: ', fstamp
    335     #exit()
    336 
    337     for k in fdict.keys():
    338         print 'fmeta: ', len(fmeta),fmeta
    339         fml = fmeta[k]
    340         fdl = fdict[k]
    341         print 'fm1: ', len(fml), fml
    342         #print 'fd1: ', fdl
    343         #print zip(fdl, fml)
    344         for fd, fm in zip(fdl, fml):
    345             print fm
    346             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 = -20
    364             lon = 20
    365             plotTS(c, fdl, fml, fsl, lat, lon,
    366                            gdict, ftitle, pname, 'png')
    367 
    368     return
    369 
    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 strings
    382             Contains some meta date for the data field to be plotted:
    383             parameter id, parameter Name, grid type, datetime, level
    384 
    385         ftimestamps: list of datetime
    386             Contains the time stamps.
    387 
    388         lat:
    389 
    390         lon:
    391 
    392         gdict:
    393 
    394         ftitle: string
    395             The title of the timeseries.
    396 
    397         filename: string
    398             The time series is stored in a file with this name.
    399 
    400         fending: string
    401             Contains the type of plot, e.g. pdf or png
    402 
    403         show: boolean
    404             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.shape
    430 
    431     #latindex = (lat + 90) * 180 / (gdict['Nj'] - 1)
    432     #lonindex = (lon + 179) * 360 / gdict['Ni']
    433     #print latindex, lonindex
    434 
    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 + '/' + filename
    445     if show == True:
    446         plt.show()
    447     fig.clf()
    448     plt.close(fig)
    449 
    450     print time.time() - t1, 's'
    451 
    452     return
    453 
    454 def plotMap(c, flist, fmetalist, gdict, ftitle, filename, fending, show=False):
    455     '''
    456     @Description:
    457189
    458190    @Input:
    … …  
    468200            documentation.
    469201
    470         flist
    471 
    472         fmetalist
    473 
    474         gdict: dict
    475             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: string
    484             The titel of the plot.
    485 
    486         filename: string
    487             The plot is stored in a file with this name.
    488 
    489         fending: string
    490             Contains the type of plot, e.g. pdf or png
    491 
    492         show: boolean
    493             Decides if the plot is shown after plotting or hidden.
    494 
    495202    @Return:
    496203        <nothing>
    497204    '''
    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')
    565325
    566326    return
    567327
    568 def getPlotArgs():
     328def plot_timeseries(c, flist, fmetalist, ftimestamps, lat, lon,
     329                    gdict, ftitle, filename, fending, show=False):
    569330    '''
    570331    @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.
    573333
    574334    @Input:
    575         <nothing>
    576 
    577     @Return:
    578         args: instance of ArgumentParser
    579             Contains the commandline arguments from script/program call.
    580 
    581335        c: instance of class ControlFile
    582336            Contains all necessary information of a CONTROL file. The parameters
    … …  
    589343            For more information about format and content of the parameter see
    590344            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
     428def 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
     567def 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.
    591590    '''
    592591    parser = ArgumentParser(description='Plot retrieved GRIB data from ' + \
    … …  
    597596    parser.add_argument("--start_date", dest="start_date",
    598597                        help="start date YYYYMMDD")
    599     parser.add_argument( "--end_date", dest="end_date",
    600                          help="end_date YYYYMMDD")
     598    parser.add_argument("--end_date", dest="end_date",
     599                        help="end_date YYYYMMDD")
    601600
    602601    parser.add_argument("--start_step", dest="start_step",
    603602                        help="start step in hours")
    604     parser.add_argument( "--end_step", dest="end_step",
    605                          help="end step in hours")
     603    parser.add_argument("--end_step", dest="end_step",
     604                        help="end step in hours")
    606605
    607606# some arguments that override the default in the CONTROL file
    … …  
    635634    except IOError:
    636635        try:
    637             c = ControlFile(localpythonpath + args.controlfile)
    638 
    639         except:
     636            c = ControlFile(LOCAL_PYTHON_PATH + args.controlfile)
     637        except IOError:
    640638            print 'Could not read CONTROL file "' + args.controlfile + '"'
    641639            print 'Either it does not exist or its syntax is wrong.'
    … …  
    682680if __name__ == "__main__":
    683681    main()
    684 
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG