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

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

Adding 'points' as paramter on 'draw_lines_time'
File size: 112.3 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];[points];[period]
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 ('|' for spaces)
1805        [leglabels]: ',' separated list of names for the legend ('None', no 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, None for automatic, single
1825          value all the same
1826        [points]: ',' list of type of points for the lines, None for automatic, single
1827          value all the same
1828        [period]: which period to plot
1829          '-1': all period
1830          [beg],[end]: beginning and end of the period in reference time-units of first file
1831      varname0= ',' list of variable names to plot (assuming only 1 variable per file)
1832      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'
1833    """
1834
1835    fname = 'draw_lines_time'
1836
1837    if values == 'h':
1838        print fname + '_____________________________________________________________'
1839        print draw_lines_time.__doc__
1840        quit()
1841
1842    expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1843    expectargs = expectargs + '[timevals];[locleg];[graphk];[collines];[points];'
1844    expectargs = expectargs + '[period]'
1845    drw.check_arguments(fname,len(expectargs.split(';')),values,';',expectargs)
1846
1847    ncfiles = ncfilens.split(',')
1848    dimvname0 = values.split(';')[0]
1849    valuesaxis = values.split(';')[1]
1850    dimtit = values.split(';')[2].replace('|',' ')
1851    leglabels = values.split(';')[3].replace('_','\_')
1852    vartit = values.split(';')[4]
1853    title = values.split(';')[5].replace('|',' ')
1854    timevals = values.split(';')[6]
1855    locleg = int(values.split(';')[7])
1856    graphk = values.split(';')[8]
1857    collines0 = values.split(';')[9]
1858    points0 = values.split(';')[10]
1859    period = values.split(';')[11]
1860
1861    Nfiles = len(ncfiles)
1862
1863# Multiple variable-dimension names?
1864    if dimvname0.find(',') != -1:
1865        dimvname = dimvname0.split(',')
1866    else:
1867        dimvname = [dimvname0]
1868
1869# Multiple variables?
1870    if varname0.find(',') != -1:
1871        varname = varname0.split(',')
1872    else:
1873        varname = [varname0]
1874
1875# Multiple color names?
1876    if collines0.find(',') != -1:
1877        collines = collines0.split(',')
1878    elif collines == 'None':
1879        collines = None
1880    else:
1881        collines = []
1882        for ip in range(Nfiles):
1883            collines.append(collines0)
1884
1885# Multiple point types?
1886    if points0.find(',') != -1:
1887        points = points0.split(',')
1888    elif points0 == 'None':
1889        points = None
1890    else:
1891        points = []
1892        for ip in range(Nfiles):
1893            points.append(points0)
1894
1895# Getting values
1896##
1897    varvalues = []
1898    dimvalues = []
1899    timvalues = []
1900    timvals0 = timvalues
1901
1902    print '  ' + fname
1903    ifn = 0
1904    mintval = 1.e20
1905    maxtval = -1.e20
1906
1907    for ifile in ncfiles:
1908        filen = ifile.split('@')[0]
1909
1910        print '    filen:',filen
1911
1912        if not os.path.isfile(filen):
1913            print errormsg
1914            print '  ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1915            quit(-1)
1916
1917        objfile = NetCDFFile(filen, 'r')
1918
1919        founddvar = False
1920        for dvar in dimvname:
1921            if objfile.variables.has_key(dvar):
1922                founddvar = True
1923                vdobj = objfile.variables[dvar]
1924                if len(vdobj.shape) != 1:
1925                    print errormsg
1926                    print '  ' + fname + ': wrong shape:',vdobj.shape," of " +       \
1927                      "variable '" + dvar +  "' !!"
1928                    quit(-1)
1929                break
1930        if not founddvar:
1931            print errormsg
1932            print '  ' + fname + ": netCDF file '" + filen +                         \
1933            "' has any variable '", dimvname, "' !!"
1934            quit(-1)
1935
1936        foundvar = False
1937        for var in varname:
1938            if objfile.variables.has_key(var):
1939                foundvar = True
1940                vvobj = objfile.variables[var]
1941                if len(vvobj.shape) != 1:
1942                    print errormsg
1943                    print '  ' + fname + ': wrong shape:',vvobj.shape," of " +       \
1944                      "variable '" + var +  "' !!"
1945                    quit(-1)
1946
1947                break
1948        if not foundvar:
1949            print errormsg
1950            print '  ' + fname + ": netCDF file '" + filen +                         \
1951              "' has any variable '", varname, "' !!"
1952            quit(-1)
1953
1954
1955# Getting period
1956        dimt = len(vdobj[:])
1957
1958        if period == '-1':
1959            varvalues.append(vvobj[:])
1960            dimvalues.append(vdobj[:])
1961            mindvals = np.min(vdobj[:])
1962            maxdvals = np.max(vdobj[:])
1963        else:
1964            ibeg=-1
1965            iend=-1
1966            tbeg = np.float(period.split(',')[0])
1967            tend = np.float(period.split(',')[1])
1968
1969            for it in range(dimt-1):
1970                if vdobj[it] <= tbeg and vdobj[it+1] > tbeg : ibeg = it
1971                if vdobj[it] <= tend and vdobj[it+1] > tend: iend = it + 1
1972                if ibeg != -1 and iend != -1: break
1973
1974            if ibeg == -1 and iend == -1:
1975                print errormsg
1976                print '  ' + fname + ': Period:',tbeg,',',tend,'not found!!'
1977                print '    ibeg:',ibeg,'iend:',iend
1978                print '    period in file:',np.min(vdobj[:]), np.max(vdobj[:])
1979                quit(-1)
1980            elif iend == -1:
1981                iend = dimt
1982                print warnmsg
1983                print '  ' + fname + ': end of Period:',tbeg,',',tend,'not found!!'
1984                print '    getting last available time instead'
1985                print '    ibeg:',ibeg,'iend:',iend
1986                print '    period in file:',np.min(vdobj[:]), np.max(vdobj[:])
1987            elif ibeg == -1:
1988                ibeg = 0
1989                print warnmsg
1990                print '  ' + fname + ': beginning of Period:',tbeg,',',tend,         \
1991                  'not found!!'
1992                print '    getting first available time instead'
1993                print '    ibeg:',ibeg,'iend:',iend
1994                print '    period in file:',np.min(vdobj[:]), np.max(vdobj[:])
1995
1996            varvalues.append(vvobj[ibeg:iend])
1997            dimvalues.append(vdobj[ibeg:iend])
1998            mindvals = np.min(vdobj[ibeg:iend])
1999            maxdvals = np.max(vdobj[ibeg:iend])
2000
2001            dimt = iend - ibeg
2002
2003        if mindvals < mintval: mintval = mindvals
2004        if maxdvals > maxtval: maxtval = maxdvals
2005
2006        if ifn == 0:
2007            varunits = drw.units_lunits(vvobj.units)
2008
2009        objfile.close()
2010
2011        ifn = ifn + 1
2012
2013# Times
2014    timename = timevals.split('|')[0]
2015    timeunit = timevals.split('|')[1].replace('!',' ')
2016    timekind = timevals.split('|')[2]
2017    timefmt = timevals.split('|')[3]
2018
2019    dtvals = (maxtval - mintval)/dimt
2020    tvals = np.arange(mintval, maxtval, dtvals/2.)
2021
2022    timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
2023
2024    if leglabels.find(',') != -1:
2025      drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit,                  \
2026        leglabels.split(','), vartit, varunits, timepos, timelabels, title, locleg,  \
2027        graphk, collines, points)
2028    else:
2029      drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit,                  \
2030        None, vartit, varunits, timepos, timelabels, title, locleg, graphk, collines)
2031    return
2032
2033def draw_Neighbourghood_evol(filen, values, variable):
2034    """ Function to draw the temporal evolution of a neighbourghood around a point
2035    draw_Neighbourghood_evol(filen, values, variable)
2036      filen= netCDF file name
2037      values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
2038       [timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
2039        [dimsval]: [dimn1]|[val1]|[dimv1],...,[dimnN]|[valN]|[dimvN] dimension names, values to get
2040          (-1, for all; no name/value pair given full length) and variable with values of the dimension
2041          NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
2042        [neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
2043        [Nneig]: Number of grid points of the full side of the box (odd value)
2044        [Ncol]: Number of columns ('auto': square final plot)
2045        [gvarname]: name of the variable to appear in the graph
2046        [timetits]: [titX],[titY] titles of the axes ('|' for spaces)
2047        [tkinds]: [tkindX]|[tkindY] kinds of time to appear in the graph
2048          'Nval': according to a given number of values as 'Nval',[Nval]
2049          'exct': according to an exact time unit as 'exct',[tunit];
2050            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2051              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2052              'l': milisecond
2053        [timefmts]: [tfmtX],[tfmtY] format of the time labels
2054        [gtitle]: title of the graphic ('|' for spaces)
2055        [shadxtrms]: Extremes for the shading
2056        [cbar]: colorbar to use
2057        [gkind]: kind of graphical output
2058        [ofile]: True/False whether the netcdf with data should be created or not
2059      variable= name of the variable
2060      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'
2061    """ 
2062
2063    fname = 'draw_Neighbourghood_evol'
2064
2065    if values == 'h':
2066        print fname + '_____________________________________________________________'
2067        print draw_Neighbourghood_evol.__doc__
2068        quit()
2069
2070    expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]',      \
2071      '[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]',   \
2072      '[gkind]', '[ofile]']
2073 
2074    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
2075
2076    gvarname = values.split(':')[0]
2077    dimsval = values.split(':')[1].split(',')
2078    neigdims = values.split(':')[2].split(',')
2079    Nneig = int(values.split(':')[3])
2080    Ncol0 = values.split(':')[4]
2081    timetits = values.split(':')[5].split(',')
2082    timekinds = values.split(':')[6].split('|')
2083    timefmts = values.split(':')[7].split(',')
2084    gtitle = values.split(':')[8].replace('|',' ')
2085    shadxtrms = values.split(':')[9].split(',')
2086    cbar = values.split(':')[10]
2087    gkind = values.split(':')[11]
2088    ofile = values.split(':')[12]
2089
2090    if Ncol0 != 'auto': 
2091        Ncol = int(Ncol0)
2092    else:
2093        Ncol = Ncol0
2094
2095    timetits[0] = timetits[0].replace('|',' ')
2096    timetits[1] = timetits[1].replace('|',' ')
2097
2098    if np.mod(Nneig,2) == 0:
2099        print errormsg
2100        print '  ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
2101        quit(-1)
2102
2103    Nneig2 = int(Nneig/2)
2104
2105# Values to slice the variable
2106    dimvslice = {}
2107    dimvvalues = {}
2108    for dimvs in dimsval:
2109        dimn = dimvs.split('|')[0]
2110        dimv = int(dimvs.split('|')[1])
2111        dimnv = dimvs.split('|')[2]
2112
2113        dimvvalues[dimn] = dimnv
2114        dimvslice[dimn] = dimv
2115
2116    ncobj = NetCDFFile(filen, 'r')
2117
2118    varobj = ncobj.variables[variable]
2119
2120    slicevar = []
2121    newdimn = []
2122    newdimsvar = {}
2123
2124    for dimn in varobj.dimensions:
2125        if not drw.searchInlist(dimvslice.keys(), dimn):
2126            dimsize = len(ncobj.dimensions[dimn])
2127            slicevar.append(slice(0, dimsize+1))
2128            newdimn.append(dimn)
2129            newdimsvar[dimn] = dimseize
2130
2131        for dimslicen in dimvslice.keys():
2132            if dimn == dimslicen:
2133                if dimvslice[dimn] != -1:
2134                    if drw.searchInlist(neigdims, dimn):
2135                        slicevar.append(slice(dimvslice[dimn]-Nneig2,                \
2136                          dimvslice[dimn]+Nneig2+1))
2137                        newdimn.append(dimn)
2138                        newdimsvar[dimn] = Nneig
2139                        break
2140                    else:
2141                        slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
2142                        break
2143                else:
2144                    dimsize = len(ncobj.dimensions[dimn])
2145                    slicevar.append(slice(0, dimsize+1))
2146                    newdimn.append(dimn)
2147                    newdimsvar[dimn] = dimsize
2148                    break
2149 
2150    varv = varobj[tuple(slicevar)]
2151
2152    if len(newdimn) != 3:
2153        print errormsg
2154        print '  ' + fname + ': sliced variable with shape=', varv.shape,            \
2155          ' must have three dimensions',len(varv.shape),'given !!'
2156        quit(-1)
2157
2158    newdims = []
2159    for nwdims in newdimn:
2160        newdims.append(newdimsvar[nwdims])
2161
2162# The dimension which is not in the neighbourhood dimensions must be time!
2163    for dim1 in newdimn:
2164        if not drw.searchInlist(neigdims, dim1):
2165            dimt = newdimsvar[dim1]
2166            dimtime = dim1
2167
2168    if Ncol == 'auto':
2169        dimtsqx = int(np.sqrt(dimt)) + 1
2170        dimtsqy = int(np.sqrt(dimt)) + 1
2171    else:
2172        dimtsqx = int(Ncol)
2173        dimtsqy = dimt/dimtsqx + 1
2174
2175    neighbourghood = np.ones((dimtsqy*Nneig,dimtsqx*Nneig), dtype=np.float)*fillValue
2176
2177    for it in range(dimt):
2178        ity = int(it/dimtsqx)
2179        itx = it-ity*dimtsqx
2180
2181        itty = (dimtsqy - ity - 1)*Nneig + Nneig2
2182        ittx = itx*Nneig + Nneig2
2183
2184        neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]=         \
2185          varv[it,::-1,:]
2186
2187    variablevals = drw.variables_values(variable)
2188    if drw.searchInlist(varobj.ncattrs(), 'units'):
2189        vunits = varobj.units
2190    else:
2191        vunits = variablevals[5]
2192
2193# Time values at the X/Y axes
2194    if ncobj.variables[dimvvalues[dimtime]].dtype == '|S1':
2195        print '    ' + fname + ': WRF time variable!'
2196        refdate = '19491201000000'
2197        tunitsval = 'hours'
2198        dimtvalues = np.zeros((dimt), dtype=np.float)
2199        tvals = ncobj.variables[dimvvalues[dimtime]]
2200        yrref=refdate[0:4]
2201        monref=refdate[4:6]
2202        dayref=refdate[6:8]
2203        horref=refdate[8:10]
2204        minref=refdate[10:12]
2205        secref=refdate[12:14]
2206
2207        refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' +        \
2208          minref + ':' + secref
2209        tunits = tunitsval + ' since ' + refdateS
2210        for it in range(dimt):
2211            wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
2212            dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
2213    else:
2214        dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
2215        tunits = ncobj.variables[newdimsvar[dimtime]].units
2216
2217    dimxv = dimtvalues[0:dimtsqx]
2218    dimyv = dimtvalues[0:dimt:dimtsqx]
2219
2220    dimn = ['time','time']
2221
2222    if ofile == 'True':
2223        ofilen = 'Neighbourghood_evol.nc'
2224        newnc = NetCDFFile(ofilen, 'w')
2225# Dimensions
2226        newdim = newnc.createDimension('time',None)
2227        newdim = newnc.createDimension('y',dimtsqy*Nneig)
2228        newdim = newnc.createDimension('x',dimtsqx*Nneig)
2229# Dimension values
2230        newvar = newnc.createVariable('time','f8',('time'))
2231        newvar[:] = np.arange(dimt)
2232        newattr = drw.basicvardef(newvar, 'time','time',tunits)
2233# Neighbourhghood variable
2234        newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'),        \
2235          fill_value=fillValue)
2236        newvar[:] = neighbourghood
2237
2238        newnc.sync()
2239        newnc.close()
2240        print fname + ": Successfull generation of file '" + ofilen + "' !!"
2241
2242# Time ticks
2243    timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
2244    timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
2245
2246    timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
2247    timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
2248
2249    for i in range(2):
2250        if shadxtrms[i][0:1] != 'S':
2251            shadxtrms[i] = np.float(shadxtrms[i])
2252
2253    drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits,   \
2254      timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
2255
2256def draw_points(filen, values):
2257    """ Function to plot a series of points
2258      [values]= [ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:[pointlabels]:
2259        [locleg]:[figureko]:[figuren]
2260        [ptasciifile]:[file],[comchar],[collon],[collat],[lab]
2261          [file]: column ASCII file with the location of the points
2262          [comchar]: '|' list of characters for commentaries
2263          [collon]: number of column with the longitude of the points
2264          [collat]: number of column with the latitude of the points
2265          [collab]: number of column with the labels of the points ('None', and will get
2266            the values from the [pointlabels] variable
2267        [gtit]: title of the figure ('|' for spaces)
2268        [mapvalues]: drawing coastaline ([proj],[res]) or None
2269          [proj]: projection
2270             * 'cyl', cilindric
2271             * 'lcc', lambert conformal
2272          [res]: resolution:
2273             * 'c', crude
2274             * 'l', low
2275             * 'i', intermediate
2276             * 'h', high
2277             * 'f', full
2278        [kindfigure]: kind of figure
2279          'legend': only points in the map with the legend with the names
2280          'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2281        [pointcolor]: color for the points ('auto' for "red")
2282        [pointlabels]: ',' of labels [only used if [collab]='None'] ('None' for no labels)
2283        [locleg]: location of the legend (0, autmoatic)
2284          1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2285          5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2286          9: 'upper center', 10: 'center'
2287        [figureko]: kind of the output file (pdf, png, ...)
2288        [figuren]: name of the figure
2289      [filen]= [ncfile],[lonvarn],[latvarn][,[varn],[dimvals],[vargn],[min],[max],[cbar],[varu]]
2290        [ncfile]: netCDF to use to geolocalize the points
2291        [lonvarn]: name of the variable with the longitudes
2292        [latvarn]: name of the variable with the latitudes
2293        [varn]: optional variable to add staff into the graph
2294        [dimval]: '@' list of [dimn]|[dimval] to get the values for [varn]
2295          [dimn]: name of the dimension
2296          [dimval]: value of the dimension (no value all range)
2297        [vargn]: name of the variable in the graph
2298        [min]: minimum value for the extra variable
2299        [max]: maximum value for the extra variable
2300        [cbar]: color bar
2301        [varu]: units of the variable
2302    """
2303    fname = 'draw_points'
2304
2305    if values == 'h':
2306        print fname + '_____________________________________________________________'
2307        print draw_points.__doc__
2308        quit()
2309
2310    expectargs = '[ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:' +     \
2311      '[pointlabels]:[locleg]:[figurek]:[figuren]'
2312 
2313    drw.check_arguments(fname,len(expectargs.split(':')),values,':',                 \
2314      expectargs.split(':'))
2315
2316    ptasciifile = values.split(':')[0]
2317    gtit = values.split(':')[1]
2318    mapvalues = values.split(':')[2]
2319    kindfigure = values.split(':')[3]
2320    pointcolor = values.split(':')[4]
2321    pointlabels = values.split(':')[5]
2322    locleg = int(values.split(':')[6])
2323    figureko = values.split(':')[7]
2324    figuren = values.split(':')[8]
2325
2326# Getting station locations
2327##
2328    filev = ptasciifile.split(',')[0]
2329    comchar = ptasciifile.split(',')[1].split('|')
2330    collon = int(ptasciifile.split(',')[2])
2331    collat = int(ptasciifile.split(',')[3])
2332    collab = ptasciifile.split(',')[4]
2333
2334    if not os.path.isfile(filev):
2335        print errormsg
2336        print '  ' + fname + ": file '" + filev + "' does not exist!!"
2337        quit(-1)
2338
2339# Getting points position and labels
2340    oascii = open(filev, 'r')
2341    xptval = []
2342    yptval = []
2343    if collab != 'None':
2344        ptlabels = []
2345        for line in oascii:
2346            if not drw.searchInlist(comchar, line[0:1]):
2347                linevals = drw.reduce_spaces(line)
2348                xptval.append(np.float(linevals[collon].replace('\n','')))
2349                yptval.append(np.float(linevals[collat].replace('\n','')))
2350                ptlabels.append(linevals[int(collab)].replace('\n',''))
2351    else:
2352        ptlabels = None
2353        for line in oascii:
2354            if  not drw.searchInlist(comchar, line[0:1]):
2355                linevals = drw.reduce_spaces(line)
2356                xptval.append(np.float(linevals[collon].replace('\n','')))
2357                yptval.append(np.float(linevals[collat].replace('\n','')))
2358
2359    oascii.close()
2360
2361    if pointlabels != 'None' and collab == 'None':
2362        ptlabels = pointlabels.split(',')
2363
2364# Getting localization of the points
2365##
2366    filev = filen.split(',')
2367    Nvals = len(filev)
2368
2369    ncfile = filev[0]
2370    lonvarn = filev[1]
2371    latvarn = filev[2]
2372    varn = None
2373    varextrav = None
2374    if Nvals == 10:
2375        varn = filev[3]
2376        dimvals = filev[4]
2377        varextrav = [filev[5], np.float(filev[6]), np.float(filev[7]), filev[8],     \
2378          filev[9]]
2379   
2380    if not os.path.isfile(ncfile):
2381        print errormsg
2382        print '  ' + fname + ": file '" + ncfile + "' does not exist!!"
2383        quit(-1)
2384
2385    onc = NetCDFFile(ncfile, 'r')
2386   
2387    lonv, latv = drw.lonlat2D(onc.variables[lonvarn], onc.variables[latvarn])
2388
2389    print 'Lluis varn',varn
2390    if varn is not None:
2391        objextra = onc.variables[varn]
2392        vard = objextra.dimensions
2393        dd = {}
2394        for dn in dimvals.split('@'):
2395            ddn = dn.split('|')[0]
2396            ddv = dn.split('|')[1]
2397            dd[ddn] = ddv
2398        slicevar = []
2399        for dv in vard:
2400            found= False
2401            for dn in dd.keys():
2402                if dn == dv:
2403                    slicevar.append(int(dd[dn]))
2404                    found = True
2405                    break
2406            if not found:
2407                slicevar.append(slice(0,len(onc.dimensions[dv])))
2408
2409        varextra = np.squeeze(objextra[tuple(slicevar)])
2410
2411    if mapvalues == 'None':
2412        mapV = None
2413    else:
2414        mapV = mapvalues
2415
2416    drw.plot_points(xptval, yptval, lonv, latv, varextra, varextrav, gtit, mapV,     \
2417      kindfigure, pointcolor, ptlabels, locleg, figureko, figuren)
2418
2419    onc.close()
2420
2421    return
2422
2423def draw_timeSeries(filen, values, variables):
2424    """ Function to draw a time-series
2425    draw_timeSeries(filen, values, variable):
2426      filen= name of the file
2427      values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]
2428      [gvarname]: name of the variable to appear in the graph
2429      [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
2430      [tkind]: kind of time to appear in the graph (assumed x-axis)
2431        'Nval': according to a given number of values as 'Nval',[Nval]
2432        'exct': according to an exact time unit as 'exct',[tunit];
2433          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2434            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2435            'l': milisecond
2436      [timefmt]: format of the time labels
2437      [title]: title of the graphic ('|' for spaces)
2438      [locleg]: location of the legend (-1, autmoatic)
2439        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2440        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2441        9: 'upper center', 10: 'center'
2442      [gkind]: kind of graphical output
2443      variables= [varname],[timename] names of variable and variable with times
2444      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')
2445    """
2446
2447    fname = 'draw_timeSeries'
2448
2449    if values == 'h':
2450        print fname + '_____________________________________________________________'
2451        print draw_timeSeries.__doc__
2452        quit()
2453
2454    expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]',      \
2455      '[locleg]', '[gkind]']
2456 
2457    drw.check_arguments(fname,len(expectargs),values,':',expectargs)
2458
2459    gvarname = values.split(':')[0]
2460    timetit = values.split(':')[1].replace('|',' ')
2461    tkind = values.split(':')[2]
2462    timefmt = values.split(':')[3]
2463    title = values.split(':')[4].replace('|',' ')
2464    locleg = int(values.split(':')[5])
2465    gkind = values.split(':')[6]
2466   
2467    ncobj = NetCDFFile(filen, 'r')
2468
2469    variable = variables.split(',')[0]
2470    timevar = variables.split(',')[1]
2471
2472    if not ncobj.variables.has_key(variable):
2473        print errormsg
2474        print '  ' + fname + ": file '" +  filen + "' does not have variable '" +    \
2475          variable + "' !!"
2476        quit(-1)
2477
2478    if not ncobj.variables.has_key(timevar):
2479        print errormsg
2480        print '  ' + fname + ": file '" +  filen + "' does not have variable time '" \
2481          + timevar + "' !!"
2482        quit(-1)
2483
2484    varobj = ncobj.variables[variable]
2485    timeobj = ncobj.variables[timevar]
2486
2487    dimt = len(timeobj[:])
2488    varvals = np.zeros((2,dimt), dtype=np.float)
2489
2490    gunits = varobj.getncattr('units')
2491    tunits = timeobj.getncattr('units')
2492
2493    varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2494    varvals[1,:] = timeobj[:]
2495
2496    tseriesvals = []
2497    tseriesvals.append(varvals)
2498
2499    drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits,        \
2500      'TimeSeries', gvarname, timetit, tkind, timefmt, title,      \
2501      gvarname.replace('_','\_'), locleg, gkind)
2502
2503    return
2504
2505#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')
2506
2507def draw_trajectories(trjfilens, values, observations):
2508    """ Function to draw different trajectories at the same time
2509    draw_trajectories(trjfilens, values, observations):
2510      trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2511         time intervals and reference maps (first one will be used to plot)
2512        [filen]: name of the file to use (lines with '#', not readed) as:
2513          [t-step] [x] [y]
2514        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2515        [map]: [file]#[lonname]#[latname]
2516          [file]; with the [lon],[lat] matrices
2517          [lonname],[latname]; names of the longitudes and latitudes variables
2518      values=[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]
2519        [leglabels]: ',' separated list of names for the legend
2520        [lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2521        [title]: title of the plot ('!' for spaces)
2522        [graphk]: kind of the graphic
2523        [mapkind]: drawing coastaline ([proj],[res]) or None
2524          [proj]: projection
2525             * 'cyl', cilindric
2526             * 'lcc', lambert conformal
2527          [res]: resolution:
2528             * 'c', crude
2529             * 'l', low
2530             * 'i', intermediate
2531             * 'h', high
2532             * 'f', full
2533      obsevations= [obsfile],[obsname],[Tint],[null]
2534        [obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2535        [obsname]: name of the observations in the graph
2536        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2537        [null]: null value for the observed trajectory
2538    """
2539
2540    fname = 'draw_trajectories'
2541
2542    if values == 'h':
2543        print fname + '_____________________________________________________________'
2544        print draw_trajectories.__doc__
2545        quit()
2546
2547    expectargs = '[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]'
2548 
2549    drw.check_arguments(fname,len(expectargs.split('|')),values,'|',expectargs)
2550
2551    trjfiles = trjfilens.split(',')
2552    leglabels = values.split('|')[0]
2553    lonlatlims = values.split('|')[1]
2554    title = values.split('|')[2].replace('!',' ')
2555    graphk = values.split('|')[3]
2556    mapkind = values.split('|')[4]
2557
2558    Nfiles = len(trjfiles)
2559
2560# Getting trajectotries
2561##
2562
2563    lontrjvalues = []
2564    lattrjvalues = []
2565
2566    print '  ' + fname
2567    ifn = 0
2568    for ifile in trjfiles:
2569        filen = ifile.split('@')[0]
2570        Tint = ifile.split('@')[1]
2571
2572        print '    trajectory:',filen
2573
2574        if Tint != '-1':
2575            Tbeg = Tint
2576            Tend = ifile.split('@')[2]
2577            mapv = ifile.split('@')[3]
2578        else:
2579            mapv = ifile.split('@')[2]
2580
2581        if not os.path.isfile(filen):
2582            print errormsg
2583            print '  ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2584            quit(-1)
2585
2586# Charging longitude and latitude values
2587##
2588        lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2589          mapv.split('#')[2])
2590
2591        if ifn == 0: mapref = mapv
2592        ifn = ifn + 1
2593
2594        objfile = open(filen, 'r')
2595        trjtimev = []
2596        trjxv = []
2597        trjyv = []
2598
2599        for line in objfile:
2600            if line[0:1] != '#':
2601                trjtimev.append(int(line.split(' ')[0]))
2602                trjxv.append(int(line.split(' ')[1]))
2603                trjyv.append(int(line.split(' ')[2]))
2604
2605        objfile.close()
2606
2607        if Tint != '-1':
2608            lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2609            lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2610        else:
2611            lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2612            lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2613
2614# lonlatlimits
2615##
2616
2617    if lonlatlims == 'None':
2618        lonlatlimsv = None
2619    else:
2620        lonlatlimsv = np.zeros((4), dtype=np.float)
2621        lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2622        lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2623        lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2624        lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2625
2626# lon/lat objects
2627##
2628    objnc = NetCDFFile(mapref.split('#')[0])
2629    lonobj = objnc.variables[mapref.split('#')[1]]
2630    latobj = objnc.variables[mapref.split('#')[2]]
2631
2632# map
2633##
2634    if mapkind == 'None':
2635        mapkindv = None
2636    else:
2637        mapkindv = mapkind
2638
2639    if observations is None:
2640        obsname = None
2641    else:
2642        obsfile = observations.split(',')[0]
2643        obsname = observations.split(',')[1]
2644        Tint = observations.split(',')[2]
2645        null = np.float(observations.split(',')[3])
2646        print '    observational trajectory:',obsfile
2647
2648        if not os.path.isfile(obsfile):
2649            print errormsg
2650            print '  ' + fname + ": observations trajectory file '" + obsfile +      \
2651              "' does not exist !!"
2652            quit(-1)
2653
2654        objfile = open(obsfile, 'r')
2655        obstrjtimev = []
2656        obstrjxv = []
2657        obstrjyv = []
2658
2659        for line in objfile:
2660            if line[0:1] != '#':
2661                lon = np.float(line.split(' ')[2])
2662                lat = np.float(line.split(' ')[1])
2663                if not lon == null and not lat == null:
2664                    obstrjtimev.append(int(line.split(' ')[0]))
2665                    obstrjxv.append(lon)
2666                    obstrjyv.append(lat)
2667                else:
2668                    obstrjtimev.append(int(line.split(' ')[0]))
2669                    obstrjxv.append(None)
2670                    obstrjyv.append(None)
2671
2672        objfile.close()
2673
2674        if Tint != '-1':
2675            Tint = int(observations.split(',')[2].split('@')[0])
2676            Tend = int(observations.split(',')[2].split('@')[1])
2677            lontrjvalues.append(obstrjxv[Tint:Tend+1])
2678            lattrjvalues.append(obstrjyv[Tint:Tend+1])
2679        else:
2680            lontrjvalues.append(obstrjxv[:])
2681            lattrjvalues.append(obstrjyv[:])
2682
2683    drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','),          \
2684      lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2685
2686    objnc.close()
2687
2688    return
2689
2690def draw_vals_trajectories(ncfile, values, variable):
2691    """ Function to draw values from the outputs from 'compute_tevolboxtraj'
2692    draw_vals_trajectories(ncfile, values, variable)
2693    ncfile= [ncfile] ',' list of files to use
2694    values= [statisticskind]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2695      [statisticskind]=[statistics][kind]
2696        [statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2697        'mean2', 'stdev'
2698        [kind]: 'box', 'circle' statistics taking the values from a box or a circle
2699        'trj': value following the trajectory
2700      [Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2701      [labels]: ',' separated list of labels for the legend
2702      [locleg]: location of the legend (-1, autmoatic)
2703        1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2704        5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2705        9: 'upper center', 10: 'center'
2706      [gvarname]: name of the variable to appear in the graph
2707      [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
2708      [tkind]: kind of time to appear in the graph (assumed x-axis)
2709        'Nval': according to a given number of values as 'Nval',[Nval]
2710        'exct': according to an exact time unit as 'exct',[tunit];
2711          tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2712            'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2713            'l': milisecond
2714      [timefmt]: format of the time labels
2715      [title]: title of the graphic ('|' for spaces)
2716      [gkind]: kind of graphical output
2717    variable= variable to use
2718    """
2719
2720    fname = 'draw_vals_trajectories'
2721
2722    if values == 'h':
2723        print fname + '_____________________________________________________________'
2724        print draw_vals_trajectories.__doc__
2725        quit()
2726
2727    sims = ncfile.split(',')
2728
2729    if len(values.split(':')) != 9:
2730        print errormsg
2731        print '  ' + fname  + ': wrong number of values!', len(values.split(':')),   \
2732          'given 9 needed!!'
2733        print '    ',values.split(':')
2734        quit(-1)
2735
2736    statisticskind = values.split(':')[0]
2737    Tint = values.split(':')[1]
2738    labels = values.split(':')[2]
2739    gvarname = values.split(':')[3]
2740    timetit = values.split(':')[4].replace('|',' ')
2741    tkind = values.split(':')[5]
2742    timefmt = values.split(':')[6]
2743    title = values.split(':')[7].replace('|',' ')
2744    gkind = values.split(':')[8]
2745
2746    leglabels = labels.split('@')[0].split(',')
2747    locleg = int(labels.split('@')[1])
2748
2749    Nsims = len(sims)
2750
2751    if Tint != '-1':
2752        tini = np.float(Tint.split('@')[0])
2753        tend = np.float(Tint.split('@')[1])
2754    else:
2755        tini = -1.
2756        tend = -1.
2757
2758    vartimetrjv = []
2759
2760    print '  ' + fname
2761    for trjfile in sims:
2762        print '    ' + trjfile + ' ...'
2763        if not os.path.isfile(trjfile):
2764            print errormsg
2765            print '  ' + fname + ": trajectory file: '" + trjfile +                  \
2766              "' does not exist !!"
2767            quit(-1)
2768
2769        trjobj = NetCDFFile(trjfile, 'r')
2770        otim = trjobj.variables['time']
2771        if not trjobj.variables.has_key(statisticskind + '_' + variable):
2772            print errormsg
2773            print '  ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2774              statisticskind + '_' + variable + "' !!"
2775            quit(-1)
2776        ovar = trjobj.variables[statisticskind + '_' + variable]
2777        dimt = otim.shape[0]
2778
2779        if trjfile == sims[0]:
2780            gunits = ovar.getncattr('units')
2781            lname = ovar.getncattr('long_name')
2782            tunits = otim.getncattr('units')
2783
2784        if tini != -1:
2785            tiniid = -1
2786            tendid = -1       
2787            for itv in range(dimt):
2788                if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2789                if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2790
2791            if tiniid == -1 or tendid == -1:
2792                print errormsg
2793                print '  ' + main + ' time interval ', tini,',',tend,' not found: ',     \
2794                  tendid, ',', tiniid, ' !!'
2795                print '    data interval [',otim[0], otim[dimt-1],']'
2796                quit(-1)
2797            dimt = tendid - tiniid + 1
2798
2799        else:
2800            dimt = otim.shape[0]
2801
2802        valsv = np.zeros((2,dimt), dtype=np.float)
2803# Checking for time consistency
2804        if otim.getncattr('units') != tunits:
2805            print warnmsg
2806            print '  ' + fname + ': different time units in the plot!!'
2807            newtimes = drw.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
2808        else:
2809            newtimes = otim[:]
2810
2811        if tini == -1:
2812            valsv[1,:] = newtimes[:]
2813            valsv[0,:] = ovar[:]
2814        else:
2815            valsv[1,:] = newtimes[tiniid:tendid+1]
2816            valsv[0,:] = ovar[tiniid:tendid+1]
2817
2818        vartimetrjv.append(valsv)
2819        trjobj.close()
2820
2821    drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits,               \
2822      'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
2823      leglabels, locleg, gkind)
2824
2825def variable_values(values):
2826    """ Function to give back values for a given variable
2827      values= [varname] name of the variable
2828    """
2829
2830    fname = 'variable_values'
2831
2832    values = drw.variables_values(values)
2833
2834    print fname,'values:',values
2835    print fname,'variable_name:',values[0]
2836    print fname,'standard_name:',values[1]
2837    print fname,'min,max:',str(values[2]) + ',' + str(values[3])
2838    print fname,'long_name:',values[4]
2839    print fname,'units:',values[5]
2840    print fname,'shad_colors:',values[6]
2841    print fname,'all_values:',drw.numVector_String(values,',')
2842
2843    return
2844
2845####### ###### ##### #### ### ## #
2846
2847ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
2848
2849### Options
2850##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
2851string_operation="""operation to make:
2852  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]
2853  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]
2854    [ckind]:
2855      'cmap': as it gets from colorbar
2856      'fixc,[colname]': fixed color [colname], all stright lines
2857      'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
2858"""
2859
2860#print string_operation
2861
2862parser = OptionParser()
2863parser.add_option("-f", "--netCDF_file", dest="ncfile", 
2864                  help="file to use", metavar="FILE")
2865parser.add_option("-o", "--operation", type='choice', dest="operation", 
2866       choices=namegraphics, 
2867                  help="operation to make: " + ngraphics, metavar="OPER")
2868parser.add_option("-S", "--valueS", dest="values", 
2869                  help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
2870parser.add_option("-v", "--variable", dest="varname",
2871                  help="[WHEN APPLICABLE] variable to check", metavar="VAR")
2872
2873(opts, args) = parser.parse_args()
2874
2875#######    #######
2876## MAIN
2877    #######
2878
2879# Not checking file operation
2880Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time',                    \
2881  'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time',    \
2882  'draw_points', 'draw_topo_geogrid_boxes', 'draw_trajectories',                     \
2883  'draw_vals_trajectories', 'variable_values']
2884
2885####### ###### ##### #### ### ## #
2886errormsg='ERROR -- error -- ERROR -- error'
2887
2888varn=opts.varname
2889oper=opts.operation
2890
2891if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and                   \
2892  not drw.searchInlist(Notcheckingfile, oper):
2893    print errormsg
2894    print '  ' + main + ': File ' + opts.ncfile + ' does not exist !!'
2895    quit(-1)
2896
2897if oper == 'create_movie':
2898    create_movie(opts.ncfile, opts.values, opts.varname)
2899elif oper == 'draw_2D_shad':
2900    draw_2D_shad(opts.ncfile, opts.values, opts.varname)
2901elif oper == 'draw_2D_shad_time':
2902    draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
2903elif oper == 'draw_2D_shad_cont':
2904    draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
2905elif oper == 'draw_2D_shad_cont_time':
2906    draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
2907elif oper == 'draw_2D_shad_line':
2908    draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
2909elif oper == 'draw_2D_shad_line_time':
2910    draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
2911elif oper == 'draw_barbs':
2912    draw_barbs(opts.ncfile, opts.values, opts.varname)
2913elif oper == 'draw_Neighbourghood_evol':
2914    draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
2915elif oper == 'draw_lines':
2916    draw_lines(opts.ncfile, opts.values, opts.varname)
2917elif oper == 'draw_lines_time':
2918    draw_lines_time(opts.ncfile, opts.values, opts.varname)
2919elif oper == 'draw_points':
2920    draw_points(opts.ncfile, opts.values)
2921elif oper == 'draw_timeSeries':
2922    draw_timeSeries(opts.ncfile, opts.values, opts.varname)
2923elif oper == 'draw_topo_geogrid':
2924    draw_topo_geogrid(opts.ncfile, opts.values)
2925elif oper == 'draw_topo_geogrid_boxes':
2926    draw_topo_geogrid_boxes(opts.ncfile, opts.values)
2927elif oper == 'draw_trajectories':
2928    draw_trajectories(opts.ncfile, opts.values, opts.varname)
2929elif oper == 'draw_vals_trajectories':
2930    draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
2931elif oper == 'list_graphics':
2932# From: http://www.diveintopython.net/power_of_introspection/all_together.html
2933    import drawing as myself
2934    object = myself
2935    for opern in namegraphics:
2936        if  opern != 'list_graphics': 
2937            print opern + '_______ ______ _____ ____ ___ __ _'
2938            print getattr(object, opern).__doc__
2939elif oper == 'variable_values':
2940    variable_values(opts.values)
2941else:
2942    print errormsg
2943    print '  ' + main + ": the graphic '" + oper + "' is not ready !!"
2944    print errormsg
2945    quit()
Note: See TracBrowser for help on using the repository browser.