Changeset ccab809 in flex_extract.git


Ignore:
Timestamp:
May 22, 2018, 10:54:43 PM (6 years ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
master, ctbto, dev
Children:
ff99eae
Parents:
812283d
Message:

further improved plotting routine for date selection

Location:
python
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • python/CONTROL.temp

    refdb01a rccab809  
    2929PREFIX EI
    3030ECSTORAGE 0
    31 ECTRANS 0
     31ECTRANS 1
    3232ECFSDIR ectmp:/${USER}/econdemand/
    3333MAILFAIL ${USER}
  • python/GribTools.py

    r991df6a rccab809  
    319319                    self.iid = grib_index_new_from_file(file, index_keys)
    320320                else:
     321                    print 'in else zweig'
    321322                    grib_index_add_file(self.iid, file)
    322323
  • python/job.ksh

    r991df6a rccab809  
    7070cwc 0
    7171date_chunk 3
    72 debug True
     72debug 1
    7373dpdeta 1
    7474dtime 3
    7575ecfsdir ectmp:/${USER}/econdemand/
    7676ecstorage 0
    77 ectrans 0
     77ectrans 1
    7878end_date 20160809
    7979eta 0
  • python/joboper.ksh

    r991df6a rccab809  
    7070cwc 0
    7171date_chunk 3
    72 debug True
     72debug 1
    7373dpdeta 1
    7474dtime 3
    7575ecfsdir ectmp:/${USER}/econdemand/
    7676ecstorage 0
    77 ectrans 0
     77ectrans 1
    7878start_date ${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}
    7979eta 0
  • python/plot_retrieved.py

    r991df6a rccab809  
    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
     22#          of the data fields from the first grib file it gets.
    2123#
    2224# @License:
     
    3032#
    3133# @Program Content:
    32 #    - plot_retrieved
    33 #    - plottimeseries
    34 #    - plotmap
    35 #    - interpret_plotargs
     34#    - main
     35#    - getBasics
     36#    - plotRetrieved
     37#    - plotTS
     38#    - plotMap
     39#    - getPlotArgs
    3640#
    3741#*******************************************************************************
     
    4044# MODULES
    4145# ------------------------------------------------------------------------------
     46import time
    4247import datetime
    43 import time
    4448import os
    4549import inspect
     
    5761import matplotlib.colors as colors
    5862#from rasotools.utils import stats
    59 from gribapi import *
     63
     64font = {'family': 'monospace', 'size': 12}
     65matplotlib.rcParams['xtick.major.pad'] = '20'
     66
     67matplotlib.rc('font', **font)
     68
     69from eccodes import *
     70import numpy as np
    6071
    6172# add path to pythonpath so that python finds its buddies
     
    6677
    6778# software specific classes and modules from flex_extract
    68 from Tools import silentremove, product
     79from Tools import silentremove
    6980from ControlFile import ControlFile
    70 from GribTools import GribTools
     81#from GribTools import GribTools
     82from UIOFiles import UIOFiles
    7183
    7284# ------------------------------------------------------------------------------
     
    7688    '''
    7789    @Description:
    78         If plot_retrieved is called from command line, this function controls
     90        If plotRetrieved is called from command line, this function controls
    7991        the program flow and calls the argumentparser function and
    80         the plot_retrieved function for plotting the retrieved GRIB data.
     92        the plotRetrieved function for plotting the retrieved GRIB data.
    8193
    8294    @Input:
     
    8698        <nothing>
    8799    '''
    88     args, c = interpret_plotargs()
    89     plot_retrieved(args, c)
     100    args, c = getPlotArgs()
     101    plotRetrieved(args, c)
    90102
    91103    return
    92104
    93 def plot_retrieved(args, c):
     105def getBasics(ifile, verb=False):
     106    """
     107    @Description:
     108        An example grib file will be opened and basic information will
     109        be extracted. These information are important for later use and the
     110        initialization of numpy arrays for data storing.
     111
     112    @Input:
     113        ifile: string
     114            Contains the full absolute path to the ECMWF grib file.
     115        verb (opt): bool
     116            Is True if there should be extra output in verbose mode.
     117            Default value is False.
     118
     119    @Return:
     120        data: dict
     121            Contains basic informations of the ECMWF grib files, e.g.
     122            'Ni', 'Nj', 'latitudeOfFirstGridPointInDegrees',
     123            'longitudeOfFirstGridPointInDegrees',
     124            'latitudeOfLastGridPointInDegrees',
     125            'longitudeOfLastGridPointInDegrees',
     126            'jDirectionIncrementInDegrees',
     127            'iDirectionIncrementInDegrees'
     128    """
     129    from eccodes import *
     130
     131    data = {}
     132
     133    # --- open file ---
     134    print("Opening file for getting information data --- %s" %
     135          os.path.basename(ifile))
     136
     137    with open(ifile) as f:
     138
     139        # load first message from file
     140        gid = codes_grib_new_from_file(f)
     141
     142        # information needed from grib message
     143        keys = [
     144                'Ni',
     145                'Nj',
     146                'latitudeOfFirstGridPointInDegrees',
     147                'longitudeOfFirstGridPointInDegrees',
     148                'latitudeOfLastGridPointInDegrees',
     149                'longitudeOfLastGridPointInDegrees',
     150                'jDirectionIncrementInDegrees',
     151                'iDirectionIncrementInDegrees'
     152               ]
     153
     154        if verb: print('\nInformations are: ')
     155        for key in keys:
     156            # 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'
     160
     161        # Free the memory for the message referred as gribid.
     162        codes_release(gid)
     163
     164    return data
     165
     166def getFilesPerDate(files, datelist):
     167    '''
     168    @Description:
     169        The filenames contain dates which are used to select a list
     170        of files for a specific time period specified in datelist.
     171
     172    @Input:
     173        files: instance of UIOFiles
     174            For description see class documentation.
     175            It contains the attribute "files" which is a list of pathes
     176            to filenames.
     177
     178        datelist: list of datetimes
     179            Contains the list of dates which should be processed for plotting.
     180
     181    @Return:
     182        filelist: list of strings
     183            Contains the selected files for the time period.
     184    '''
     185
     186    filelist = []
     187    for file in files:
     188        fdate = file[-8:]
     189        ddate = datetime.datetime.strptime(fdate, '%y%m%d%H')
     190        if ddate in datelist:
     191            filelist.append(file)
     192
     193    return filelist
     194
     195def plotRetrieved(args, c):
    94196    '''
    95197    @Description:
     
    118220    end = datetime.datetime.strptime(c.end_date, '%Y%m%d%H')
    119221
     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
    120231    c.paramIds = asarray(c.paramIds, dtype='int')
    121232    c.levels = asarray(c.levels, dtype='int')
    122233    c.area = asarray(c.area)
    123234
    124     index_keys = ["date", "time", "step"]
    125     indexfile = c.inputdir + "/date_time_stepRange.idx"
    126     silentremove(indexfile)
    127     files = glob.glob(c.inputdir + '/' + c.prefix + '*')
    128     grib = GribTools(files)
    129     iid = grib.index(index_keys=index_keys, index_file = indexfile)
    130 
    131     gdict = dict(Ni = 360, Nj = 181,
    132                 iScansNegatively = 0,  jScansPositively = 0,
    133                 jPointsAreConsecutive = 0,  alternativeRowScanning = 0,
    134                 latitudeOfFirstGridPointInDegrees = 90,
    135                 longitudeOfFirstGridPointInDegrees = 181,
    136                 latitudeOfLastGridPointInDegrees = -90,
    137                 longitudeOfLastGridPointInDegrees = 180,
    138                 iDirectionIncrementInDegrees = 1,
    139                 jDirectionIncrementInDegrees = 1
    140                 )
    141 
    142     index_vals = []
    143     for key in index_keys:
    144         key_vals = grib_index_get(iid, key)
    145         print key_vals
    146 
    147         index_vals.append(key_vals)
     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)
    148265
    149266    fdict = dict()
     
    156273            fmeta[key] = []
    157274            fstamp[key] = []
    158     for prod in product(*index_vals):
    159         for i in range(len(index_keys)):
    160             grib_index_select(iid, index_keys[i], prod[i])
    161 
    162         gid = grib_new_from_index(iid)
    163 
     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)
    164283        while(gid is not None):
    165             date = grib_get(gid, 'date')
    166             fdate = datetime.datetime(date/10000, mod(date,10000)/100,
    167                                       mod(date,100))
    168             gtime = grib_get(gid, 'time')
    169             step = grib_get(gid, 'step')
    170             fdatetime = fdate + datetime.timedelta(hours=gtime/100)
    171             gtype = grib_get(gid, 'type')
    172             paramId = grib_get(gid, 'paramId')
    173             parameterName = grib_get(gid, 'parameterName')
    174             level = grib_get(gid, 'level')
    175             if step >= c.start_step and step <= c.end_step and \
    176                fdatetime >= start and fdatetime <= end and \
    177                paramId in c.paramIds and level in c.levels:
     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:
    178290                key = '{:0>3}_{:0>3}'.format(paramId, level)
    179                 print key
    180                 fdatetimestep = fdatetime + datetime.timedelta(hours=step)
    181                 if len(fstamp) == 0:
    182                     fstamp[key].append(fdatetimestamp)
     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)
    183323                    fmeta[key].append((paramId, parameterName, gtype,
    184                                        fdatetime, gtime, step, level))
     324                                       fdate, level))
    185325                    fdict[key].append(flipud(reshape(
    186                             grib_get_values(gid), [gdict['Nj'], gdict['Ni']])))
    187                 else:
    188                     i = 0
    189                     inserted = False
    190                     for i in range(len(fstamp[key])):
    191                         if fdatetimestep < fstamp[key][i]:
    192                             fstamp[key][i:i] = [fdatetimestep]
    193                             fmeta[key][i:i] = [(paramId, parameterName, gtype,
    194                                                 fdatetime, gtime, step, level)]
    195                             fdict[key][i:i] = [flipud(reshape(
    196                                                 grib_get_values(gid),
    197                                                 [gdict['Nj'], gdict['Ni']]))]
    198                             inserted = True
    199                             break
    200                     if not inserted:
    201                         fstamp[key].append(fdatetimestep)
    202                         fmeta[key].append((paramId, parameterName, gtype,
    203                                            fdatetime, gtime, step, level))
    204                         fdict[key].append(flipud(reshape(
    205                             grib_get_values(gid), [gdict['Nj'], gdict['Ni']])))
    206 
    207             grib_release(gid)
    208             gid = grib_new_from_index(iid)
     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()
    209336
    210337    for k in fdict.keys():
     338        print 'fmeta: ', len(fmeta),fmeta
    211339        fml = fmeta[k]
    212340        fdl = fdict[k]
    213 
     341        print 'fm1: ', len(fml), fml
     342        #print 'fd1: ', fdl
     343        #print zip(fdl, fml)
    214344        for fd, fm in zip(fdl, fml):
     345            print fm
    215346            ftitle = fm[1] + ' {} '.format(fm[-1]) + \
    216347                datetime.datetime.strftime(fm[3], '%Y%m%d%H') #+ ' ' + stats(fd)
    217348            pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \
    218                 datetime.datetime.strftime(fm[3], '%Y%m%d%H') + \
    219                 '.{:0>3}'.format(fm[5])
    220             plotmap(fd, fm, gdict, ftitle, pname + '.eps')
     349                datetime.datetime.strftime(fm[3], '%Y%m%d%H')
     350            plotMap(c, fd, fm, gdict, ftitle, pname, 'png')
    221351
    222352    for k in fdict.keys():
     
    230360                datetime.datetime.strftime(fm[3], '%Y%m%d%H') #+ ' ' + stats(fd)
    231361            pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \
    232                 datetime.datetime.strftime(fm[3], '%Y%m%d%H') + \
    233                 '.{:0>3}'.format(fm[5])
     362                datetime.datetime.strftime(fm[3], '%Y%m%d%H')
    234363            lat = -20
    235364            lon = 20
    236             plottimeseries(fdl, fml, fsl, lat, lon,
    237                            gdict, ftitle, pname + '.eps')
     365            plotTS(c, fdl, fml, fsl, lat, lon,
     366                           gdict, ftitle, pname, 'png')
    238367
    239368    return
    240369
    241 def plottimeseries(flist, fmetalist, ftimestamps, lat, lon,
    242                    gdict, ftitle, filename):
     370def plotTS(c, flist, fmetalist, ftimestamps, lat, lon,
     371                   gdict, ftitle, filename, fending, show=False):
    243372    '''
    244373    @Description:
    245374
    246375    @Input:
     376        c:
     377
    247378        flist:
    248379            The actual data values to be plotted from the grib messages.
     
    250381        fmetalist: list of strings
    251382            Contains some meta date for the data field to be plotted:
    252             parameter id, parameter Name, grid type, date and time,
    253             time, forecast step, level
     383            parameter id, parameter Name, grid type, datetime, level
    254384
    255385        ftimestamps: list of datetime
    256             Contains the time stamps in a datetime format, e.g.
     386            Contains the time stamps.
    257387
    258388        lat:
     
    268398            The time series is stored in a file with this name.
    269399
     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
    270406    @Return:
    271407        <nothing>
    272408    '''
     409    print 'plotting timeseries'
     410
    273411    t1 = time.time()
    274     latindex = (lat + 90) * 180 / (gdict['Nj'] - 1)
    275     lonindex = (lon + 179) * 360 / gdict['Ni']
     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
    276424    farr = asarray(flist)
    277     ts = farr[:, latindex, lonindex]
    278     f = plt.figure(figsize=(12,6.7))
     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
    279440    plt.plot(ftimestamps, ts)
    280441    plt.title(ftitle)
    281     savefig(c.outputdir + '/' + filename)
     442
     443    plt.savefig(c.outputdir+'/'+filename+'_TS.'+fending, facecolor=fig.get_facecolor(), edgecolor='none',format=fending)
    282444    print 'created ', c.outputdir + '/' + filename
    283     plt.close(f)
     445    if show == True:
     446        plt.show()
     447    fig.clf()
     448    plt.close(fig)
     449
    284450    print time.time() - t1, 's'
    285451
    286452    return
    287453
    288 def plotmap(flist, fmetalist, gdict, ftitle, filename):
     454def plotMap(c, flist, fmetalist, gdict, ftitle, filename, fending, show=False):
    289455    '''
    290456    @Description:
    291457
    292458    @Input:
    293         flist
    294         fmetalist
    295         gdict
    296         ftitle
    297         filename
    298 
    299     @Return:
    300         <nothing>
    301     '''
    302     t1 = time.time()
    303     f = plt.figure(figsize=(12, 6.7))
    304     mbaxes = f.add_axes([0.05, 0.15, 0.8, 0.7])
    305     m = Basemap(llcrnrlon=-180., llcrnrlat=-90., urcrnrlon=180, urcrnrlat=90.)
    306     #if bw==0 :
    307         #fill_color=rgb(0.6,0.8,1)
    308     #else:
    309         #fill_color=rgb(0.85,0.85,0.85)
    310 
    311     lw = 0.3
    312     m.drawmapboundary()
    313     parallels = arange(-90., 91, 90.)
    314     # labels = [left, right, top, bottom]
    315     m.drawparallels(parallels, labels=[True, True, True, True], linewidth=lw)
    316     meridians = arange(-180., 181., 60.)
    317     m.drawmeridians(meridians, labels=[True, True, True, True], linewidth=lw)
    318     m.drawcoastlines(linewidth=lw)
    319     xleft = gdict['longitudeOfFirstGridPointInDegrees']
    320     if xleft > 180.0:
    321         xleft -= 360.
    322     x = linspace(xleft, gdict['longitudeOfLastGridPointInDegrees'], gdict['Ni'])
    323     y = linspace(gdict['latitudeOfLastGridPointInDegrees'],
    324                  gdict['latitudeOfFirstGridPointInDegrees'], gdict['Nj'])
    325     xx, yy = m(*meshgrid(x, y))
    326 
    327     s = m.contourf(xx, yy, flist)
    328     title(ftitle, y=1.1)
    329     cbaxes = f.add_axes([0.9, 0.2, 0.04, 0.6])
    330     cb = colorbar(cax=cbaxes)
    331 
    332     savefig(c.outputdir + '/' + filename)
    333     print 'created ', c.outputdir + '/' + filename
    334     plt.close(f)
    335     print time.time() - t1, 's'
    336 
    337     return
    338 
    339 def interpret_plotargs():
    340     '''
    341     @Description:
    342         Assigns the command line arguments and reads CONTROL file
    343         content. Apply default values for non mentioned arguments.
    344 
    345     @Input:
    346         <nothing>
    347 
    348     @Return:
    349         args: instance of ArgumentParser
    350             Contains the commandline arguments from script/program call.
    351 
    352459        c: instance of class ControlFile
    353460            Contains all necessary information of a CONTROL file. The parameters
     
    360467            For more information about format and content of the parameter see
    361468            documentation.
    362     '''
    363     parser = ArgumentParser(description='Retrieve FLEXPART input from ' + \
     469
     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
     495    @Return:
     496        <nothing>
     497    '''
     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'
     565
     566    return
     567
     568def getPlotArgs():
     569    '''
     570    @Description:
     571        Assigns the command line arguments and reads CONTROL file
     572        content. Apply default values for non mentioned arguments.
     573
     574    @Input:
     575        <nothing>
     576
     577    @Return:
     578        args: instance of ArgumentParser
     579            Contains the commandline arguments from script/program call.
     580
     581        c: instance of class ControlFile
     582            Contains all necessary information of a CONTROL file. The parameters
     583            are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
     584            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
     585            RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
     586            SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
     587            ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR,
     588            OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
     589            For more information about format and content of the parameter see
     590            documentation.
     591    '''
     592    parser = ArgumentParser(description='Plot retrieved GRIB data from ' + \
    364593                            'ECMWF MARS archive',
    365594                            formatter_class=ArgumentDefaultsHelpFormatter)
     
    374603                        help="start step in hours")
    375604    parser.add_argument( "--end_step", dest="end_step",
    376                          help="end_step in hours")
     605                         help="end step in hours")
    377606
    378607# some arguments that override the default in the CONTROL file
    379608    parser.add_argument("--levelist", dest="levelist",
    380                         help="Vertical levels to be retrieved, e.g. 30/to/60")
     609                        help="vertical levels to be retrieved, e.g. 30/to/60")
    381610    parser.add_argument("--area", dest="area",
    382611                        help="area defined as north/west/south/east")
     
    398627
    399628    parser.add_argument("--controlfile", dest="controlfile",
    400                         default='CONTROL.temp', help="file with CONTROL parameters")
     629                        default='CONTROL.temp',
     630                        help="file with CONTROL parameters")
    401631    args = parser.parse_args()
    402632
     
    418648    else:
    419649        c.levels = [0]
     650
    420651    if args.area:
    421652        c.area = args.area.split('/')
     
    424655
    425656    c.paramIds = args.paramIds.split('/')
     657
    426658    if args.start_step:
    427659        c.start_step = int(args.start_step)
    428660    else:
    429661        c.start_step = 0
     662
    430663    if args.end_step:
    431664        c.end_step = int(args.end_step)
     
    435668    c.start_date = args.start_date
    436669    c.end_date = args.end_date
     670
    437671    c.prefix = args.prefix
     672
    438673    c.inputdir = args.inputdir
     674
    439675    if args.outputdir:
    440676        c.outputdir = args.outputdir
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG