source: flex_extract.git/python/plot_retrieved.py @ 64cf353

ctbtodev
Last change on this file since 64cf353 was 64cf353, checked in by Anne Philipp <bscannephilipp@…>, 6 years ago

pep8 changes + documentation added + minor code style changes

  • Property mode set to 100644
File size: 10.5 KB
Line 
1#!/usr/bin/env python
2# This software is licensed under the terms of the Apache Licence Version 2.0
3# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
4#
5# Functionality provided: Simple tool for creating maps and time series of retrieved fields.
6#
7# Requirements:
8# in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed
9# ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/
10# dateutils
11# matplotlib (optional, for debugging)
12
13import datetime
14import time
15import os,inspect,sys,glob
16import socket
17# add path to submit.py to pythonpath so that python finds its buddies
18localpythonpath=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
19if localpythonpath not in sys.path:
20    sys.path.append(localpythonpath)
21
22from matplotlib.pylab import *
23import matplotlib.patches as mpatches
24from mpl_toolkits.basemap import Basemap,addcyclic
25import matplotlib.colors as mcolors
26from matplotlib.font_manager import FontProperties
27from matplotlib.patches import Polygon
28import matplotlib.cm as cmx
29import matplotlib.colors as colors
30from argparse import ArgumentParser,ArgumentDefaultsHelpFormatter
31
32from FlexpartTools import interpret_args_and_control,silentremove,product,Control
33from GribTools import GribTools
34from gribapi import *
35from rasotools.utils import stats
36
37def plot_retrieved(args,c):
38
39    start = datetime.datetime.strptime(c.start_date,'%Y%m%d%H')
40    end = datetime.datetime.strptime(c.end_date,'%Y%m%d%H')
41
42    c.paramIds=asarray(c.paramIds,dtype='int')
43    c.levels=asarray(c.levels,dtype='int')
44    c.area=asarray(c.area)
45
46    index_keys=["date","time","step"]
47    indexfile=c.inputdir+"/date_time_stepRange.idx"
48    silentremove(indexfile)
49    files=glob.glob(c.inputdir+'/'+c.prefix+'*')
50    grib=GribTools(files)
51    iid=grib.index(index_keys=index_keys, index_file = indexfile)
52
53    gdict=dict(Ni = 360, Nj = 181,  iScansNegatively = 0,  jScansPositively = 0,
54               jPointsAreConsecutive = 0,  alternativeRowScanning = 0,
55               latitudeOfFirstGridPointInDegrees = 90,
56               longitudeOfFirstGridPointInDegrees = 181,
57               latitudeOfLastGridPointInDegrees = -90,
58               longitudeOfLastGridPointInDegrees = 180,
59               iDirectionIncrementInDegrees = 1,
60               jDirectionIncrementInDegrees = 1
61               )
62
63    index_vals = []
64    for key in index_keys:
65        key_vals = grib_index_get(iid,key)
66        print key_vals
67
68        index_vals.append(key_vals)
69
70    fdict=dict()
71    fmeta=dict()
72    fstamp=dict()
73    for p in c.paramIds:
74        for l in c.levels:
75            key='{:0>3}_{:0>3}'.format(p,l)
76            fdict[key]=[]
77            fmeta[key]=[]
78            fstamp[key]=[]
79    for prod in product(*index_vals):
80        for i in range(len(index_keys)):
81            grib_index_select(iid,index_keys[i],prod[i])
82
83        gid = grib_new_from_index(iid)
84#        if gid is not None:
85
86        while(gid is not None):
87            date = grib_get(gid, 'date')
88            fdate= datetime.datetime(date/10000,mod(date,10000)/100,mod(date,100))
89            gtime = grib_get(gid, 'time')
90            step = grib_get(gid, 'step')
91            fdatetime=fdate+datetime.timedelta(hours=gtime/100)
92#            if fdatetime<start or fdatetime>end:
93#                break
94            gtype = grib_get(gid, 'type')
95            paramId = grib_get(gid, 'paramId')
96            parameterName = grib_get(gid, 'parameterName')
97            level=grib_get(gid,'level')
98            if step>=c.start_step and step <=c.end_step and fdatetime>=start and fdatetime<=end and paramId in c.paramIds and level in c.levels:
99                key='{:0>3}_{:0>3}'.format(paramId,level)
100                print key
101                fdatetimestep=fdatetime+datetime.timedelta(hours=step)
102                if len(fstamp)==0:
103                    fstamp[key].append(fdatetimestamp)
104                    fmeta[key].append((paramId,parameterName,gtype,fdatetime,gtime,step,level))
105                    fdict[key].append(flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']])))
106                else:
107                    i=0
108                    inserted=False
109                    for i in range(len(fstamp[key])):
110                        if fdatetimestep<fstamp[key][i]:
111                            fstamp[key][i:i]=[fdatetimestep]
112                            fmeta[key][i:i]=[(paramId,parameterName,gtype,fdatetime,gtime,step,level)]
113                            fdict[key][i:i]=[flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']]))]
114                            inserted=True
115                            break
116                    if not inserted:
117                        fstamp[key].append(fdatetimestep)
118                        fmeta[key].append((paramId,parameterName,gtype,fdatetime,gtime,step,level))
119                        fdict[key].append(flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']])))
120
121
122            grib_release(gid)
123            gid = grib_new_from_index(iid)
124
125    for k in fdict.keys():
126        fml=fmeta[k]
127        fdl=fdict[k]
128
129        for fd,fm in zip(fdl,fml):
130            ftitle=fm[1]+' {} '.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+' '+stats(fd)
131            pname='_'.join(fm[1].split())+'_{}_'.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+'.{:0>3}'.format(fm[5])
132            plotmap(fd, fm,gdict,ftitle,pname+'.eps')
133
134    for k in fdict.keys():
135        fml=fmeta[k]
136        fdl=fdict[k]
137        fsl=fstamp[k]
138        if fdl:
139            fm=fml[0]
140            fd=fdl[0]
141            ftitle=fm[1]+' {} '.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+' '+stats(fd)
142            pname='_'.join(fm[1].split())+'_{}_'.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+'.{:0>3}'.format(fm[5])
143            lat=-20
144            lon=20
145            plottimeseries(fdl,fml,fsl,lat,lon, gdict, ftitle, pname+'.eps')
146
147def plottimeseries(flist,fmetalist,ftimestamps,lat,lon,gdict,ftitle,filename):
148    t1=time.time()
149    latindex=(lat+90)*180/(gdict['Nj']-1)
150    lonindex=(lon+179)*360/gdict['Ni']
151    farr=asarray(flist)
152    ts=farr[:,latindex,lonindex]
153    f=plt.figure(figsize=(12,6.7))
154    plt.plot(ftimestamps,ts)
155    plt.title(ftitle)
156    savefig(c.outputdir+'/'+filename)
157    print 'created ',c.outputdir+'/'+filename
158    plt.close(f)
159    print time.time()-t1,'s'
160
161def plotmap(flist,fmetalist,gdict,ftitle,filename):
162    t1=time.time()
163    f=plt.figure(figsize=(12,6.7))
164    mbaxes = f.add_axes([0.05, 0.15, 0.8, 0.7])
165    m =Basemap(llcrnrlon=-180.,llcrnrlat=-90.,urcrnrlon=180,urcrnrlat=90.)
166    #if bw==0 :
167        #fill_color=rgb(0.6,0.8,1)
168    #else:
169        #fill_color=rgb(0.85,0.85,0.85)
170
171    lw=0.3
172    m.drawmapboundary()
173    parallels = arange(-90.,91,90.)
174    # labels = [left,right,top,bottom]
175    m.drawparallels(parallels,labels=[True,True,True,True],linewidth=lw)
176    meridians = arange(-180.,181.,60.)
177    m.drawmeridians(meridians,labels=[True,True,True,True],linewidth=lw)
178    m.drawcoastlines(linewidth=lw)
179    xleft=gdict['longitudeOfFirstGridPointInDegrees']
180    if xleft>180.0:
181        xleft-=360.
182    x=linspace(xleft,gdict['longitudeOfLastGridPointInDegrees'],gdict['Ni'])
183    y=linspace(gdict['latitudeOfLastGridPointInDegrees'],gdict['latitudeOfFirstGridPointInDegrees'],gdict['Nj'])
184    xx, yy = m(*meshgrid(x,y))
185
186    s=m.contourf(xx,yy, flist)
187    title(ftitle,y=1.1)
188    cbaxes = f.add_axes([0.9, 0.2, 0.04, 0.6])
189    cb=colorbar(cax=cbaxes)
190
191    savefig(c.outputdir+'/'+filename)
192    print 'created ',c.outputdir+'/'+filename
193    plt.close(f)
194    print time.time()-t1,'s'
195
196
197def interpret_plotargs(*args,**kwargs):
198
199    parser = ArgumentParser(description='Retrieve FLEXPART input from ECMWF MARS archive',
200                            formatter_class=ArgumentDefaultsHelpFormatter)
201
202# the most important arguments
203    parser.add_argument("--start_date", dest="start_date",
204                        help="start date YYYYMMDD")
205    parser.add_argument( "--end_date", dest="end_date",
206                         help="end_date YYYYMMDD")
207
208    parser.add_argument("--start_step", dest="start_step",
209                        help="start step in hours")
210    parser.add_argument( "--end_step", dest="end_step",
211                         help="end_step in hours")
212
213# some arguments that override the default in the control file
214    parser.add_argument("--levelist", dest="levelist",help="Vertical levels to be retrieved, e.g. 30/to/60")
215    parser.add_argument("--area", dest="area", help="area defined as north/west/south/east")
216    parser.add_argument("--paramIds", dest="paramIds", help="parameter IDs")
217    parser.add_argument("--prefix", dest="prefix",default='EN', help="output file name prefix")
218
219# set the working directories
220    parser.add_argument("--inputdir", dest="inputdir",default=None,
221                        help="root directory for storing intermediate files")
222    parser.add_argument("--outputdir", dest="outputdir",default=None,
223                        help="root directory for storing output files")
224    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
225                        help="FLEXPART root directory (to find grib2flexpart and COMMAND file)\n\
226                              Normally ECMWFDATA resides in the scripts directory of the FLEXPART distribution")
227
228    parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp',
229                        help="file with control parameters")
230    args = parser.parse_args()
231
232    try:
233        c=Control(args.controlfile)
234    except IOError:
235        try:
236            c=Control(localpythonpath+args.controlfile)
237
238        except:
239            print 'Could not read control file "'+args.controlfile+'"'
240            print 'Either it does not exist or its syntax is wrong.'
241            print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information'
242            exit(1)
243
244    if args.levelist:
245        c.levels=args.levelist.split('/')
246    else:
247        c.levels=[0]
248    if args.area:
249        c.area=args.area.split('/')
250    else:
251        c.area='[0,0]'
252
253    c.paramIds=args.paramIds.split('/')
254    if args.start_step:
255        c.start_step=int(args.start_step)
256    else:
257        c.start_step=0
258    if args.end_step:
259        c.end_step=int(args.end_step)
260    else:
261        c.end_step=0
262
263    c.start_date=args.start_date
264    c.end_date=args.end_date
265    c.prefix=args.prefix
266    c.inputdir=args.inputdir
267    if args.outputdir:
268        c.outputdir=args.outputdir
269    else:
270        c.outputdir=c.inputdir
271
272        return args,c
273
274if __name__ == "__main__":
275
276    args,c=interpret_plotargs()
277    plot_retrieved(args,c)
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG