source: lmdz_wrf/trunk/tools/casefigures.py @ 2607

Last change on this file since 2607 was 1463, checked in by lfita, 9 years ago

Creation of:

Plotting case figures tacking GRIB files as reference
File size: 10.6 KB
Line 
1# Plotting case figures tacking ERA-Interim files as reference
2import os
3import numpy as np
4import generic_tools as gen
5import nc_var_tools as ncvar
6import drawing_tools as drw
7import subprocess as sub
8
9main = 'casefigures.py'
10
11# Starting files from scratch
12filescratch = False
13
14# Starting figures from scratch
15figscratch = False
16
17# PyNCplot folder
18pyHOME = '/home/lluis/PyNCplot'
19
20# CDO folder
21cdoHOME = '/usr/bin/cdo'
22
23# GRIB folder
24gribfold='/media/lluis/ExtDiskC_ext4/DATA/ECMWF/ERA-Interim/'
25
26# head, middle and tail of grib name files
27gribheadf = 'ERAI_'
28gribmidf = '200512'
29gribtailf = '.grib'
30
31# Does the the code of the variable appears in the name of the grib file
32gribcodeInfile = True
33
34# GRIB table
35gribtable = '128ECMWF'
36
37# Folder where to keep generated netcdf files and figures
38outfold = 'out'
39
40# Dates to plot
41idate='20051212120000'
42edate='20051216120000'
43
44# Slice netcdf
45slicevals = 'lev|0,lat|-1,lon|-1'
46
47# Map values
48mapvalues='cyl,l'
49
50# kind of figure output
51kindfig = 'png'
52
53# Reverse of matrices in figure
54remap = 'flip@y'
55
56# Map range as ['SWlon','SWlat','NElon','NElat']
57maprange=[-10., 25., 40., 55.]
58
59# Axis values
60# [dimxyfmt]=[dxs],[dxf],[Ndx],[ordx],[dys],[dyf],[Ndy],[ordx]
61dimxyfmt = ' auto'
62
63# Characteristics variables to make the files
64#  [CFvarn]@[level (hpa)] level='sfc' for surface values
65varlevs = ['z|100000', 'ta|100000', 'ua|100000', 'va|100000', 'hur|100000',          \
66  'z|50000', 'ta|50000', 'ua|50000', 'va|50000', 'hur|50000', 'z|30000', 'ta|30000', \
67  'ua|30000', 'va|30000', 'hur|30000', 'ps|sfc', 'psl|sfc']
68
69# Figures
70#   shdcont: shadow-contour plots
71#     [shadvnl]@[contvnl]@[minshad],[maxshad]@[mincnt],[maxcnt],[Nlev]@[colbar]@[cntype]
72#   shd2cont: shadow, 2-contour plots
73#     [shadvnl]@[contvnl1]@[contvnl2]@[minshad],[maxshad]@[mincnt1],[maxcnt1],[Nlev1]@[mincnt2],[maxcnt2],[Nlev2]@[colbar]@[cntype1]@[cntype2]
74#   shdvec: shadow-vector plots
75#     [shadvnl]@[xvecnl],[yvecnl]@[minshad],[maxshad]@[colbar]@[vectype]
76#   shdcontvec: shadow-contour-vector plots
77#     [shadvnl]@[contvnl]@[xvecnl],[yvecnl]@[minshad],[maxshad]@[mincnt],[maxcnt],[Nlev]@[colbar]@[cntype]@[vectype]
78#   shdcontbarb: shadow-contour-barb plots
79#     [shadvnl]@[contvnl]@[xnbarbnl],[ybarbnl]@[minshad],[maxshad]@[mincnt],[maxcnt],[Nlev]@[colbar]@[cntype]@[barbtype]
80
81#### ### ## #
82#    [XXXvnl]: [vn]|[lev], [vn] CF-variable name at level [lev]
83#    [minXXX],[maxXXX]: minimun and maximum values to plot
84#    [colbar]: name of the colormap to use
85#    [cntype]=[ckind]|[clabfmt] values for the contour plot
86#      ckind]: kind of contours
87#        'cmap': as it gets from colorbar
88#        'fixc,[colname]': fixed color [colname], all stright lines
89#        'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
90#      clabfmt: format of the labels in the contour (None, also possible)
91#    [vectype]= [frequency],[color],[length],[windname],[windunits]  values for the vector plot
92#    [barbtype]: values for the barb plot
93
94shd2cont = ['hur|100000@z|100000@ta|100000@60.,100.@-700.,4000.,8@265.,300.,8@BuPu@fixc,green|None@fixc,red|None',
95'hur|100000@z|30000@psl|sfc@60.,100.@82000.,93000.,8@99000.,104000.,8@BuPu@fixc,green|None@fixc,red|None']
96#shdcontbarb = [[shadvnl]@[contvnl]@[xnbarbnl],[ybarbnl]@[minshad],[maxshad]@[mincnt],[maxcnt],[Nlev]@[colbar]@[cntype]@[barbtype]]
97
98availfigk = ['shad2cont']
99
100# dictionary of figures
101figures = {}
102figures['shad2cont'] = shd2cont
103
104####### ###### ##### #### ### ## #
105
106
107#######    #######
108## MAIN
109    #######
110
111cdolonlatS = str(maprange[0]) + ',' + str(maprange[2]) + ',' + str(maprange[1]) +    \
112  ',' + str(maprange[3])
113
114cdoidate = gen.datetimeStr_conversion(idate, 'YmdHMS', 'Y-m-d_H:M:S').replace('_','T')
115cdoedate = gen.datetimeStr_conversion(edate, 'YmdHMS', 'Y-m-d_H:M:S').replace('_','T')
116
117# Files in folder
118gribfiles = gen.files_folder_HMT(folder=gribfold, head=gribheadf, middle=gribmidf,   \
119  tail=gribtailf)
120print 'grib files to use:', gribfiles
121
122# Mapping variables in file
123gribvarfiles = {}
124for gribf in gribfiles:
125    ins = sub.Popen([cdoHOME,"showcode",gribfold + '/' + gribf], stdout=sub.PIPE)
126    codes = ins.communicate()
127    codels = codes[0]
128    codelist = codels.replace('\n','').split(' ')
129    while gen.searchInlist(codelist,''): codelist.remove('')
130    gribvarfiles[gribf] = list(np.array(codelist, dtype=int))
131
132# creation of netcdf variable & level file from grib
133ncfiles = []
134# Time ranges for each file
135dimts = {}
136for vln in varlevs:
137    varn = vln.split('|')[0]
138    varlev = vln.split('|')[1]
139    varcode = gen.CF_gribequiv(varn, gribtable)[0]
140    varGRIBname = gen.CF_gribequiv(varn, gribtable)[1]
141
142    varcodeS = str(varcode)
143
144    ofilen = outfold + '/' + varn + '_' + varcodeS + '_' + varlev + '_' + idate +    \
145      '-' + edate + '.nc'
146    ncfiles.append(ofilen)
147
148    if filescratch: 
149        ins = 'rm -f ' + ofilen #+ ' >& /dev/null'
150        print 'ins:', ins
151        sub.call(ins, shell=True)
152
153    if not os.path.isfile(ofilen):
154        # Looking in which file there is the variable
155        for gribfilen in gribvarfiles.keys():
156            codesfile = gribvarfiles[gribfilen]
157            if gen.searchInlist(codesfile, varcode):
158                foundcode = True
159                if varlev != 'sfc':
160                    ins = cdoHOME + ' -f nc selcode,' + varcodeS + ' -sellevel,' +   \
161                      varlev + ' -sellonlatbox,' + cdolonlatS + ' -seldate,' +       \
162                      cdoidate + ',' + cdoedate
163                else:
164                    ins = cdoHOME + ' -f nc selcode,' + varcodeS + ' -sellonlatbox,'+\
165                      cdolonlatS + ' -seldate,' + cdoidate + ',' + cdoedate
166                try:
167                    with gen.Capturing() as output:
168                        sout = sub.call(ins + ' ' + gribfold+'/' + gribfilen + ' ' + \
169                          ofilen, shell=True)
170                except:
171                    print gen.errormsg
172                    print ins+' '+gribfold+'/' + gribfilen + ' ' + ofilen
173                    print sout
174                    for s1out in output: print s1out
175                    quit(-1)
176
177                ncvar.chvarname(varn, ofilen, varGRIBname)
178                print "    created file '" + ofilen + "'"
179                break
180
181        if not foundcode:
182            print gen.errormsg
183            print "  variable: '"+varn+"' with code:",varcode, " not found !!"
184            quit(-1)
185
186# Checking time consitency
187for ofilen in ncfiles:
188    ins = sub.Popen([cdoHOME,"showtimestamp",ofilen], stdout=sub.PIPE)
189    dts = ins.communicate()
190    Ts = dts[0]
191    dt = Ts.replace('\n','').split(' ')
192    while gen.searchInlist(dt,''): 
193        dt.remove('')
194    dimts[ofilen] = dt
195
196Nncfiles = len(ncfiles)
197dt0 = len(dimts[ncfiles[0]])
198for ifn in range(Nncfiles-1):
199    gfilen = ncfiles[ifn+1]
200    dt1 = len(dimts[gfilen])
201    if dt1 != dt0:
202        print gen.warnmsg
203        print "  file '" + gfilen + "' has different time-steps:", dt1, " than '" +  \
204          ncfiles[0] + "' with:", dt0
205
206# Plotting
207for figk in figures.keys():
208    print figk
209    figplot = figures[figk]
210
211    for figp in figplot:
212        print "  '" + figp + "' ..."
213        figvalues = figp.split('@')
214
215        # shad2cont
216        if figk == 'shad2cont':
217            shdvarn = figvalues[0].split('|')[0]
218            shdlev = figvalues[0].split('|')[1]
219            cnt1varn = figvalues[1].split('|')[0]
220            cnt1lev = figvalues[1].split('|')[1]
221            cnt2varn = figvalues[2].split('|')[0]
222            cnt2lev = figvalues[2].split('|')[1]
223            shdrange = figvalues[3]
224            cnt1range = figvalues[4]
225            cnt2range = figvalues[5]
226            shdcolbar = figvalues[6]
227            cnt1type = figvalues[7].replace('|',':')
228            cnt2type = figvalues[8].replace('|',':')
229           
230            shdcode = gen.CF_gribequiv(shdvarn, gribtable)[0]
231            cnt1code = gen.CF_gribequiv(cnt1varn, gribtable)[0]
232            cnt2code = gen.CF_gribequiv(cnt2varn, gribtable)[0]
233
234            shdfile= outfold + '/' + shdvarn + '_' + str(shdcode) + '_' + shdlev +   \
235              '_' + idate+'-'+edate+'.nc'
236            cnt1file= outfold + '/' + cnt1varn + '_' + str(cnt1code) + '_' + cnt1lev+\
237              '_' + idate+'-'+edate+'.nc'
238            cnt2file= outfold + '/' + cnt2varn + '_' + str(cnt2code) + '_' + cnt2lev+\
239               '_' + idate+'-'+edate+'.nc'
240
241            ncfile= shdfile + ',' + cnt1file + ',' + cnt2file
242            varfignames = shdvarn + '@' + shdlev + ',' + cnt1varn + '@' + cnt1lev +  \
243              ',' +  cnt2varn + '@' + cnt2lev
244            varnames = shdvarn + ',' + cnt1varn + ',' +  cnt2varn
245
246            timevals = dimts[shdfile]
247            for it in range(len(timevals)):
248                slicev = 'time|' + str(it) + ',' + slicevals
249                slicesv = slicev + ':' + slicev + ':' + slicev
250                timev = timevals[it].replace(':','').replace('-','').replace('T','')
251 
252                if shdlev == cnt1lev == cnt2lev:
253                   titlefig= varnames.replace(',','|')+'|@'+ shdlev + '|on|' + timev
254                else:
255                   titlefig = varfignames.replace(',', '|') + '|on|' + timev
256
257                values = varfignames + ':' + slicesv + ':lon:lat:auto:' + shdcolbar+ \
258                  ',auto,auto:' + cnt1type + ':' + cnt2type + ':' + shdrange + ':' + \
259                  cnt1range + ':' + cnt2range + ':' + titlefig + ':' + kindfig +     \
260                  ':' + remap + ':' + mapvalues + ':yes'
261
262                ins = pyHOME + '/drawing.py -o draw_2D_shad_2cont -f ' + ncfile +    \
263                  " -S '" + values + "' -v " + varnames
264
265                oimgn = outfold + '/' + figk + '_' + varfignames.replace('@','_p').replace(',','_') +\
266                  '_' + timev + '.' + kindfig
267
268                print 'Lluis ins:', ins
269
270                if figscratch: sub.call('rm ' + oimgn, shell=True)
271
272                if not os.path.isfile(oimgn):
273                    try:
274                        with gen.Capturing() as output:
275                            sout = sub.call('python ' + ins, shell=True)
276                    except:
277                        print gen.errormsg
278                        print 'python ' + ins
279                        for s1out in output: print s1out
280                        quit(-1)
281                    if sout != 0:
282                        print gen.errormsg
283                        print '  : $ python ' + ins
284                        sub.call('python ' + ins, shell=True)
285                        quit(-1)
286
287                    ins = 'mv 2Dfields_shadow-2contour.' + kindfig + ' ' + oimgn
288                    sub.call(ins, shell=True)
289
290                    #quit()
291
292        else:
293            print gen.errormsg
294            print "  figure kind '" + figk + "' not found !!"
295            print '    available ones:',  availfigk
296
297
Note: See TracBrowser for help on using the repository browser.