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

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

Adding in `draw_points':
[kindfigure]: kind of figure

'legend': only points in the map with the legend with the names
'labelled',[txtsize],[txtcol]: points with the names and size, color of text
File size: 109.4 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## e.g. # 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## e.g. # 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## e.g. # 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## e.g. # 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## e.g. # 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## e.g. # drawing.py -f 'geo_em.d02.nc' -o draw_topo_geogrid -S '0.,3000.:None:FF_3dom d02:png:cyl,i'
15## e.g. # 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## e.g. # 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## e.g. # 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## e.g. # 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## e.g. # drawing.py -o variable_values -S PSFC
20## e.g. # 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## e.g. # 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## e.g. # drawing.py -o draw_lines_time -f wrfout_d01_1980-03-01_00:00:00_Time_B0-E48-I1_south_north_B15-E15-I1_west_east_B15-E15-I1.nc -S 'time;y;time ([DD]${[HH]}$);file1;tas;evolution;time|hours!since!1949-12-01_00:00:00|exct,12,h|%d$^{%H}$;pdf' -v T2
23## e.g. # drawing.py -o draw_barbs -f ERAI_pl199501_131-132.nc -S 'X|lon|lon|-1,Y|lat|lat|-1,Z|lev|lev|4,T|time|time|0:auto,auto,auto:wind,ms-1:cyl,c:ERA-Interim|winds|at|1000|hPa|on|1996|January|1st|at|00|UTC:pdf:ERAI_pl199501_131-132' -v var131,var132
24## e.g. # ~/etudes/WRF_LMDZ/svn/LMDZ_WRF/tools/drawing.py -o draw_points -S 'tslist.dat,#,3,2,1:SuperStorm|sfc|stations:cyl,i:legend:auto:None:0:png:stations_loc' -f $HOME'/etudes/domains/WRFdynamicoSENS_SuperStorm/geo_em.d02.nc,XLONG_M,XLAT_M,HGT_M,Time|0,height,0.,3000.,terrain,m'
25## e.g. # drawing.py -o draw_points -S 'tslist.dat,#,3,2,1:SuperStorm|sfc|stations:cyl,i:labelled,8,black:auto:None:0:png:stations_loc' -f $HOME'/etudes/domains/WRFdynamicoSENS_SuperStorm/geo_em_west_east_B25-E180-I1_south_north_B160-E262-I1.nc,XLONG_M,XLAT_M,HGT_M,Time|0,height,0.,1500.,terrain,m'
26main = 'drawing.py'
27
28errormsg = 'ERROR -- error -- ERROR -- error'
29warnmsg = 'WARNING -- waring -- WARNING -- warning'
30fillValue=1.e20
31
32namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time',                 \
33  'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line',                \
34  'draw_2D_shad_line_time', 'draw_barbs', 'draw_points', 'draw_timeSeries',          \
35   'draw_topo_geogrid',                                                              \
36  'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
37  'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', 'list_graphics',      \
38  'variable_values']
39
40def draw_2D_shad(ncfile, values, varn):
41    """ plotting a fields with shading
42    draw_2D_shad(ncfile, values, varn)
43      ncfile= file to use
44      values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
45       [kindfig]:[reverse]:[mapv]:[close]
46        [vnamefs]: Name in the figure of the variable to be shaded
47        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
48          variable a given value is required (-1, all the length)
49        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
50        [colorbar]: name of the color bar
51        [smin/axv]: minimum and maximum value for the shading or:
52          'Srange': for full range
53          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
54          'Saroundminmax@val': for min*val,max*val
55          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
56            percentile_(100-val)-median)
57          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
58          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
59          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
60             percentile_(100-val)-median)
61        [figt]: title of the figure ('|' for spaces)
62        [kindfig]: kind of figure
63        [reverse]: Transformation of the values
64          * 'transpose': reverse the axes (x-->y, y-->x)
65          * 'flip'@[x/y]: flip the axis x or y
66        [mapv]: map characteristics: [proj],[res]
67          see full documentation: http://matplotlib.org/basemap/
68          [proj]: projection
69            * 'cyl', cilindric
70            * 'lcc', lamvbert conformal
71          [res]: resolution:
72            * 'c', crude
73            * 'l', low
74            * 'i', intermediate
75            * 'h', high
76            * 'f', full
77      valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
78      varn= [varsn] name of the variable to plot with shading
79    """
80
81    fname = 'draw_2D_shad'
82    if values == 'h':
83        print fname + '_____________________________________________________________'
84        print draw_2D_shad.__doc__
85        quit()
86
87    expectargs = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]',    \
88      '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[mapv]', '[close]']
89 
90    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
91
92    vnamesfig = values.split(':')[0]
93    dimvals= values.split(':')[1].replace('|',':')
94    vdimxn = values.split(':')[2]
95    vdimyn = values.split(':')[3]
96    colbarn = values.split(':')[4]
97    shadminmax = values.split(':')[5]
98    figtitle = values.split(':')[6].replace('|',' ')
99    figkind = values.split(':')[7]
100    revals = values.split(':')[8]
101    mapvalue = values.split(':')[9]
102#    varn = values.split(':')[10]
103
104    ncfiles = ncfile
105   
106    if not os.path.isfile(ncfiles):
107        print errormsg
108        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
109        quit(-1)   
110
111    objsf = NetCDFFile(ncfiles, 'r')
112   
113    varns = varn.split(',')[0]
114
115    if  not objsf.variables.has_key(varns):
116        print errormsg
117        print '  ' + fname + ': shading file "' + ncfiles +                          \
118          '" does not have variable "' +  varns + '" !!'
119        quit(-1)
120
121# Variables' values
122    objvars = objsf.variables[varns]
123
124    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
125
126# Dimensions names
127##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
128##    dimnamesv = []
129##    for idd in range(len(objvars.dimensions)):
130##        cutdim = False
131##        for idc in range(len(dimvals.split(','))):
132##            dimcutn = dimvals.split(',')[idc].split(':')[0]
133##            print objvars.dimensions[idd], dimcutn
134##            if objvars.dimensions[idd] == dimcutn:
135##                cutdim = True
136##                break
137##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
138    dimnamesv = [vdimyn, vdimxn]
139
140    if drw.searchInlist(objvars.ncattrs(),'units'):
141        varunits = objvars.getncattr('units')
142    else:
143        print warnmsg
144        print '  ' + fname + ": variable '" + varn + "' without units!!"
145        varunits = '-'
146
147    if  not objsf.variables.has_key(vdimxn):
148        print errormsg
149        print '  ' + fname + ': shading file "' + ncfiles +                          \
150          '" does not have dimension variable "' +  vdimxn + '" !!'
151        quit(-1)
152    if  not objsf.variables.has_key(vdimyn):
153        print errormsg
154        print '  ' + fname + ': shading file "' + ncfiles +                          \
155          '" does not have dimension variable "' +  vdimyn + '" !!'
156        quit(-1)
157
158    objdimx = objsf.variables[vdimxn]
159    objdimy = objsf.variables[vdimyn]
160    if drw.searchInlist(objdimx.ncattrs(),'units'):
161        odimxu = objdimx.getncattr('units')
162    else:
163        print warnmsg
164        print '  ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
165        odimxu = '-'
166
167    if drw.searchInlist(objdimy.ncattrs(),'units'):
168        odimyu = objdimy.getncattr('units')
169    else:
170        print warnmsg
171        print '  ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
172        odimyu = '-'
173
174    odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions,     \
175      objdimy.dimensions, dimvals.replace(':','|').split(','))
176
177
178#    if len(objdimx.shape) <= 2:
179##        odimxv = objdimx[valshad.shape]
180##        odimyv = objdimy[valshad.shape]
181#        odimxv = objdimx[:]
182#        odimyv = objdimy[:]
183
184#    elif len(objdimx.shape) == 3:
185##        dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
186##        odimxv = objdimx[tuple(dimcut)]
187##        odimyv = objdimy[tuple(dimcut)]
188#        odimxv = objdimx[0,:]
189#        odimyv = objdimy[0,:]
190#    else:
191#        print errormsg
192#        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
193#          ' not ready!!'
194#        quit(-1)
195
196    shading_nx = []
197    if shadminmax.split(',')[0][0:1] != 'S':
198            shading_nx.append(np.float(shadminmax.split(',')[0]))
199    else:
200        shading_nx.append(shadminmax.split(',')[0])
201
202    if shadminmax.split(',')[1][0:1] != 'S':
203        shading_nx.append(np.float(shadminmax.split(',')[1]))
204    else:
205        shading_nx.append(shadminmax.split(',')[1])
206
207    if mapvalue == 'None': mapvalue = None
208
209    drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, dimnamesv,\
210      colbarn, shading_nx, varunits, figtitle, figkind, revals, mapvalue, True)
211
212    return
213
214def draw_2D_shad_time(ncfile, values, varn):
215    """ plotting a fields with shading with time values
216    draw_2D_shad(ncfile, values, varn)
217      ncfile= file to use
218      values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
219       [kindfig]:[reverse]:[timevals]:[close]
220        [vnamefs]: Name in the figure of the variable to be shaded
221        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
222          variable a given value is required (-1, all the length)
223        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
224        [colorbar]: name of the color bar
225        [smin/axv]: minimum and maximum value for the shading or:
226          'Srange': for full range
227          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
228          'Saroundminmax@val': for min*val,max*val
229          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
230            percentile_(100-val)-median)
231          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
232          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
233          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
234             percentile_(100-val)-median)
235        [figt]: title of the figure ('|' for spaces)
236        [kindfig]: kind of figure
237        [reverse]: Transformation of the values
238          * 'transpose': reverse the axes (x-->y, y-->x)
239          * 'flip'@[x/y]: flip the axis x or y
240        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
241           [timen]; name of the time variable
242           [units]; units string according to CF conventions ([tunits] since
243             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
244           [kind]; kind of output
245             'Nval': according to a given number of values as 'Nval',[Nval]
246             'exct': according to an exact time unit as 'exct',[tunit];
247               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
248                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
249                'l': milisecond
250           [tfmt]; desired format
251           [label]; label at the graph ('!' for spaces)
252        [close]: should figure be closed (finished)
253      values='dtcon:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|'
254        'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True
255      varn= [varsn] name of the variable to plot with shading
256    """
257    fname = 'draw_2D_shad_time'
258
259    if values == 'h':
260        print fname + '_____________________________________________________________'
261        print draw_2D_shad_time.__doc__
262        quit()
263
264    farguments = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]',     \
265      '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', '[close]']
266    drw.check_arguments(fname,11,values,':',farguments)
267
268    vnamesfig = values.split(':')[0]
269    dimvals= values.split(':')[1].replace('|',':')
270    vdimxn = values.split(':')[2]
271    vdimyn = values.split(':')[3]
272    colbarn = values.split(':')[4]
273    shadminmax = values.split(':')[5]
274    figtitle = values.split(':')[6].replace('|',' ')
275    figkind = values.split(':')[7]
276    revals = values.split(':')[8]
277    timevals = values.split(':')[9]
278    close = values.split(':')[10]
279
280    ncfiles = ncfile
281   
282    if not os.path.isfile(ncfiles):
283        print errormsg
284        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
285        quit(-1)   
286
287    objsf = NetCDFFile(ncfiles, 'r')
288   
289    varns = varn.split(',')[0]
290
291    if  not objsf.variables.has_key(varns):
292        print errormsg
293        print '  ' + fname + ': shading file "' + ncfiles +                          \
294          '" does not have variable "' +  varns + '" !!'
295        quit(-1)
296
297# Variables' values
298    objvars = objsf.variables[varns]
299
300    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
301
302    dimnamesv = [vdimyn, vdimxn]
303
304    varunits = objvars.getncattr('units')
305
306    if  not objsf.variables.has_key(vdimxn):
307        print errormsg
308        print '  ' + fname + ': shading file "' + ncfiles +                          \
309          '" does not have dimension variable "' +  vdimxn + '" !!'
310        quit(-1)
311    if  not objsf.variables.has_key(vdimyn):
312        print errormsg
313        print '  ' + fname + ': shading file "' + ncfiles +                          \
314          '" does not have dimensino variable "' +  vdimyn + '" !!'
315        quit(-1)
316
317    objdimx = objsf.variables[vdimxn]
318    objdimy = objsf.variables[vdimyn]
319    odimxu = objdimx.getncattr('units')
320    odimyu = objdimy.getncattr('units')
321
322    if len(objdimx.shape) <= 2:
323        odimxv = objdimx[:]
324        odimyv = objdimy[:]
325
326    elif len(objdimx.shape) == 3:
327        odimxv = objdimx[0,:]
328        odimyv = objdimy[0,:]
329    else:
330        print errormsg
331        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
332          ' not ready!!'
333        quit(-1)
334
335    timename = timevals.split('|')[0]
336    timeunit = timevals.split('|')[1].replace('!',' ')
337    timekind = timevals.split('|')[2]
338    timefmt = timevals.split('|')[3]
339    timelabel = timevals.split('|')[4].replace('!',' ')
340
341    if vdimxn == timename:
342        odimxv = objsf.variables[vdimxn][:]
343        odimxu = timelabel
344        timeaxis = 'x'
345        odimyv = objsf.variables[vdimyn]
346        odimyu = odimyv.getncattr('units')
347        timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
348    elif vdimyn == timename:
349        odimyv = objsf.variables[vdimyn][:]
350        odimyu = timelabel
351        timeaxis = 'y'
352        odimxv = objsf.variables[vdimxn]
353        odimxu = odimxv.getncattr('units')
354        timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
355    else:
356        print errormsg
357        print '  ' + fname + ": time variable '" + timename + "' not found!!"
358        quit(-1)
359
360    shading_nx = []
361    if shadminmax.split(',')[0][0:1] != 'S':
362        shading_nx.append(np.float(shadminmax.split(',')[0]))
363    else:
364        shading_nx.append(shadminmax.split(',')[0])
365
366    if shadminmax.split(',')[1][0:1] != 'S':
367        shading_nx.append(np.float(shadminmax.split(',')[1]))
368    else:
369        shading_nx.append(shadminmax.split(',')[1])
370
371    closeval = drw.Str_Bool(close)
372
373    drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu,      \
374      dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
375      timepos, timelabels, closeval)
376
377    return
378
379def draw_2D_shad_cont(ncfile, values, varn):
380    """ plotting two fields, one with shading and the other with contour lines
381    draw_2D_shad_cont(ncfile, values, varn)
382      ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
383      values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]
384        [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
385        [dimvals/c]: list of [dimname]|[value] telling at which dimension of the
386          variable a given value is required (no dimension name, all the length)
387        [dimx/yvn]: names of the variables with the values of the dimensions for the plot
388        [colorbar]: name of the color bar
389        [ckind]: kind of contours
390          'cmap': as it gets from colorbar
391          'fixc,[colname]': fixed color [colname], all stright lines
392          'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
393        [clabfmt]: format of the labels in the contour (None, also possible)
394        [smin/axv]: minimum and maximum value for the shading
395        [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
396        [figt]: title of the figure ('|' for spaces)
397        [kindfig]: kind of figure
398        [reverse]: does the values be transposed? 'True/False',
399        [mapv]: map characteristics: [proj],[res]
400          see full documentation: http://matplotlib.org/basemap/
401          [proj]: projection
402            * 'cyl', cilindric
403            * 'lcc', lamvbert conformal
404          [res]: resolution:
405            * 'c', crude
406            * 'l', low
407            * 'i', intermediate
408            * 'h', high
409            * 'f', full
410      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'
411      varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
412    """
413
414    fname = 'draw_2D_shad_cont'
415    if values == 'h':
416        print fname + '_____________________________________________________________'
417        print draw_2D_shad_cont.__doc__
418        quit()
419
420    expectargs = '[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:'
421    expectargs = expectargs + '[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],'
422    expectargs = expectargs + '[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]'
423 
424    drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
425
426    vnamesfig = values.split(':')[0].split(',')
427    dimvals= values.split(':')[1].replace('|',':')
428    dimvalc= values.split(':')[2].replace('|',':')
429    vdimxn = values.split(':')[3]
430    vdimyn = values.split(':')[4]
431    colbarn = values.split(':')[5]
432    countkind = values.split(':')[6]
433    countlabelfmt = values.split(':')[7]
434    shadminmax = values.split(':')[8]
435    contlevels = values.split(':')[9]
436    figtitle = values.split(':')[10].replace('|',' ')
437    figkind = values.split(':')[11]
438    revals = values.split(':')[12]
439    mapvalue = values.split(':')[13]
440
441    if2filenames = ncfile.find(',')
442
443    if if2filenames != -1:
444        ncfiles = ncfile.split(',')[0]
445        ncfilec = ncfile.split(',')[1]
446    else:
447        ncfiles = ncfile
448        ncfilec = ncfile
449
450    if not os.path.isfile(ncfiles):
451        print errormsg
452        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
453        quit(-1)   
454
455    if not os.path.isfile(ncfilec):
456        print errormsg
457        print '  ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
458        quit(-1)   
459
460    objsf = NetCDFFile(ncfiles, 'r')
461    objcf = NetCDFFile(ncfilec, 'r')
462   
463    varns = varn.split(',')[0]
464    varnc = varn.split(',')[1]
465
466    if  not objsf.variables.has_key(varns):
467        print errormsg
468        print '  ' + fname + ': shading file "' + ncfiles +                          \
469          '" does not have variable "' +  varns + '" !!'
470        quit(-1)
471
472    if  not objcf.variables.has_key(varnc):
473        print errormsg
474        print '  ' + fname + ': contour file "' + ncfilec +                          \
475          '" does not have variable "' +  varnc + '" !!'
476        quit(-1)
477
478# Variables' values
479    objvars = objsf.variables[varns]
480    objvarc = objcf.variables[varnc]
481
482    if len(objvars.shape) != len(objvarc.shape):
483        print errormsg
484        print '  ' + fname + ': shading variable "' + varns + '" has a shape: ',     \
485          objvars.shape,  'different than contour variable "' +  varnc + '": ',      \
486          objvarc.shape,' !!!'
487        quit(-1)
488
489    for idim in range(len(objvars.shape)):
490        if objvars.shape[idim] != objvarc.shape[idim]:
491            print errormsg
492            print '  ' + fname + ': shading variable "' + varns + '" has a shape: ', \
493              objvars.shape,  'different than contour variable "' +  varnc + '": ',  \
494              objvarc.shape,' !!!'
495            quit(-1)
496
497    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
498    valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|'))
499
500# Dimensions names
501##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
502##    dimnamesv = []
503##    for idd in range(len(objvars.dimensions)):
504##        cutdim = False
505##        for idc in range(len(dimvals.split(','))):
506##            dimcutn = dimvals.split(',')[idc].split(':')[0]
507##            print objvars.dimensions[idd], dimcutn
508##            if objvars.dimensions[idd] == dimcutn:
509##                cutdim = True
510##                break
511##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
512    dimnamesv = [vdimyn, vdimxn]
513
514    varunits = []
515    varunits.append(objvars.getncattr('units'))
516    varunits.append(objvarc.getncattr('units'))
517
518    if  not objsf.variables.has_key(vdimxn):
519        print errormsg
520        print '  ' + fname + ': shading file "' + ncfiles +                          \
521          '" does not have dimension variable "' +  vdimxn + '" !!'
522        quit(-1)
523    if  not objsf.variables.has_key(vdimyn):
524        print errormsg
525        print '  ' + fname + ': shading file "' + ncfiles +                          \
526          '" does not have dimensino variable "' +  vdimyn + '" !!'
527        quit(-1)
528
529    objdimx = objsf.variables[vdimxn]
530    objdimy = objsf.variables[vdimyn]
531    odimxu = objdimx.getncattr('units')
532    odimyu = objdimy.getncattr('units')
533
534# Getting only that dimensions with coincident names
535    odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions,     \
536      objdimy.dimensions, dimvals.replace(':','|').split(','))
537
538#    dimnvx = objdimx.dimensions
539#    cutslice = []
540#    for idimn in objdimx.dimensions:
541#        found = False
542#        for dimsn in dimsshad:
543#            if idimn == dimsn:
544#                cutslice.append(slice(0,len(objsf.dimensions[idimn])))
545#                found = True
546#        if not found: cutslice.append(0)
547#
548#    odimxv = objdimx[tuple(cutslice)]
549#
550#    dimnvy = objdimy.dimensions
551#    cutslice = []
552#    for idimn in objdimy.dimensions:
553#        found = False
554#        for dimsn in dimsshad:
555#            if idimn == dimsn:
556#                cutslice.append(slice(0,len(objsf.dimensions[idimn])))
557#                found = True
558#        if not found: cutslice.append(0)
559#
560#    odimyv = objdimy[tuple(cutslice)]
561
562#    if len(objdimx.shape) <= 2:
563#        odimxv = objdimx[:]
564#        odimyv = objdimy[:]
565#    elif len(objdimx.shape) == 3:
566#        odimxv = objdimx[0,:]
567#        odimyv = objdimy[0,:]
568#    else:
569#        print errormsg
570#        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
571#          ' not ready!!'
572#        quit(-1)
573
574    if countlabelfmt == 'None': 
575        countlfmt = None
576    else:
577        countlfmt = countlabelfmt
578
579    shading_nx = np.zeros((2), dtype=np.float)
580    shading_nx[0] = np.float(shadminmax.split(',')[0])
581    shading_nx[1] = np.float(shadminmax.split(',')[1])
582
583    clevmin = np.float(contlevels.split(',')[0])
584    clevmax = np.float(contlevels.split(',')[1])
585    Nclevels = int(contlevels.split(',')[2])
586
587    levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
588
589    if len(levels_cont) <= 1: 
590        print warnmsg
591        print '  ' + fname + ': wrong contour levels:', levels_cont,' !!'
592        del(levels_cont)
593        levels_cont = np.zeros((Nclevels), dtype=np.float)
594        levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
595        print '    generating default ones: ',levels_cont
596
597    if mapvalue == 'None': mapvalue = None
598
599    drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu,  \
600      odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont,     \
601      varunits, figtitle, figkind, revals, mapvalue)
602
603    return
604
605def draw_2D_shad_cont_time(ncfile, values, varn):
606    """ plotting two fields, one with shading and the other with contour lines being
607    one of the dimensions of time characteristics
608    draw_2D_shad_cont(ncfile, values, varn)
609      ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
610      values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv]
611        [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
612        [dimvals/c]: list of [dimname]|[value] telling at which dimension of the
613          variable a given value is required (no dimension name, all the length)
614        [dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
615        [colorbar]: name of the color bar
616        [ckind]: kind of contours
617          'cmap': as it gets from colorbar
618          'fixc,[colname]': fixed color [colname], all stright lines
619          'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
620        [clabfmt]: format of the labels in the contour (None, also possible)
621        [smin/axv]: minimum and maximum value for the shading
622        [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
623        [figt]: title of the figure ('|' for spaces)
624        [kindfig]: kind of figure
625        [reverse]: modification to the dimensions:
626          'transposed': transpose matrices
627          'flip',[x/y]: flip only the dimension [x] or [y]
628        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label] time labels characteristics
629           [timen]; name of the time variable
630           [units]; units string according to CF conventions ([tunits] since
631             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
632           [kind]; kind of output
633             'Nval': according to a given number of values as 'Nval',[Nval]
634             'exct': according to an exact time unit as 'exct',[tunit];
635               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
636                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
637                'l': milisecond
638           [tfmt]; desired format
639           [label]; label at the graph ('!' for spaces)
640        [mapv]: map characteristics: [proj],[res]
641          see full documentation: http://matplotlib.org/basemap/
642          [proj]: projection
643            * 'cyl', cilindric
644            * 'lcc', lamvbert conformal
645          [res]: resolution:
646            * 'c', crude
647            * 'l', low
648            * 'i', intermediate
649            * 'h', high
650            * 'f', full
651      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'
652      varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
653    """
654
655    fname = 'draw_2D_shad_cont_time'
656    if values == 'h':
657        print fname + '_____________________________________________________________'
658        print draw_2D_shad_cont_time.__doc__
659        quit()
660
661    expectargs = ['[vnamefs]', '[dimvals]', '[dimvalc]', '[dimxvn]', '[dimyvn]',     \
662      '[colorbar]', '[ckind]', '[clabfmt]', '[sminv],[smaxv]',                       \
663      '[sminc],[smaxv],[Nlev]', '[figt]', '[kindfig]', '[reverse]', '[timevals]',    \
664      '[mapv]']
665 
666    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
667
668    vnamesfig = values.split(':')[0].split(',')
669    dimvals= values.split(':')[1].replace('|',':')
670    dimvalc= values.split(':')[2].replace('|',':')
671    vdimxn = values.split(':')[3]
672    vdimyn = values.split(':')[4]
673    colbarn = values.split(':')[5]
674    countkind = values.split(':')[6]
675    countlabelfmt = values.split(':')[7]
676    shadminmax = values.split(':')[8]
677    contlevels = values.split(':')[9]
678    figtitle = values.split(':')[10].replace('|',' ')
679    figkind = values.split(':')[11]
680    revals = values.split(':')[12]
681    timevals = values.split(':')[13]
682    mapvalue = values.split(':')[14]
683
684    if2filenames = ncfile.find(',')
685
686    if if2filenames != -1:
687        ncfiles = ncfile.split(',')[0]
688        ncfilec = ncfile.split(',')[1]
689    else:
690        ncfiles = ncfile
691        ncfilec = ncfile
692
693    if not os.path.isfile(ncfiles):
694        print errormsg
695        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
696        quit(-1)   
697
698    if not os.path.isfile(ncfilec):
699        print errormsg
700        print '  ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
701        quit(-1)   
702
703    objsf = NetCDFFile(ncfiles, 'r')
704    objcf = NetCDFFile(ncfilec, 'r')
705   
706    varns = varn.split(',')[0]
707    varnc = varn.split(',')[1]
708
709    if  not objsf.variables.has_key(varns):
710        print errormsg
711        print '  ' + fname + ': shading file "' + ncfiles +                          \
712          '" does not have variable "' +  varns + '" !!'
713        quit(-1)
714
715    if  not objcf.variables.has_key(varnc):
716        print errormsg
717        print '  ' + fname + ': contour file "' + ncfilec +                          \
718          '" does not have variable "' +  varnc + '" !!'
719        quit(-1)
720
721# Variables' values
722    objvars = objsf.variables[varns]
723    objvarc = objcf.variables[varnc]
724
725    if len(objvars.shape) != len(objvarc.shape):
726        print errormsg
727        print '  ' + fname + ': shading variable "' + varns + '" has a shape: ',     \
728          objvars.shape,  'different than contour variable "' +  varnc + '": ',      \
729          objvarc.shape,' !!!'
730        quit(-1)
731
732    for idim in range(len(objvars.shape)):
733        if objvars.shape[idim] != objvarc.shape[idim]:
734            print errormsg
735            print '  ' + fname + ': shading variable "' + varns + '" has a shape: ', \
736              objvars.shape,  'different than contour variable "' +  varnc + '": ',  \
737              objvarc.shape,' !!!'
738            quit(-1)
739
740    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
741    valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|'))
742
743# Dimensions names
744##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
745##    dimnamesv = []
746##    for idd in range(len(objvars.dimensions)):
747##        cutdim = False
748##        for idc in range(len(dimvals.split(','))):
749##            dimcutn = dimvals.split(',')[idc].split(':')[0]
750##            print objvars.dimensions[idd], dimcutn
751##            if objvars.dimensions[idd] == dimcutn:
752##                cutdim = True
753##                break
754##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
755    dimnamesv = [vdimyn, vdimxn]
756
757    varunits = []
758    varunits.append(objvars.getncattr('units'))
759    varunits.append(objvarc.getncattr('units'))
760
761    if  not objsf.variables.has_key(vdimxn):
762        print errormsg
763        print '  ' + fname + ': shading file "' + ncfiles +                          \
764          '" does not have dimension variable "' +  vdimxn + '" !!'
765        quit(-1)
766    if  not objsf.variables.has_key(vdimyn):
767        print errormsg
768        print '  ' + fname + ': shading file "' + ncfiles +                          \
769          '" does not have dimension variable "' +  vdimyn + '" !!'
770        quit(-1)
771
772    timename = timevals.split('|')[0]
773    timeunit = timevals.split('|')[1].replace('!',' ')
774    timekind = timevals.split('|')[2]
775    timefmt = timevals.split('|')[3]
776    timelabel = timevals.split('|')[4].replace('!',' ')
777
778    if vdimxn == timename:
779        timevals = objsf.variables[vdimxn][:]
780        timedims = objsf.variables[vdimxn].dimensions
781        dimt = 'x'
782        ovalaxis = objsf.variables[vdimyn]
783        ovalu = ovalaxis.getncattr('units')
784    elif vdimyn == timename:
785        timevals = objsf.variables[vdimyn][:]
786        timedims = objsf.variables[vdimyn].dimensions
787        dimt = 'y'
788        ovalaxis = objsf.variables[vdimxn]
789        ovalu = ovalaxis.getncattr('units')
790    else:
791        print errormsg
792        print '  ' + fname + ": time variable '" + timename + "' not found!!"
793        quit(-1)
794
795    timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt)
796
797# Getting only that dimensions with coincident names
798    dimnvx = ovalaxis.dimensions
799
800    cutslice = []
801    for idimn in dimsshad:
802        found = False
803        for dimsn in dimnvx:
804            if idimn == dimsn:
805                cutslice.append(slice(0,len(objsf.dimensions[idimn])))
806                found = True
807        if not found: cutslice.append(0)
808
809    ovalaxisv = ovalaxis[tuple(cutslice)]
810
811##    if len(ovalaxis.shape) <= 2:
812##        ovalaxisv = ovalaxis[:]
813
814##    elif len(ovalaxis.shape) == 3:
815##        ovalaxisv = ovalaxis[0,:]
816##    else:
817##        print errormsg
818##        print '  ' + fname + ': shape of dimension variable:', ovalaxis.shape,       \
819##          ' not ready!!'
820##        quit(-1)
821
822    if countlabelfmt == 'None': 
823        countlfmt = None
824    else:
825        countlfmt = countlabelfmt
826
827    shading_nx = np.zeros((2), dtype=np.float)
828    shading_nx[0] = np.float(shadminmax.split(',')[0])
829    shading_nx[1] = np.float(shadminmax.split(',')[1])
830
831    clevmin = np.float(contlevels.split(',')[0])
832    clevmax = np.float(contlevels.split(',')[1])
833    Nclevels = int(contlevels.split(',')[2])
834
835    levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
836
837    if len(levels_cont) <= 1: 
838        print warnmsg
839        print '  ' + fname + ': wrong contour levels:', levels_cont,' !!'
840        del(levels_cont)
841        levels_cont = np.zeros((Nclevels), dtype=np.float)
842        levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
843        print '    generating default ones: ',levels_cont
844
845    if mapvalue == 'None': mapvalue = None
846
847    drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv,         \
848      timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn,    \
849      countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind,   \
850      revals, mapvalue)
851
852    return
853
854def draw_2D_shad_line(ncfile, values, varn):
855    """ plotting a fields with shading and another with line
856    draw_2D_shad_line(ncfile, values, varn)
857      ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
858      values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
859       [kindfig]:[reverse]:[mapv]:[close]
860        [vnamefs]: Name in the figure of the variable to be shaded
861        [vnamefl]: Name in the figure of the variable to be lined
862        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
863          variable a given value is required (-1, all the length)
864        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
865        [colorbar]: name of the color bar
866        [colline]: name of the color for the line
867        [smin/axv]: minimum and maximum value for the shading
868        [figt]: title of the figure ('|' for spaces)
869        [kindfig]: kind of figure
870        [reverse]: Transformation of the values
871          * 'transpose': reverse the axes (x-->y, y-->x)
872          * 'flip'@[x/y]: flip the axis x or y
873        [mapv]: map characteristics: [proj],[res]
874          see full documentation: http://matplotlib.org/basemap/
875          [proj]: projection
876            * 'cyl', cilindric
877            * 'lcc', lamvbert conformal
878          [res]: resolution:
879            * 'c', crude
880            * 'l', low
881            * 'i', intermediate
882            * 'h', high
883            * 'f', full
884      valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
885      varn= [varsn],[varnl] name of the variable to plot with shading and with line
886    """
887
888    fname = 'draw_2D_shad_line'
889    if values == 'h':
890        print fname + '_____________________________________________________________'
891        print draw_2D_shad_line.__doc__
892        quit()
893
894    farguments = ['[vnamefs],[vnamefl]', '[dimvals]', '[dimxvn]', '[dimyvn]',        \
895      '[colorbar],[colline]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
896      '[mapv]', '[close]']
897    drw.check_arguments(fname,11,values,':',farguments)
898
899    vnamesfig = values.split(':')[0].split(',')[0]
900    dimvals= values.split(':')[1].replace('|',':')
901    vdimxn = values.split(':')[2]
902    vdimyn = values.split(':')[3]
903    colbarn = values.split(':')[4].split(',')[0]
904    shadminmax = values.split(':')[5]
905    figtitle = values.split(':')[6].replace('|',' ')
906    figkind = values.split(':')[7]
907    revals = values.split(':')[8]
908    mapvalue = values.split(':')[9]
909#    varn = values.split(':')[10]
910
911    ncfiles = ncfile.split(',')[0]
912   
913    if not os.path.isfile(ncfiles):
914        print errormsg
915        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
916        quit(-1)   
917
918    objsf = NetCDFFile(ncfiles, 'r')
919   
920    varns = varn.split(',')[0]
921
922    if  not objsf.variables.has_key(varns):
923        print errormsg
924        print '  ' + fname + ': shading file "' + ncfiles +                          \
925          '" does not have variable "' +  varns + '" !!'
926        quit(-1)
927
928# Variables' values
929    objvars = objsf.variables[varns]
930
931    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
932
933# Dimensions names
934##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
935##    dimnamesv = []
936##    for idd in range(len(objvars.dimensions)):
937##        cutdim = False
938##        for idc in range(len(dimvals.split(','))):
939##            dimcutn = dimvals.split(',')[idc].split(':')[0]
940##            print objvars.dimensions[idd], dimcutn
941##            if objvars.dimensions[idd] == dimcutn:
942##                cutdim = True
943##                break
944##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
945    dimnamesv = [vdimyn, vdimxn]
946
947    varunits = objvars.getncattr('units')
948
949    if  not objsf.variables.has_key(vdimxn):
950        print errormsg
951        print '  ' + fname + ': shading file "' + ncfiles +                          \
952          '" does not have dimension variable "' +  vdimxn + '" !!'
953        quit(-1)
954    if  not objsf.variables.has_key(vdimyn):
955        print errormsg
956        print '  ' + fname + ': shading file "' + ncfiles +                          \
957          '" does not have dimensino variable "' +  vdimyn + '" !!'
958        quit(-1)
959
960    objdimx = objsf.variables[vdimxn]
961    objdimy = objsf.variables[vdimyn]
962    odimxu = objdimx.getncattr('units')
963    odimyu = objdimy.getncattr('units')
964
965    if len(objdimx.shape) <= 2:
966#        odimxv = objdimx[valshad.shape]
967#        odimyv = objdimy[valshad.shape]
968        odimxv = objdimx[:]
969        odimyv = objdimy[:]
970
971    elif len(objdimx.shape) == 3:
972#        dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
973#        odimxv = objdimx[tuple(dimcut)]
974#        odimyv = objdimy[tuple(dimcut)]
975        odimxv = objdimx[0,:]
976        odimyv = objdimy[0,:]
977    else:
978        print errormsg
979        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
980          ' not ready!!'
981        quit(-1)
982
983    shading_nx = np.zeros((2), dtype=np.float)
984    shading_nx[0] = np.float(shadminmax.split(',')[0])
985    shading_nx[1] = np.float(shadminmax.split(',')[1])
986
987    if mapvalue == 'None': mapvalue = None
988
989# line plot
990##
991    ncfilel = ncfile.split(',')[1]
992    vnamelfig = values.split(':')[0].split(',')[1]
993    varnl = varn.split(',')[1]
994    colline = values.split(':')[4].split(',')[1]
995
996    objlf = NetCDFFile(ncfilel,'r')
997    objlvar = objlf.variables[varnl]
998
999    linevals = objlvar[:]
1000    varlunits = objlvar.units
1001
1002    drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
1003      odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits,  \
1004      figtitle, figkind, revals, mapvalue, True)
1005
1006    objsf.close()
1007    objlf.close()
1008
1009    return
1010
1011def draw_2D_shad_line_time(ncfile, values, varn):
1012    """ plotting a fields with shading and a line with time values
1013    draw_2D_shad_line(ncfile, values, varn)
1014      ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
1015      values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
1016       [kindfig]:[reverse]:[timevals]:[close]
1017        [vnamefs]: Name in the figure of the variable to be shaded
1018        [vnamefl]: Name in the figure of the variable to be lined
1019        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
1020          variable a given value is required (-1, all the length)
1021        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
1022        [colorbar]: name of the color bar
1023        [smin/axv]: minimum and maximum value for the shading
1024        [figt]: title of the figure ('|' for spaces)
1025        [kindfig]: kind of figure
1026        [reverse]: Transformation of the values
1027          * 'transpose': reverse the axes (x-->y, y-->x)
1028          * 'flip'@[x/y]: flip the axis x or y
1029        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
1030           [timen]; name of the time variable
1031           [units]; units string according to CF conventions ([tunits] since
1032             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1033           [kind]; kind of output
1034             'Nval': according to a given number of values as 'Nval',[Nval]
1035             'exct': according to an exact time unit as 'exct',[tunit];
1036               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1037                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1038                'l': milisecond
1039           [tfmt]; desired format
1040           [label]; label at the graph ('!' for spaces)
1041        [close]: should figure be closed (finished)
1042      values='dtcon,prc:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|'
1043        'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True
1044      varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
1045    """
1046    fname = 'draw_2D_shad_line_time'
1047    if values == 'h':
1048        print fname + '_____________________________________________________________'
1049        print draw_2D_shad__line_time.__doc__
1050        quit()
1051
1052    farguments = ['[vnamefs],[vanemefl]', '[dimvals]', '[dimxvn]', '[dimyvn]',       \
1053      '[colorbar]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]',           \
1054      '[timevals]', '[close]']
1055    drw.check_arguments(fname,11,values,':',farguments)
1056
1057    vnamesfig = values.split(':')[0].split(',')[0]
1058    dimvals= values.split(':')[1].replace('|',':')
1059    vdimxn = values.split(':')[2]
1060    vdimyn = values.split(':')[3]
1061    colbarn = values.split(':')[4]
1062    shadminmax = values.split(':')[5]
1063    figtitle = values.split(':')[6].replace('|',' ')
1064    figkind = values.split(':')[7]
1065    revals = values.split(':')[8]
1066    timevals = values.split(':')[9]
1067    close = values.split(':')[10]
1068
1069    ncfiles = ncfile.split(',')[0]
1070   
1071    if not os.path.isfile(ncfiles):
1072        print errormsg
1073        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1074        quit(-1)   
1075
1076    objsf = NetCDFFile(ncfiles, 'r')
1077   
1078    varns = varn.split(',')[0]
1079
1080    if  not objsf.variables.has_key(varns):
1081        print errormsg
1082        print '  ' + fname + ': shading file "' + ncfiles +                          \
1083          '" does not have variable "' +  varns + '" !!'
1084        quit(-1)
1085
1086# Variables' values
1087    objvars = objsf.variables[varns]
1088
1089    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
1090
1091    dimnamesv = [vdimyn, vdimxn]
1092
1093    varunits = objvars.getncattr('units')
1094
1095    if  not objsf.variables.has_key(vdimxn):
1096        print errormsg
1097        print '  ' + fname + ': shading file "' + ncfiles +                          \
1098          '" does not have dimension variable "' +  vdimxn + '" !!'
1099        quit(-1)
1100    if  not objsf.variables.has_key(vdimyn):
1101        print errormsg
1102        print '  ' + fname + ': shading file "' + ncfiles +                          \
1103          '" does not have dimensino variable "' +  vdimyn + '" !!'
1104        quit(-1)
1105
1106    objdimx = objsf.variables[vdimxn]
1107    objdimy = objsf.variables[vdimyn]
1108    odimxu = objdimx.getncattr('units')
1109    odimyu = objdimy.getncattr('units')
1110
1111    if len(objdimx.shape) <= 2:
1112        odimxv = objdimx[:]
1113        odimyv = objdimy[:]
1114
1115    elif len(objdimx.shape) == 3:
1116        odimxv = objdimx[0,:]
1117        odimyv = objdimy[0,:]
1118    else:
1119        print errormsg
1120        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
1121          ' not ready!!'
1122        quit(-1)
1123
1124    timename = timevals.split('|')[0]
1125    timeunit = timevals.split('|')[1].replace('!',' ')
1126    timekind = timevals.split('|')[2]
1127    timefmt = timevals.split('|')[3]
1128    timelabel = timevals.split('|')[4].replace('!',' ')
1129
1130    if vdimxn == timename:
1131        odimxv = objsf.variables[vdimxn][:]
1132        odimxu = timelabel
1133        timeaxis = 'x'
1134        odimyv = objsf.variables[vdimyn]
1135        odimyu = odimyv.getncattr('units')
1136        timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1137    elif vdimyn == timename:
1138        odimyv = objsf.variables[vdimyn][:]
1139        odimyu = timelabel
1140        timeaxis = 'y'
1141        odimxv = objsf.variables[vdimxn]
1142        odimxu = odimxv.getncattr('units')
1143        timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1144    else:
1145        print errormsg
1146        print '  ' + fname + ": time variable '" + timename + "' not found!!"
1147        quit(-1)
1148
1149    shading_nx = np.zeros((2), dtype=np.float)
1150    shading_nx[0] = np.float(shadminmax.split(',')[0])
1151    shading_nx[1] = np.float(shadminmax.split(',')[1])
1152
1153    closeval = drw.Str_Bool(close)
1154
1155    drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu,      \
1156      dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1157      timepos, timelabels, False)
1158
1159# Line values
1160##
1161    ncfilel = ncfile.split(',')[1]
1162
1163    vnamelfig = values.split(':')[0].split(',')[1]
1164    varnl = varn.split(',')[1]
1165
1166    objlf = NetCDFFile(ncfilel,'r')
1167    objlvar = objlf.variables[varnl]
1168
1169    linevals = objlvar[:]
1170    if reva0 == 'tranpose':
1171        plt.plot (linevals, odimxv, '-', color='k')
1172    else:
1173        plt.plot (odimxv, linevals, '-', color='k')
1174
1175    objsf.close()
1176    objsl.close()
1177
1178    return
1179
1180def draw_barbs(ncfile, values, varns):
1181    """ Function to plot wind barbs
1182      values= [dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:[mapvalues]:
1183        [gtit]:[kindfig]:[figuren]
1184        'X/Y/Z/T'|[dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
1185          [dimname]: name of the dimension in the file
1186          [vardimname]: name of the variable with the values for the dimension in the file
1187          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
1188          No value takes all the range of the dimension
1189        [vecvals]= [frequency],[color],[length]
1190          [frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
1191            'auto', computed automatically to have 20 vectors along each axis)
1192          [color]: color of the vectors ('auto', for 'red')
1193          [length]: length of the wind barbs ('auto', for 9)
1194        [windlabs]= [windname],[windunits]
1195          [windname]: name of the wind variable in the graph
1196          [windunits]: units of the wind variable in the graph ('None', for the value in the file)
1197        [mapvalues]= map characteristics: [proj],[res]
1198          see full documentation: http://matplotlib.org/basemap/
1199          [proj]: projection
1200            * 'cyl', cilindric
1201            * 'lcc', lambert conformal
1202          [res]: resolution:
1203            * 'c', crude
1204            * 'l', low
1205            * 'i', intermediate
1206            * 'h', high
1207            * 'f', full
1208        gtit= title of the graph ('|', for spaces)
1209        kindfig= kind of figure
1210        figuren= name of the figure
1211      ncfile= file to use
1212      varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
1213    """
1214    fname = 'draw_barbs'
1215
1216    if values == 'h':
1217        print fname + '_____________________________________________________________'
1218        print draw_barbs.__doc__
1219        quit()
1220
1221    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:' +  \
1222      '[mapvalues]:[gtit]:[kindfig]:[figuren]'
1223 
1224    drw.check_arguments(fname,len(expectargs.split(':')),values,':',                 \
1225      expectargs.split(':'))
1226
1227    dimvals = values.split(':')[0]
1228    vecvals = values.split(':')[1]
1229    windlabels = values.split(':')[2]
1230    mapvalues = values.split(':')[3]
1231    gtit = values.split(':')[4]
1232    kindfig = values.split(':')[5]
1233    figuren = values.split(':')[6]
1234
1235    of = NetCDFFile(ncfile,'r')
1236
1237    dims = {}
1238    for dimv in dimvals.split(','):
1239        dns = dimv.split('|')
1240        dims[dns[0]] = [dns[1], dns[2], dns[3]]
1241
1242    varNs = []
1243    for dn in dims.keys():
1244        if dn == 'X':
1245            varNs.append(dims[dn][1])
1246            dimx = len(of.dimensions[dims[dn][0]])
1247        elif dn == 'Y':
1248            varNs.append(dims[dn][1])
1249            dimy = len(of.dimensions[dims[dn][0]])
1250
1251    ivar = 0
1252    for wvar in varns.split(','):
1253        if not drw.searchInlist(of.variables.keys(), wvar):
1254            print errormsg
1255            print '  ' + fname + ": file does not have variable '" + wvar + "' !!"
1256            quit(-1)
1257        if ivar == 0:
1258            varNs.append(wvar)
1259        else:
1260            varNs.append(wvar)
1261
1262    ivar = 0
1263    for varN in varNs:
1264        varslice = []
1265
1266        ovarN = of.variables[varN]
1267        vard = ovarN.dimensions
1268        for vdn in vard:
1269            found = False
1270            for dd in dims.keys():
1271                if dims[dd][0] == vdn:
1272                    if dims[dd][2].find('@') != -1:
1273                        rvals = dims[dd][2].split('@')
1274                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
1275                    elif dims[dd][2] == '-1':
1276                        varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
1277                    else:
1278                        varslice.append(int(dims[dd][2]))
1279
1280                    found = True
1281                    break
1282            if not found:
1283                varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
1284
1285        if varN == dims['X'][1]:
1286            lonvals0 = np.squeeze(ovarN[tuple(varslice)])
1287        elif varN == dims['Y'][1]:
1288            latvals0 = np.squeeze(ovarN[tuple(varslice)])
1289        elif ivar == 2:
1290            uwvals = np.squeeze(np.array(ovarN[tuple(varslice)]))
1291        elif ivar == 3:
1292            vwvals = np.squeeze(ovarN[tuple(varslice)])
1293
1294        ivar = ivar + 1
1295
1296#    print 'Final shapes:',lonvals0.shape,':',latvals0.shape,':',uwvals.shape,':',
1297#      vwvals.shape
1298
1299    if len(uwvals.shape) != 2 or len(vwvals.shape) != 2:
1300        print errormsg
1301        print '  ' + fname + ': wrong size of the wind fields! they must be ' +      \
1302          '2-dimensional!'
1303        print '    u-winds shape:',uwvals.shape,'dims:',of.variables[varNs[2]]
1304        print '    v-winds shape:',vwvals.shape,'dims:',of.variables[varNs[3]]
1305        print '      provide more values for their dimensions!!'
1306        quit(-1)
1307
1308    if len(lonvals0.shape) == 1:
1309        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
1310    else:
1311        lonvals = lonvals0
1312        latvals = latvals0
1313
1314# Vecor values
1315    if vecvals.split(',')[0] == 'None':
1316        freqv = None
1317    else:
1318        freqv = vecvals.split(',')[0] 
1319    colorv = vecvals.split(',')[1]
1320    lengthv = vecvals.split(',')[2]
1321
1322# Vector labels
1323    windname = windlabels.split(',')[0]
1324    windunits = windlabels.split(',')[1]
1325
1326    drw.plot_barbs(lonvals, latvals, uwvals, vwvals, freqv, colorv, lengthv, 
1327      windname, windunits, mapvalues, gtit, kindfig, figuren)
1328
1329    return
1330 
1331def draw_topo_geogrid(ncfile, values):
1332    """ plotting geo_em.d[nn].nc topography from WPS files
1333    draw_topo_geogrid(ncfile, values)
1334      ncfile= geo_em.d[nn].nc file to use
1335      values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1336        [min/max]Topo: minimum and maximum values of topography to draw
1337        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1338        title: title of the graph
1339        graphic_kind: kind of figure (jpg, pdf, png)
1340        mapvalues: map characteristics [proj],[res]
1341          see full documentation: http://matplotlib.org/basemap/
1342          [proj]: projection
1343            * 'cyl', cilindric
1344            * 'lcc', lambert conformal
1345          [res]: resolution:
1346            * 'c', crude
1347            * 'l', low
1348            * 'i', intermediate
1349            * 'h', high
1350            * 'f', full
1351    """
1352    fname = 'draw_topo_geogrid'
1353
1354    if values == 'h':
1355        print fname + '_____________________________________________________________'
1356        print draw_topo_geogrid.__doc__
1357        quit()
1358
1359    expectargs = ['[minTopo]','[maxTopo]', '[lonlatL]', '[title]', '[graphic_kind]', \
1360      '[mapvalues]']
1361 
1362    drw.check_arguments(fname,5,values,':',expectargs)
1363
1364    mintopo = values.split(':')[0].split(',')[0]
1365    maxtopo = values.split(':')[0].split(',')[1]
1366
1367    lonlatLS = values.split(':')[1]
1368    lonlatLv = lonlatLS.split(',')[0]
1369
1370    if lonlatLv == 'None':
1371        lonlatL = None
1372    else:
1373        lonlatL = np.zeros((4), dtype=np.float)
1374        lonlatL[0] = np.float(lonlatLS.split(',')[0])
1375        lonlatL[1] = np.float(lonlatLS.split(',')[1])
1376        lonlatL[2] = np.float(lonlatLS.split(',')[2])
1377        lonlatL[3] = np.float(lonlatLS.split(',')[3])
1378
1379    grtit = values.split(':')[2]
1380    kindfig = values.split(':')[3]
1381    mapvalues = values.split(':')[4]
1382
1383    if not os.path.isfile(ncfile):
1384        print errormsg
1385        print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1386        quit(-1)   
1387
1388    objdomf = NetCDFFile(ncfile, 'r')
1389   
1390    objhgt = objdomf.variables['HGT_M']
1391    objlon = objdomf.variables['XLONG_M']
1392    objlat = objdomf.variables['XLAT_M']
1393
1394    topography = objhgt[0,:,:]
1395
1396    drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL,     \
1397      grtit, kindfig, mapvalues, True)
1398
1399    objdomf.close()
1400
1401    return
1402
1403def draw_topo_geogrid_boxes(ncfiles, values):
1404    """ plotting different geo_em.d[nn].nc topography from WPS files
1405    draw_topo_geogrid_boxes(ncfile, values)
1406      ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference)
1407      values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels]
1408        [min/max]Topo: minimum and maximum values of topography to draw
1409        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1410        title: title of the graph
1411        graphic_kind: kind of figure (jpg, pdf, png)
1412        mapvalues: map characteristics [proj],[res]
1413          see full documentation: http://matplotlib.org/basemap/
1414          [proj]: projection
1415            * 'cyl', cilindric
1416            * 'lcc', lambert conformal
1417          [res]: resolution:
1418            * 'c', crude
1419            * 'l', low
1420            * 'i', intermediate
1421            * 'h', high
1422            * 'f', full
1423        labels= labels to write in the graph
1424    """
1425#    import matplotlib as mpl
1426#    mpl.use('Agg')
1427    import matplotlib.pyplot as plt
1428
1429    fname = 'draw_topo_geogrid_boxes'
1430
1431    if values == 'h':
1432        print fname + '_____________________________________________________________'
1433        print draw_topo_geogrid_boxes.__doc__
1434        quit()
1435
1436    mintopo = values.split(':')[0].split(',')[0]
1437    maxtopo = values.split(':')[0].split(',')[1]
1438
1439    lonlatLS = values.split(':')[1]
1440    lonlatLv = lonlatLS.split(',')[0]
1441
1442    if lonlatLv == 'None':
1443        lonlatL = None
1444    else:
1445        lonlatL = np.zeros((4), dtype=np.float)
1446        lonlatL[0] = np.float(lonlatLS.split(',')[0])
1447        lonlatL[1] = np.float(lonlatLS.split(',')[1])
1448        lonlatL[2] = np.float(lonlatLS.split(',')[2])
1449        lonlatL[3] = np.float(lonlatLS.split(',')[3])
1450
1451    grtit = values.split(':')[2]
1452    kindfig = values.split(':')[3]
1453    mapvalues = values.split(':')[4]
1454    labels = values.split(':')[5]
1455
1456    ncfile = ncfiles.split(',')[0]
1457    if not os.path.isfile(ncfile):
1458        print errormsg
1459        print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1460        quit(-1)   
1461
1462    objdomf = NetCDFFile(ncfile, 'r')
1463   
1464    objhgt = objdomf.variables['HGT_M']
1465    objlon0 = objdomf.variables['XLONG_M']
1466    objlat0 = objdomf.variables['XLAT_M']
1467
1468    topography = objhgt[0,:,:]
1469
1470    Nfiles = len(ncfiles.split(','))
1471    boxlabels = labels.split(',')
1472
1473    Xboxlines = []
1474    Yboxlines = []
1475
1476    for ifile in range(Nfiles):
1477        ncfile = ncfiles.split(',')[ifile]
1478#        print ifile, ncfile
1479        if not os.path.isfile(ncfile):
1480            print errormsg
1481            print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1482            quit(-1)   
1483
1484        objdomfi = NetCDFFile(ncfile, 'r')
1485   
1486        objlon = objdomfi.variables['XLONG_M']
1487        objlat = objdomfi.variables['XLAT_M']
1488
1489        dx = objlon.shape[2]
1490        dy = objlon.shape[1]
1491
1492        Xboxlines.append(objlon[0,0,:])
1493        Yboxlines.append(objlat[0,0,:])
1494        Xboxlines.append(objlon[0,dy-1,:])
1495        Yboxlines.append(objlat[0,dy-1,:])
1496        Xboxlines.append(objlon[0,:,0])
1497        Yboxlines.append(objlat[0,:,0])
1498        Xboxlines.append(objlon[0,:,dx-1])
1499        Yboxlines.append(objlat[0,:,dx-1])
1500
1501        objdomfi.close()
1502
1503    drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels,         \
1504      objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, True)
1505
1506    objdomf.close()
1507
1508    return
1509
1510def movievalslice(origslice, dimmovien, framenum):
1511    """ Function to provide variable slice according to a geneation of a movie
1512    movievals(origslice, dimmovien, framenum)
1513      [origslice]= slice original as [dimname1]|[val1],[...,[dimnameN]|[valN]]
1514        ([val] = -1, full length)
1515      [dimmovien]= name of the dimension to produce the movie
1516      [framenum]= value of the frame to substitue in [origslice] as
1517        [dimmovien]|[framenum]
1518    >>> movievalslice('East_West|-1,North_South|-1,Time|2','Time',0)
1519    East_West|-1,North_South|-1,Time|0
1520    """
1521
1522    fname = 'movievalslice'
1523
1524    if origslice == 'h':
1525        print fname + '_____________________________________________________________'
1526        print movievalslice.__doc__
1527        quit()
1528   
1529    dims = origslice.split(',')
1530
1531    movieslice = ''
1532    idim = 0
1533
1534    for dimn in dims:
1535        dn = dimn.split('|')[0]
1536        if dn == dimmovien:
1537            movieslice = movieslice + dn + '|' + str(framenum)
1538        else:
1539            movieslice = movieslice + dimn
1540        if idim < len(dims)-1: movieslice = movieslice + ','
1541
1542        idim = idim + 1
1543
1544    return movieslice
1545
1546class Capturing(list):
1547    """ Class to capture function output as a list
1548    from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call
1549    """
1550#    from cStringIO import StringIO
1551
1552    def __enter__(self):
1553        self._stdout = sys.stdout
1554        sys.stdout = self._stringio = StringIO()
1555        return self
1556    def __exit__(self, *args):
1557        self.extend(self._stringio.getvalue().splitlines())
1558        sys.stdout = self._stdout
1559
1560def create_movie(netcdfile, values, variable):
1561    """ Function to create a movie assuming ImageMagick installed!
1562      values= [graph]#[movie_dimension]#[graph_values]
1563        [graph]: which graphic
1564        [movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval]
1565          [dimnmovie]; name of the dimension from which make the movie
1566          [dimvmovie]; name of the variable with the values of the dimension
1567          [moviedelay]; delay between frames
1568          [interval]; [beg]@[end]@[freq] or -1 (all)
1569        [graph_values]: values to generate the graphic
1570      netcdfile= netCDF file
1571      variable= variable to use (when applicable)
1572    """ 
1573    fname = 'create_movie'
1574
1575    if values == 'h':
1576        print fname + '_____________________________________________________________'
1577        print create_movie.__doc__
1578        quit()
1579
1580    graph = values.split('#')[0]
1581    movie_dim = values.split('#')[1]
1582    graph_vals = values.split('#')[2]
1583
1584    ncobj = NetCDFFile(netcdfile, 'r')
1585
1586# Movie dimension
1587##
1588    dimnmovie = movie_dim.split('@')[0]
1589    dimvmovie = movie_dim.split('@')[1]
1590    moviedelay = movie_dim.split('@')[2]
1591    moviebeg = int(movie_dim.split('@')[3])
1592
1593    if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie):
1594        print errormsg
1595        print '  ' + fname + ": file '" + netcdfile + "' has not dimension '" +      \
1596          dimnmovie + "' !!!"
1597        quit(-1)
1598
1599    objdmovie = ncobj.dimensions[dimnmovie]
1600    dmovie = len(objdmovie)
1601    if moviebeg != -1:
1602        moviend = int(movie_dim.split('@')[4])
1603        moviefreq = int(movie_dim.split('@')[5])
1604    else:
1605        moviebeg = 0
1606        moviend = dmovie
1607        moviefreq = 1
1608
1609    if dimvmovie == 'WRFTimes':
1610        objvdmovie = ncobj.variables['Times']
1611        vdmovieunits = ''
1612        valsdmovie = []
1613        for it in range(objvdmovie.shape[0]):
1614            valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:],           \
1615              'WRFdatetime', 'Y/m/d H-M-S'))
1616    elif dimvmovie == 'CFtime':
1617        objvdmovie = ncobj.variables['time']
1618        vdmovieunits = ''
1619        print objvdmovie.units
1620        valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard',  \
1621          objvdmovie[:])
1622        valsdmovie = []
1623        for it in range(objvdmovie.shape[0]):
1624            valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:],          \
1625              'matYmdHMS', 'Y/m/d H-M-S'))
1626    else:
1627        if  not drw.searchInlist(ncobj.variables.keys(),dimvmovie):
1628            print errormsg
1629            print '  ' + fname + ": file '" + netcdfile + "' has not variable '" +   \
1630              dimvmovie + "' !!!"
1631            quit(-1)
1632        vdmovieunits = objvdmovie.getncattr('units')
1633        objvdmovie = ncobj.variables[dimvmovie]
1634        if len(objvdmovie.shape) == 1:
1635            vasldmovie = objvdmovie[:]
1636        else:
1637            print errormsg
1638            print '  ' + fname + ': shape', objvdmovie.shape, 'of variable with ' +  \
1639              'dimension movie values not ready!!!'
1640            quit(-1)
1641
1642    ncobj.close()
1643    os.system('rm frame_*.png > /dev/null')
1644
1645# graphic
1646##
1647    if graph == 'draw_2D_shad':
1648        graphvals = graph_vals.split(':')
1649
1650        for iframe in range(moviebeg,moviend,moviefreq):
1651            iframeS = str(iframe).zfill(4)
1652
1653            drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \
1654              5, 'frames')
1655            titgraph = dimnmovie + '|=|' + str(valsdmovie[iframe]) + '|' +           \
1656              vdmovieunits
1657
1658            graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe)
1659            graphvals[6] = titgraph
1660            graphvals[7] = 'png'
1661
1662            graphv = drw.numVector_String(graphvals, ":")
1663
1664            with Capturing() as output:
1665                draw_2D_shad(netcdfile, graphv, variable)
1666
1667            os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png')
1668    else:
1669        print errormsg
1670        print '  ' + fname + ": graphic '" +  graph + "' not defined !!!"
1671        quit(-1)
1672
1673    os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif')
1674
1675    print "Succesfuly creation of movie file 'create_movie.gif' !!!"
1676
1677    return
1678
1679def draw_lines(ncfilens, values, varname):
1680    """ Function to draw different lines at the same time from different files
1681    draw_lines(ncfilens, values, varname):
1682      ncfilens= [filen] ',' separated list of netCDF files
1683      values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]
1684        [dimvname]: ',' list of names of the variable with he values of the common dimension
1685        [valuesaxis]: which axis will be used for the values ('x', or 'y')
1686        [dimtit]: title for the common dimension
1687        [leglabels]: ',' separated list of names for the legend
1688        [vartit]: name of the variable in the graph
1689        [title]: title of the plot ('|' for spaces)
1690        [locleg]: location of the legend (-1, autmoatic)
1691          1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1692          5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1693          9: 'upper center', 10: 'center'
1694        [graphk]: kind of the graphic
1695      varname= variable to plot
1696      values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1697    """
1698
1699    fname = 'draw_lines'
1700
1701    if values == 'h':
1702        print fname + '_____________________________________________________________'
1703        print draw_lines.__doc__
1704        quit()
1705
1706    expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]'
1707    drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
1708
1709    ncfiles = ncfilens.split(',')
1710    dimvnames = values.split(':')[0]
1711    valuesaxis = values.split(':')[1]
1712    dimtit = values.split(':')[2]
1713    leglabels = values.split(':')[3].replace('_','\_')
1714    vartit = values.split(':')[4]
1715    title = values.split(':')[5].replace('|',' ')
1716    locleg = values.split(':')[6]
1717    graphk = values.split(':')[7]
1718
1719    Nfiles = len(ncfiles)
1720
1721# Getting trajectotries
1722##
1723
1724    varvalues = []
1725    dimvalues = []
1726
1727    print '  ' + fname
1728    ifn = 0
1729    for ifile in ncfiles:
1730        filen = ifile.split('@')[0]
1731
1732        print '    filen:',filen
1733
1734        if not os.path.isfile(filen):
1735            print errormsg
1736            print '  ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1737            quit(-1)
1738
1739        objfile = NetCDFFile(filen, 'r')
1740
1741        if dimvnames.find(',') != -1:
1742            dimvname = dimvnames.split(',')
1743        else:
1744            dimvname = [dimvnames]
1745   
1746        found = False
1747        for dvn in dimvname:
1748            if objfile.variables.has_key(dvn):
1749                found = True
1750                break
1751   
1752        if not found:
1753            print errormsg
1754            print '  ' + fname + ": netCDF file '" + filen +                         \
1755              "' does not have variables '" + dimvnames + "' !!"
1756            quit(-1)
1757
1758        if not objfile.variables.has_key(varname):
1759            print errormsg
1760            print '  ' + fname + ": netCDF file '" + filen +                         \
1761              "' does not have variable '" + varname + "' !!"
1762            quit(-1)
1763
1764        vvobj = objfile.variables[varname]
1765        if len(vvobj.shape) != 1:
1766            print errormsg
1767            print '  ' + fname + ': wrong shape:',vvobj.shape," of variable '" +     \
1768              varname +  "' !!"
1769            quit(-1)
1770
1771        for dimvn in dimvname:
1772            if drw.searchInlist(objfile.variables, dimvn):
1773                vdobj = objfile.variables[dimvn]
1774                if len(vdobj.shape) != 1:
1775                    print errormsg
1776                    print '  ' + fname + ': wrong shape:',vdobj.shape,               \
1777                      " of variable '" + dimvn +  "' !!"
1778                    quit(-1)
1779                break
1780
1781        varvalues.append(vvobj[:])
1782        dimvalues.append(vdobj[:])
1783
1784        if ifn == 0:
1785            varunits = vvobj.units
1786
1787        objfile.close()
1788
1789        ifn = ifn + 1
1790
1791    drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','),   \
1792      vartit, varunits, title, locleg, graphk)
1793
1794    return
1795
1796def draw_lines_time(ncfilens, values, varname0):
1797    """ Function to draw different lines at the same time from different files with times
1798    draw_lines_time(ncfilens, values, varname):
1799      ncfilens= [filen] ',' separated list of netCDF files
1800      values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[timevals];[locleg];
1801        [graphk];[collines]
1802        [dimvname]: ',' list of names of the variables with he values of the common dimension
1803        [valuesaxis]: which axis will be used for the values ('x', or 'y')
1804        [dimtit]: title for the common dimension
1805        [leglabels]: ',' separated list of names for the legend
1806        [vartit]: name of the variable in the graph
1807        [title]: title of the plot ('|' for spaces)
1808        [timevals]: [timen]|[units]|[kind]|[tfmt] time labels characteristics
1809           [timen]; name of the time variable
1810           [units]; units string according to CF conventions ([tunits] since
1811             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1812           [kind]; kind of output
1813             'Nval': according to a given number of values as 'Nval',[Nval]
1814             'exct': according to an exact time unit as 'exct',[tunit];
1815               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1816                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1817                'l': milisecond
1818           [tfmt]; desired format
1819        [locleg]: location of the legend (-1, autmoatic)
1820          1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1821          5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1822          9: 'upper center', 10: 'center'
1823        [graphk]: kind of the graphic
1824        [collines]: ',' list of colors for the lines or None for automatic
1825      varname0= ',' list of variable names to plot (assuming only 1 variable per file)
1826      values= 'time;y;time ([DD]${[HH]}$);32x32;$wss^{*}$;wss Taylor's turbulence term;time|hours!since!1949-12-01_00:00:00;exct,12,h|%d$^{%H}$;2;pdf'
1827    """
1828
1829    fname = 'draw_lines_time'
1830
1831    if values == 'h':
1832        print fname + '_____________________________________________________________'
1833        print draw_lines_time.__doc__
1834        quit()
1835
1836    expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1837    expectargs = expectargs + '[timevals];[locleg];[graphk];[collines]'
1838    drw.check_arguments(fname,len(expectargs.split(';')),values,';',expectargs)
1839
1840    ncfiles = ncfilens.split(',')
1841    dimvname0 = values.split(';')[0]
1842    valuesaxis = values.split(';')[1]
1843    dimtit = values.split(';')[2]
1844    leglabels = values.split(';')[3].replace('_','\_')
1845    vartit = values.split(';')[4]
1846    title = values.split(';')[5].replace('|',' ')
1847    timevals = values.split(';')[6]
1848    locleg = int(values.split(';')[7])
1849    graphk = values.split(';')[8]
1850    collines0 = values.split(';')[9]
1851
1852    Nfiles = len(ncfiles)
1853
1854# Multiple variable-dimension names?
1855    if dimvname0.find(',') != -1:
1856        dimvname = dimvname0.split(',')
1857    else:
1858        dimvname = [dimvname0]
1859
1860# Multiple variables?
1861    if varname0.find(',') != -1:
1862        varname = varname0.split(',')
1863    else:
1864        varname = [varname0]
1865
1866# Multiple color names?
1867    if collines0.find(',') != -1:
1868        collines = collines0.split(',')
1869    else:
1870        collines = None
1871
1872# Getting values
1873##
1874    varvalues = []
1875    dimvalues = []
1876    timvalues = []
1877    timvals0 = timvalues
1878
1879    print '  ' + fname
1880    ifn = 0
1881    mintval = 1.e20
1882    maxtval = -1.e20
1883
1884    for ifile in ncfiles:
1885        filen = ifile.split('@')[0]
1886
1887        print '    filen:',filen
1888
1889        if not os.path.isfile(filen):
1890            print errormsg
1891            print '  ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1892            quit(-1)
1893
1894        objfile = NetCDFFile(filen, 'r')
1895
1896        founddvar = False
1897        for dvar in dimvname:
1898            if objfile.variables.has_key(dvar):
1899                founddvar = True
1900                vdobj = objfile.variables[dvar]
1901                if len(vdobj.shape) != 1:
1902                    print errormsg
1903                    print '  ' + fname + ': wrong shape:',vdobj.shape," of " +       \
1904                      "variable '" + dvar +  "' !!"
1905                    quit(-1)
1906                break
1907        if not founddvar:
1908            print errormsg
1909            print '  ' + fname + ": netCDF file '" + filen +                         \
1910            "' has any variable '", dimvname, "' !!"
1911            quit(-1)
1912
1913        foundvar = False
1914        for var in varname:
1915            if objfile.variables.has_key(var):
1916                foundvar = True
1917                vvobj = objfile.variables[var]
1918                if len(vvobj.shape) != 1:
1919                    print errormsg
1920                    print '  ' + fname + ': wrong shape:',vvobj.shape," of " +       \
1921                      "variable '" + var +  "' !!"
1922                    quit(-1)
1923
1924                break
1925        if not foundvar:
1926            print errormsg
1927            print '  ' + fname + ": netCDF file '" + filen +                         \
1928              "' has any variable '", varname, "' !!"
1929            quit(-1)
1930
1931        varvalues.append(vvobj[:])
1932        dimvalues.append(vdobj[:])
1933
1934        mindvals = np.min(vdobj[:])
1935        maxdvals = np.max(vdobj[:])
1936
1937        if mindvals < mintval: mintval = mindvals
1938        if maxdvals > maxtval: maxtval = maxdvals
1939
1940        if ifn == 0:
1941            varunits = drw.units_lunits(vvobj.units)
1942
1943        objfile.close()
1944
1945        ifn = ifn + 1
1946
1947# Times
1948    timename = timevals.split('|')[0]
1949    timeunit = timevals.split('|')[1].replace('!',' ')
1950    timekind = timevals.split('|')[2]
1951    timefmt = timevals.split('|')[3]
1952
1953    dtvals = (maxtval - mintval)/5.
1954    tvals = np.arange(mintval, maxtval, dtvals/2.)
1955
1956    timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
1957
1958    drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','),   \
1959      vartit, varunits, timepos, timelabels, title, locleg, graphk, collines)
1960
1961    return
1962
1963def draw_Neighbourghood_evol(filen, values, variable):
1964    """ Function to draw the temporal evolution of a neighbourghood around a point
1965    draw_Neighbourghood_evol(filen, values, variable)
1966      filen= netCDF file name
1967      values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
1968       [timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
1969        [dimsval]: [dimn1]|[val1]|[dimv1],...,[dimnN]|[valN]|[dimvN] dimension names, values to get
1970          (-1, for all; no name/value pair given full length) and variable with values of the dimension
1971          NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
1972        [neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
1973        [Nneig]: Number of grid points of the full side of the box (odd value)
1974        [Ncol]: Number of columns ('auto': square final plot)
1975        [gvarname]: name of the variable to appear in the graph
1976        [timetits]: [titX],[titY] titles of the axes ('|' for spaces)
1977        [tkinds]: [tkindX]|[tkindY] kinds of time to appear in the graph
1978          'Nval': according to a given number of values as 'Nval',[Nval]
1979          'exct': according to an exact time unit as 'exct',[tunit];
1980            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1981              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1982              'l': milisecond
1983        [timefmts]: [tfmtX],[tfmtY] format of the time labels
1984        [gtitle]: title of the graphic ('|' for spaces)
1985        [shadxtrms]: Extremes for the shading
1986        [cbar]: colorbar to use
1987        [gkind]: kind of graphical output
1988        [ofile]: True/False whether the netcdf with data should be created or not
1989      variable= name of the variable
1990      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'
1991    """ 
1992
1993    fname = 'draw_Neighbourghood_evol'
1994
1995    if values == 'h':
1996        print fname + '_____________________________________________________________'
1997        print draw_Neighbourghood_evol.__doc__
1998        quit()
1999
2000    expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]',      \
2001      '[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]',   \
2002      '[gkind]', '[ofile]']
2003 
2004    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
2005
2006    gvarname = values.split(':')[0]
2007    dimsval = values.split(':')[1].split(',')
2008    neigdims = values.split(':')[2].split(',')
2009    Nneig = int(values.split(':')[3])
2010    Ncol0 = values.split(':')[4]
2011    timetits = values.split(':')[5].split(',')
2012    timekinds = values.split(':')[6].split('|')
2013    timefmts = values.split(':')[7].split(',')
2014    gtitle = values.split(':')[8].replace('|',' ')
2015    shadxtrms = values.split(':')[9].split(',')
2016    cbar = values.split(':')[10]
2017    gkind = values.split(':')[11]
2018    ofile = values.split(':')[12]
2019
2020    if Ncol0 != 'auto': 
2021        Ncol = int(Ncol0)
2022    else:
2023        Ncol = Ncol0
2024
2025    timetits[0] = timetits[0].replace('|',' ')
2026    timetits[1] = timetits[1].replace('|',' ')
2027
2028    if np.mod(Nneig,2) == 0:
2029        print errormsg
2030        print '  ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
2031        quit(-1)
2032
2033    Nneig2 = int(Nneig/2)
2034
2035# Values to slice the variable
2036    dimvslice = {}
2037    dimvvalues = {}
2038    for dimvs in dimsval:
2039        dimn = dimvs.split('|')[0]
2040        dimv = int(dimvs.split('|')[1])
2041        dimnv = dimvs.split('|')[2]
2042
2043        dimvvalues[dimn] = dimnv
2044        dimvslice[dimn] = dimv
2045
2046    ncobj = NetCDFFile(filen, 'r')
2047
2048    varobj = ncobj.variables[variable]
2049
2050    slicevar = []
2051    newdimn = []
2052    newdimsvar = {}
2053
2054    for dimn in varobj.dimensions:
2055        if not drw.searchInlist(dimvslice.keys(), dimn):
2056            dimsize = len(ncobj.dimensions[dimn])
2057            slicevar.append(slice(0, dimsize+1))
2058            newdimn.append(dimn)
2059            newdimsvar[dimn] = dimseize
2060
2061        for dimslicen in dimvslice.keys():
2062            if dimn == dimslicen:
2063                if dimvslice[dimn] != -1:
2064                    if drw.searchInlist(neigdims, dimn):
2065                        slicevar.append(slice(dimvslice[dimn]-Nneig2,                \
2066                          dimvslice[dimn]+Nneig2+1))
2067                        newdimn.append(dimn)
2068                        newdimsvar[dimn] = Nneig
2069                        break
2070                    else:
2071                        slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
2072                        break
2073                else:
2074                    dimsize = len(ncobj.dimensions[dimn])
2075                    slicevar.append(slice(0, dimsize+1))
2076                    newdimn.append(dimn)
2077                    newdimsvar[dimn] = dimsize
2078                    break
2079 
2080    varv = varobj[tuple(slicevar)]
2081
2082    if len(newdimn) != 3:
2083        print errormsg
2084        print '  ' + fname + ': sliced variable with shape=', varv.shape,            \
2085          ' must have three dimensions',len(varv.shape),'given !!'
2086        quit(-1)
2087
2088    newdims = []
2089    for nwdims in newdimn:
2090        newdims.append(newdimsvar[nwdims])
2091
2092# The dimension which is not in the neighbourhood dimensions must be time!
2093    for dim1 in newdimn:
2094        if not drw.searchInlist(neigdims, dim1):
2095            dimt = newdimsvar[dim1]
2096            dimtime = dim1
2097
2098    if Ncol == 'auto':
2099        dimtsqx = int(np.sqrt(dimt)) + 1
2100        dimtsqy = int(np.sqrt(dimt)) + 1
2101    else:
2102        dimtsqx = int(Ncol)
2103        dimtsqy = dimt/dimtsqx + 1
2104
2105    neighbourghood = np.ones((dimtsqy*Nneig,dimtsqx*Nneig), dtype=np.float)*fillValue
2106
2107    for it in range(dimt):
2108        ity = int(it/dimtsqx)
2109        itx = it-ity*dimtsqx
2110
2111        itty = (dimtsqy - ity - 1)*Nneig + Nneig2
2112        ittx = itx*Nneig + Nneig2
2113
2114        neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]=         \
2115          varv[it,::-1,:]
2116
2117    variablevals = drw.variables_values(variable)
2118    if drw.searchInlist(varobj.ncattrs(), 'units'):
2119        vunits = varobj.units
2120    else:
2121        vunits = variablevals[5]
2122
2123# Time values at the X/Y axes
2124    if ncobj.variables[dimvvalues[dimtime]].dtype == '|S1':
2125        print '    ' + fname + ': WRF time variable!'
2126        refdate = '19491201000000'
2127        tunitsval = 'hours'
2128        dimtvalues = np.zeros((dimt), dtype=np.float)
2129        tvals = ncobj.variables[dimvvalues[dimtime]]
2130        yrref=refdate[0:4]
2131        monref=refdate[4:6]
2132        dayref=refdate[6:8]
2133        horref=refdate[8:10]
2134        minref=refdate[10:12]
2135        secref=refdate[12:14]
2136
2137        refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' +        \
2138          minref + ':' + secref
2139        tunits = tunitsval + ' since ' + refdateS
2140        for it in range(dimt):
2141            wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
2142            dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
2143    else:
2144        dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
2145        tunits = ncobj.variables[newdimsvar[dimtime]].units
2146
2147    dimxv = dimtvalues[0:dimtsqx]
2148    dimyv = dimtvalues[0:dimt:dimtsqx]
2149
2150    dimn = ['time','time']
2151
2152    if ofile == 'True':
2153        ofilen = 'Neighbourghood_evol.nc'
2154        newnc = NetCDFFile(ofilen, 'w')
2155# Dimensions
2156        newdim = newnc.createDimension('time',None)
2157        newdim = newnc.createDimension('y',dimtsqy*Nneig)
2158        newdim = newnc.createDimension('x',dimtsqx*Nneig)
2159# Dimension values
2160        newvar = newnc.createVariable('time','f8',('time'))
2161        newvar[:] = np.arange(dimt)
2162        newattr = drw.basicvardef(newvar, 'time','time',tunits)
2163# Neighbourhghood variable
2164        newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'),        \
2165          fill_value=fillValue)
2166        newvar[:] = neighbourghood
2167
2168        newnc.sync()
2169        newnc.close()
2170        print fname + ": Successfull generation of file '" + ofilen + "' !!"
2171
2172# Time ticks
2173    timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
2174    timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
2175
2176    timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
2177    timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
2178
2179    for i in range(2):
2180        if shadxtrms[i][0:1] != 'S':
2181            shadxtrms[i] = np.float(shadxtrms[i])
2182
2183    drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits,   \
2184      timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
2185
2186def draw_points(filen, values):
2187    """ Function to plot a series of points
2188      [values]= [ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:[pointlabels]:
2189        [locleg]:[figureko]:[figuren]
2190        [ptasciifile]:[file],[comchar],[collon],[collat],[lab]
2191          [file]: column ASCII file with the location of the points
2192          [comchar]: '|' list of characters for commentaries
2193          [collon]: number of column with the longitude of the points
2194          [collat]: number of column with the latitude of the points
2195          [collab]: number of column with the labels of the points ('None', and will get
2196            the values from the [pointlabels] variable
2197        [gtit]: title of the figure ('|' for spaces)
2198        [mapvalues]: drawing coastaline ([proj],[res]) or None
2199          [proj]: projection
2200             * 'cyl', cilindric
2201             * 'lcc', lambert conformal
2202          [res]: resolution:
2203             * 'c', crude
2204             * 'l', low
2205             * 'i', intermediate
2206             * 'h', high
2207             * 'f', full
2208        [kindfigure]: kind of figure
2209          'legend': only points in the map with the legend with the names
2210          'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2211        [pointcolor]: color for the points ('auto' for "red")
2212        [pointlabels]: ',' of labels [only used if [collab]='None'] ('None' for no labels)
2213        [locleg]: location of the legend (0, autmoatic)
2214          1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2215          5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2216          9: 'upper center', 10: 'center'
2217        [figureko]: kind of the output file (pdf, png, ...)
2218        [figuren]: name of the figure
2219      [filen]= [ncfile],[lonvarn],[latvarn][,[varn],[dimvals],[vargn],[min],[max],[cbar],[varu]]
2220        [ncfile]: netCDF to use to geolocalize the points
2221        [lonvarn]: name of the variable with the longitudes
2222        [latvarn]: name of the variable with the latitudes
2223        [varn]: optional variable to add staff into the graph
2224        [dimval]: '@' list of [dimn]|[dimval] to get the values for [varn]
2225          [dimn]: name of the dimension
2226          [dimval]: value of the dimension (no value all range)
2227        [vargn]: name of the variable in the graph
2228        [min]: minimum value for the extra variable
2229        [max]: maximum value for the extra variable
2230        [cbar]: color bar
2231        [varu]: units of the variable
2232    """
2233    fname = 'draw_points'
2234
2235    if values == 'h':
2236        print fname + '_____________________________________________________________'
2237        print draw_points.__doc__
2238        quit()
2239
2240    expectargs = '[ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:' +     \
2241      '[pointlabels]:[locleg]:[figurek]:[figuren]'
2242 
2243    drw.check_arguments(fname,len(expectargs.split(':')),values,':',                 \
2244      expectargs.split(':'))
2245
2246    ptasciifile = values.split(':')[0]
2247    gtit = values.split(':')[1]
2248    mapvalues = values.split(':')[2]
2249    kindfigure = values.split(':')[3]
2250    pointcolor = values.split(':')[4]
2251    pointlabels = values.split(':')[5]
2252    locleg = int(values.split(':')[6])
2253    figureko = values.split(':')[7]
2254    figuren = values.split(':')[8]
2255
2256# Getting station locations
2257##
2258    filev = ptasciifile.split(',')[0]
2259    comchar = ptasciifile.split(',')[1].split('|')
2260    collon = int(ptasciifile.split(',')[2])
2261    collat = int(ptasciifile.split(',')[3])
2262    collab = ptasciifile.split(',')[4]
2263
2264    if not os.path.isfile(filev):
2265        print errormsg
2266        print '  ' + fname + ": file '" + filev + "' does not exist!!"
2267        quit(-1)
2268
2269# Getting points position and labels
2270    oascii = open(filev, 'r')
2271    xptval = []
2272    yptval = []
2273    if collab != 'None':
2274        ptlabels = []
2275        for line in oascii:
2276            if not drw.searchInlist(comchar, line[0:1]):
2277                linevals = drw.reduce_spaces(line)
2278                xptval.append(np.float(linevals[collon].replace('\n','')))
2279                yptval.append(np.float(linevals[collat].replace('\n','')))
2280                ptlabels.append(linevals[int(collab)].replace('\n',''))
2281    else:
2282        ptlabels = None
2283        for line in oascii:
2284            if  not drw.searchInlist(comchar, line[0:1]):
2285                linevals = drw.reduce_spaces(line)
2286                xptval.append(np.float(linevals[collon].replace('\n','')))
2287                yptval.append(np.float(linevals[collat].replace('\n','')))
2288
2289    oascii.close()
2290
2291    if pointlabels != 'None' and collab == 'None':
2292        ptlabels = pointlabels.split(',')
2293
2294# Getting localization of the points
2295##
2296    filev = filen.split(',')
2297    Nvals = len(filev)
2298
2299    ncfile = filev[0]
2300    lonvarn = filev[1]
2301    latvarn = filev[2]
2302    varn = None
2303    varextrav = None
2304    if Nvals == 10:
2305        varn = filev[3]
2306        dimvals = filev[4]
2307        varextrav = [filev[5], np.float(filev[6]), np.float(filev[7]), filev[8],     \
2308          filev[9]]
2309   
2310    if not os.path.isfile(ncfile):
2311        print errormsg
2312        print '  ' + fname + ": file '" + ncfile + "' does not exist!!"
2313        quit(-1)
2314
2315    onc = NetCDFFile(ncfile, 'r')
2316   
2317    lonv, latv = drw.lonlat2D(onc.variables[lonvarn], onc.variables[latvarn])
2318
2319    print 'Lluis varn',varn
2320    if varn is not None:
2321        objextra = onc.variables[varn]
2322        vard = objextra.dimensions
2323        dd = {}
2324        for dn in dimvals.split('@'):
2325            ddn = dn.split('|')[0]
2326            ddv = dn.split('|')[1]
2327            dd[ddn] = ddv
2328        slicevar = []
2329        for dv in vard:
2330            found= False
2331            for dn in dd.keys():
2332                if dn == dv:
2333                    slicevar.append(int(dd[dn]))
2334                    found = True
2335                    break
2336            if not found:
2337                slicevar.append(slice(0,len(onc.dimensions[dv])))
2338
2339        varextra = np.squeeze(objextra[tuple(slicevar)])
2340
2341    drw.plot_points(xptval, yptval, lonv, latv, varextra, varextrav, gtit, mapvalues,\
2342      kindfigure, pointcolor, ptlabels, locleg, figureko, figuren)
2343
2344    onc.close()
2345
2346    return
2347
2348def draw_timeSeries(filen, values, variables):
2349    """ Function to draw a time-series
2350    draw_timeSeries(filen, values, variable):
2351      filen= name of the file
2352      values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]
2353      [gvarname]: name of the variable to appear in the graph
2354      [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
2355      [tkind]: kind of time to appear in the graph (assumed x-axis)
2356        'Nval': according to a given number of values as 'Nval',[Nval]
2357        'exct': according to an exact time unit as 'exct',[tunit];
2358          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2359            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2360            'l': milisecond
2361      [timefmt]: format of the time labels
2362      [title]: title of the graphic ('|' for spaces)
2363      [locleg]: location of the legend (-1, autmoatic)
2364        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2365        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2366        9: 'upper center', 10: 'center'
2367      [gkind]: kind of graphical output
2368      variables= [varname],[timename] names of variable and variable with times
2369      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')
2370    """
2371
2372    fname = 'draw_timeSeries'
2373
2374    if values == 'h':
2375        print fname + '_____________________________________________________________'
2376        print draw_timeSeries.__doc__
2377        quit()
2378
2379    expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]',      \
2380      '[locleg]', '[gkind]']
2381 
2382    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
2383
2384    gvarname = values.split(':')[0]
2385    timetit = values.split(':')[1].replace('|',' ')
2386    tkind = values.split(':')[2]
2387    timefmt = values.split(':')[3]
2388    title = values.split(':')[4].replace('|',' ')
2389    locleg = int(values.split(':')[5])
2390    gkind = values.split(':')[6]
2391   
2392    ncobj = NetCDFFile(filen, 'r')
2393
2394    variable = variables.split(',')[0]
2395    timevar = variables.split(',')[1]
2396
2397    if not ncobj.variables.has_key(variable):
2398        print errormsg
2399        print '  ' + fname + ": file '" +  filen + "' does not have variable '" +    \
2400          variable + "' !!"
2401        quit(-1)
2402
2403    if not ncobj.variables.has_key(timevar):
2404        print errormsg
2405        print '  ' + fname + ": file '" +  filen + "' does not have variable time '" \
2406          + timevar + "' !!"
2407        quit(-1)
2408
2409    varobj = ncobj.variables[variable]
2410    timeobj = ncobj.variables[timevar]
2411
2412    dimt = len(timeobj[:])
2413    varvals = np.zeros((2,dimt), dtype=np.float)
2414
2415    gunits = varobj.getncattr('units')
2416    tunits = timeobj.getncattr('units')
2417
2418    varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2419    varvals[1,:] = timeobj[:]
2420
2421    tseriesvals = []
2422    tseriesvals.append(varvals)
2423
2424    drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits,        \
2425      'TimeSeries', gvarname, timetit, tkind, timefmt, title,      \
2426      gvarname.replace('_','\_'), locleg, gkind)
2427
2428    return
2429
2430#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')
2431
2432def draw_trajectories(trjfilens, values, observations):
2433    """ Function to draw different trajectories at the same time
2434    draw_trajectories(trjfilens, values, observations):
2435      trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2436         time intervals and reference maps (first one will be used to plot)
2437        [filen]: name of the file to use (lines with '#', not readed) as:
2438          [t-step] [x] [y]
2439        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2440        [map]: [file]#[lonname]#[latname]
2441          [file]; with the [lon],[lat] matrices
2442          [lonname],[latname]; names of the longitudes and latitudes variables
2443      values=[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]
2444        [leglabels]: ',' separated list of names for the legend
2445        [lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2446        [title]: title of the plot ('!' for spaces)
2447        [graphk]: kind of the graphic
2448        [mapkind]: drawing coastaline ([proj],[res]) or None
2449          [proj]: projection
2450             * 'cyl', cilindric
2451             * 'lcc', lambert conformal
2452          [res]: resolution:
2453             * 'c', crude
2454             * 'l', low
2455             * 'i', intermediate
2456             * 'h', high
2457             * 'f', full
2458      obsevations= [obsfile],[obsname],[Tint],[null]
2459        [obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2460        [obsname]: name of the observations in the graph
2461        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2462        [null]: null value for the observed trajectory
2463    """
2464
2465    fname = 'draw_trajectories'
2466
2467    if values == 'h':
2468        print fname + '_____________________________________________________________'
2469        print draw_trajectories.__doc__
2470        quit()
2471
2472    expectargs = '[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]'
2473 
2474    drw.check_arguments(fname,len(expectargs.split('|')),values,'|',expectargs)
2475
2476    trjfiles = trjfilens.split(',')
2477    leglabels = values.split('|')[0]
2478    lonlatlims = values.split('|')[1]
2479    title = values.split('|')[2].replace('!',' ')
2480    graphk = values.split('|')[3]
2481    mapkind = values.split('|')[4]
2482
2483    Nfiles = len(trjfiles)
2484
2485# Getting trajectotries
2486##
2487
2488    lontrjvalues = []
2489    lattrjvalues = []
2490
2491    print '  ' + fname
2492    ifn = 0
2493    for ifile in trjfiles:
2494        filen = ifile.split('@')[0]
2495        Tint = ifile.split('@')[1]
2496
2497        print '    trajectory:',filen
2498
2499        if Tint != '-1':
2500            Tbeg = Tint
2501            Tend = ifile.split('@')[2]
2502            mapv = ifile.split('@')[3]
2503        else:
2504            mapv = ifile.split('@')[2]
2505
2506        if not os.path.isfile(filen):
2507            print errormsg
2508            print '  ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2509            quit(-1)
2510
2511# Charging longitude and latitude values
2512##
2513        lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2514          mapv.split('#')[2])
2515
2516        if ifn == 0: mapref = mapv
2517        ifn = ifn + 1
2518
2519        objfile = open(filen, 'r')
2520        trjtimev = []
2521        trjxv = []
2522        trjyv = []
2523
2524        for line in objfile:
2525            if line[0:1] != '#':
2526                trjtimev.append(int(line.split(' ')[0]))
2527                trjxv.append(int(line.split(' ')[1]))
2528                trjyv.append(int(line.split(' ')[2]))
2529
2530        objfile.close()
2531
2532        if Tint != '-1':
2533            lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2534            lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2535        else:
2536            lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2537            lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2538
2539# lonlatlimits
2540##
2541
2542    if lonlatlims == 'None':
2543        lonlatlimsv = None
2544    else:
2545        lonlatlimsv = np.zeros((4), dtype=np.float)
2546        lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2547        lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2548        lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2549        lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2550
2551# lon/lat objects
2552##
2553    objnc = NetCDFFile(mapref.split('#')[0])
2554    lonobj = objnc.variables[mapref.split('#')[1]]
2555    latobj = objnc.variables[mapref.split('#')[2]]
2556
2557# map
2558##
2559    if mapkind == 'None':
2560        mapkindv = None
2561    else:
2562        mapkindv = mapkind
2563
2564    if observations is None:
2565        obsname = None
2566    else:
2567        obsfile = observations.split(',')[0]
2568        obsname = observations.split(',')[1]
2569        Tint = observations.split(',')[2]
2570        null = np.float(observations.split(',')[3])
2571        print '    observational trajectory:',obsfile
2572
2573        if not os.path.isfile(obsfile):
2574            print errormsg
2575            print '  ' + fname + ": observations trajectory file '" + obsfile +      \
2576              "' does not exist !!"
2577            quit(-1)
2578
2579        objfile = open(obsfile, 'r')
2580        obstrjtimev = []
2581        obstrjxv = []
2582        obstrjyv = []
2583
2584        for line in objfile:
2585            if line[0:1] != '#':
2586                lon = np.float(line.split(' ')[2])
2587                lat = np.float(line.split(' ')[1])
2588                if not lon == null and not lat == null:
2589                    obstrjtimev.append(int(line.split(' ')[0]))
2590                    obstrjxv.append(lon)
2591                    obstrjyv.append(lat)
2592                else:
2593                    obstrjtimev.append(int(line.split(' ')[0]))
2594                    obstrjxv.append(None)
2595                    obstrjyv.append(None)
2596
2597        objfile.close()
2598
2599        if Tint != '-1':
2600            Tint = int(observations.split(',')[2].split('@')[0])
2601            Tend = int(observations.split(',')[2].split('@')[1])
2602            lontrjvalues.append(obstrjxv[Tint:Tend+1])
2603            lattrjvalues.append(obstrjyv[Tint:Tend+1])
2604        else:
2605            lontrjvalues.append(obstrjxv[:])
2606            lattrjvalues.append(obstrjyv[:])
2607
2608    drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','),          \
2609      lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2610
2611    objnc.close()
2612
2613    return
2614
2615def draw_vals_trajectories(ncfile, values, variable):
2616    """ Function to draw values from the outputs from 'compute_tevolboxtraj'
2617    draw_vals_trajectories(ncfile, values, variable)
2618    ncfile= [ncfile] ',' list of files to use
2619    values= [statisticskind]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2620      [statisticskind]=[statistics][kind]
2621        [statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2622        'mean2', 'stdev'
2623        [kind]: 'box', 'circle' statistics taking the values from a box or a circle
2624        'trj': value following the trajectory
2625      [Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2626      [labels]: ',' separated list of labels for the legend
2627      [locleg]: location of the legend (-1, autmoatic)
2628        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2629        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2630        9: 'upper center', 10: 'center'
2631      [gvarname]: name of the variable to appear in the graph
2632      [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
2633      [tkind]: kind of time to appear in the graph (assumed x-axis)
2634        'Nval': according to a given number of values as 'Nval',[Nval]
2635        'exct': according to an exact time unit as 'exct',[tunit];
2636          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2637            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2638            'l': milisecond
2639      [timefmt]: format of the time labels
2640      [title]: title of the graphic ('|' for spaces)
2641      [gkind]: kind of graphical output
2642    variable= variable to use
2643    """
2644
2645    fname = 'draw_vals_trajectories'
2646
2647    if values == 'h':
2648        print fname + '_____________________________________________________________'
2649        print draw_vals_trajectories.__doc__
2650        quit()
2651
2652    sims = ncfile.split(',')
2653
2654    if len(values.split(':')) != 9:
2655        print errormsg
2656        print '  ' + fname  + ': wrong number of values!', len(values.split(':')),   \
2657          'given 9 needed!!'
2658        print '    ',values.split(':')
2659        quit(-1)
2660
2661    statisticskind = values.split(':')[0]
2662    Tint = values.split(':')[1]
2663    labels = values.split(':')[2]
2664    gvarname = values.split(':')[3]
2665    timetit = values.split(':')[4].replace('|',' ')
2666    tkind = values.split(':')[5]
2667    timefmt = values.split(':')[6]
2668    title = values.split(':')[7].replace('|',' ')
2669    gkind = values.split(':')[8]
2670
2671    leglabels = labels.split('@')[0].split(',')
2672    locleg = int(labels.split('@')[1])
2673
2674    Nsims = len(sims)
2675
2676    if Tint != '-1':
2677        tini = np.float(Tint.split('@')[0])
2678        tend = np.float(Tint.split('@')[1])
2679    else:
2680        tini = -1.
2681        tend = -1.
2682
2683    vartimetrjv = []
2684
2685    print '  ' + fname
2686    for trjfile in sims:
2687        print '    ' + trjfile + ' ...'
2688        if not os.path.isfile(trjfile):
2689            print errormsg
2690            print '  ' + fname + ": trajectory file: '" + trjfile +                  \
2691              "' does not exist !!"
2692            quit(-1)
2693
2694        trjobj = NetCDFFile(trjfile, 'r')
2695        otim = trjobj.variables['time']
2696        if not trjobj.variables.has_key(statisticskind + '_' + variable):
2697            print errormsg
2698            print '  ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2699              statisticskind + '_' + variable + "' !!"
2700            quit(-1)
2701        ovar = trjobj.variables[statisticskind + '_' + variable]
2702        dimt = otim.shape[0]
2703
2704        if trjfile == sims[0]:
2705            gunits = ovar.getncattr('units')
2706            lname = ovar.getncattr('long_name')
2707            tunits = otim.getncattr('units')
2708
2709        if tini != -1:
2710            tiniid = -1
2711            tendid = -1       
2712            for itv in range(dimt):
2713                if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2714                if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2715
2716            if tiniid == -1 or tendid == -1:
2717                print errormsg
2718                print '  ' + main + ' time interval ', tini,',',tend,' not found: ',     \
2719                  tendid, ',', tiniid, ' !!'
2720                print '    data interval [',otim[0], otim[dimt-1],']'
2721                quit(-1)
2722            dimt = tendid - tiniid + 1
2723
2724        else:
2725            dimt = otim.shape[0]
2726
2727        valsv = np.zeros((2,dimt), dtype=np.float)
2728# Checking for time consistency
2729        if otim.getncattr('units') != tunits:
2730            print warnmsg
2731            print '  ' + fname + ': different time units in the plot!!'
2732            newtimes = drw.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
2733        else:
2734            newtimes = otim[:]
2735
2736        if tini == -1:
2737            valsv[1,:] = newtimes[:]
2738            valsv[0,:] = ovar[:]
2739        else:
2740            valsv[1,:] = newtimes[tiniid:tendid+1]
2741            valsv[0,:] = ovar[tiniid:tendid+1]
2742
2743        vartimetrjv.append(valsv)
2744        trjobj.close()
2745
2746    drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits,               \
2747      'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
2748      leglabels, locleg, gkind)
2749
2750def variable_values(values):
2751    """ Function to give back values for a given variable
2752      values= [varname] name of the variable
2753    """
2754
2755    fname = 'variable_values'
2756
2757    values = drw.variables_values(values)
2758
2759    print fname,'values:',values
2760    print fname,'variable_name:',values[0]
2761    print fname,'standard_name:',values[1]
2762    print fname,'min,max:',str(values[2]) + ',' + str(values[3])
2763    print fname,'long_name:',values[4]
2764    print fname,'units:',values[5]
2765    print fname,'shad_colors:',values[6]
2766    print fname,'all_values:',drw.numVector_String(values,',')
2767
2768    return
2769
2770####### ###### ##### #### ### ## #
2771
2772ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
2773
2774### Options
2775##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
2776string_operation="""operation to make:
2777  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]
2778  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]
2779    [ckind]:
2780      'cmap': as it gets from colorbar
2781      'fixc,[colname]': fixed color [colname], all stright lines
2782      'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
2783"""
2784
2785#print string_operation
2786
2787parser = OptionParser()
2788parser.add_option("-f", "--netCDF_file", dest="ncfile", 
2789                  help="file to use", metavar="FILE")
2790parser.add_option("-o", "--operation", type='choice', dest="operation", 
2791       choices=namegraphics, 
2792                  help="operation to make: " + ngraphics, metavar="OPER")
2793parser.add_option("-S", "--valueS", dest="values", 
2794                  help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
2795parser.add_option("-v", "--variable", dest="varname",
2796                  help="[WHEN APPLICABLE] variable to check", metavar="VAR")
2797
2798(opts, args) = parser.parse_args()
2799
2800#######    #######
2801## MAIN
2802    #######
2803
2804# Not checking file operation
2805Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time',                    \
2806  'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time',    \
2807  'draw_points', 'draw_topo_geogrid_boxes', 'draw_trajectories',                     \
2808  'draw_vals_trajectories', 'variable_values']
2809
2810####### ###### ##### #### ### ## #
2811errormsg='ERROR -- error -- ERROR -- error'
2812
2813varn=opts.varname
2814oper=opts.operation
2815
2816if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and                   \
2817  not drw.searchInlist(Notcheckingfile, oper):
2818    print errormsg
2819    print '  ' + main + ': File ' + opts.ncfile + ' does not exist !!'
2820    quit(-1)
2821
2822if oper == 'create_movie':
2823    create_movie(opts.ncfile, opts.values, opts.varname)
2824elif oper == 'draw_2D_shad':
2825    draw_2D_shad(opts.ncfile, opts.values, opts.varname)
2826elif oper == 'draw_2D_shad_time':
2827    draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
2828elif oper == 'draw_2D_shad_cont':
2829    draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
2830elif oper == 'draw_2D_shad_cont_time':
2831    draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
2832elif oper == 'draw_2D_shad_line':
2833    draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
2834elif oper == 'draw_2D_shad_line_time':
2835    draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
2836elif oper == 'draw_barbs':
2837    draw_barbs(opts.ncfile, opts.values, opts.varname)
2838elif oper == 'draw_Neighbourghood_evol':
2839    draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
2840elif oper == 'draw_lines':
2841    draw_lines(opts.ncfile, opts.values, opts.varname)
2842elif oper == 'draw_lines_time':
2843    draw_lines_time(opts.ncfile, opts.values, opts.varname)
2844elif oper == 'draw_points':
2845    draw_points(opts.ncfile, opts.values)
2846elif oper == 'draw_timeSeries':
2847    draw_timeSeries(opts.ncfile, opts.values, opts.varname)
2848elif oper == 'draw_topo_geogrid':
2849    draw_topo_geogrid(opts.ncfile, opts.values)
2850elif oper == 'draw_topo_geogrid_boxes':
2851    draw_topo_geogrid_boxes(opts.ncfile, opts.values)
2852elif oper == 'draw_trajectories':
2853    draw_trajectories(opts.ncfile, opts.values, opts.varname)
2854elif oper == 'draw_vals_trajectories':
2855    draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
2856elif oper == 'list_graphics':
2857# From: http://www.diveintopython.net/power_of_introspection/all_together.html
2858    import drawing as myself
2859    object = myself
2860    for opern in namegraphics:
2861        if  opern != 'list_graphics': 
2862            print opern + '_______ ______ _____ ____ ___ __ _'
2863            print getattr(object, opern).__doc__
2864elif oper == 'variable_values':
2865    variable_values(opts.values)
2866else:
2867    print errormsg
2868    print '  ' + main + ": the graphic '" + oper + "' is not ready !!"
2869    print errormsg
2870    quit()
Note: See TracBrowser for help on using the repository browser.