source: flex_extract.git/python/plot_retrieved.py @ 991df6a

ctbtodev
Last change on this file since 991df6a was 991df6a, checked in by Anne Philipp <anne.philipp@…>, 6 years ago

finished documentation (except plot_retrieved)

  • Property mode set to 100755
File size: 15.6 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3#************************************************************************
4# TODO AP
5# - documentation der Funktionen
6# - docu der progam functionality
7# - apply pep8
8#************************************************************************
9#*******************************************************************************
10# @Author: Leopold Haimberger (University of Vienna)
11#
12# @Date: November 2015
13#
14# @Change History:
15#
16#    February 2018 - Anne Philipp (University of Vienna):
17#        - applied PEP8 style guide
18#        - added documentation
19#        - created function main and moved the two function calls for
20#          arguments and plotting into it
21#
22# @License:
23#    (C) Copyright 2015-2018.
24#
25#    This software is licensed under the terms of the Apache Licence Version 2.0
26#    which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
27#
28# @Program Functionality:
29#    Simple tool for creating maps and time series of retrieved fields.
30#
31# @Program Content:
32#    - plot_retrieved
33#    - plottimeseries
34#    - plotmap
35#    - interpret_plotargs
36#
37#*******************************************************************************
38
39# ------------------------------------------------------------------------------
40# MODULES
41# ------------------------------------------------------------------------------
42import datetime
43import time
44import os
45import inspect
46import sys
47import glob
48from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
49
50from matplotlib.pylab import *
51import matplotlib.patches as mpatches
52from mpl_toolkits.basemap import Basemap, addcyclic
53import matplotlib.colors as mcolors
54from matplotlib.font_manager import FontProperties
55from matplotlib.patches import Polygon
56import matplotlib.cm as cmx
57import matplotlib.colors as colors
58#from rasotools.utils import stats
59from gribapi import *
60
61# add path to pythonpath so that python finds its buddies
62localpythonpath = os.path.dirname(os.path.abspath(
63    inspect.getfile(inspect.currentframe())))
64if localpythonpath not in sys.path:
65    sys.path.append(localpythonpath)
66
67# software specific classes and modules from flex_extract
68from Tools import silentremove, product
69from ControlFile import ControlFile
70from GribTools import GribTools
71
72# ------------------------------------------------------------------------------
73# FUNCTION
74# ------------------------------------------------------------------------------
75def main():
76    '''
77    @Description:
78        If plot_retrieved is called from command line, this function controls
79        the program flow and calls the argumentparser function and
80        the plot_retrieved function for plotting the retrieved GRIB data.
81
82    @Input:
83        <nothing>
84
85    @Return:
86        <nothing>
87    '''
88    args, c = interpret_plotargs()
89    plot_retrieved(args, c)
90
91    return
92
93def plot_retrieved(args, c):
94    '''
95    @Description:
96        Reads GRIB data from a specified time period, a list of levels
97        and a specified list of parameter.
98
99    @Input:
100        args: instance of ArgumentParser
101            Contains the commandline arguments from script/program call.
102
103        c: instance of class ControlFile
104            Contains all necessary information of a CONTROL file. The parameters
105            are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
106            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
107            RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
108            SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
109            ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR,
110            OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
111            For more information about format and content of the parameter see
112            documentation.
113
114    @Return:
115        <nothing>
116    '''
117    start = datetime.datetime.strptime(c.start_date, '%Y%m%d%H')
118    end = datetime.datetime.strptime(c.end_date, '%Y%m%d%H')
119
120    c.paramIds = asarray(c.paramIds, dtype='int')
121    c.levels = asarray(c.levels, dtype='int')
122    c.area = asarray(c.area)
123
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)
148
149    fdict = dict()
150    fmeta = dict()
151    fstamp = dict()
152    for p in c.paramIds:
153        for l in c.levels:
154            key = '{:0>3}_{:0>3}'.format(p, l)
155            fdict[key] = []
156            fmeta[key] = []
157            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
164        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:
178                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)
183                    fmeta[key].append((paramId, parameterName, gtype,
184                                       fdatetime, gtime, step, level))
185                    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)
209
210    for k in fdict.keys():
211        fml = fmeta[k]
212        fdl = fdict[k]
213
214        for fd, fm in zip(fdl, fml):
215            ftitle = fm[1] + ' {} '.format(fm[-1]) + \
216                datetime.datetime.strftime(fm[3], '%Y%m%d%H') #+ ' ' + stats(fd)
217            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')
221
222    for k in fdict.keys():
223        fml = fmeta[k]
224        fdl = fdict[k]
225        fsl = fstamp[k]
226        if fdl:
227            fm = fml[0]
228            fd = fdl[0]
229            ftitle = fm[1] + ' {} '.format(fm[-1]) + \
230                datetime.datetime.strftime(fm[3], '%Y%m%d%H') #+ ' ' + stats(fd)
231            pname = '_'.join(fm[1].split()) + '_{}_'.format(fm[-1]) + \
232                datetime.datetime.strftime(fm[3], '%Y%m%d%H') + \
233                '.{:0>3}'.format(fm[5])
234            lat = -20
235            lon = 20
236            plottimeseries(fdl, fml, fsl, lat, lon,
237                           gdict, ftitle, pname + '.eps')
238
239    return
240
241def plottimeseries(flist, fmetalist, ftimestamps, lat, lon,
242                   gdict, ftitle, filename):
243    '''
244    @Description:
245
246    @Input:
247        flist:
248            The actual data values to be plotted from the grib messages.
249
250        fmetalist: list of strings
251            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
254
255        ftimestamps: list of datetime
256            Contains the time stamps in a datetime format, e.g.
257
258        lat:
259
260        lon:
261
262        gdict:
263
264        ftitle: string
265            The title of the timeseries.
266
267        filename: string
268            The time series is stored in a file with this name.
269
270    @Return:
271        <nothing>
272    '''
273    t1 = time.time()
274    latindex = (lat + 90) * 180 / (gdict['Nj'] - 1)
275    lonindex = (lon + 179) * 360 / gdict['Ni']
276    farr = asarray(flist)
277    ts = farr[:, latindex, lonindex]
278    f = plt.figure(figsize=(12,6.7))
279    plt.plot(ftimestamps, ts)
280    plt.title(ftitle)
281    savefig(c.outputdir + '/' + filename)
282    print 'created ', c.outputdir + '/' + filename
283    plt.close(f)
284    print time.time() - t1, 's'
285
286    return
287
288def plotmap(flist, fmetalist, gdict, ftitle, filename):
289    '''
290    @Description:
291
292    @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
339def 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
352        c: instance of class ControlFile
353            Contains all necessary information of a CONTROL file. The parameters
354            are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM,
355            NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST,
356            RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA,
357            SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS,
358            ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, DEBUG, INPUTDIR,
359            OUTPUTDIR, FLEXPART_ROOT_SCRIPTS
360            For more information about format and content of the parameter see
361            documentation.
362    '''
363    parser = ArgumentParser(description='Retrieve FLEXPART input from ' + \
364                            'ECMWF MARS archive',
365                            formatter_class=ArgumentDefaultsHelpFormatter)
366
367# the most important arguments
368    parser.add_argument("--start_date", dest="start_date",
369                        help="start date YYYYMMDD")
370    parser.add_argument( "--end_date", dest="end_date",
371                         help="end_date YYYYMMDD")
372
373    parser.add_argument("--start_step", dest="start_step",
374                        help="start step in hours")
375    parser.add_argument( "--end_step", dest="end_step",
376                         help="end_step in hours")
377
378# some arguments that override the default in the CONTROL file
379    parser.add_argument("--levelist", dest="levelist",
380                        help="Vertical levels to be retrieved, e.g. 30/to/60")
381    parser.add_argument("--area", dest="area",
382                        help="area defined as north/west/south/east")
383    parser.add_argument("--paramIds", dest="paramIds",
384                        help="parameter IDs")
385    parser.add_argument("--prefix", dest="prefix", default='EN',
386                        help="output file name prefix")
387
388# set the working directories
389    parser.add_argument("--inputdir", dest="inputdir", default=None,
390                        help="root directory for storing intermediate files")
391    parser.add_argument("--outputdir", dest="outputdir", default=None,
392                        help="root directory for storing output files")
393    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
394                        help="FLEXPART root directory (to find \
395                        'grib2flexpart and COMMAND file)\n \
396                        Normally ECMWFDATA resides in the scripts directory \
397                        of the FLEXPART distribution")
398
399    parser.add_argument("--controlfile", dest="controlfile",
400                        default='CONTROL.temp', help="file with CONTROL parameters")
401    args = parser.parse_args()
402
403    try:
404        c = ControlFile(args.controlfile)
405    except IOError:
406        try:
407            c = ControlFile(localpythonpath + args.controlfile)
408
409        except:
410            print 'Could not read CONTROL file "' + args.controlfile + '"'
411            print 'Either it does not exist or its syntax is wrong.'
412            print 'Try "' + sys.argv[0].split('/')[-1] + \
413                  ' -h" to print usage information'
414            exit(1)
415
416    if args.levelist:
417        c.levels = args.levelist.split('/')
418    else:
419        c.levels = [0]
420    if args.area:
421        c.area = args.area.split('/')
422    else:
423        c.area = '[0,0]'
424
425    c.paramIds = args.paramIds.split('/')
426    if args.start_step:
427        c.start_step = int(args.start_step)
428    else:
429        c.start_step = 0
430    if args.end_step:
431        c.end_step = int(args.end_step)
432    else:
433        c.end_step = 0
434
435    c.start_date = args.start_date
436    c.end_date = args.end_date
437    c.prefix = args.prefix
438    c.inputdir = args.inputdir
439    if args.outputdir:
440        c.outputdir = args.outputdir
441    else:
442        c.outputdir = c.inputdir
443
444    return args, c
445
446if __name__ == "__main__":
447    main()
448
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG