source: lmdz_wrf/trunk/tools/drawing.py @ 193

Last change on this file since 193 was 193, checked in by lfita, 10 years ago

Adding non-checking of files in 'draw_lines'
Adding '|' character for spaces in title of 'draw_lines'
File size: 87.6 KB
Line 
1import numpy as np
2import os
3from netCDF4 import Dataset as NetCDFFile
4import drawing_tools as drw
5from optparse import OptionParser
6import sys
7from cStringIO import StringIO
8
9## g.e. # drawing.py -f /media/data2/etudes/WRF_LMDZ/WL_HyMeX/IIphase/medic950116/wlmdza/wrfout/wrfout_d01_1995-01-13_00:00:00 -o create_movie -S 'draw_2D_shad#Time@WRFTimes@10@95@191@1#tas:East_West|-1,North_South|-1,Time|2:longitude:latitude:Summer:270.,300.:tas|at|t=0:pdf:None:None' -v T2
10## g.e. # drawing.py -f wrfout_d01_1980-03-01_00\:00\:00_Time_B0-E48-I1.nc -o draw_2D_shad -S 'tas:East_West|-1,North_South|-1,Time|2:longitude:latitude:Summer:270.,300.:tas|at|t=0:pdf:None:cyl,i' -v T2
11## g.e. # drawing.py -f ~/etudes/domains/MEDCORDEX/geo_em.d01.nc -o draw_2D_shad_cont -S 'landmask,height:Time|0:Time|0:XLONG_M:XLAT_M:rainbow:fixc,k:%.2f:0,1:0.,3000.,10:landmask & height:pdf:False:lcc,i' -v LANDMASK,HGT_M
12## g.e. # drawing.py -f ~/etudes/domains/MEDCORDEX/geo_em.d01.nc -o draw_2D_shad_cont -S 'height,landmask:Time|0:Time|0:XLONG_M:XLAT_M:terrain:fixc,k:None:0.,3000.:0,1,10:MEDCORDEX height & landmask:pdf:False:lcc,i' -v HGT_M,LANDMASK
13## g.e. # drawing.py -o draw_2D_shad_line -f 'mean_dtcon-pluc-pres_lat.nc,mean_dtcon-pluc-pres_lat.nc' -S 'dtcon,prc:bottom_top|-1,south_north|-1:latmean:presmean:seismic,k:-5.,5.:monthly|dtcon|&|prc:pdf:flip@y:None:True' -v 'dtconmean,prcmean'
14## g.e. # drawing.py -f 'geo_em.d02.nc' -o draw_topo_geogrid -S '0.,3000.:None:FF_3dom d02:png:cyl,i'
15## g.e. # drawing.py -o draw_topo_geogrid_boxes -S '0.,3000.:None:FF domains:pdf:lcc,i:d01,d02,d03' -f geo_em.d01.nc,geo_em.d02.nc,geo_em.d03.nc
16## g.e. # drawing.py -o draw_trajectories -f 'WRF/control/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_HighRes_C/geo_em.d03.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdza/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdzb/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdzb_ii/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M' -S '$WRF_{CRM}$,$LMDZ_{AR4.0}$,$LMDZ_{NPv3.1}$,$LMDZ_{NPv3.1b}$|None|medicane trajectories|pdf|cyl,i' -v obs/trajectory.dat,satellite,-1
17## g.e. # drawing.py -o draw_vals_trajectories -f WRF_LMDZ/wlmdza/tevolboxtraj_T2.nc,WRF_LMDZ/wlmdzb/tevolboxtraj_T2.nc,WRF/control/tevolboxtraj_T2.nc -S 'mean:-1:$WRF_{CRM}$,$LMDZ_{AR4.0}$,$LMDZ_{NPv3.1}$@4:tas:time|($[DD]^[HH]$):exct,6,h:$%d^{%H}$:trajectory|following|mean:pdf' -v T2
18## g.e. # drawing.py -o draw_2D_shad_time -f 'netcdf_concatenated.nc' -S 'dtcon:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True' -v 'dtconmean'
19## g.e. # drawing.py -o variable_values -S PSFC
20## g.e. # drawing.py -o draw_timeSeries -f wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc -S 'dt_con:time|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time|evolution|at|-1|6|3|26:1:pdf' -v 'LDQCON,time'
21## g.e. # drawing.py -f wrfout_d01_1979-12-01_00:00:00 -o draw_Neighbourghood_evol -S 'q:Time|-1|Times,bottom_top|6|ZNU,south_north|3|XLAT,west_east|26|XLONG:south_north,west_east:5:auto:time|($[DD]^{[HH]}$),time|($[DD]^{[HH]}$):exct,2,h|exct,1,d:$%d^{%H}$,$%d^{%H}$:5|pts|neighbourghood|temporal|evolution:0.0,0.004:BuPu:pdf:True' -v QVAPOR
22
23main = 'drawing.py'
24
25errormsg = 'ERROR -- error -- ERROR -- error'
26warnmsg = 'WARNING -- waring -- WARNING -- warning'
27fillValue=1.e20
28
29namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time',                 \
30  'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line',                \
31  'draw_2D_shad_line_time', 'draw_timeSeries', 'draw_topo_geogrid',                  \
32  'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
33  'draw_lines', 'draw_Neighbourghood_evol', 'list_graphics', 'variable_values']
34
35def draw_2D_shad(ncfile, values, varn):
36    """ plotting a fields with shading
37    draw_2D_shad(ncfile, values, varn)
38      ncfile= file to use
39      values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
40       [kindfig]:[reverse]:[mapv]:[close]
41        [vnamefs]: Name in the figure of the variable to be shaded
42        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
43          variable a given value is required (-1, all the length)
44        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
45        [colorbar]: name of the color bar
46        [smin/axv]: minimum and maximum value for the shading or:
47          'Srange': for full range
48          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
49          'Saroundminmax@val': for min*val,max*val
50          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
51            percentile_(100-val)-median)
52          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
53          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
54          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
55             percentile_(100-val)-median)
56        [figt]: title of the figure ('|' for spaces)
57        [kindfig]: kind of figure
58        [reverse]: Transformation of the values
59          * 'transpose': reverse the axes (x-->y, y-->x)
60          * 'flip'@[x/y]: flip the axis x or y
61        [mapv]: map characteristics: [proj],[res]
62          see full documentation: http://matplotlib.org/basemap/
63          [proj]: projection
64            * 'cyl', cilindric
65            * 'lcc', lamvbert conformal
66          [res]: resolution:
67            * 'c', crude
68            * 'l', low
69            * 'i', intermediate
70            * 'h', high
71            * 'f', full
72      valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
73      varn= [varsn] name of the variable to plot with shading
74    """
75
76    fname = 'draw_2D_shad'
77    if values == 'h':
78        print fname + '_____________________________________________________________'
79        print draw_2D_shad.__doc__
80        quit()
81
82    expectargs = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]',    \
83      '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[mapv]', '[close]']
84 
85    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
86
87    vnamesfig = values.split(':')[0]
88    dimvals= values.split(':')[1].replace('|',':')
89    vdimxn = values.split(':')[2]
90    vdimyn = values.split(':')[3]
91    colbarn = values.split(':')[4]
92    shadminmax = values.split(':')[5]
93    figtitle = values.split(':')[6].replace('|',' ')
94    figkind = values.split(':')[7]
95    revals = values.split(':')[8]
96    mapvalue = values.split(':')[9]
97#    varn = values.split(':')[10]
98
99    ncfiles = ncfile
100   
101    if not os.path.isfile(ncfiles):
102        print errormsg
103        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
104        quit(-1)   
105
106    objsf = NetCDFFile(ncfiles, 'r')
107   
108    varns = varn.split(',')[0]
109
110    if  not objsf.variables.has_key(varns):
111        print errormsg
112        print '  ' + fname + ': shading file "' + ncfiles +                          \
113          '" does not have variable "' +  varns + '" !!'
114        quit(-1)
115
116# Variables' values
117    objvars = objsf.variables[varns]
118
119    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
120
121# Dimensions names
122##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
123##    dimnamesv = []
124##    for idd in range(len(objvars.dimensions)):
125##        cutdim = False
126##        for idc in range(len(dimvals.split(','))):
127##            dimcutn = dimvals.split(',')[idc].split(':')[0]
128##            print objvars.dimensions[idd], dimcutn
129##            if objvars.dimensions[idd] == dimcutn:
130##                cutdim = True
131##                break
132##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
133    dimnamesv = [vdimyn, vdimxn]
134
135    if drw.searchInlist(objvars.ncattrs(),'units'):
136        varunits = objvars.getncattr('units')
137    else:
138        print warnmsg
139        print '  ' + fname + ": variable '" + varn + "' without units!!"
140        varunits = '-'
141
142    if  not objsf.variables.has_key(vdimxn):
143        print errormsg
144        print '  ' + fname + ': shading file "' + ncfiles +                          \
145          '" does not have dimension variable "' +  vdimxn + '" !!'
146        quit(-1)
147    if  not objsf.variables.has_key(vdimyn):
148        print errormsg
149        print '  ' + fname + ': shading file "' + ncfiles +                          \
150          '" does not have dimension variable "' +  vdimyn + '" !!'
151        quit(-1)
152
153    objdimx = objsf.variables[vdimxn]
154    objdimy = objsf.variables[vdimyn]
155    if drw.searchInlist(objdimx.ncattrs(),'units'):
156        odimxu = objdimx.getncattr('units')
157    else:
158        print warnmsg
159        print '  ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
160        odimxu = '-'
161
162    if drw.searchInlist(objdimy.ncattrs(),'units'):
163        odimyu = objdimy.getncattr('units')
164    else:
165        print warnmsg
166        print '  ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
167        odimyu = '-'
168
169    if len(objdimx.shape) <= 2:
170#        odimxv = objdimx[valshad.shape]
171#        odimyv = objdimy[valshad.shape]
172        odimxv = objdimx[:]
173        odimyv = objdimy[:]
174
175    elif len(objdimx.shape) == 3:
176#        dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
177#        odimxv = objdimx[tuple(dimcut)]
178#        odimyv = objdimy[tuple(dimcut)]
179        odimxv = objdimx[0,:]
180        odimyv = objdimy[0,:]
181    else:
182        print errormsg
183        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
184          ' not ready!!'
185        quit(-1)
186
187    shading_nx = []
188    if shadminmax.split(',')[0][0:1] != 'S':
189            shading_nx.append(np.float(shadminmax.split(',')[0]))
190    else:
191        shading_nx.append(shadminmax.split(',')[0])
192
193    if shadminmax.split(',')[1][0:1] != 'S':
194        shading_nx.append(np.float(shadminmax.split(',')[1]))
195    else:
196        shading_nx.append(shadminmax.split(',')[1])
197
198    if mapvalue == 'None': mapvalue = None
199
200    drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, dimnamesv,\
201      colbarn, shading_nx, varunits, figtitle, figkind, revals, mapvalue, True)
202
203    return
204
205def draw_2D_shad_time(ncfile, values, varn):
206    """ plotting a fields with shading with time values
207    draw_2D_shad(ncfile, values, varn)
208      ncfile= file to use
209      values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
210       [kindfig]:[reverse]:[timevals]:[close]
211        [vnamefs]: Name in the figure of the variable to be shaded
212        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
213          variable a given value is required (-1, all the length)
214        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
215        [colorbar]: name of the color bar
216        [smin/axv]: minimum and maximum value for the shading or:
217          'Srange': for full range
218          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
219          'Saroundminmax@val': for min*val,max*val
220          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
221            percentile_(100-val)-median)
222          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
223          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
224          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
225             percentile_(100-val)-median)
226        [figt]: title of the figure ('|' for spaces)
227        [kindfig]: kind of figure
228        [reverse]: Transformation of the values
229          * 'transpose': reverse the axes (x-->y, y-->x)
230          * 'flip'@[x/y]: flip the axis x or y
231        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
232           [timen]; name of the time variable
233           [units]; units string according to CF conventions ([tunits] since
234             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
235           [kind]; kind of output
236             'Nval': according to a given number of values as 'Nval',[Nval]
237             'exct': according to an exact time unit as 'exct',[tunit];
238               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
239                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
240                'l': milisecond
241           [tfmt]; desired format
242           [label]; label at the graph ('!' for spaces)
243        [close]: should figure be closed (finished)
244      values='dtcon:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|'
245        'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True
246      varn= [varsn] name of the variable to plot with shading
247    """
248    fname = 'draw_2D_shad_time'
249
250    if values == 'h':
251        print fname + '_____________________________________________________________'
252        print draw_2D_shad_time.__doc__
253        quit()
254
255    farguments = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]',     \
256      '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', '[close]']
257    drw.check_arguments(fname,11,values,':',farguments)
258
259    vnamesfig = values.split(':')[0]
260    dimvals= values.split(':')[1].replace('|',':')
261    vdimxn = values.split(':')[2]
262    vdimyn = values.split(':')[3]
263    colbarn = values.split(':')[4]
264    shadminmax = values.split(':')[5]
265    figtitle = values.split(':')[6].replace('|',' ')
266    figkind = values.split(':')[7]
267    revals = values.split(':')[8]
268    timevals = values.split(':')[9]
269    close = values.split(':')[10]
270
271    ncfiles = ncfile
272   
273    if not os.path.isfile(ncfiles):
274        print errormsg
275        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
276        quit(-1)   
277
278    objsf = NetCDFFile(ncfiles, 'r')
279   
280    varns = varn.split(',')[0]
281
282    if  not objsf.variables.has_key(varns):
283        print errormsg
284        print '  ' + fname + ': shading file "' + ncfiles +                          \
285          '" does not have variable "' +  varns + '" !!'
286        quit(-1)
287
288# Variables' values
289    objvars = objsf.variables[varns]
290
291    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
292
293    dimnamesv = [vdimyn, vdimxn]
294
295    varunits = objvars.getncattr('units')
296
297    if  not objsf.variables.has_key(vdimxn):
298        print errormsg
299        print '  ' + fname + ': shading file "' + ncfiles +                          \
300          '" does not have dimension variable "' +  vdimxn + '" !!'
301        quit(-1)
302    if  not objsf.variables.has_key(vdimyn):
303        print errormsg
304        print '  ' + fname + ': shading file "' + ncfiles +                          \
305          '" does not have dimensino variable "' +  vdimyn + '" !!'
306        quit(-1)
307
308    objdimx = objsf.variables[vdimxn]
309    objdimy = objsf.variables[vdimyn]
310    odimxu = objdimx.getncattr('units')
311    odimyu = objdimy.getncattr('units')
312
313    if len(objdimx.shape) <= 2:
314        odimxv = objdimx[:]
315        odimyv = objdimy[:]
316
317    elif len(objdimx.shape) == 3:
318        odimxv = objdimx[0,:]
319        odimyv = objdimy[0,:]
320    else:
321        print errormsg
322        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
323          ' not ready!!'
324        quit(-1)
325
326    timename = timevals.split('|')[0]
327    timeunit = timevals.split('|')[1].replace('!',' ')
328    timekind = timevals.split('|')[2]
329    timefmt = timevals.split('|')[3]
330    timelabel = timevals.split('|')[4].replace('!',' ')
331
332    if vdimxn == timename:
333        odimxv = objsf.variables[vdimxn][:]
334        odimxu = timelabel
335        timeaxis = 'x'
336        odimyv = objsf.variables[vdimyn]
337        odimyu = odimyv.getncattr('units')
338        timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
339    elif vdimyn == timename:
340        odimyv = objsf.variables[vdimyn][:]
341        odimyu = timelabel
342        timeaxis = 'y'
343        odimxv = objsf.variables[vdimxn]
344        odimxu = odimxv.getncattr('units')
345        timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
346    else:
347        print errormsg
348        print '  ' + fname + ": time variable '" + timename + "' not found!!"
349        quit(-1)
350
351    shading_nx = []
352    if shadminmax.split(',')[0][0:1] != 'S':
353        shading_nx.append(np.float(shadminmax.split(',')[0]))
354    else:
355        shading_nx.append(shadminmax.split(',')[0])
356
357    if shadminmax.split(',')[1][0:1] != 'S':
358        shading_nx.append(np.float(shadminmax.split(',')[1]))
359    else:
360        shading_nx.append(shadminmax.split(',')[1])
361
362    closeval = drw.Str_Bool(close)
363
364    drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu,      \
365      dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
366      timepos, timelabels, closeval)
367
368    return
369
370def draw_2D_shad_cont(ncfile, values, varn):
371    """ plotting two fields, one with shading and the other with contour lines
372    draw_2D_shad_cont(ncfile, values, varn)
373      ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
374      values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]
375        [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
376        [dimvals/c]: list of [dimname]|[value] telling at which dimension of the
377          variable a given value is required (no dimension name, all the length)
378        [dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
379        [colorbar]: name of the color bar
380        [ckind]: kind of contours
381          'cmap': as it gets from colorbar
382          'fixc,[colname]': fixed color [colname], all stright lines
383          'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
384        [clabfmt]: format of the labels in the contour (None, also possible)
385        [smin/axv]: minimum and maximum value for the shading
386        [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
387        [figt]: title of the figure ('|' for spaces)
388        [kindfig]: kind of figure
389        [reverse]: does the values be transposed? 'True/False',
390        [mapv]: map characteristics: [proj],[res]
391          see full documentation: http://matplotlib.org/basemap/
392          [proj]: projection
393            * 'cyl', cilindric
394            * 'lcc', lamvbert conformal
395          [res]: resolution:
396            * 'c', crude
397            * 'l', low
398            * 'i', intermediate
399            * 'h', high
400            * 'f', full
401      valules= 'rh,ta:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal|average|of|rh|&|ta:pdf:flip@y:None'
402      varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
403    """
404
405    fname = 'draw_2D_shad_cont'
406    if values == 'h':
407        print fname + '_____________________________________________________________'
408        print draw_2D_shad_cont.__doc__
409        quit()
410
411    vnamesfig = values.split(':')[0].split(',')
412    dimvals= values.split(':')[1].replace('|',':')
413    dimvalc= values.split(':')[2].replace('|',':')
414    vdimxn = values.split(':')[3]
415    vdimyn = values.split(':')[4]
416    colbarn = values.split(':')[5]
417    countkind = values.split(':')[6]
418    countlabelfmt = values.split(':')[7]
419    shadminmax = values.split(':')[8]
420    contlevels = values.split(':')[9]
421    figtitle = values.split(':')[10].replace('|',' ')
422    figkind = values.split(':')[11]
423    revals = values.split(':')[12]
424    mapvalue = values.split(':')[13]
425
426    if2filenames = ncfile.find(',')
427
428    if if2filenames != -1:
429        ncfiles = ncfile.split(',')[0]
430        ncfilec = ncfile.split(',')[1]
431    else:
432        ncfiles = ncfile
433        ncfilec = ncfile
434
435    if not os.path.isfile(ncfiles):
436        print errormsg
437        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
438        quit(-1)   
439
440    if not os.path.isfile(ncfilec):
441        print errormsg
442        print '  ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
443        quit(-1)   
444
445    objsf = NetCDFFile(ncfiles, 'r')
446    objcf = NetCDFFile(ncfilec, 'r')
447   
448    varns = varn.split(',')[0]
449    varnc = varn.split(',')[1]
450
451    if  not objsf.variables.has_key(varns):
452        print errormsg
453        print '  ' + fname + ': shading file "' + ncfiles +                          \
454          '" does not have variable "' +  varns + '" !!'
455        quit(-1)
456
457    if  not objcf.variables.has_key(varnc):
458        print errormsg
459        print '  ' + fname + ': contour file "' + ncfilec +                          \
460          '" does not have variable "' +  varnc + '" !!'
461        quit(-1)
462
463# Variables' values
464    objvars = objsf.variables[varns]
465    objvarc = objcf.variables[varnc]
466
467    if len(objvars.shape) != len(objvarc.shape):
468        print errormsg
469        print '  ' + fname + ': shading variable "' + varns + '" has a shape: ',     \
470          objvars.shape,  'different than contour variable "' +  varnc + '": ',      \
471          objvarc.shape,' !!!'
472        quit(-1)
473
474    for idim in range(len(objvars.shape)):
475        if objvars.shape[idim] != objvarc.shape[idim]:
476            print errormsg
477            print '  ' + fname + ': shading variable "' + varns + '" has a shape: ', \
478              objvars.shape,  'different than contour variable "' +  varnc + '": ',  \
479              objvarc.shape,' !!!'
480            quit(-1)
481
482    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
483    valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|'))
484
485# Dimensions names
486##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
487##    dimnamesv = []
488##    for idd in range(len(objvars.dimensions)):
489##        cutdim = False
490##        for idc in range(len(dimvals.split(','))):
491##            dimcutn = dimvals.split(',')[idc].split(':')[0]
492##            print objvars.dimensions[idd], dimcutn
493##            if objvars.dimensions[idd] == dimcutn:
494##                cutdim = True
495##                break
496##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
497    dimnamesv = [vdimyn, vdimxn]
498
499    varunits = []
500    varunits.append(objvars.getncattr('units'))
501    varunits.append(objvarc.getncattr('units'))
502
503    if  not objsf.variables.has_key(vdimxn):
504        print errormsg
505        print '  ' + fname + ': shading file "' + ncfiles +                          \
506          '" does not have dimension variable "' +  vdimxn + '" !!'
507        quit(-1)
508    if  not objsf.variables.has_key(vdimyn):
509        print errormsg
510        print '  ' + fname + ': shading file "' + ncfiles +                          \
511          '" does not have dimensino variable "' +  vdimyn + '" !!'
512        quit(-1)
513
514    objdimx = objsf.variables[vdimxn]
515    objdimy = objsf.variables[vdimyn]
516    odimxu = objdimx.getncattr('units')
517    odimyu = objdimy.getncattr('units')
518
519    if len(objdimx.shape) <= 2:
520#        odimxv = objdimx[valshad.shape]
521#        odimyv = objdimy[valshad.shape]
522        odimxv = objdimx[:]
523        odimyv = objdimy[:]
524
525    elif len(objdimx.shape) == 3:
526#        dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
527#        odimxv = objdimx[tuple(dimcut)]
528#        odimyv = objdimy[tuple(dimcut)]
529        odimxv = objdimx[0,:]
530        odimyv = objdimy[0,:]
531    else:
532        print errormsg
533        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
534          ' not ready!!'
535        quit(-1)
536
537    if countlabelfmt == 'None': 
538        countlfmt = None
539    else:
540        countlfmt = countlabelfmt
541
542    shading_nx = np.zeros((2), dtype=np.float)
543    shading_nx[0] = np.float(shadminmax.split(',')[0])
544    shading_nx[1] = np.float(shadminmax.split(',')[1])
545
546    clevmin = np.float(contlevels.split(',')[0])
547    clevmax = np.float(contlevels.split(',')[1])
548    Nclevels = int(contlevels.split(',')[2])
549
550    levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
551
552    if len(levels_cont) <= 1: 
553        print warnmsg
554        print '  ' + fname + ': wrong contour levels:', levels_cont,' !!'
555        del(levels_cont)
556        levels_cont = np.zeros((Nclevels), dtype=np.float)
557        levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
558        print '    generating default ones: ',levels_cont
559
560    if mapvalue == 'None': mapvalue = None
561
562    drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu,  \
563      odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont,     \
564      varunits, figtitle, figkind, revals, mapvalue)
565
566    return
567
568def draw_2D_shad_cont_time(ncfile, values, varn):
569    """ plotting two fields, one with shading and the other with contour lines being
570    one of the dimensions of time characteristics
571    draw_2D_shad_cont(ncfile, values, varn)
572      ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
573      values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv]
574        [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
575        [dimvals/c]: list of [dimname]|[value] telling at which dimension of the
576          variable a given value is required (no dimension name, all the length)
577        [dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
578        [colorbar]: name of the color bar
579        [ckind]: kind of contours
580          'cmap': as it gets from colorbar
581          'fixc,[colname]': fixed color [colname], all stright lines
582          'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
583        [clabfmt]: format of the labels in the contour (None, also possible)
584        [smin/axv]: minimum and maximum value for the shading
585        [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
586        [figt]: title of the figure ('|' for spaces)
587        [kindfig]: kind of figure
588        [reverse]: modification to the dimensions:
589          'transposed': transpose matrices
590          'flip',[x/y]: flip only the dimension [x] or [y]
591        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label] time labels characteristics
592           [timen]; name of the time variable
593           [units]; units string according to CF conventions ([tunits] since
594             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
595           [kind]; kind of output
596             'Nval': according to a given number of values as 'Nval',[Nval]
597             'exct': according to an exact time unit as 'exct',[tunit];
598               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
599                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
600                'l': milisecond
601           [tfmt]; desired format
602           [label]; label at the graph ('!' for spaces)
603        [mapv]: map characteristics: [proj],[res]
604          see full documentation: http://matplotlib.org/basemap/
605          [proj]: projection
606            * 'cyl', cilindric
607            * 'lcc', lamvbert conformal
608          [res]: resolution:
609            * 'c', crude
610            * 'l', low
611            * 'i', intermediate
612            * 'h', high
613            * 'f', full
614      valules= 'rh,ta:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal|average|of|rh|&|ta:pdf:flip@y:time!hours!since!1949/12/01|exct,5d|%d|date!([DD]):None'
615      varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
616    """
617
618    fname = 'draw_2D_shad_cont_time'
619    if values == 'h':
620        print fname + '_____________________________________________________________'
621        print draw_2D_shad_cont_time.__doc__
622        quit()
623
624    expectargs = ['[vnamefs]', '[dimvals]', '[dimvalc]', '[dimxvn]', '[dimyvn]',     \
625      '[colorbar]', '[ckind]', '[clabfmt]', '[sminv],[smaxv]',                       \
626      '[sminc],[smaxv],[Nlev]', '[figt]', '[kindfig]', '[reverse]', '[timevals]',    \
627      '[mapv]']
628 
629    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
630
631    vnamesfig = values.split(':')[0].split(',')
632    dimvals= values.split(':')[1].replace('|',':')
633    dimvalc= values.split(':')[2].replace('|',':')
634    vdimxn = values.split(':')[3]
635    vdimyn = values.split(':')[4]
636    colbarn = values.split(':')[5]
637    countkind = values.split(':')[6]
638    countlabelfmt = values.split(':')[7]
639    shadminmax = values.split(':')[8]
640    contlevels = values.split(':')[9]
641    figtitle = values.split(':')[10].replace('|',' ')
642    figkind = values.split(':')[11]
643    revals = values.split(':')[12]
644    timevals = values.split(':')[13]
645    mapvalue = values.split(':')[14]
646
647    if2filenames = ncfile.find(',')
648
649    if if2filenames != -1:
650        ncfiles = ncfile.split(',')[0]
651        ncfilec = ncfile.split(',')[1]
652    else:
653        ncfiles = ncfile
654        ncfilec = ncfile
655
656    if not os.path.isfile(ncfiles):
657        print errormsg
658        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
659        quit(-1)   
660
661    if not os.path.isfile(ncfilec):
662        print errormsg
663        print '  ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
664        quit(-1)   
665
666    objsf = NetCDFFile(ncfiles, 'r')
667    objcf = NetCDFFile(ncfilec, 'r')
668   
669    varns = varn.split(',')[0]
670    varnc = varn.split(',')[1]
671
672    if  not objsf.variables.has_key(varns):
673        print errormsg
674        print '  ' + fname + ': shading file "' + ncfiles +                          \
675          '" does not have variable "' +  varns + '" !!'
676        quit(-1)
677
678    if  not objcf.variables.has_key(varnc):
679        print errormsg
680        print '  ' + fname + ': contour file "' + ncfilec +                          \
681          '" does not have variable "' +  varnc + '" !!'
682        quit(-1)
683
684# Variables' values
685    objvars = objsf.variables[varns]
686    objvarc = objcf.variables[varnc]
687
688    if len(objvars.shape) != len(objvarc.shape):
689        print errormsg
690        print '  ' + fname + ': shading variable "' + varns + '" has a shape: ',     \
691          objvars.shape,  'different than contour variable "' +  varnc + '": ',      \
692          objvarc.shape,' !!!'
693        quit(-1)
694
695    for idim in range(len(objvars.shape)):
696        if objvars.shape[idim] != objvarc.shape[idim]:
697            print errormsg
698            print '  ' + fname + ': shading variable "' + varns + '" has a shape: ', \
699              objvars.shape,  'different than contour variable "' +  varnc + '": ',  \
700              objvarc.shape,' !!!'
701            quit(-1)
702
703    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
704    valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|'))
705
706# Dimensions names
707##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
708##    dimnamesv = []
709##    for idd in range(len(objvars.dimensions)):
710##        cutdim = False
711##        for idc in range(len(dimvals.split(','))):
712##            dimcutn = dimvals.split(',')[idc].split(':')[0]
713##            print objvars.dimensions[idd], dimcutn
714##            if objvars.dimensions[idd] == dimcutn:
715##                cutdim = True
716##                break
717##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
718    dimnamesv = [vdimyn, vdimxn]
719
720    varunits = []
721    varunits.append(objvars.getncattr('units'))
722    varunits.append(objvarc.getncattr('units'))
723
724    if  not objsf.variables.has_key(vdimxn):
725        print errormsg
726        print '  ' + fname + ': shading file "' + ncfiles +                          \
727          '" does not have dimension variable "' +  vdimxn + '" !!'
728        quit(-1)
729    if  not objsf.variables.has_key(vdimyn):
730        print errormsg
731        print '  ' + fname + ': shading file "' + ncfiles +                          \
732          '" does not have dimension variable "' +  vdimyn + '" !!'
733        quit(-1)
734
735    timename = timevals.split('|')[0]
736    timeunit = timevals.split('|')[1].replace('!',' ')
737    timekind = timevals.split('|')[2]
738    timefmt = timevals.split('|')[3]
739    timelabel = timevals.split('|')[4].replace('!',' ')
740
741    if vdimxn == timename:
742        timevals = objsf.variables[vdimxn][:]
743        dimt = 'x'
744        ovalaxis = objsf.variables[vdimyn]
745        ovalu = ovalaxis.getncattr('units')
746    elif vdimyn == timename:
747        timevals = objsf.variables[vdimyn][:]
748        dimt = 'y'
749        ovalaxis = objsf.variables[vdimxn]
750        ovalu = ovalaxis.getncattr('units')
751    else:
752        print errormsg
753        print '  ' + fname + ": time variable '" + timename + "' not found!!"
754        quit(-1)
755
756    timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt)
757
758    if len(ovalaxis.shape) <= 2:
759        ovalaxisv = ovalaxis[:]
760
761    elif len(ovalaxis.shape) == 3:
762        ovalaxisv = ovalaxis[0,:]
763    else:
764        print errormsg
765        print '  ' + fname + ': shape of dimension variable:', ovalaxis.shape,       \
766          ' not ready!!'
767        quit(-1)
768
769    if countlabelfmt == 'None': 
770        countlfmt = None
771    else:
772        countlfmt = countlabelfmt
773
774    shading_nx = np.zeros((2), dtype=np.float)
775    shading_nx[0] = np.float(shadminmax.split(',')[0])
776    shading_nx[1] = np.float(shadminmax.split(',')[1])
777
778    clevmin = np.float(contlevels.split(',')[0])
779    clevmax = np.float(contlevels.split(',')[1])
780    Nclevels = int(contlevels.split(',')[2])
781
782    levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
783
784    if len(levels_cont) <= 1: 
785        print warnmsg
786        print '  ' + fname + ': wrong contour levels:', levels_cont,' !!'
787        del(levels_cont)
788        levels_cont = np.zeros((Nclevels), dtype=np.float)
789        levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
790        print '    generating default ones: ',levels_cont
791
792    if mapvalue == 'None': mapvalue = None
793
794    drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv,         \
795      timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn,    \
796      countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind,   \
797      revals, mapvalue)
798
799    return
800
801def draw_2D_shad_line(ncfile, values, varn):
802    """ plotting a fields with shading and another with line
803    draw_2D_shad_line(ncfile, values, varn)
804      ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
805      values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
806       [kindfig]:[reverse]:[mapv]:[close]
807        [vnamefs]: Name in the figure of the variable to be shaded
808        [vnamefl]: Name in the figure of the variable to be lined
809        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
810          variable a given value is required (-1, all the length)
811        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
812        [colorbar]: name of the color bar
813        [colline]: name of the color for the line
814        [smin/axv]: minimum and maximum value for the shading
815        [figt]: title of the figure ('|' for spaces)
816        [kindfig]: kind of figure
817        [reverse]: Transformation of the values
818          * 'transpose': reverse the axes (x-->y, y-->x)
819          * 'flip'@[x/y]: flip the axis x or y
820        [mapv]: map characteristics: [proj],[res]
821          see full documentation: http://matplotlib.org/basemap/
822          [proj]: projection
823            * 'cyl', cilindric
824            * 'lcc', lamvbert conformal
825          [res]: resolution:
826            * 'c', crude
827            * 'l', low
828            * 'i', intermediate
829            * 'h', high
830            * 'f', full
831      valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
832      varn= [varsn],[varnl] name of the variable to plot with shading and with line
833    """
834
835    fname = 'draw_2D_shad_line'
836    if values == 'h':
837        print fname + '_____________________________________________________________'
838        print draw_2D_shad_line.__doc__
839        quit()
840
841    farguments = ['[vnamefs],[vnamefl]', '[dimvals]', '[dimxvn]', '[dimyvn]',        \
842      '[colorbar],[colline]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
843      '[mapv]', '[close]']
844    drw.check_arguments(fname,11,values,':',farguments)
845
846    vnamesfig = values.split(':')[0].split(',')[0]
847    dimvals= values.split(':')[1].replace('|',':')
848    vdimxn = values.split(':')[2]
849    vdimyn = values.split(':')[3]
850    colbarn = values.split(':')[4].split(',')[0]
851    shadminmax = values.split(':')[5]
852    figtitle = values.split(':')[6].replace('|',' ')
853    figkind = values.split(':')[7]
854    revals = values.split(':')[8]
855    mapvalue = values.split(':')[9]
856#    varn = values.split(':')[10]
857
858    ncfiles = ncfile.split(',')[0]
859   
860    if not os.path.isfile(ncfiles):
861        print errormsg
862        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
863        quit(-1)   
864
865    objsf = NetCDFFile(ncfiles, 'r')
866   
867    varns = varn.split(',')[0]
868
869    if  not objsf.variables.has_key(varns):
870        print errormsg
871        print '  ' + fname + ': shading file "' + ncfiles +                          \
872          '" does not have variable "' +  varns + '" !!'
873        quit(-1)
874
875# Variables' values
876    objvars = objsf.variables[varns]
877
878    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
879
880# Dimensions names
881##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
882##    dimnamesv = []
883##    for idd in range(len(objvars.dimensions)):
884##        cutdim = False
885##        for idc in range(len(dimvals.split(','))):
886##            dimcutn = dimvals.split(',')[idc].split(':')[0]
887##            print objvars.dimensions[idd], dimcutn
888##            if objvars.dimensions[idd] == dimcutn:
889##                cutdim = True
890##                break
891##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
892    dimnamesv = [vdimyn, vdimxn]
893
894    varunits = objvars.getncattr('units')
895
896    if  not objsf.variables.has_key(vdimxn):
897        print errormsg
898        print '  ' + fname + ': shading file "' + ncfiles +                          \
899          '" does not have dimension variable "' +  vdimxn + '" !!'
900        quit(-1)
901    if  not objsf.variables.has_key(vdimyn):
902        print errormsg
903        print '  ' + fname + ': shading file "' + ncfiles +                          \
904          '" does not have dimensino variable "' +  vdimyn + '" !!'
905        quit(-1)
906
907    objdimx = objsf.variables[vdimxn]
908    objdimy = objsf.variables[vdimyn]
909    odimxu = objdimx.getncattr('units')
910    odimyu = objdimy.getncattr('units')
911
912    if len(objdimx.shape) <= 2:
913#        odimxv = objdimx[valshad.shape]
914#        odimyv = objdimy[valshad.shape]
915        odimxv = objdimx[:]
916        odimyv = objdimy[:]
917
918    elif len(objdimx.shape) == 3:
919#        dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
920#        odimxv = objdimx[tuple(dimcut)]
921#        odimyv = objdimy[tuple(dimcut)]
922        odimxv = objdimx[0,:]
923        odimyv = objdimy[0,:]
924    else:
925        print errormsg
926        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
927          ' not ready!!'
928        quit(-1)
929
930    shading_nx = np.zeros((2), dtype=np.float)
931    shading_nx[0] = np.float(shadminmax.split(',')[0])
932    shading_nx[1] = np.float(shadminmax.split(',')[1])
933
934    if mapvalue == 'None': mapvalue = None
935
936# line plot
937##
938    ncfilel = ncfile.split(',')[1]
939    vnamelfig = values.split(':')[0].split(',')[1]
940    varnl = varn.split(',')[1]
941    colline = values.split(':')[4].split(',')[1]
942
943    objlf = NetCDFFile(ncfilel,'r')
944    objlvar = objlf.variables[varnl]
945
946    linevals = objlvar[:]
947    varlunits = objlvar.units
948
949    drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
950      odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits,  \
951      figtitle, figkind, revals, mapvalue, True)
952
953    objsf.close()
954    objlf.close()
955
956    return
957
958def draw_2D_shad_line_time(ncfile, values, varn):
959    """ plotting a fields with shading and a line with time values
960    draw_2D_shad_line(ncfile, values, varn)
961      ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
962      values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
963       [kindfig]:[reverse]:[timevals]:[close]
964        [vnamefs]: Name in the figure of the variable to be shaded
965        [vnamefl]: Name in the figure of the variable to be lined
966        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
967          variable a given value is required (-1, all the length)
968        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
969        [colorbar]: name of the color bar
970        [smin/axv]: minimum and maximum value for the shading
971        [figt]: title of the figure ('|' for spaces)
972        [kindfig]: kind of figure
973        [reverse]: Transformation of the values
974          * 'transpose': reverse the axes (x-->y, y-->x)
975          * 'flip'@[x/y]: flip the axis x or y
976        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
977           [timen]; name of the time variable
978           [units]; units string according to CF conventions ([tunits] since
979             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
980           [kind]; kind of output
981             'Nval': according to a given number of values as 'Nval',[Nval]
982             'exct': according to an exact time unit as 'exct',[tunit];
983               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
984                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
985                'l': milisecond
986           [tfmt]; desired format
987           [label]; label at the graph ('!' for spaces)
988        [close]: should figure be closed (finished)
989      values='dtcon,prc:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|'
990        'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True
991      varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
992    """
993    fname = 'draw_2D_shad_line_time'
994    if values == 'h':
995        print fname + '_____________________________________________________________'
996        print draw_2D_shad__line_time.__doc__
997        quit()
998
999    farguments = ['[vnamefs],[vanemefl]', '[dimvals]', '[dimxvn]', '[dimyvn]',       \
1000      '[colorbar]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]',           \
1001      '[timevals]', '[close]']
1002    drw.check_arguments(fname,11,values,':',farguments)
1003
1004    vnamesfig = values.split(':')[0].split(',')[0]
1005    dimvals= values.split(':')[1].replace('|',':')
1006    vdimxn = values.split(':')[2]
1007    vdimyn = values.split(':')[3]
1008    colbarn = values.split(':')[4]
1009    shadminmax = values.split(':')[5]
1010    figtitle = values.split(':')[6].replace('|',' ')
1011    figkind = values.split(':')[7]
1012    revals = values.split(':')[8]
1013    timevals = values.split(':')[9]
1014    close = values.split(':')[10]
1015
1016    ncfiles = ncfile.split(',')[0]
1017   
1018    if not os.path.isfile(ncfiles):
1019        print errormsg
1020        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1021        quit(-1)   
1022
1023    objsf = NetCDFFile(ncfiles, 'r')
1024   
1025    varns = varn.split(',')[0]
1026
1027    if  not objsf.variables.has_key(varns):
1028        print errormsg
1029        print '  ' + fname + ': shading file "' + ncfiles +                          \
1030          '" does not have variable "' +  varns + '" !!'
1031        quit(-1)
1032
1033# Variables' values
1034    objvars = objsf.variables[varns]
1035
1036    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
1037
1038    dimnamesv = [vdimyn, vdimxn]
1039
1040    varunits = objvars.getncattr('units')
1041
1042    if  not objsf.variables.has_key(vdimxn):
1043        print errormsg
1044        print '  ' + fname + ': shading file "' + ncfiles +                          \
1045          '" does not have dimension variable "' +  vdimxn + '" !!'
1046        quit(-1)
1047    if  not objsf.variables.has_key(vdimyn):
1048        print errormsg
1049        print '  ' + fname + ': shading file "' + ncfiles +                          \
1050          '" does not have dimensino variable "' +  vdimyn + '" !!'
1051        quit(-1)
1052
1053    objdimx = objsf.variables[vdimxn]
1054    objdimy = objsf.variables[vdimyn]
1055    odimxu = objdimx.getncattr('units')
1056    odimyu = objdimy.getncattr('units')
1057
1058    if len(objdimx.shape) <= 2:
1059        odimxv = objdimx[:]
1060        odimyv = objdimy[:]
1061
1062    elif len(objdimx.shape) == 3:
1063        odimxv = objdimx[0,:]
1064        odimyv = objdimy[0,:]
1065    else:
1066        print errormsg
1067        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
1068          ' not ready!!'
1069        quit(-1)
1070
1071    timename = timevals.split('|')[0]
1072    timeunit = timevals.split('|')[1].replace('!',' ')
1073    timekind = timevals.split('|')[2]
1074    timefmt = timevals.split('|')[3]
1075    timelabel = timevals.split('|')[4].replace('!',' ')
1076
1077    if vdimxn == timename:
1078        odimxv = objsf.variables[vdimxn][:]
1079        odimxu = timelabel
1080        timeaxis = 'x'
1081        odimyv = objsf.variables[vdimyn]
1082        odimyu = odimyv.getncattr('units')
1083        timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1084    elif vdimyn == timename:
1085        odimyv = objsf.variables[vdimyn][:]
1086        odimyu = timelabel
1087        timeaxis = 'y'
1088        odimxv = objsf.variables[vdimxn]
1089        odimxu = odimxv.getncattr('units')
1090        timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1091    else:
1092        print errormsg
1093        print '  ' + fname + ": time variable '" + timename + "' not found!!"
1094        quit(-1)
1095
1096    shading_nx = np.zeros((2), dtype=np.float)
1097    shading_nx[0] = np.float(shadminmax.split(',')[0])
1098    shading_nx[1] = np.float(shadminmax.split(',')[1])
1099
1100    closeval = drw.Str_Bool(close)
1101
1102    drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu,      \
1103      dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1104      timepos, timelabels, False)
1105
1106# Line values
1107##
1108    ncfilel = ncfile.split(',')[1]
1109
1110    vnamelfig = values.split(':')[0].split(',')[1]
1111    varnl = varn.split(',')[1]
1112
1113    objlf = NetCDFFile(ncfilel,'r')
1114    objlvar = objlf.variables[varnl]
1115
1116    linevals = objlvar[:]
1117    if reva0 == 'tranpose':
1118        plt.plot (linevals, odimxv, '-', color='k')
1119    else:
1120        plt.plot (odimxv, linevals, '-', color='k')
1121
1122    objsf.close()
1123    objsl.close()
1124
1125    return
1126
1127def draw_topo_geogrid(ncfile, values):
1128    """ plotting geo_em.d[nn].nc topography from WPS files
1129    draw_topo_geogrid(ncfile, values)
1130      ncfile= geo_em.d[nn].nc file to use
1131      values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1132        [min/max]Topo: minimum and maximum values of topography to draw
1133        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1134        title: title of the graph
1135        graphic_kind: kind of figure (jpg, pdf, png)
1136        mapvalues: map characteristics [proj],[res]
1137          see full documentation: http://matplotlib.org/basemap/
1138          [proj]: projection
1139            * 'cyl', cilindric
1140            * 'lcc', lambert conformal
1141          [res]: resolution:
1142            * 'c', crude
1143            * 'l', low
1144            * 'i', intermediate
1145            * 'h', high
1146            * 'f', full
1147    """
1148    fname = 'draw_topo_geogrid'
1149
1150    if values == 'h':
1151        print fname + '_____________________________________________________________'
1152        print draw_topo_geogrid.__doc__
1153        quit()
1154
1155    expectargs = ['[minTopo]','[maxTopo]', '[lonlatL]', '[title]', '[graphic_kind]', \
1156      '[mapvalues]']
1157 
1158    drw.check_arguments(fname,5,values,':',expectargs)
1159
1160    mintopo = values.split(':')[0].split(',')[0]
1161    maxtopo = values.split(':')[0].split(',')[1]
1162
1163    lonlatLS = values.split(':')[1]
1164    lonlatLv = lonlatLS.split(',')[0]
1165
1166    if lonlatLv == 'None':
1167        lonlatL = None
1168    else:
1169        lonlatL = np.zeros((4), dtype=np.float)
1170        lonlatL[0] = np.float(lonlatLS.split(',')[0])
1171        lonlatL[1] = np.float(lonlatLS.split(',')[1])
1172        lonlatL[2] = np.float(lonlatLS.split(',')[2])
1173        lonlatL[3] = np.float(lonlatLS.split(',')[3])
1174
1175    grtit = values.split(':')[2]
1176    kindfig = values.split(':')[3]
1177    mapvalues = values.split(':')[4]
1178
1179    if not os.path.isfile(ncfile):
1180        print errormsg
1181        print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1182        quit(-1)   
1183
1184    objdomf = NetCDFFile(ncfile, 'r')
1185   
1186    objhgt = objdomf.variables['HGT_M']
1187    objlon = objdomf.variables['XLONG_M']
1188    objlat = objdomf.variables['XLAT_M']
1189
1190    topography = objhgt[0,:,:]
1191
1192    drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL,     \
1193      grtit, kindfig, mapvalues, True)
1194
1195    objdomf.close()
1196
1197    return
1198
1199def draw_topo_geogrid_boxes(ncfiles, values):
1200    """ plotting different geo_em.d[nn].nc topography from WPS files
1201    draw_topo_geogrid_boxes(ncfile, values)
1202      ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference)
1203      values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels]
1204        [min/max]Topo: minimum and maximum values of topography to draw
1205        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1206        title: title of the graph
1207        graphic_kind: kind of figure (jpg, pdf, png)
1208        mapvalues: map characteristics [proj],[res]
1209          see full documentation: http://matplotlib.org/basemap/
1210          [proj]: projection
1211            * 'cyl', cilindric
1212            * 'lcc', lambert conformal
1213          [res]: resolution:
1214            * 'c', crude
1215            * 'l', low
1216            * 'i', intermediate
1217            * 'h', high
1218            * 'f', full
1219        labels= labels to write in the graph
1220    """
1221#    import matplotlib as mpl
1222#    mpl.use('Agg')
1223    import matplotlib.pyplot as plt
1224
1225    fname = 'draw_topo_geogrid_boxes'
1226
1227    if values == 'h':
1228        print fname + '_____________________________________________________________'
1229        print draw_topo_geogrid_boxes.__doc__
1230        quit()
1231
1232    mintopo = values.split(':')[0].split(',')[0]
1233    maxtopo = values.split(':')[0].split(',')[1]
1234
1235    lonlatLS = values.split(':')[1]
1236    lonlatLv = lonlatLS.split(',')[0]
1237
1238    if lonlatLv == 'None':
1239        lonlatL = None
1240    else:
1241        lonlatL = np.zeros((4), dtype=np.float)
1242        lonlatL[0] = np.float(lonlatLS.split(',')[0])
1243        lonlatL[1] = np.float(lonlatLS.split(',')[1])
1244        lonlatL[2] = np.float(lonlatLS.split(',')[2])
1245        lonlatL[3] = np.float(lonlatLS.split(',')[3])
1246
1247    grtit = values.split(':')[2]
1248    kindfig = values.split(':')[3]
1249    mapvalues = values.split(':')[4]
1250    labels = values.split(':')[5]
1251
1252    ncfile = ncfiles.split(',')[0]
1253    if not os.path.isfile(ncfile):
1254        print errormsg
1255        print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1256        quit(-1)   
1257
1258    objdomf = NetCDFFile(ncfile, 'r')
1259   
1260    objhgt = objdomf.variables['HGT_M']
1261    objlon0 = objdomf.variables['XLONG_M']
1262    objlat0 = objdomf.variables['XLAT_M']
1263
1264    topography = objhgt[0,:,:]
1265
1266    Nfiles = len(ncfiles.split(','))
1267    boxlabels = labels.split(',')
1268
1269    Xboxlines = []
1270    Yboxlines = []
1271
1272    for ifile in range(Nfiles):
1273        ncfile = ncfiles.split(',')[ifile]
1274#        print ifile, ncfile
1275        if not os.path.isfile(ncfile):
1276            print errormsg
1277            print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1278            quit(-1)   
1279
1280        objdomfi = NetCDFFile(ncfile, 'r')
1281   
1282        objlon = objdomfi.variables['XLONG_M']
1283        objlat = objdomfi.variables['XLAT_M']
1284
1285        dx = objlon.shape[2]
1286        dy = objlon.shape[1]
1287
1288        Xboxlines.append(objlon[0,0,:])
1289        Yboxlines.append(objlat[0,0,:])
1290        Xboxlines.append(objlon[0,dy-1,:])
1291        Yboxlines.append(objlat[0,dy-1,:])
1292        Xboxlines.append(objlon[0,:,0])
1293        Yboxlines.append(objlat[0,:,0])
1294        Xboxlines.append(objlon[0,:,dx-1])
1295        Yboxlines.append(objlat[0,:,dx-1])
1296
1297        objdomfi.close()
1298
1299    drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels,         \
1300      objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, True)
1301
1302    objdomf.close()
1303
1304    return
1305
1306def movievalslice(origslice, dimmovien, framenum):
1307    """ Function to provide variable slice according to a geneation of a movie
1308    movievals(origslice, dimmovien, framenum)
1309      [origslice]= slice original as [dimname1]|[val1],[...,[dimnameN]|[valN]]
1310        ([val] = -1, full length)
1311      [dimmovien]= name of the dimension to produce the movie
1312      [framenum]= value of the frame to substitue in [origslice] as
1313        [dimmovien]|[framenum]
1314    >>> movievalslice('East_West|-1,North_South|-1,Time|2','Time',0)
1315    East_West|-1,North_South|-1,Time|0
1316    """
1317
1318    fname = 'movievalslice'
1319
1320    if origslice == 'h':
1321        print fname + '_____________________________________________________________'
1322        print movievalslice.__doc__
1323        quit()
1324   
1325    dims = origslice.split(',')
1326
1327    movieslice = ''
1328    idim = 0
1329
1330    for dimn in dims:
1331        dn = dimn.split('|')[0]
1332        if dn == dimmovien:
1333            movieslice = movieslice + dn + '|' + str(framenum)
1334        else:
1335            movieslice = movieslice + dimn
1336        if idim < len(dims)-1: movieslice = movieslice + ','
1337
1338        idim = idim + 1
1339
1340    return movieslice
1341
1342class Capturing(list):
1343    """ Class to capture function output as a list
1344    from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call
1345    """
1346#    from cStringIO import StringIO
1347
1348    def __enter__(self):
1349        self._stdout = sys.stdout
1350        sys.stdout = self._stringio = StringIO()
1351        return self
1352    def __exit__(self, *args):
1353        self.extend(self._stringio.getvalue().splitlines())
1354        sys.stdout = self._stdout
1355
1356def create_movie(netcdfile, values, variable):
1357    """ Function to create a movie assuming ImageMagick installed!
1358      values= [graph]#[movie_dimension]#[graph_values]
1359        [graph]: which graphic
1360        [movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval]
1361          [dimnmovie]; name of the dimension from which make the movie
1362          [dimvmovie]; name of the variable with the values of the dimension
1363          [moviedelay]; delay between frames
1364          [interval]; [beg]@[end]@[freq] or -1 (all)
1365        [graph_values]: values to generate the graphic
1366      netcdfile= netCDF file
1367      variable= variable to use (when applicable)
1368    """ 
1369    fname = 'create_movie'
1370
1371    if values == 'h':
1372        print fname + '_____________________________________________________________'
1373        print create_movie.__doc__
1374        quit()
1375
1376    graph = values.split('#')[0]
1377    movie_dim = values.split('#')[1]
1378    graph_vals = values.split('#')[2]
1379
1380    ncobj = NetCDFFile(netcdfile, 'r')
1381
1382# Movie dimension
1383##
1384    dimnmovie = movie_dim.split('@')[0]
1385    dimvmovie = movie_dim.split('@')[1]
1386    moviedelay = movie_dim.split('@')[2]
1387    moviebeg = int(movie_dim.split('@')[3])
1388
1389    if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie):
1390        print errormsg
1391        print '  ' + fname + ": file '" + netcdfile + "' has not dimension '" +      \
1392          dimnmovie + "' !!!"
1393        quit(-1)
1394
1395    objdmovie = ncobj.dimensions[dimnmovie]
1396    dmovie = len(objdmovie)
1397    if moviebeg != -1:
1398        moviend = int(movie_dim.split('@')[4])
1399        moviefreq = int(movie_dim.split('@')[5])
1400    else:
1401        moviebeg = 0
1402        moviend = dmovie
1403        moviefreq = 1
1404
1405    if dimvmovie == 'WRFTimes':
1406        objvdmovie = ncobj.variables['Times']
1407        vdmovieunits = ''
1408        valsdmovie = []
1409        for it in range(objvdmovie.shape[0]):
1410            valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:],           \
1411              'WRFdatetime', 'Y/m/d H-M-S'))
1412    elif dimvmovie == 'CFtime':
1413        objvdmovie = ncobj.variables['time']
1414        vdmovieunits = ''
1415        print objvdmovie.units
1416        valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard',  \
1417          objvdmovie[:])
1418        valsdmovie = []
1419        for it in range(objvdmovie.shape[0]):
1420            valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:],          \
1421              'matYmdHMS', 'Y/m/d H-M-S'))
1422    else:
1423        if  not drw.searchInlist(ncobj.variables.keys(),dimvmovie):
1424            print errormsg
1425            print '  ' + fname + ": file '" + netcdfile + "' has not variable '" +   \
1426              dimvmovie + "' !!!"
1427            quit(-1)
1428        vdmovieunits = objvdmovie.getncattr('units')
1429        objvdmovie = ncobj.variables[dimvmovie]
1430        if len(objvdmovie.shape) == 1:
1431            vasldmovie = objvdmovie[:]
1432        else:
1433            print errormsg
1434            print '  ' + fname + ': shape', objvdmovie.shape, 'of variable with ' +  \
1435              'dimension movie values not ready!!!'
1436            quit(-1)
1437
1438    ncobj.close()
1439    os.system('rm frame_*.png > /dev/null')
1440
1441# graphic
1442##
1443    if graph == 'draw_2D_shad':
1444        graphvals = graph_vals.split(':')
1445
1446        for iframe in range(moviebeg,moviend,moviefreq):
1447            iframeS = str(iframe).zfill(4)
1448
1449            drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \
1450              5, 'frames')
1451            titgraph = dimnmovie + '|=|' + str(valsdmovie[iframe]) + '|' +           \
1452              vdmovieunits
1453
1454            graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe)
1455            graphvals[6] = titgraph
1456            graphvals[7] = 'png'
1457
1458            graphv = drw.numVector_String(graphvals, ":")
1459
1460            with Capturing() as output:
1461                draw_2D_shad(netcdfile, graphv, variable)
1462
1463            os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png')
1464    else:
1465        print errormsg
1466        print '  ' + fname + ": graphic '" +  graph + "' not defined !!!"
1467        quit(-1)
1468
1469    os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif')
1470
1471    print "Succesfuly creation of movie file 'create_movie.gif' !!!"
1472
1473    return
1474
1475def draw_lines(ncfilens, values, varname):
1476    """ Function to draw different lines at the same time from different files
1477    draw_trajectories(ncfilens, values, varname):
1478      ncfilens= [filen] ',' separated list of netCDF files
1479      values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[graphk]
1480        [dimvname]: name of the variable with he values of the common dimension
1481        [valuesaxis]: which axis will be used for the values ('x', or 'y')
1482        [dimtit]: title for the common dimension
1483        [leglabels]: ',' separated list of names for the legend
1484        [vartit]: name of the variable in the graph
1485        [title]: title of the plot ('|' for spaces)
1486        [graphk]: kind of the graphic
1487      varname= variable to plot
1488      values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1489    """
1490
1491    fname = 'draw_lines'
1492
1493    if values == 'h':
1494        print fname + '_____________________________________________________________'
1495        print draw_lines.__doc__
1496        quit()
1497
1498    expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[graphk]'
1499    drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
1500
1501    ncfiles = ncfilens.split(',')
1502    dimvname = values.split(':')[0]
1503    valuesaxis = values.split(':')[1]
1504    dimtit = values.split(':')[2]
1505    leglabels = values.split(':')[3].replace('_','\_')
1506    vartit = values.split(':')[4]
1507    title = values.split(':')[5].replace('|',' ')
1508    graphk = values.split(':')[6]
1509
1510    Nfiles = len(ncfiles)
1511
1512# Getting trajectotries
1513##
1514
1515    varvalues = []
1516    dimvalues = []
1517
1518    print '  ' + fname
1519    ifn = 0
1520    for ifile in ncfiles:
1521        filen = ifile.split('@')[0]
1522
1523        print '    filen:',filen
1524
1525        if not os.path.isfile(filen):
1526            print errormsg
1527            print '  ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1528            quit(-1)
1529
1530        objfile = NetCDFFile(filen, 'r')
1531
1532        if not objfile.variables.has_key(dimvname):
1533            print errormsg
1534            print '  ' + fname + ": netCDF file '" + filen +                         \
1535              "' does not have variable '" + dimvname + "' !!"
1536            quit(-1)
1537
1538        if not objfile.variables.has_key(varname):
1539            print errormsg
1540            print '  ' + fname + ": netCDF file '" + filen +                         \
1541              "' does not have variable '" + varname + "' !!"
1542            quit(-1)
1543
1544        vvobj = objfile.variables[varname]
1545        if len(vvobj.shape) != 1:
1546            print errormsg
1547            print '  ' + fname + ': wrong shape:',vvobj.shape," of variable '" +     \
1548              varname +  "' !!"
1549            quit(-1)
1550
1551        vdobj = objfile.variables[dimvname]
1552        if len(vdobj.shape) != 1:
1553            print errormsg
1554            print '  ' + fname + ': wrong shape:',vdobj.shape," of variable '" +     \
1555              dimvname +  "' !!"
1556            quit(-1)
1557
1558        varvalues.append(vvobj[:])
1559        dimvalues.append(vdobj[:])
1560
1561        if ifn == 0:
1562            varunits = vvobj.units
1563
1564        objfile.close()
1565
1566        ifn = ifn + 1
1567
1568
1569    drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','),   \
1570      vartit, varunits, title, graphk)
1571
1572    return
1573
1574def draw_Neighbourghood_evol(filen, values, variable):
1575    """ Function to draw the temporal evolution of a neighbourghood around a point
1576    draw_Neighbourghood_evol(filen, values, variable)
1577      filen= netCDF file name
1578      values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
1579       [timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
1580        [dimsval]: [dimn1]|[val1]|[dimv1],...,[dimnN]|[valN]|[dimvN] dimension names, values to get
1581          (-1, for all; no name/value pair given full length) and variable with values of the dimension
1582          NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
1583        [neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
1584        [Nneig]: Number of grid points of the full side of the box (odd value)
1585        [Ncol]: Number of columns ('auto': square final plot)
1586        [gvarname]: name of the variable to appear in the graph
1587        [timetits]: [titX],[titY] titles of the axes ('|' for spaces)
1588        [tkinds]: [tkindX]|[tkindY] kinds of time to appear in the graph
1589          'Nval': according to a given number of values as 'Nval',[Nval]
1590          'exct': according to an exact time unit as 'exct',[tunit];
1591            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1592              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1593              'l': milisecond
1594        [timefmts]: [tfmtX],[tfmtY] format of the time labels
1595        [gtitle]: title of the graphic ('|' for spaces)
1596        [shadxtrms]: Extremes for the shading
1597        [cbar]: colorbar to use
1598        [gkind]: kind of graphical output
1599        [ofile]: True/False whether the netcdf with data should be created or not
1600      variable= name of the variable
1601      values = 'q:Time|-1|Times,bottom_top|6|ZNU,south_north|3|XLAT,west_east|26|XLONG:south_north,west_east:5:auto:time|($[DD]^{[HH]}$),time|($[DD]^{[HH]}$):exct,2,h|exct,1,d:$%d^{%H}$,$%d^{%H}$:5|pts|neighbourghood|temporal|evolution:0.0,0.004:BuPu:pdf:True'
1602    """ 
1603
1604    fname = 'draw_Neighbourghood_evol'
1605
1606    if values == 'h':
1607        print fname + '_____________________________________________________________'
1608        print draw_Neighbourghood_evol.__doc__
1609        quit()
1610
1611    expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]',      \
1612      '[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]',   \
1613      '[gkind]', '[ofile]']
1614 
1615    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
1616
1617    gvarname = values.split(':')[0]
1618    dimsval = values.split(':')[1].split(',')
1619    neigdims = values.split(':')[2].split(',')
1620    Nneig = int(values.split(':')[3])
1621    Ncol0 = values.split(':')[4]
1622    timetits = values.split(':')[5].split(',')
1623    timekinds = values.split(':')[6].split('|')
1624    timefmts = values.split(':')[7].split(',')
1625    gtitle = values.split(':')[8].replace('|',' ')
1626    shadxtrms = values.split(':')[9].split(',')
1627    cbar = values.split(':')[10]
1628    gkind = values.split(':')[11]
1629    ofile = values.split(':')[12]
1630
1631    if Ncol0 != 'auto': 
1632        Ncol = int(Ncol0)
1633    else:
1634        Ncol = Ncol0
1635
1636    timetits[0] = timetits[0].replace('|',' ')
1637    timetits[1] = timetits[1].replace('|',' ')
1638
1639    if np.mod(Nneig,2) == 0:
1640        print errormsg
1641        print '  ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
1642        quit(-1)
1643
1644    Nneig2 = int(Nneig/2)
1645
1646# Values to slice the variable
1647    dimvslice = {}
1648    dimvvalues = {}
1649    for dimvs in dimsval:
1650        dimn = dimvs.split('|')[0]
1651        dimv = int(dimvs.split('|')[1])
1652        dimnv = dimvs.split('|')[2]
1653
1654        dimvvalues[dimn] = dimnv
1655        dimvslice[dimn] = dimv
1656
1657    ncobj = NetCDFFile(filen, 'r')
1658
1659    varobj = ncobj.variables[variable]
1660
1661    slicevar = []
1662    newdimn = []
1663    newdimsvar = {}
1664
1665    for dimn in varobj.dimensions:
1666        if not drw.searchInlist(dimvslice.keys(), dimn):
1667            dimsize = len(ncobj.dimensions[dimn])
1668            slicevar.append(slice(0, dimsize+1))
1669            newdimn.append(dimn)
1670            newdimsvar[dimn] = dimseize
1671
1672        for dimslicen in dimvslice.keys():
1673            if dimn == dimslicen:
1674                if dimvslice[dimn] != -1:
1675                    if drw.searchInlist(neigdims, dimn):
1676                        slicevar.append(slice(dimvslice[dimn]-Nneig2,                \
1677                          dimvslice[dimn]+Nneig2+1))
1678                        newdimn.append(dimn)
1679                        newdimsvar[dimn] = Nneig
1680                        break
1681                    else:
1682                        slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
1683                        break
1684                else:
1685                    dimsize = len(ncobj.dimensions[dimn])
1686                    slicevar.append(slice(0, dimsize+1))
1687                    newdimn.append(dimn)
1688                    newdimsvar[dimn] = dimsize
1689                    break
1690 
1691    varv = varobj[tuple(slicevar)]
1692
1693    if len(newdimn) != 3:
1694        print errormsg
1695        print '  ' + fname + ': sliced variable with shape=', varv.shape,            \
1696          ' must have three dimensions',len(varv.shape),'given !!'
1697        quit(-1)
1698
1699    newdims = []
1700    for nwdims in newdimn:
1701        newdims.append(newdimsvar[nwdims])
1702
1703# The dimension which is not in the neighbourhood dimensions must be time!
1704    for dim1 in newdimn:
1705        if not drw.searchInlist(neigdims, dim1):
1706            dimt = newdimsvar[dim1]
1707            dimtime = dim1
1708
1709    if Ncol == 'auto':
1710        dimtsqx = int(np.sqrt(dimt)) + 1
1711        dimtsqy = int(np.sqrt(dimt)) + 1
1712    else:
1713        dimtsqx = int(Ncol)
1714        dimtsqy = dimt/dimtsqx + 1
1715
1716    neighbourghood = np.ones((dimtsqy*Nneig,dimtsqx*Nneig), dtype=np.float)*fillValue
1717
1718    for it in range(dimt):
1719        ity = int(it/dimtsqx)
1720        itx = it-ity*dimtsqx
1721
1722        itty = (dimtsqy - ity - 1)*Nneig + Nneig2
1723        ittx = itx*Nneig + Nneig2
1724
1725        neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]=         \
1726          varv[it,::-1,:]
1727
1728    variablevals = drw.variables_values(variable)
1729    if drw.searchInlist(varobj.ncattrs(), 'units'):
1730        vunits = varobj.units
1731    else:
1732        vunits = variablevals[5]
1733
1734# Time values at the X/Y axes
1735    if ncobj.variables[dimvvalues[dimtime]].dtype == '|S1':
1736        print '    ' + fname + ': WRF time variable!'
1737        refdate = '19491201000000'
1738        tunitsval = 'hours'
1739        dimtvalues = np.zeros((dimt), dtype=np.float)
1740        tvals = ncobj.variables[dimvvalues[dimtime]]
1741        yrref=refdate[0:4]
1742        monref=refdate[4:6]
1743        dayref=refdate[6:8]
1744        horref=refdate[8:10]
1745        minref=refdate[10:12]
1746        secref=refdate[12:14]
1747
1748        refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' +        \
1749          minref + ':' + secref
1750        tunits = tunitsval + ' since ' + refdateS
1751        for it in range(dimt):
1752            wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
1753            dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
1754    else:
1755        dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
1756        tunits = ncobj.variables[newdimsvar[dimtime]].units
1757
1758    dimxv = dimtvalues[0:dimtsqx]
1759    dimyv = dimtvalues[0:dimt:dimtsqx]
1760
1761    print 'Lluis lens: dimxv, dimyv:',len(dimxv), len(dimyv)
1762
1763    dimn = ['time','time']
1764
1765    if ofile == 'True':
1766        ofilen = 'Neighbourghood_evol.nc'
1767        newnc = NetCDFFile(ofilen, 'w')
1768# Dimensions
1769        newdim = newnc.createDimension('time',None)
1770        newdim = newnc.createDimension('y',dimtsqy*Nneig)
1771        newdim = newnc.createDimension('x',dimtsqx*Nneig)
1772# Dimension values
1773        newvar = newnc.createVariable('time','f8',('time'))
1774        newvar[:] = np.arange(dimt)
1775        newattr = drw.basicvardef(newvar, 'time','time',tunits)
1776# Neighbourhghood variable
1777        newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'),        \
1778          fill_value=fillValue)
1779        newvar[:] = neighbourghood
1780
1781        newnc.sync()
1782        newnc.close()
1783        print fname + ": Successfull generation of file '" + ofilen + "' !!"
1784
1785# Time ticks
1786    timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
1787    timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
1788
1789    timepos = [timeposX, timeposY[::-1]]
1790    timelabels = [timelabelsX, timelabelsY[::-1]]
1791
1792    for i in range(2):
1793        if shadxtrms[i][0:1] != 'S':
1794            shadxtrms[i] = np.float(shadxtrms[i])
1795
1796    drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits,   \
1797      timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
1798
1799def draw_timeSeries(filen, values, variables):
1800    """ Function to draw a time-series
1801    draw_timeSeries(filen, values, variable):
1802      filen= name of the file
1803      values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]
1804      [gvarname]: name of the variable to appear in the graph
1805      [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
1806      [tkind]: kind of time to appear in the graph (assumed x-axis)
1807        'Nval': according to a given number of values as 'Nval',[Nval]
1808        'exct': according to an exact time unit as 'exct',[tunit];
1809          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1810            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1811            'l': milisecond
1812      [timefmt]: format of the time labels
1813      [title]: title of the graphic ('|' for spaces)
1814      [locleg]: location of the legend (-1, autmoatic)
1815        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1816        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1817        9: 'upper center', 10: 'center'
1818      [gkind]: kind of graphical output
1819      variables= [varname],[timename] names of variable and variable with times
1820      draw_timeSeries('wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc', 'dt_con:time|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time|evolution|at|-1|6|3|26:1:pdf', 'LDQCON,time')
1821    """
1822
1823    fname = 'draw_timeSeries'
1824
1825    if values == 'h':
1826        print fname + '_____________________________________________________________'
1827        print draw_timeSeries.__doc__
1828        quit()
1829
1830    expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]',      \
1831      '[locleg]', '[gkind]']
1832 
1833    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
1834
1835    gvarname = values.split(':')[0]
1836    timetit = values.split(':')[1].replace('|',' ')
1837    tkind = values.split(':')[2]
1838    timefmt = values.split(':')[3]
1839    title = values.split(':')[4].replace('|',' ')
1840    locleg = int(values.split(':')[5])
1841    gkind = values.split(':')[6]
1842   
1843    ncobj = NetCDFFile(filen, 'r')
1844
1845    variable = variables.split(',')[0]
1846    timevar = variables.split(',')[1]
1847
1848    if not ncobj.variables.has_key(variable):
1849        print errormsg
1850        print '  ' + fname + ": file '" +  filen + "' does not have variable '" +    \
1851          variable + "' !!"
1852        quit(-1)
1853
1854    if not ncobj.variables.has_key(timevar):
1855        print errormsg
1856        print '  ' + fname + ": file '" +  filen + "' does not have variable time '" \
1857          + timevar + "' !!"
1858        quit(-1)
1859
1860    varobj = ncobj.variables[variable]
1861    timeobj = ncobj.variables[timevar]
1862
1863    dimt = len(timeobj[:])
1864    varvals = np.zeros((2,dimt), dtype=np.float)
1865
1866    gunits = varobj.getncattr('units')
1867    tunits = timeobj.getncattr('units')
1868
1869    varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
1870    varvals[1,:] = timeobj[:]
1871
1872    tseriesvals = []
1873    tseriesvals.append(varvals)
1874
1875    drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits,        \
1876      'TimeSeries', gvarname, timetit, tkind, timefmt, title,                        \
1877      [gvarname.replace('_','\_')], locleg, gkind)
1878
1879    return
1880
1881#draw_timeSeries('wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc', 'dt_con:time|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time|evolution|at|-1|6|3|26:1:pdf', 'LDQCON,time')
1882
1883def draw_trajectories(trjfilens, values, observations):
1884    """ Function to draw different trajectories at the same time
1885    draw_trajectories(trjfilens, values, observations):
1886      trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
1887         time intervals and reference maps (first one will be used to plot)
1888        [filen]: name of the file to use (lines with '#', not readed) as:
1889          [t-step] [x] [y]
1890        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
1891        [map]: [file]#[lonname]#[latname]
1892          [file]; with the [lon],[lat] matrices
1893          [lonname],[latname]; names of the longitudes and latitudes variables
1894      values=[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]
1895        [leglabels]: ',' separated list of names for the legend
1896        [lonlatlims]: limits of the map [lonmin, latmin, lonmax, latmax] or None
1897        [title]: title of the plot
1898        [graphk]: kind of the graphic
1899        [mapkind]: drawing coastaline ([proj],[res]) or None
1900          [proj]: projection
1901             * 'cyl', cilindric
1902             * 'lcc', lambert conformal
1903          [res]: resolution:
1904             * 'c', crude
1905             * 'l', low
1906             * 'i', intermediate
1907             * 'h', high
1908             * 'f', full
1909      obsevations= [obsfile],[obsname],[Tint],[null]
1910        [obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
1911        [obsname]: name of the observations in the graph
1912        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
1913        [null]: null value for the observed trajectory
1914    """
1915
1916    fname = 'draw_trajectories'
1917
1918    if values == 'h':
1919        print fname + '_____________________________________________________________'
1920        print draw_trajectories.__doc__
1921        quit()
1922
1923    trjfiles = trjfilens.split(',')
1924    leglabels = values.split('|')[0]
1925    lonlatlims = values.split('|')[1]
1926    title = values.split('|')[2]
1927    graphk = values.split('|')[3]
1928    mapkind = values.split('|')[4]
1929
1930    Nfiles = len(trjfiles)
1931
1932# Getting trajectotries
1933##
1934
1935    lontrjvalues = []
1936    lattrjvalues = []
1937
1938    print '  ' + fname
1939    ifn = 0
1940    for ifile in trjfiles:
1941        filen = ifile.split('@')[0]
1942        Tint = ifile.split('@')[1]
1943
1944        print '    trajectory:',filen
1945
1946        if Tint != '-1':
1947            Tbeg = Tint
1948            Tend = ifile.split('@')[2]
1949            mapv = ifile.split('@')[3]
1950        else:
1951            mapv = ifile.split('@')[2]
1952
1953        if not os.path.isfile(filen):
1954            print errormsg
1955            print '  ' + fname + ": trajectory file '" + filen + "' does not exist !!"
1956            quit(-1)
1957
1958# Charging longitude and latitude values
1959##
1960        lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
1961          mapv.split('#')[2])
1962
1963        if ifn == 0: mapref = mapv
1964        ifn = ifn + 1
1965
1966        objfile = open(filen, 'r')
1967        trjtimev = []
1968        trjxv = []
1969        trjyv = []
1970
1971        for line in objfile:
1972            if line[0:1] != '#':
1973                trjtimev.append(int(line.split(' ')[0]))
1974                trjxv.append(int(line.split(' ')[1]))
1975                trjyv.append(int(line.split(' ')[2]))
1976
1977        objfile.close()
1978
1979        if Tint != '-1':
1980            lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
1981            lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
1982        else:
1983            lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
1984            lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
1985
1986# lonlatlimits
1987##
1988
1989    if lonlatlims == 'None':
1990        lonlatlimsv = None
1991    else:
1992        lonlatlimsv = np.zeros((4), dtype=np.float)
1993        lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
1994        lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
1995        lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
1996        lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
1997
1998# lon/lat objects
1999##
2000    objnc = NetCDFFile(mapref.split('#')[0])
2001    lonobj = objnc.variables[mapref.split('#')[1]]
2002    latobj = objnc.variables[mapref.split('#')[2]]
2003
2004# map
2005##
2006    if mapkind == 'None':
2007        mapkindv = None
2008    else:
2009        mapkindv = mapkind
2010
2011    if observations is None:
2012        obsname = None
2013    else:
2014        obsfile = observations.split(',')[0]
2015        obsname = observations.split(',')[1]
2016        Tint = observations.split(',')[2]
2017        null = np.float(observations.split(',')[3])
2018        print '    observational trajectory:',obsfile
2019
2020        if not os.path.isfile(obsfile):
2021            print errormsg
2022            print '  ' + fname + ": observations trajectory file '" + obsfile +      \
2023              "' does not exist !!"
2024            quit(-1)
2025
2026        objfile = open(obsfile, 'r')
2027        obstrjtimev = []
2028        obstrjxv = []
2029        obstrjyv = []
2030
2031        for line in objfile:
2032            if line[0:1] != '#':
2033                lon = np.float(line.split(' ')[2])
2034                lat = np.float(line.split(' ')[1])
2035                if not lon == null and not lat == null:
2036                    obstrjtimev.append(int(line.split(' ')[0]))
2037                    obstrjxv.append(lon)
2038                    obstrjyv.append(lat)
2039                else:
2040                    obstrjtimev.append(int(line.split(' ')[0]))
2041                    obstrjxv.append(None)
2042                    obstrjyv.append(None)
2043
2044        objfile.close()
2045
2046        if Tint != '-1':
2047            Tint = int(observations.split(',')[2].split('@')[0])
2048            Tend = int(observations.split(',')[2].split('@')[1])
2049            lontrjvalues.append(obstrjxv[Tint:Tend+1])
2050            lattrjvalues.append(obstrjyv[Tint:Tend+1])
2051        else:
2052            lontrjvalues.append(obstrjxv[:])
2053            lattrjvalues.append(obstrjyv[:])
2054
2055    drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','),          \
2056      lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2057
2058    objnc.close()
2059
2060    return
2061
2062def draw_vals_trajectories(ncfile, values, variable):
2063    """ Function to draw values from the outputs from 'compute_tevolboxtraj'
2064    draw_vals_trajectories(ncfile, values, variable)
2065    ncfile= [ncfile] ',' list of files to use
2066    values= [statistics]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2067      [statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2068        'mean2', 'stdev'
2069      [Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2070      [labels]: ',' separated list of labels for the legend
2071      [locleg]: location of the legend (-1, autmoatic)
2072        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2073        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2074        9: 'upper center', 10: 'center'
2075      [gvarname]: name of the variable to appear in the graph
2076      [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
2077      [tkind]: kind of time to appear in the graph (assumed x-axis)
2078        'Nval': according to a given number of values as 'Nval',[Nval]
2079        'exct': according to an exact time unit as 'exct',[tunit];
2080          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2081            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2082            'l': milisecond
2083      [timefmt]: format of the time labels
2084      [title]: title of the graphic ('|' for spaces)
2085      [gkind]: kind of graphical output
2086    variable= variable to use
2087    """
2088
2089    fname = 'draw_vals_trajectories'
2090
2091    if values == 'h':
2092        print fname + '_____________________________________________________________'
2093        print draw_vals_trajectories.__doc__
2094        quit()
2095
2096    sims = ncfile.split(',')
2097
2098    if len(values.split(':')) != 9:
2099        print errormsg
2100        print '  ' + fname  + ': wrong number of values!', len(values.split(':')),   \
2101          'given 9 needed!!'
2102        print '    ',values.split(':')
2103        quit(-1)
2104
2105    statistics = values.split(':')[0]
2106    Tint = values.split(':')[1]
2107    labels = values.split(':')[2]
2108    gvarname = values.split(':')[3]
2109    timetit = values.split(':')[4].replace('|',' ')
2110    tkind = values.split(':')[5]
2111    timefmt = values.split(':')[6]
2112    title = values.split(':')[7].replace('|',' ')
2113    gkind = values.split(':')[8]
2114
2115    leglabels = labels.split('@')[0].split(',')
2116    locleg = int(labels.split('@')[1])
2117
2118    Nsims = len(sims)
2119
2120    if Tint != '-1':
2121        tini = np.float(Tint.split('@')[0])
2122        tend = np.float(Tint.split('@')[1])
2123    else:
2124        tini = -1.
2125        tend = -1.
2126
2127    vartimetrjv = []
2128
2129    print '  ' + fname
2130    for trjfile in sims:
2131        print '    ' + trjfile + ' ...'
2132        if not os.path.isfile(trjfile):
2133            print errormsg
2134            print '  ' + fname + ": trajectory file: '" + trjfile +                  \
2135              "' does not exist !!"
2136            quit(-1)
2137
2138        trjobj = NetCDFFile(trjfile, 'r')
2139        otim = trjobj.variables['time']
2140        if not trjobj.variables.has_key(statistics + 'box_' + variable):
2141            print errormsg
2142            print '  ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2143              statistics + 'box_' + variable + "' !!"
2144            quit(-1)
2145        ovar = trjobj.variables[statistics + 'box_' + variable]
2146        dimt = otim.shape[0]
2147
2148        if trjfile == sims[0]:
2149            gunits = ovar.getncattr('units')
2150            lname = ovar.getncattr('long_name')
2151            tunits = otim.getncattr('units')
2152
2153        if tini != -1:
2154            tiniid = -1
2155            tendid = -1       
2156            for itv in range(dimt):
2157                if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2158                if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2159
2160            if tiniid == -1 or tendid == -1:
2161                print errormsg
2162                print '  ' + main + ' time interval ', tini,',',tend,' not found: ',     \
2163                  tendid, ',', tiniid, ' !!'
2164                print '    data interval [',otim[0], otim[dimt-1],']'
2165                quit(-1)
2166            dimt = tendid - tiniid + 1
2167
2168        else:
2169            dimt = otim.shape[0]
2170
2171        valsv = np.zeros((2,dimt), dtype=np.float)
2172
2173        if tini == -1:
2174            valsv[1,:] = otim[:]
2175            valsv[0,:] = ovar[:]
2176        else:
2177            valsv[1,:] = otim[tiniid:tendid+1]
2178            valsv[0,:] = ovar[tiniid:tendid+1]
2179
2180        vartimetrjv.append(valsv)
2181        trjobj.close()
2182
2183    drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits,               \
2184      'val_trajectories_' + statistics, gvarname, timetit, tkind, timefmt, title,    \
2185      leglabels, locleg, gkind)
2186
2187def variable_values(values):
2188    """ Function to give back values for a given variable
2189      values= [varname] name of the variable
2190    """
2191
2192    fname = 'variable_values'
2193
2194    values = drw.variables_values(values)
2195
2196    print fname,'values:',values
2197    print fname,'variable_name:',values[0]
2198    print fname,'standard_name:',values[1]
2199    print fname,'min,max:',str(values[2]) + ',' + str(values[3])
2200    print fname,'long_name:',values[4]
2201    print fname,'units:',values[5]
2202    print fname,'shad_colors:',values[6]
2203    print fname,'all_values:',drw.numVector_String(values,',')
2204
2205    return
2206
2207####### ###### ##### #### ### ## #
2208
2209ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
2210
2211### Options
2212##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
2213string_operation="""operation to make:
2214  draw_topo_geogrid, draws topography from a WPS geo_em.d[nn].nc: -S [minTopo],[maxTopo]:[SW_lon],[SW_lat],[NE_lon],[NE_lat]:[title]:[graphic_kind]:[projection],[res_coastline]
2215  draw_2D_shad_cont, draws two 2D fields, first with shading second with contour lines: -v [varns],[varnc] -S [vnamefs],[vnamefc],[dimxvn],[dimyvn],[colorbar],[ckind],[clabfmt],[sminv]:[smaxv],[sminc]:[smaxv]:[Nlev],[figt],[kindfig],[reverse]
2216    [ckind]:
2217      'cmap': as it gets from colorbar
2218      'fixc,[colname]': fixed color [colname], all stright lines
2219      'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
2220"""
2221
2222#print string_operation
2223
2224parser = OptionParser()
2225parser.add_option("-f", "--netCDF_file", dest="ncfile", 
2226                  help="file to use", metavar="FILE")
2227parser.add_option("-o", "--operation", type='choice', dest="operation", 
2228       choices=namegraphics, 
2229                  help="operation to make: " + ngraphics, metavar="OPER")
2230parser.add_option("-S", "--valueS", dest="values", 
2231                  help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
2232parser.add_option("-v", "--variable", dest="varname",
2233                  help="[WHEN APPLICABLE] variable to check", metavar="VAR")
2234
2235(opts, args) = parser.parse_args()
2236
2237#######    #######
2238## MAIN
2239    #######
2240
2241# Not checking file operation
2242Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time',                    \
2243  'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines',                       \
2244  'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
2245  'variable_values']
2246
2247####### ###### ##### #### ### ## #
2248errormsg='ERROR -- error -- ERROR -- error'
2249
2250varn=opts.varname
2251oper=opts.operation
2252
2253if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and                   \
2254  not drw.searchInlist(Notcheckingfile, oper):
2255    print errormsg
2256    print '  ' + main + ': File ' + opts.ncfile + ' does not exist !!'
2257    quit(-1)
2258
2259if oper == 'create_movie':
2260    create_movie(opts.ncfile, opts.values, opts.varname)
2261elif oper == 'draw_2D_shad':
2262    draw_2D_shad(opts.ncfile, opts.values, opts.varname)
2263elif oper == 'draw_2D_shad_time':
2264    draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
2265elif oper == 'draw_2D_shad_cont':
2266    draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
2267elif oper == 'draw_2D_shad_cont_time':
2268    draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
2269elif oper == 'draw_2D_shad_line':
2270    draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
2271elif oper == 'draw_2D_shad_line_time':
2272    draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
2273elif oper == 'draw_Neighbourghood_evol':
2274    draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
2275elif oper == 'draw_lines':
2276    draw_lines(opts.ncfile, opts.values, opts.varname)
2277elif oper == 'draw_timeSeries':
2278    draw_timeSeries(opts.ncfile, opts.values, opts.varname)
2279elif oper == 'draw_topo_geogrid':
2280    draw_topo_geogrid(opts.ncfile, opts.values)
2281elif oper == 'draw_topo_geogrid_boxes':
2282    draw_topo_geogrid_boxes(opts.ncfile, opts.values)
2283elif oper == 'draw_trajectories':
2284    draw_trajectories(opts.ncfile, opts.values, opts.varname)
2285elif oper == 'draw_vals_trajectories':
2286    draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
2287elif oper == 'list_graphics':
2288# From: http://www.diveintopython.net/power_of_introspection/all_together.html
2289    import drawing as myself
2290    object = myself
2291    for opern in namegraphics:
2292        if  opern != 'list_graphics': 
2293            print opern + '_______ ______ _____ ____ ___ __ _'
2294            print getattr(object, opern).__doc__
2295elif oper == 'variable_values':
2296    variable_values(opts.values)
2297else:
2298    print errormsg
2299    print '  ' + main + ": the graphic '" + oper + "' is not ready !!"
2300    print errormsg
2301    quit()
Note: See TracBrowser for help on using the repository browser.