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

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

Fixing to the final standard format the: shad_time', shad_cont', `shad_cont_time'
File size: 187.6 KB
Line 
1import numpy as np
2import os
3from netCDF4 import Dataset as NetCDFFile
4import drawing_tools as drw
5import generic_tools as gen
6from optparse import OptionParser
7import sys
8from cStringIO import StringIO
9
10## 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
11## e.g. # drawing.py -o draw_2D_shad -f wrfout_d01_2001-11-11_00:00:00 -S 'tas:west_east|-1,south_north|-1,Time|2:XLONG:XLAT:auto:rainbow,auto,auto:Srange,Srange:tas|at|2001-11-11|02|UTC:png:None:cyl,l:True' -v T2
12## e.g. # drawing.py -o draw_2D_shad_time -f wrfout_d01_2001-11-11_00:00:00 -S 'hurs~west_east|-1,south_north|27,Time|-1~XLONG~WRFtime~auto~Blues,auto,auto~Srange,Srange~tas|evolution|at|south_north=27~png~None~WRFtime|hours!since!1949-12-01|exct,1,h|$%d^{%H}$|date!($[DD]^{[HH]}$)~True' -v Q2
13## e.g. # drawing.py -o draw_2D_shad_cont -f wrfout_d01_2001-11-11_00:00:00 -S 'huss,tas:west_east|-1,south_north|-1,Time|2:Time|2:XLONG:XLAT:auto:Blues,auto,auto:fixc,r:%3g:Srange,Srange:260,300,9:huss|&|tas|at|2001-11-11|02|UTC:png:None:cyl,c:True' -v Q2,T2
14## e.g. # drawing.py -o draw_2D_shad_cont_time -f wrfout_d01_2001-11-11_00:00:00 -S 'hfls,tas;west_east|-1,south_north|27,Time|-1;south_north|27;XLONG;WRFtime;auto;BuPu,auto,auto;fixc,y;%3g;Srange,Srange;260,300,9;huss|&|tas|evolution|at|south_north=27;png;None;WRFtime|hours!since!1949-12-01|exct,1,h|$%d^{%H}$|date!($[DD]^{[HH]}$);True' -v LH,T2
15## 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
16## 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'
17## e.g. # drawing.py -f 'geo_em.d02.nc' -o draw_topo_geogrid -S '0.,3000.:None:FF_3dom d02:png:cyl,i'
18## e.g. # drawing.py -o draw_topo_geogrid_boxes -S '0.,3000.:None:FF domains:pdf:lcc,i:d01,d02,d03:0' -f geo_em.d01.nc,geo_em.d02.nc,geo_em.d03.nc
19## 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
20## 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
21## 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'
22## e.g. # drawing.py -o variable_values -S PSFC
23## 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'
24## 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
25## 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
26## 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
27## 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'
28## 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'
29## e.g. # drawing.py -o draw_ptZvals -f geo_v2_2012102123_RR1.nc -S 'pracc:lon,lat:o:80:2,42,7,47,:values!of!values:Blues:cyl,l:pdf' -v pr
30## e.g. # drawing.py -f carteveg5km.nc -o draw_points_lonlat -S 'longitude:latitude:pdf:points!veget|type:green:.:0.5:None:0:legend'
31## e.g. # drawing.py -o draw_vectors -f wrfout_d01_2001-11-11_00:00:00 -S 'T|Time|Times|2,Y|south_north|XLAT|-1,X|west_east|XLONG|-1:3@3,wind@rainbow,9:10m wind,ms-1:cyl,l:WRF 10 m winds:pdf:winds' -v U10,V10
32## e.g. # drawing.py -o draw_basins -f routing.py -S 'Y|y|nav_lat|-1,X|x|nav_lon|-1:1@1,rainbow,9:basins,-:cyl,l:ORCDHIEE river-basins:pdf:basins_named' -v nav_lon,nav_lat,trip,basins
33## e.g. # drawing.py -o draw_river_desc -f diver_desc.nc -S 'Y|lat|lat|-1,X|lon|lon|-1:red,green:Blues:cyl,l:ORCDHIEE rivers:pdf:0:or_rivers -v Amazon
34## e.g. # drawing.py -o draw_vertical_levels -f wrfout_d01_2001-11-11_00:00:00 -S 'true:false:wrfout!vertical!levels!(standard!40):png:4' -v WRFz
35## e.g. # drawing.py -o draw_subbasin -f Caceres_subbasin.nc -S 'Caceres:None:cyl,l:2,True:Caceres:pdf:0:Caceres_subbasin'
36## e.g. # drawing.py -o draw_2lines -f /home/lluis/etudes/WRF_LMDZ/WaquaL_highres/tests/model_graphics/WRF/current/wss_wrfout_tvar_xmean.nc,/home/lluis/etudes/WRF_LMDZ/WaquaL_highres/tests/model_graphics/WRF/current/tas_wrfout_tvar_xmean.nc -v wssvarmean,tasvarmean -S 'lat:0.,20.:0.,4.:-90.,90.:y:wss,tas:red,blue:2.,2.:-:2.,2.:,:wss!tas!mean!meridional!tvar:lon:0:wss_tas_wrfout_tvar_xmean:pdf'
37## e.g. # drawing.py -o draw_2lines_time -f /home/lluis/etudes/WRF_LMDZ/WaquaL_highres/tests/model_graphics/WRF/current/wss_wrfout_xvar_ymean.nc,/home/lluis/etudes/WRF_LMDZ/WaquaL_highres/tests/model_graphics/WRF/current/tas_wrfout_xvar_ymean.nc -v wssvarmean,tasvarmean -S 'time:0.,20.:0.,4.:exct,1,d;%d:x:wss,tas:red,blue:2.,2.:-:2.,2.:,:wss!tas!mean!longitudinal!xvar:time!($[DD]$):0:wss_tas_wrfout_xvar_ymean:pdf'
38
39mainn = 'drawing.py'
40
41errormsg = 'ERROR -- error -- ERROR -- error'
42warnmsg = 'WARNING -- waring -- WARNING -- warning'
43fillValue=1.e20
44
45namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time',                 \
46  'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line',                \
47  'draw_2D_shad_line_time', 'draw_barbs', 'draw_basins',                             \
48  'draw_2lines', 'draw_2lines_time', 'draw_lines', 'draw_lines_time',                \
49  'draw_Neighbourghood_evol',                                                        \
50  'draw_points', 'draw_points_lonlat',                                               \
51  'draw_ptZvals', 'draw_river_desc', 'draw_subbasin', 'draw_timeSeries',             \
52  'draw_topo_geogrid',                                                               \
53  'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories',          \
54  'draw_vectors',  'draw_vertical_levels', 'list_graphics', 'variable_values']
55
56def draw_2D_shad(ncfile, values, varn, axfig=None, fig=None):
57    """ plotting a fields with shading
58    draw_2D_shad(ncfile, values, varn)
59      ncfile= file to use
60      values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[dimxyfmt]:[colorbarvals]:[sminv],[smaxv]:[figt]:
61       [kindfig]:[reverse]:[mapv]:[close]
62        [vnamefs]: Name in the figure of the variable to be shaded
63        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
64          variable a given value is required (-1, all the length)
65        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
66        [dimxyfmt]=[dxs],[dxf],[Ndx],[ordx],[dys],[dyf],[Ndy],[ordx]: format of the values at each axis (or 'auto')
67          [dxs]: style of x-axis ('auto' for 'pretty')
68            'Nfix', values computed at even 'Ndx'
69            'Vfix', values computed at even 'Ndx' increments
70            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
71          [dxf]: format of the labels at the x-axis ('auto' for '%5g')
72          [Ndx]: Number of ticks at the x-axis ('auto' for 5)
73          [ordx]: angle of orientation of ticks at the x-axis ('auto' for horizontal)
74          [dys]: style of y-axis ('auto' for 'pretty')
75          [dyf]: format of the labels at the y-axis ('auto' for '%5g')
76          [Ndy]: Number of ticks at the y-axis ('auto' for 5)
77          [ordy]: angle of orientation of ticks at the y-axis ('auto' for horizontal)
78        [colorbarvals]=[colbarn],[fmtcolorbar],[orientation]
79          [colorbarn]: name of the color bar
80          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
81          [orientation]: orientation of the colorbar ('vertical' (default, by 'auto'), 'horizontal')
82        [smin/axv]: minimum and maximum value for the shading or:
83          'Srange': for full range
84          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
85          'Saroundminmax@val': for min*val,max*val
86          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
87            percentile_(100-val)-median)
88          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
89          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
90          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
91             percentile_(100-val)-median)
92        [figt]: title of the figure ('|' for spaces)
93        [kindfig]: kind of figure
94        [reverse]: Transformation of the values
95          * 'transpose': reverse the axes (x-->y, y-->x)
96          * 'flip'@[x/y]: flip the axis x or y
97        [mapv]: map characteristics: [proj],[res]
98          see full documentation: http://matplotlib.org/basemap/
99          [proj]: projection
100            * 'cyl', cilindric
101            * 'lcc', lambert conformal
102          [res]: resolution:
103            * 'c', crude
104            * 'l', low
105            * 'i', intermediate
106            * 'h', high
107            * 'f', full
108        [close]: Whether figure should be finished or not
109      valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
110      varn= [varsn] name of the variable to plot with shading
111    """
112
113    fname = 'draw_2D_shad'
114    if values == 'h':
115        print fname + '_____________________________________________________________'
116        print draw_2D_shad.__doc__
117        quit()
118
119    expectargs = '[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[dimxyf]:[colbarvals]:' +    \
120      '[sminv],[smaxv]:[figt]:[kindfig]:[reverse]:[mapv]:[close]'
121 
122    drw.check_arguments(fname,values,expectargs,':')
123
124    vnamesfig = values.split(':')[0]
125    dimvals= values.split(':')[1].replace('|',':')
126    vdimxn = values.split(':')[2]
127    vdimyn = values.split(':')[3]
128    dimxyf = values.split(':')[4]
129    colorbarvals = values.split(':')[5]
130    shadminmax = values.split(':')[6]
131    figtitle = values.split(':')[7].replace('|',' ')
132    figkind = values.split(':')[8]
133    revals = values.split(':')[9]
134    mapvalue = values.split(':')[10]
135    close = gen.Str_Bool(values.split(':')[11])
136
137    ncfiles = ncfile
138   
139    if not os.path.isfile(ncfiles):
140        print errormsg
141        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
142        quit(-1)   
143
144    objsf = NetCDFFile(ncfiles, 'r')
145   
146    varns = varn.split(',')[0]
147
148    if  not objsf.variables.has_key(varns):
149        print errormsg
150        print '  ' + fname + ': shading file "' + ncfiles +                          \
151          '" does not have variable "' +  varns + '" !!'
152        quit(-1)
153
154# Variables' values
155    objvars = objsf.variables[varns]
156
157    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
158
159# Dimensions names
160##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
161##    dimnamesv = []
162##    for idd in range(len(objvars.dimensions)):
163##        cutdim = False
164##        for idc in range(len(dimvals.split(','))):
165##            dimcutn = dimvals.split(',')[idc].split(':')[0]
166##            print objvars.dimensions[idd], dimcutn
167##            if objvars.dimensions[idd] == dimcutn:
168##                cutdim = True
169##                break
170##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
171    dimnamesv = [vdimxn, vdimyn]
172
173    if drw.searchInlist(objvars.ncattrs(),'units'):
174        varunits = objvars.getncattr('units')
175    else:
176        print warnmsg
177        print '  ' + fname + ": variable '" + varn + "' without units!!"
178        varunits = '-'
179
180    if  not objsf.variables.has_key(vdimxn):
181        print errormsg
182        print '  ' + fname + ': shading file "' + ncfiles +                          \
183          '" does not have dimension variable "' +  vdimxn + '" !!'
184        quit(-1)
185    if  not objsf.variables.has_key(vdimyn):
186        print errormsg
187        print '  ' + fname + ': shading file "' + ncfiles +                          \
188          '" does not have dimension variable "' +  vdimyn + '" !!'
189        quit(-1)
190
191    objdimx = objsf.variables[vdimxn]
192    objdimy = objsf.variables[vdimyn]
193    if drw.searchInlist(objdimx.ncattrs(),'units'):
194        odimxu = objdimx.getncattr('units')
195    else:
196        print warnmsg
197        print '  ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
198        odimxu = '-'
199
200    if drw.searchInlist(objdimy.ncattrs(),'units'):
201        odimyu = objdimy.getncattr('units')
202    else:
203        print warnmsg
204        print '  ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
205        odimyu = '-'
206
207    odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions,     \
208      objdimy.dimensions, dimvals.replace(':','|').split(','))
209
210    shading_nx = []
211    if shadminmax.split(',')[0][0:1] != 'S':
212            shading_nx.append(np.float(shadminmax.split(',')[0]))
213    else:
214        shading_nx.append(shadminmax.split(',')[0])
215
216    if shadminmax.split(',')[1][0:1] != 'S':
217        shading_nx.append(np.float(shadminmax.split(',')[1]))
218    else:
219        shading_nx.append(shadminmax.split(',')[1])
220
221    if mapvalue == 'None': mapvalue = None
222
223    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvals,',')
224    colormapv = [colbarn, fmtcolbar, colbaror]
225
226    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimxyf,',')
227    xaxis = [xstyl, xaxf, Nxax, xaxor]
228    yaxis = [ystyl, yaxf, Nyax, yaxor]
229
230    if revals == 'None':
231        revals = None
232
233    drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, xaxis,    \
234      yaxis, dimnamesv, colormapv, shading_nx, varunits, figtitle, figkind, revals,  \
235      mapvalue, close)
236
237    return
238
239def draw_2D_shad_time(ncfile, values, varn):
240    """ plotting a fields with shading with time values
241    draw_2D_shad_time(ncfile, values, varn)
242      ncfile= file to use
243      values=[vnamefs]~[dimvals]~[dimxvn]~[dimyvn]~[dimvfmt]~[colorbarvals]~[sminv],[smaxv]~[figt]~
244       [kindfig]~[reverse]~[timevals]~[close]
245        [vnamefs]: Name in the figure of the variable to be shaded
246        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
247          variable a given value is required (-1, all the length, [beg]@[end] for an interval)
248        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)  ('WRFtime' for WRF times)
249        [dimvfmt]=[dvs],[dvf],[Ndv],[ordv]: format of the values for the non-temporal axis (or 'auto')
250          [dvs]: style of non-temporal axis ('auto' for 'pretty')
251            'Nfix', values computed at even 'Ndx'
252            'Vfix', values computed at even 'Ndx' increments
253            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
254          [dvf]: format of the labels at the x-axis ('auto' for '%5g')
255          [Ndv]: Number of ticks at the x-axis ('auto' for 5)
256          [ordv]: angle of orientation of ticks at the x-axis ('auto' for horizontal)
257        [colorbarvals]=[colbarn],[fmtcolorbar],[orientation]
258          [colorbarn]: name of the color bar
259          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
260          [orientation]: orientation of the colorbar ('vertical' (default, by 'auto'), 'horizontal')
261        [smin/axv]: minimum and maximum value for the shading or:
262          'Srange': for full range
263          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
264          'Saroundminmax@val': for min*val,max*val
265          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
266            percentile_(100-val)-median)
267          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
268          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
269          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
270             percentile_(100-val)-median)
271        [figt]: title of the figure ('|' for spaces)
272        [kindfig]: kind of figure
273        [reverse]: Transformation of the values
274          * 'transpose': reverse the axes (x-->y, y-->x)
275          * 'flip'@[x/y]: flip the axis x or y
276        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
277           [timen]; name of the time variable ('WRFtime' for WRF times)
278           [units]; units string according to CF conventions ([tunits] since
279             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
280           [kind]; kind of output
281             'Nval': according to a given number of values as 'Nval',[Nval]
282             'exct': according to an exact time unit as 'exct',[tunit];
283               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
284                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
285                'l': milisecond
286           [tfmt]; desired format
287           [label]; label at the graph ('!' for spaces)
288        [close]: should figure be closed (finished)
289      values='dtcon~Time|-1,bottom_top|-1~presmean~time~auto~seismic,auto,auto~-3.e-6,3.e-6~monthly|'
290        'dtcon~pdf~transpose~time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])~True
291      varn= [varsn] name of the variable to plot with shading
292    """
293    fname = 'draw_2D_shad_time'
294
295    if values == 'h':
296        print fname + '_____________________________________________________________'
297        print draw_2D_shad_time.__doc__
298        quit()
299
300    farguments = '[vnamefs]~[dimvals]~[dimxvn]~[dimyvn]~[dimvfmt]~[colorbarvals]~' + \
301      '[sminv],[smaxv]~[figt]~[kindfig]~[reverse]~[timevals]~[close]'
302    drw.check_arguments(fname,values,farguments,'~')
303
304    vnamesfig = values.split('~')[0]
305    dimvals= values.split('~')[1].replace('|',':')
306    vdimxn = values.split('~')[2]
307    vdimyn = values.split('~')[3]
308    dimvfmt = values.split('~')[4]   
309    colorbarvals = values.split('~')[5]
310    shadminmax = values.split('~')[6]
311    figtitle = values.split('~')[7].replace('|',' ')
312    figkind = values.split('~')[8]
313    revals = values.split('~')[9]
314    timevals = values.split('~')[10]
315    close = gen.Str_Bool(values.split('~')[11])
316
317    ncfiles = ncfile
318   
319    if not os.path.isfile(ncfiles):
320        print errormsg
321        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
322        quit(-1)   
323
324    objsf = NetCDFFile(ncfiles, 'r')
325   
326    varns = varn.split(',')[0]
327
328    if  not objsf.variables.has_key(varns):
329        print errormsg
330        print '  ' + fname + ': shading file "' + ncfiles +                          \
331          '" does not have variable "' +  varns + '" !!'
332        quit(-1)
333
334# Variables' values
335    objvars = objsf.variables[varns]
336
337    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
338
339    dimnamesv = [vdimxn, vdimyn]
340
341    varunits = objvars.getncattr('units')
342
343    if vdimxn != 'WRFtime' and not objsf.variables.has_key(vdimxn):
344        print errormsg
345        print '  ' + fname + ": shading file '" + ncfiles +                          \
346          "' does not have dimension variable '" +  vdimxn + "' !!"
347        quit(-1)
348    if vdimyn != 'WRFtime' and not objsf.variables.has_key(vdimyn):
349        print errormsg
350        print '  ' + fname + ": shading file '" + ncfiles +                          \
351          "' does not have dimension variable '" +  vdimyn + " !!'"
352        quit(-1)
353
354    timename = timevals.split('|')[0]
355    timeunit = timevals.split('|')[1].replace('!',' ')
356    timekind = timevals.split('|')[2]
357    timefmt = timevals.split('|')[3]
358    timelabel = timevals.split('|')[4].replace('!',' ')
359
360    if vdimxn == 'WRFtime' or vdimyn == 'WRFtime':
361        tunitsec = timeunit.split(' ')
362        if len(tunitsec) == 4:
363            refdate = tunitsec[2][0:4]+tunitsec[2][5:7]+tunitsec[2][8:10] +          \
364              tunitsec[3][0:2] + tunitsec[3][3:5] + tunitsec[3][6:8]
365        else:
366            refdate = tunitsec[2][0:4]+tunitsec[2][5:7]+tunitsec[2][8:10]+ '000000'
367        tunitsval = tunitsec[0]
368           
369        timewrfv = objsf.variables['Times']
370        dt = timewrfv.shape[0]
371        cftimes = np.zeros((dt), dtype=np.float)
372
373        for it in range(dt):
374            wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime',      \
375              'matYmdHMS')
376            cftimes[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
377
378    if vdimxn != 'WRFtime':
379        objdimx = objsf.variables[vdimxn]
380        objdimxv = objdimx[:]
381        odimxu = objdimx.getncattr('units')
382    else:
383        objdimxv = cftimes
384        odimxu = timeunit
385
386    if vdimyn != 'WRFtime':
387        objdimy = objsf.variables[vdimyn]
388        objdimyv = objdimy[:]
389        odimyu = objdimy.getncattr('units')
390    else:
391        objdimyv = cftimes
392        odimyu = timeunit
393
394    if len(objdimxv.shape) <= 2:
395        odimxv0 = objdimxv[:]
396    elif len(objdimxv.shape) == 3:
397        odimxv0 = objdimxv[0,:]
398    else:
399        print errormsg
400        print '  ' + fname + ': shape of dimension variable:', objdimxv.shape,       \
401          ' not ready!!'
402        quit(-1)
403    if len(objdimyv.shape) <= 2:
404        odimyv0 = objdimyv[:]
405    elif len(objdimyv.shape) == 3:
406        odimyv0 = objdimyv[0,:]
407    else:
408        print errormsg
409        print '  ' + fname + ': shape of dimension variable:', objdimyv.shape,       \
410          ' not ready!!'
411        quit(-1)
412
413# Dimensional values
414    if vdimxn != 'WRFtime':
415        odxv, dimsdxv = drw.slice_variable(objsf.variables[vdimxn], dimvals.replace(',','|'))
416    else:
417        odxv = cftimes
418        dimsdxv = ['Time']
419
420    if vdimyn != 'WRFtime':
421        odyv, dimsdyv = drw.slice_variable(objsf.variables[vdimyn], dimvals.replace(',','|'))
422    else:
423        odyv = cftimes
424        dimsdyv = ['Time']
425
426    if vdimxn == timename:
427        odimxv = odxv
428        odimxu = timelabel
429        timeaxis = 'x'
430        odimyv = odyv
431        odimyu = objdimy.getncattr('units')
432        timepos, timelabels = drw.CFtimes_plot(odxv, timeunit, timekind, timefmt)
433    elif vdimyn == timename:
434        odimyv = odyv
435        odimyu = timelabel
436        timeaxis = 'y'
437        odimxv = odxv
438        odimxu = objdimx.getncattr('units')
439        timepos, timelabels = drw.CFtimes_plot(odyv, timeunit, timekind, timefmt)
440    else:
441        print errormsg
442        print '  ' + fname + ": time variable '" + timename + "' not found!!"
443        quit(-1)
444
445    shading_nx = []
446    if shadminmax.split(',')[0][0:1] != 'S':
447        shading_nx.append(np.float(shadminmax.split(',')[0]))
448    else:
449        shading_nx.append(shadminmax.split(',')[0])
450
451    if shadminmax.split(',')[1][0:1] != 'S':
452        shading_nx.append(np.float(shadminmax.split(',')[1]))
453    else:
454        shading_nx.append(shadminmax.split(',')[1])
455
456    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvals,',')
457    colormapv = [colbarn, fmtcolbar, colbaror]
458
459    if dimvfmt != 'auto': dimvfmt = dimvfmt + 'auto,auto,auto,auto'
460    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimvfmt,',')
461    vaxis = [xstyl, xaxf, Nxax, xaxor]
462
463    if revals == 'None':
464        revals = None
465
466    drw.plot_2D_shadow_time(valshad, vnamesfig, odxv, odyv, odimxu, odimyu,          \
467      dimnamesv, vaxis, colormapv, shading_nx, varunits, figtitle, figkind, revals,  \
468      timeaxis, timepos, timelabels, close)
469
470    return
471
472def draw_2D_shad_cont(ncfile, values, varn):
473    """ plotting two fields, one with shading and the other with contour lines
474    draw_2D_shad_cont(ncfile, values, varn)
475      ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
476      values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[dimxyfmt]:[colorbarvals]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]:[close]
477        [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
478        [dimvals/c]: list of [dimname]|[value] telling at which dimension of the
479          variable a given value is required (no dimension name, all the length)
480        [dimx/yvn]: names of the variables with the values of the dimensions for the plot
481        [dimxyfmt]=[dxf],[Ndx],[dyf],[Ndy]: format of the values at each axis
482          [dxs]: style of x-axis ('auto' for 'pretty')
483            'Nfix', values computed at even 'Ndx'
484            'Vfix', values computed at even 'Ndx' increments
485            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
486          [dxf]: format of the labels at the x-axis
487          [Ndx]: Number of ticks at the x-axis
488          [ordx]: angle of orientation of ticks at the x-axis ('auto' for horizontal)
489          [dys]: style of y-axis ('auto' for 'pretty')
490          [dyf]: format of the labels at the y-axis
491          [Ndy]: Number of ticks at the y-axis
492          [ordy]: angle of orientation of ticks at the y-axis ('auto' for horizontal)
493        [colorbarvals]=[colbarn],[fmtcolorbar],[orientation]
494          [colorbarn]: name of the color bar
495          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
496          [orientation]: orientation of the colorbar ('vertical' (default, by 'auto'), 'horizontal')
497        [ckind]: kind of contours
498          'cmap': as it gets from colorbar
499          'fixc,[colname]': fixed color [colname], all stright lines
500          'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
501        [clabfmt]: format of the labels in the contour (None, also possible)
502        [smin/axv]: minimum and maximum value for the shading or
503          'Srange': for full range
504          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
505          'Saroundminmax@val': for min*val,max*val
506          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
507            percentile_(100-val)-median)
508          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
509          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
510          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
511             percentile_(100-val)-median)
512        [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
513        [figt]: title of the figure ('|' for spaces)
514        [kindfig]: kind of figure
515        [reverse]: does the values be transposed? 'True/False',
516        [mapv]: map characteristics: [proj],[res]
517          see full documentation: http://matplotlib.org/basemap/
518          [proj]: projection
519            * 'cyl', cilindric
520            * 'lcc', lamvbert conformal
521          [res]: resolution:
522            * 'c', crude
523            * 'l', low
524            * 'i', intermediate
525            * 'h', high
526            * 'f', full
527        [close]: Whether figure should be finished or not
528      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'
529      varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
530    """
531
532    fname = 'draw_2D_shad_cont'
533    if values == 'h':
534        print fname + '_____________________________________________________________'
535        print draw_2D_shad_cont.__doc__
536        quit()
537
538    expectargs = '[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[dimxyf]:' +       \
539      '[colorbarvals]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:' +   \
540      '[figt]:[kindfig]:[reverse]:[mapv]:[close]'
541 
542    drw.check_arguments(fname,values,expectargs,':')
543
544    vnamesfig = values.split(':')[0].split(',')
545    dimvals= values.split(':')[1].replace('|',':')
546    dimvalc= values.split(':')[2].replace('|',':')
547    vdimxn = values.split(':')[3]
548    vdimyn = values.split(':')[4]
549    dimxyf = values.split(':')[5]
550    colorbarvals = values.split(':')[6]
551    countkind = values.split(':')[7]
552    countlabelfmt = values.split(':')[8]
553    shadminmax = values.split(':')[9].split(',')
554    contlevels = values.split(':')[10]
555    figtitle = values.split(':')[11].replace('|',' ')
556    figkind = values.split(':')[12]
557    revals = values.split(':')[13]
558    mapvalue = values.split(':')[14]
559    close = gen.Str_Bool(values.split(':')[15])
560
561    if2filenames = ncfile.find(',')
562
563    if if2filenames != -1:
564        ncfiles = ncfile.split(',')[0]
565        ncfilec = ncfile.split(',')[1]
566    else:
567        ncfiles = ncfile
568        ncfilec = ncfile
569
570    if not os.path.isfile(ncfiles):
571        print errormsg
572        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
573        quit(-1)   
574
575    if not os.path.isfile(ncfilec):
576        print errormsg
577        print '  ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
578        quit(-1)   
579
580    objsf = NetCDFFile(ncfiles, 'r')
581    objcf = NetCDFFile(ncfilec, 'r')
582   
583    varns = varn.split(',')[0]
584    varnc = varn.split(',')[1]
585
586    if  not objsf.variables.has_key(varns):
587        print errormsg
588        print '  ' + fname + ': shading file "' + ncfiles +                          \
589          '" does not have variable "' +  varns + '" !!'
590        quit(-1)
591
592    if  not objcf.variables.has_key(varnc):
593        print errormsg
594        print '  ' + fname + ': contour file "' + ncfilec +                          \
595          '" does not have variable "' +  varnc + '" !!'
596        quit(-1)
597
598# Variables' values
599    objvars = objsf.variables[varns]
600    objvarc = objcf.variables[varnc]
601
602    if len(objvars.shape) != len(objvarc.shape):
603        print errormsg
604        print '  ' + fname + ': shading variable "' + varns + '" has a shape: ',     \
605          objvars.shape,  'different than contour variable "' +  varnc + '": ',      \
606          objvarc.shape,' !!!'
607        quit(-1)
608
609    for idim in range(len(objvars.shape)):
610        if objvars.shape[idim] != objvarc.shape[idim]:
611            print errormsg
612            print '  ' + fname + ': shading variable "' + varns + '" has a shape: ', \
613              objvars.shape,  'different than contour variable "' +  varnc + '": ',  \
614              objvarc.shape,' !!!'
615            quit(-1)
616
617    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
618    valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|'))
619
620# Dimensions names
621##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
622##    dimnamesv = []
623##    for idd in range(len(objvars.dimensions)):
624##        cutdim = False
625##        for idc in range(len(dimvals.split(','))):
626##            dimcutn = dimvals.split(',')[idc].split(':')[0]
627##            print objvars.dimensions[idd], dimcutn
628##            if objvars.dimensions[idd] == dimcutn:
629##                cutdim = True
630##                break
631##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
632    dimnamesv = [vdimxn, vdimyn]
633
634    varunits = []
635    varunits.append(objvars.getncattr('units'))
636    varunits.append(objvarc.getncattr('units'))
637
638    if  not objsf.variables.has_key(vdimxn):
639        print errormsg
640        print '  ' + fname + ': shading file "' + ncfiles +                          \
641          '" does not have dimension variable "' +  vdimxn + '" !!'
642        quit(-1)
643    if  not objsf.variables.has_key(vdimyn):
644        print errormsg
645        print '  ' + fname + ': shading file "' + ncfiles +                          \
646          '" does not have dimension variable "' +  vdimyn + '" !!'
647        quit(-1)
648
649    objdimx = objsf.variables[vdimxn]
650    objdimy = objsf.variables[vdimyn]
651    odimxu = objdimx.getncattr('units')
652    odimyu = objdimy.getncattr('units')
653
654# Getting only that dimensions with coincident names
655    odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions,     \
656      objdimy.dimensions, dimvals.replace(':','|').split(','))
657
658#    dimnvx = objdimx.dimensions
659#    cutslice = []
660#    for idimn in objdimx.dimensions:
661#        found = False
662#        for dimsn in dimsshad:
663#            if idimn == dimsn:
664#                cutslice.append(slice(0,len(objsf.dimensions[idimn])))
665#                found = True
666#        if not found: cutslice.append(0)
667#
668#    odimxv = objdimx[tuple(cutslice)]
669#
670#    dimnvy = objdimy.dimensions
671#    cutslice = []
672#    for idimn in objdimy.dimensions:
673#        found = False
674#        for dimsn in dimsshad:
675#            if idimn == dimsn:
676#                cutslice.append(slice(0,len(objsf.dimensions[idimn])))
677#                found = True
678#        if not found: cutslice.append(0)
679#
680#    odimyv = objdimy[tuple(cutslice)]
681
682#    if len(objdimx.shape) <= 2:
683#        odimxv = objdimx[:]
684#        odimyv = objdimy[:]
685#    elif len(objdimx.shape) == 3:
686#        odimxv = objdimx[0,:]
687#        odimyv = objdimy[0,:]
688#    else:
689#        print errormsg
690#        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
691#          ' not ready!!'
692#        quit(-1)
693
694    if countlabelfmt == 'None': 
695        countlfmt = None
696    else:
697        countlfmt = countlabelfmt
698
699    # Shading limits
700    shading_nx = drw.graphic_range(shadminmax,valshad)
701
702    # Contour limits
703    clevmin = np.float(contlevels.split(',')[0])
704    clevmax = np.float(contlevels.split(',')[1])
705    Nclevels = int(contlevels.split(',')[2])
706
707    levels_cont = gen.pretty_int(clevmin, clevmax, Nclevels)
708
709    if len(levels_cont) <= 1: 
710        print warnmsg
711        print '  ' + fname + ': wrong contour levels:', levels_cont,' !!'
712        del(levels_cont)
713        levels_cont = np.zeros((Nclevels), dtype=np.float)
714        levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
715        print '    generating default ones: ',levels_cont
716
717    if mapvalue == 'None': mapvalue = None
718
719    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvals,',')
720
721    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimxyf,',')
722    xaxis = [xstyl, xaxf, Nxax, xaxor]
723    yaxis = [ystyl, yaxf, Nyax, yaxor]
724
725    if revals == 'None': revals = None
726
727    drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu,  \
728      odimyu, xaxis, yaxis, dimnamesv, [colbarn, fmtcolbar, colbaror], countkind,    \
729      countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind, revals,       \
730      mapvalue, close)
731
732    return
733
734def draw_2D_shad_cont_time(ncfile, values, varn):
735    """ plotting two fields, one with shading and the other with contour lines being
736    one of the dimensions of time characteristics
737    draw_2D_shad_cont(ncfile, values, varn)
738      ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
739      values=[vnamefs];[dimvals];[dimvalc];[dimxvn];[dimyvn];[dimxyf];[colorbarvals];[ckind];[clabfmt];[sminv],[smaxv];[sminc],[smaxv],[Nlev];[figt];[kindfig];[reverse];[timevals];[close]
740        [vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
741        [dimvals/c]: list of [dimname]|[value] telling at which dimension of the
742          variable a given value is required (no dimension name, all the length)
743        [dimxvn]: name of the variables with the values of the dimension of the x-axis ('WRFtime' for WRF times)
744        [dimyvn]: name of the variables with the values of the dimension of the y-axis ('WRFtime' for WRF times)
745        [dimxyf]=[dxf],[Ndx],[dyf],[Ndy]: format of the values at each axis
746          [dxs]: style of x-axis ('auto' for 'pretty')
747            'Nfix', values computed at even 'Ndx' (unique map plotted with constant pixel size)
748            'Vfix', values computed at even 'Ndx' increments
749            'pretty', values computed following aprox. 'Ndx' at 'pretty' intervals (2.,2.5,4,5,10)
750          [dxf]: format of the labels at the x-axis
751          [Ndx]: Number of ticks at the x-axis
752          [ordx]: angle of orientation of ticks at the x-axis ('auto' for horizontal)
753          [dys]: style of y-axis ('auto' for 'pretty')
754          [dyf]: format of the labels at the y-axis
755          [Ndy]: Number of ticks at the y-axis
756          [ordy]: angle of orientation of ticks at the y-axis ('auto' for horizontal)
757        [colorbarvals]=[colbarn],[fmtcolorbar],[orientation]
758          [colorbarn]: name of the color bar
759          [fmtcolorbar]: format of the numbers in the color bar 'C'-like ('auto' for %6g)
760          [orientation]: orientation of the colorbar ('vertical' (default, by 'auto'), 'horizontal')
761        [ckind]: kind of contours
762          'cmap': as it gets from colorbar
763          'fixc,[colname]': fixed color [colname], all stright lines
764          'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
765        [clabfmt]: format of the labels in the contour (None, also possible)
766        [smin/axv]: minimum and maximum value for the shading or:
767          'Srange': for full range
768          'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
769          'Saroundminmax@val': for min*val,max*val
770          'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
771            percentile_(100-val)-median)
772          'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min*@val,max*@val-mean)
773          'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
774          'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
775             percentile_(100-val)-median)
776        [sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
777        [figt]: title of the figure ('|' for spaces)
778        [kindfig]: kind of figure
779        [reverse]: modification to the dimensions:
780          'transposed': transpose matrices
781          'flip',[x/y]: flip only the dimension [x] or [y]
782        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label] time labels characteristics
783          [timen]; name of the time variable ('WRFtime' for WRF times)
784          [units]; units string according to CF conventions ([tunits] since
785            [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
786          [kind]; kind of output
787            'Nval': according to a given number of values as 'Nval',[Nval]
788            'exct': according to an exact time unit as 'exct',[tunit];
789              tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
790               'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
791               'l': milisecond
792          [tfmt]; desired format
793          [label]; label at the graph ('!' for spaces)
794        [close]: Whether figure should be finished or not
795      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])'
796      varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
797    """
798
799    fname = 'draw_2D_shad_cont_time'
800    if values == 'h':
801        print fname + '_____________________________________________________________'
802        print draw_2D_shad_cont_time.__doc__
803        quit()
804
805    expectargs = '[vnamefs];[dimvals];[dimvalc];[dimxvn];[dimyvn];[dimxyf];' +       \
806      '[colorbarvals];[ckind];[clabfmt];[sminv],[smaxv];[sminc],[smaxv],[Nlev];' +   \
807      '[figt];[kindfig];[reverse];[timevals];[close]'
808 
809    drw.check_arguments(fname,values,expectargs,';')
810
811    vnamesfig = values.split(';')[0].split(',')
812    dimvals= values.split(';')[1].replace('|',':')
813    dimvalc= values.split(';')[2].replace('|',':')
814    vdimxn = values.split(';')[3]
815    vdimyn = values.split(';')[4]
816    dimxyf = values.split(';')[5]
817    colorbarvals = values.split(';')[6]
818    countkind = values.split(';')[7]
819    countlabelfmt = values.split(';')[8]
820    shadminmax = values.split(';')[9]
821    contlevels = values.split(';')[10]
822    figtitle = values.split(';')[11].replace('|',' ')
823    figkind = values.split(';')[12]
824    revals = values.split(';')[13]
825    timevals = values.split(';')[14]
826    close = gen.Str_Bool(values.split(';')[15])
827
828    if2filenames = ncfile.find(',')
829
830    if if2filenames != -1:
831        ncfiles = ncfile.split(',')[0]
832        ncfilec = ncfile.split(',')[1]
833    else:
834        ncfiles = ncfile
835        ncfilec = ncfile
836
837    if not os.path.isfile(ncfiles):
838        print errormsg
839        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
840        quit(-1)   
841
842    if not os.path.isfile(ncfilec):
843        print errormsg
844        print '  ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
845        quit(-1)   
846
847    objsf = NetCDFFile(ncfiles, 'r')
848    objcf = NetCDFFile(ncfilec, 'r')
849   
850    varns = varn.split(',')[0]
851    varnc = varn.split(',')[1]
852
853    if  not objsf.variables.has_key(varns):
854        print errormsg
855        print '  ' + fname + ': shading file "' + ncfiles +                          \
856          '" does not have variable "' +  varns + '" !!'
857        quit(-1)
858
859    if  not objcf.variables.has_key(varnc):
860        print errormsg
861        print '  ' + fname + ': contour file "' + ncfilec +                          \
862          '" does not have variable "' +  varnc + '" !!'
863        quit(-1)
864
865# Variables' values
866    objvars = objsf.variables[varns]
867    objvarc = objcf.variables[varnc]
868
869    if len(objvars.shape) != len(objvarc.shape):
870        print errormsg
871        print '  ' + fname + ': shading variable "' + varns + '" has a shape: ',     \
872          objvars.shape,  'different than contour variable "' +  varnc + '": ',      \
873          objvarc.shape,' !!!'
874        quit(-1)
875
876    for idim in range(len(objvars.shape)):
877        if objvars.shape[idim] != objvarc.shape[idim]:
878            print errormsg
879            print '  ' + fname + ': shading variable "' + varns + '" has a shape: ', \
880              objvars.shape,  'different than contour variable "' +  varnc + '": ',  \
881              objvarc.shape,' !!!'
882            quit(-1)
883
884    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
885    valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','|'))
886
887    dimnamesv = [vdimxn, vdimyn]
888
889    varunits = []
890    varunits.append(objvars.getncattr('units'))
891    varunits.append(objvarc.getncattr('units'))
892
893    if vdimxn != 'WRFtime' and not objsf.variables.has_key(vdimxn):
894        print errormsg
895        print '  ' + fname + ': shading file "' + ncfiles +                          \
896          '" does not have dimension variable "' +  vdimxn + '" !!'
897        quit(-1)
898    if vdimyn != 'WRFtime' and not objsf.variables.has_key(vdimyn):
899        print errormsg
900        print '  ' + fname + ': shading file "' + ncfiles +                          \
901          '" does not have dimension variable "' +  vdimyn + '" !!'
902        quit(-1)
903
904    timename = timevals.split('|')[0]
905    timeunit = timevals.split('|')[1].replace('!',' ')
906    timekind = timevals.split('|')[2]
907    timefmt = timevals.split('|')[3]
908    timelabel = timevals.split('|')[4].replace('!',' ')
909
910    if vdimxn == 'WRFtime' or vdimyn == 'WRFtime':
911        tunitsec = timeunit.split(' ')
912        if len(tunitsec) == 4:
913            refdate = tunitsec[2][0:4]+tunitsec[2][5:7]+tunitsec[2][8:10] +          \
914              tunitsec[3][0:2] + tunitsec[3][3:5] + tunitsec[3][6:8]
915        else:
916            refdate = tunitsec[2][0:4]+tunitsec[2][5:7]+tunitsec[2][8:10]+ '000000'
917        tunitsval = tunitsec[0]
918
919        timewrfv = objsf.variables['Times']
920        dt = timewrfv.shape[0]
921        cftimes = np.zeros((dt), dtype=np.float)
922
923        for it in range(dt):
924            wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],'WRFdatetime',      \
925              'matYmdHMS')
926            cftimes[it] = gen.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
927
928    if vdimxn == timename:
929        if vdimxn == 'WRFtime':
930            timevals = cftimes
931            timedims = ['Time']
932        else:
933            timevals = objsf.variables[vdimxn][:]
934            timedims = objsf.variables[vdimxn].dimensions
935        dimt = 'x'
936        ovalaxis = objsf.variables[vdimyn]
937        ovalu = ovalaxis.getncattr('units')
938    elif vdimyn == timename:
939        if vdimyn == 'WRFtime':
940            timevals = cftimes
941            timedims = ['Time']
942        else:
943            timevals = objsf.variables[vdimyn][:]
944            timedims = objsf.variables[vdimyn].dimensions
945        dimt = 'y'
946        ovalaxis = objsf.variables[vdimxn]
947        ovalu = ovalaxis.getncattr('units')
948    else:
949        print errormsg
950        print '  ' + fname + ": time variable '" + timename + "' not found!!"
951        quit(-1)
952
953    timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt)
954
955# Getting only that dimensions with coincident names
956    dimnvx = ovalaxis.dimensions
957    dimsliceaxis = []
958    for dimv in dimvals.split(','):
959        adimvn = dimv.split(':')[0]
960        adimval = dimv.split(':')[1]
961        found = False
962        for dimtn in timedims:
963            if adimvn == dimtn:
964                dimsliceaxis.append(adimvn + ':0')
965                found = True
966                break
967        if not found:
968            dimsliceaxis.append(dimv)
969
970    dimsliceaxisS = '|'.join(dimsliceaxis)
971
972    ovalaxisv, odimaxisv = drw.slice_variable(ovalaxis,dimsliceaxisS)
973    print fname + '; Lluis odimaxisv:', odimaxisv, timedims[0]+':0'
974
975    if countlabelfmt == 'None': 
976        countlfmt = None
977    else:
978        countlfmt = countlabelfmt
979
980    shading_nx = []
981    if shadminmax.split(',')[0][0:1] != 'S':
982            shading_nx.append(np.float(shadminmax.split(',')[0]))
983    else:
984        shading_nx.append(shadminmax.split(',')[0])
985
986    if shadminmax.split(',')[1][0:1] != 'S':
987        shading_nx.append(np.float(shadminmax.split(',')[1]))
988    else:
989        shading_nx.append(shadminmax.split(',')[1])
990
991    clevmin = np.float(contlevels.split(',')[0])
992    clevmax = np.float(contlevels.split(',')[1])
993    Nclevels = int(contlevels.split(',')[2])
994
995    levels_cont = gen.pretty_int(clevmin, clevmax, Nclevels)
996
997    if len(levels_cont) <= 1: 
998        print warnmsg
999        print '  ' + fname + ': wrong contour levels:', levels_cont,' !!'
1000        del(levels_cont)
1001        levels_cont = np.zeros((Nclevels), dtype=np.float)
1002        levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
1003        print '    generating default ones: ',levels_cont
1004
1005    colbarn, fmtcolbar, colbaror = drw.colorbar_vals(colorbarvals,',')
1006
1007    xstyl, xaxf, Nxax, xaxor, ystyl, yaxf, Nyax, yaxor = drw.format_axes(dimxyf,',')
1008    xaxis = [xstyl, xaxf, Nxax, xaxor]
1009    yaxis = [ystyl, yaxf, Nyax, yaxor]
1010
1011    if revals == 'None': revals = None
1012
1013    drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv,         \
1014      timevals, timepos, timelabels, ovalu, timelabel, dimt, xaxis, yaxis,          \
1015      dimnamesv, [colbarn, fmtcolbar, colbaror], countkind, countlfmt, shading_nx,  \
1016      levels_cont, varunits, figtitle, figkind, revals, close)
1017
1018    return
1019
1020def draw_2D_shad_line(ncfile, values, varn):
1021    """ plotting a fields with shading and another with line
1022    draw_2D_shad_line(ncfile, values, varn)
1023      ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
1024      values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
1025       [kindfig]:[reverse]:[mapv]:[close]
1026        [vnamefs]: Name in the figure of the variable to be shaded
1027        [vnamefl]: Name in the figure of the variable to be lined
1028        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
1029          variable a given value is required (-1, all the length)
1030        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
1031        [colorbar]: name of the color bar
1032        [colline]: name of the color for the line
1033        [smin/axv]: minimum and maximum value for the shading
1034        [figt]: title of the figure ('|' for spaces)
1035        [kindfig]: kind of figure
1036        [reverse]: Transformation of the values
1037          * 'transpose': reverse the axes (x-->y, y-->x)
1038          * 'flip'@[x/y]: flip the axis x or y
1039        [mapv]: map characteristics: [proj],[res]
1040          see full documentation: http://matplotlib.org/basemap/
1041          [proj]: projection
1042            * 'cyl', cilindric
1043            * 'lcc', lamvbert conformal
1044          [res]: resolution:
1045            * 'c', crude
1046            * 'l', low
1047            * 'i', intermediate
1048            * 'h', high
1049            * 'f', full
1050      valules= 'rh:z|-1,x|-1:z|-1,x|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
1051      varn= [varsn],[varnl] name of the variable to plot with shading and with line
1052    """
1053
1054    fname = 'draw_2D_shad_line'
1055    if values == 'h':
1056        print fname + '_____________________________________________________________'
1057        print draw_2D_shad_line.__doc__
1058        quit()
1059
1060    farguments = '[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:' +                \
1061      '[colorbar],[colline]:[sminv],[smaxv]:[figt]:[kindfig]:[reverse]:' +           \
1062      '[mapv]:[close]'
1063    drw.check_arguments(fname,values,farguments,':')
1064
1065    vnamesfig = values.split(':')[0].split(',')[0]
1066    dimvals= values.split(':')[1].replace('|',':')
1067    vdimxn = values.split(':')[2]
1068    vdimyn = values.split(':')[3]
1069    colbarn = values.split(':')[4].split(',')[0]
1070    shadminmax = values.split(':')[5]
1071    figtitle = values.split(':')[6].replace('|',' ')
1072    figkind = values.split(':')[7]
1073    revals = values.split(':')[8]
1074    mapvalue = values.split(':')[9]
1075#    varn = values.split(':')[10]
1076
1077    ncfiles = ncfile.split(',')[0]
1078   
1079    if not os.path.isfile(ncfiles):
1080        print errormsg
1081        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1082        quit(-1)   
1083
1084    objsf = NetCDFFile(ncfiles, 'r')
1085   
1086    varns = varn.split(',')[0]
1087
1088    if  not objsf.variables.has_key(varns):
1089        print errormsg
1090        print '  ' + fname + ': shading file "' + ncfiles +                          \
1091          '" does not have variable "' +  varns + '" !!'
1092        quit(-1)
1093
1094# Variables' values
1095    objvars = objsf.variables[varns]
1096
1097    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
1098
1099# Dimensions names
1100##    print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
1101##    dimnamesv = []
1102##    for idd in range(len(objvars.dimensions)):
1103##        cutdim = False
1104##        for idc in range(len(dimvals.split(','))):
1105##            dimcutn = dimvals.split(',')[idc].split(':')[0]
1106##            print objvars.dimensions[idd], dimcutn
1107##            if objvars.dimensions[idd] == dimcutn:
1108##                cutdim = True
1109##                break
1110##        if not cutdim: dimnamesv.append(objvars.dimensions[idd])
1111    dimnamesv = [vdimxn, vdimyn]
1112
1113    varunits = objvars.getncattr('units')
1114
1115    if  not objsf.variables.has_key(vdimxn):
1116        print errormsg
1117        print '  ' + fname + ': shading file "' + ncfiles +                          \
1118          '" does not have dimension variable "' +  vdimxn + '" !!'
1119        quit(-1)
1120    if  not objsf.variables.has_key(vdimyn):
1121        print errormsg
1122        print '  ' + fname + ': shading file "' + ncfiles +                          \
1123          '" does not have dimension variable "' +  vdimyn + '" !!'
1124        quit(-1)
1125
1126    objdimx = objsf.variables[vdimxn]
1127    objdimy = objsf.variables[vdimyn]
1128    odimxu = objdimx.getncattr('units')
1129    odimyu = objdimy.getncattr('units')
1130
1131    if len(objdimx.shape) <= 2:
1132#        odimxv = objdimx[valshad.shape]
1133#        odimyv = objdimy[valshad.shape]
1134        odimxv = objdimx[:]
1135        odimyv = objdimy[:]
1136
1137    elif len(objdimx.shape) == 3:
1138#        dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
1139#        odimxv = objdimx[tuple(dimcut)]
1140#        odimyv = objdimy[tuple(dimcut)]
1141        odimxv = objdimx[0,:]
1142        odimyv = objdimy[0,:]
1143    else:
1144        print errormsg
1145        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
1146          ' not ready!!'
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    if mapvalue == 'None': mapvalue = None
1154
1155# line plot
1156##
1157    ncfilel = ncfile.split(',')[1]
1158    vnamelfig = values.split(':')[0].split(',')[1]
1159    varnl = varn.split(',')[1]
1160    colline = values.split(':')[4].split(',')[1]
1161
1162    objlf = NetCDFFile(ncfilel,'r')
1163    objlvar = objlf.variables[varnl]
1164
1165    linevals = objlvar[:]
1166    varlunits = objlvar.units
1167
1168    drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
1169      odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits,  \
1170      figtitle, figkind, revals, mapvalue, True)
1171
1172    objsf.close()
1173    objlf.close()
1174
1175    return
1176
1177def draw_2D_shad_line_time(ncfile, values, varn):
1178    """ plotting a fields with shading and a line with time values
1179    draw_2D_shad_line(ncfile, values, varn)
1180      ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
1181      values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
1182       [kindfig]:[reverse]:[timevals]:[close]
1183        [vnamefs]: Name in the figure of the variable to be shaded
1184        [vnamefl]: Name in the figure of the variable to be lined
1185        [dimvals]: ',' list of [dimname]|[value] telling at which dimension of the
1186          variable a given value is required (-1, all the length)
1187        [dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
1188        [colorbar]: name of the color bar
1189        [smin/axv]: minimum and maximum value for the shading
1190        [figt]: title of the figure ('|' for spaces)
1191        [kindfig]: kind of figure
1192        [reverse]: Transformation of the values
1193          * 'transpose': reverse the axes (x-->y, y-->x)
1194          * 'flip'@[x/y]: flip the axis x or y
1195        [timevals]: [timen]|[units]|[kind]|[tfmt]|[label]|[timeaxis] time labels characteristics
1196           [timen]; name of the time variable
1197           [units]; units string according to CF conventions ([tunits] since
1198             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1199           [kind]; kind of output
1200             'Nval': according to a given number of values as 'Nval',[Nval]
1201             'exct': according to an exact time unit as 'exct',[tunit];
1202               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1203                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1204                'l': milisecond
1205           [tfmt]; desired format
1206           [label]; label at the graph ('!' for spaces)
1207        [close]: should figure be closed (finished)
1208      values='dtcon,prc:Time|-1,bottom_top|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly|'
1209        'dtcon:pdf:transpose:time|hours!since!1949-12-01|exct,2,d|%d|date!([DD])|x:True
1210      varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
1211    """
1212    fname = 'draw_2D_shad_line_time'
1213    if values == 'h':
1214        print fname + '_____________________________________________________________'
1215        print draw_2D_shad__line_time.__doc__
1216        quit()
1217
1218    farguments = '[vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:' +               \
1219      '[colorbar]:[sminv],[smaxv]:[figt]:[kindfig]:[reverse]:[timevals]:[close]'
1220    drw.check_arguments(fname,values,farguments,':')
1221
1222    vnamesfig = values.split(':')[0].split(',')[0]
1223    dimvals= values.split(':')[1].replace('|',':')
1224    vdimxn = values.split(':')[2]
1225    vdimyn = values.split(':')[3]
1226    colbarn = values.split(':')[4]
1227    shadminmax = values.split(':')[5]
1228    figtitle = values.split(':')[6].replace('|',' ')
1229    figkind = values.split(':')[7]
1230    revals = values.split(':')[8]
1231    timevals = values.split(':')[9]
1232    close = values.split(':')[10]
1233
1234    ncfiles = ncfile.split(',')[0]
1235   
1236    if not os.path.isfile(ncfiles):
1237        print errormsg
1238        print '  ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1239        quit(-1)   
1240
1241    objsf = NetCDFFile(ncfiles, 'r')
1242   
1243    varns = varn.split(',')[0]
1244
1245    if  not objsf.variables.has_key(varns):
1246        print errormsg
1247        print '  ' + fname + ': shading file "' + ncfiles +                          \
1248          '" does not have variable "' +  varns + '" !!'
1249        quit(-1)
1250
1251# Variables' values
1252    objvars = objsf.variables[varns]
1253
1254    valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','|'))
1255
1256    dimnamesv = [vdimxn, vdimyn]
1257
1258    varunits = objvars.getncattr('units')
1259
1260    if  not objsf.variables.has_key(vdimxn):
1261        print errormsg
1262        print '  ' + fname + ': shading file "' + ncfiles +                          \
1263          '" does not have dimension variable "' +  vdimxn + '" !!'
1264        quit(-1)
1265    if  not objsf.variables.has_key(vdimyn):
1266        print errormsg
1267        print '  ' + fname + ': shading file "' + ncfiles +                          \
1268          '" does not have dimension variable "' +  vdimyn + '" !!'
1269        quit(-1)
1270
1271    objdimx = objsf.variables[vdimxn]
1272    objdimy = objsf.variables[vdimyn]
1273    odimxu = objdimx.getncattr('units')
1274    odimyu = objdimy.getncattr('units')
1275
1276    if len(objdimx.shape) <= 2:
1277        odimxv = objdimx[:]
1278        odimyv = objdimy[:]
1279
1280    elif len(objdimx.shape) == 3:
1281        odimxv = objdimx[0,:]
1282        odimyv = objdimy[0,:]
1283    else:
1284        print errormsg
1285        print '  ' + fname + ': shape of dimension variable:', objdimx.shape,        \
1286          ' not ready!!'
1287        quit(-1)
1288
1289    timename = timevals.split('|')[0]
1290    timeunit = timevals.split('|')[1].replace('!',' ')
1291    timekind = timevals.split('|')[2]
1292    timefmt = timevals.split('|')[3]
1293    timelabel = timevals.split('|')[4].replace('!',' ')
1294
1295    if vdimxn == timename:
1296        odimxv = objsf.variables[vdimxn][:]
1297        odimxu = timelabel
1298        timeaxis = 'x'
1299        odimyv = objsf.variables[vdimyn]
1300        odimyu = odimyv.getncattr('units')
1301        timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1302    elif vdimyn == timename:
1303        odimyv = objsf.variables[vdimyn][:]
1304        odimyu = timelabel
1305        timeaxis = 'y'
1306        odimxv = objsf.variables[vdimxn]
1307        odimxu = odimxv.getncattr('units')
1308        timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1309    else:
1310        print errormsg
1311        print '  ' + fname + ": time variable '" + timename + "' not found!!"
1312        quit(-1)
1313
1314    shading_nx = np.zeros((2), dtype=np.float)
1315    shading_nx[0] = np.float(shadminmax.split(',')[0])
1316    shading_nx[1] = np.float(shadminmax.split(',')[1])
1317
1318    closeval = drw.Str_Bool(close)
1319
1320    drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu,      \
1321      dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1322      timepos, timelabels, False)
1323
1324# Line values
1325##
1326    ncfilel = ncfile.split(',')[1]
1327
1328    vnamelfig = values.split(':')[0].split(',')[1]
1329    varnl = varn.split(',')[1]
1330
1331    objlf = NetCDFFile(ncfilel,'r')
1332    objlvar = objlf.variables[varnl]
1333
1334    linevals = objlvar[:]
1335    if reva0 == 'tranpose':
1336        plt.plot (linevals, odimxv, '-', color='k')
1337    else:
1338        plt.plot (odimxv, linevals, '-', color='k')
1339
1340    objsf.close()
1341    objsl.close()
1342
1343    return
1344
1345def draw_barbs(ncfile, values, varns):
1346    """ Function to plot wind barbs
1347      values= [dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:[mapvalues]:
1348        [gtit]:[kindfig]:[figuren]
1349        'X/Y/Z/T'|[dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
1350          [dimname]: name of the dimension in the file
1351          [vardimname]: name of the variable with the values for the dimension in the file
1352          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
1353          No value takes all the range of the dimension
1354        [vecvals]= [frequency],[color],[length]
1355          [frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
1356            'auto', computed automatically to have 20 vectors along each axis)
1357          [color]: color of the vectors ('auto', for 'red')
1358          [length]: length of the wind barbs ('auto', for 9)
1359        [windlabs]= [windname],[windunits]
1360          [windname]: name of the wind variable in the graph
1361          [windunits]: units of the wind variable in the graph ('None', for the value in the file)
1362        [mapvalues]= map characteristics: [proj],[res]
1363          see full documentation: http://matplotlib.org/basemap/
1364          [proj]: projection
1365            * 'cyl', cilindric
1366            * 'lcc', lambert conformal
1367          [res]: resolution:
1368            * 'c', crude
1369            * 'l', low
1370            * 'i', intermediate
1371            * 'h', high
1372            * 'f', full
1373        gtit= title of the graph ('|', for spaces)
1374        kindfig= kind of figure
1375        figuren= name of the figure
1376      ncfile= file to use
1377      varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
1378    """
1379    fname = 'draw_barbs'
1380
1381    if values == 'h':
1382        print fname + '_____________________________________________________________'
1383        print draw_barbs.__doc__
1384        quit()
1385
1386    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:' +  \
1387      '[mapvalues]:[gtit]:[kindfig]:[figuren]'
1388 
1389    drw.check_arguments(fname,values,expectargs,':')
1390
1391    dimvals = values.split(':')[0]
1392    vecvals = values.split(':')[1]
1393    windlabels = values.split(':')[2]
1394    mapvalues = values.split(':')[3]
1395    gtit = values.split(':')[4]
1396    kindfig = values.split(':')[5]
1397    figuren = values.split(':')[6]
1398
1399    of = NetCDFFile(ncfile,'r')
1400
1401    dims = {}
1402    for dimv in dimvals.split(','):
1403        dns = dimv.split('|')
1404        dims[dns[0]] = [dns[1], dns[2], dns[3]]
1405
1406    varNs = []
1407    for dn in dims.keys():
1408        if dn == 'X':
1409            varNs.append(dims[dn][1])
1410            dimx = len(of.dimensions[dims[dn][0]])
1411        elif dn == 'Y':
1412            varNs.append(dims[dn][1])
1413            dimy = len(of.dimensions[dims[dn][0]])
1414
1415    ivar = 0
1416    for wvar in varns.split(','):
1417        if not drw.searchInlist(of.variables.keys(), wvar):
1418            print errormsg
1419            print '  ' + fname + ": file does not have variable '" + wvar + "' !!"
1420            quit(-1)
1421        if ivar == 0:
1422            varNs.append(wvar)
1423        else:
1424            varNs.append(wvar)
1425
1426    ivar = 0
1427    for varN in varNs:
1428        varslice = []
1429
1430        ovarN = of.variables[varN]
1431        vard = ovarN.dimensions
1432        for vdn in vard:
1433            found = False
1434            for dd in dims.keys():
1435                if dims[dd][0] == vdn:
1436                    if dims[dd][2].find('@') != -1:
1437                        rvals = dims[dd][2].split('@')
1438                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
1439                    elif dims[dd][2] == '-1':
1440                        varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
1441                    else:
1442                        varslice.append(int(dims[dd][2]))
1443
1444                    found = True
1445                    break
1446            if not found:
1447                varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
1448
1449        if varN == dims['X'][1]:
1450            lonvals0 = np.squeeze(ovarN[tuple(varslice)])
1451        elif varN == dims['Y'][1]:
1452            latvals0 = np.squeeze(ovarN[tuple(varslice)])
1453        elif ivar == 2:
1454            uwvals = np.squeeze(np.array(ovarN[tuple(varslice)]))
1455        elif ivar == 3:
1456            vwvals = np.squeeze(ovarN[tuple(varslice)])           
1457
1458        ivar = ivar + 1
1459
1460#    print 'Final shapes:',lonvals0.shape,':',latvals0.shape,':',uwvals.shape,':',
1461#      vwvals.shape
1462
1463    if len(uwvals.shape) != 2 or len(vwvals.shape) != 2:
1464        print errormsg
1465        print '  ' + fname + ': wrong size of the wind fields! they must be ' +      \
1466          '2-dimensional!'
1467        print '    u-winds shape:',uwvals.shape,'dims:',of.variables[varNs[2]]
1468        print '    v-winds shape:',vwvals.shape,'dims:',of.variables[varNs[3]]
1469        print '      provide more values for their dimensions!!'
1470        quit(-1)
1471
1472    if len(lonvals0.shape) == 1:
1473        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
1474    else:
1475        lonvals = lonvals0
1476        latvals = latvals0
1477
1478# Vecor values
1479    if vecvals.split(',')[0] == 'None':
1480        freqv = None
1481    else:
1482        freqv = vecvals.split(',')[0] 
1483    colorv = vecvals.split(',')[1]
1484    lengthv = vecvals.split(',')[2]
1485
1486# Vector labels
1487    windname = windlabels.split(',')[0]
1488    windunits = windlabels.split(',')[1]
1489
1490    drw.plot_barbs(lonvals, latvals, uwvals, vwvals, freqv, colorv, lengthv, 
1491      windname, windunits, mapvalues, gtit, kindfig, figuren)
1492
1493    return
1494 
1495def draw_topo_geogrid(ncfile, values):
1496    """ plotting geo_em.d[nn].nc topography from WPS files
1497    draw_topo_geogrid(ncfile, values)
1498      ncfile= geo_em.d[nn].nc file to use
1499      values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1500        [min/max]Topo: minimum and maximum values of topography to draw
1501        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1502        title: title of the graph ('!' for spaces)
1503        graphic_kind: kind of figure (jpg, pdf, png)
1504        mapvalues: map characteristics [proj],[res]
1505          see full documentation: http://matplotlib.org/basemap/
1506          [proj]: projection
1507            * 'cyl', cilindric
1508            * 'lcc', lambert conformal
1509          [res]: resolution:
1510            * 'c', crude
1511            * 'l', low
1512            * 'i', intermediate
1513            * 'h', high
1514            * 'f', full
1515    """
1516    fname = 'draw_topo_geogrid'
1517
1518    if values == 'h':
1519        print fname + '_____________________________________________________________'
1520        print draw_topo_geogrid.__doc__
1521        quit()
1522
1523    expectargs = '[minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]'
1524 
1525    drw.check_arguments(fname,values,expectargs,':')
1526
1527    mintopo = np.float(values.split(':')[0].split(',')[0])
1528    maxtopo = np.float(values.split(':')[0].split(',')[1])
1529
1530    lonlatLS = values.split(':')[1]
1531    lonlatLv = lonlatLS.split(',')[0]
1532
1533    if lonlatLv == 'None':
1534        lonlatL = None
1535    else:
1536        lonlatL = np.zeros((4), dtype=np.float)
1537        lonlatL[0] = np.float(lonlatLS.split(',')[0])
1538        lonlatL[1] = np.float(lonlatLS.split(',')[1])
1539        lonlatL[2] = np.float(lonlatLS.split(',')[2])
1540        lonlatL[3] = np.float(lonlatLS.split(',')[3])
1541
1542    grtit = values.split(':')[2].replace('!',' ')
1543    kindfig = values.split(':')[3]
1544    mapvalues = values.split(':')[4]
1545
1546    if not os.path.isfile(ncfile):
1547        print errormsg
1548        print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1549        quit(-1)   
1550
1551    objdomf = NetCDFFile(ncfile, 'r')
1552   
1553    objhgt = objdomf.variables['HGT_M']
1554    objlon = objdomf.variables['XLONG_M']
1555    objlat = objdomf.variables['XLAT_M']
1556
1557    topography = objhgt[0,:,:]
1558
1559    drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL,     \
1560      grtit, kindfig, mapvalues, True)
1561
1562    objdomf.close()
1563
1564    return
1565
1566def draw_topo_geogrid_boxes(ncfiles, values):
1567    """ plotting different geo_em.d[nn].nc topography from WPS files
1568    draw_topo_geogrid_boxes(ncfile, values)
1569      ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference)
1570      values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels]:[legvals]
1571        [min/max]Topo: minimum and maximum values of topography to draw
1572        lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1573        title: title of the graph ('!' for spaces)
1574        graphic_kind: kind of figure (jpg, pdf, png)
1575        mapvalues: map characteristics [proj],[res]
1576          see full documentation: http://matplotlib.org/basemap/
1577          [proj]: projection
1578            * 'cyl', cilindric
1579            * 'lcc', lambert conformal
1580          [res]: resolution:
1581            * 'c', crude
1582            * 'l', low
1583            * 'i', intermediate
1584            * 'h', high
1585            * 'f', full
1586        [legvals]=[locleg]|[fontsize]:
1587          [locleg]: location of the legend (0, autmoatic)
1588            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1589            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1590            9: 'upper center', 10: 'center'
1591          [fontsize]: font size for the legend (auto for 12)
1592        labels= labels to write in the graph
1593    """
1594#    import matplotlib as mpl
1595#    mpl.use('Agg')
1596    import matplotlib.pyplot as plt
1597
1598    fname = 'draw_topo_geogrid_boxes'
1599
1600    if values == 'h':
1601        print fname + '_____________________________________________________________'
1602        print draw_topo_geogrid_boxes.__doc__
1603        quit()
1604
1605    mintopo = np.float(values.split(':')[0].split(',')[0])
1606    maxtopo = np.float(values.split(':')[0].split(',')[1])
1607
1608    lonlatLS = values.split(':')[1]
1609    lonlatLv = lonlatLS.split(',')[0]
1610
1611    if lonlatLv == 'None':
1612        lonlatL = None
1613    else:
1614        lonlatL = np.zeros((4), dtype=np.float)
1615        lonlatL[0] = np.float(lonlatLS.split(',')[0])
1616        lonlatL[1] = np.float(lonlatLS.split(',')[1])
1617        lonlatL[2] = np.float(lonlatLS.split(',')[2])
1618        lonlatL[3] = np.float(lonlatLS.split(',')[3])
1619
1620    grtit = values.split(':')[2].replace('!', ' ')
1621    kindfig = values.split(':')[3]
1622    mapvalues = values.split(':')[4]
1623    labels = values.split(':')[5]
1624    legloc = int(values.split(':')[6])
1625
1626    ncfile = ncfiles.split(',')[0]
1627    if not os.path.isfile(ncfile):
1628        print errormsg
1629        print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1630        quit(-1)   
1631
1632    objdomf = NetCDFFile(ncfile, 'r')
1633   
1634    objhgt = objdomf.variables['HGT_M']
1635    objlon0 = objdomf.variables['XLONG_M']
1636    objlat0 = objdomf.variables['XLAT_M']
1637
1638    topography = objhgt[0,:,:]
1639
1640    Nfiles = len(ncfiles.split(','))
1641    boxlabels = labels.split(',')
1642
1643    Xboxlines = []
1644    Yboxlines = []
1645
1646    for ifile in range(Nfiles):
1647        ncfile = ncfiles.split(',')[ifile]
1648#        print ifile, ncfile
1649        if not os.path.isfile(ncfile):
1650            print errormsg
1651            print '  ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1652            quit(-1)   
1653
1654        objdomfi = NetCDFFile(ncfile, 'r')
1655   
1656        objlon = objdomfi.variables['XLONG_M']
1657        objlat = objdomfi.variables['XLAT_M']
1658
1659        dx = objlon.shape[2]
1660        dy = objlon.shape[1]
1661
1662        Xboxlines.append(objlon[0,0,:])
1663        Yboxlines.append(objlat[0,0,:])
1664        Xboxlines.append(objlon[0,dy-1,:])
1665        Yboxlines.append(objlat[0,dy-1,:])
1666        Xboxlines.append(objlon[0,:,0])
1667        Yboxlines.append(objlat[0,:,0])
1668        Xboxlines.append(objlon[0,:,dx-1])
1669        Yboxlines.append(objlat[0,:,dx-1])
1670
1671        objdomfi.close()
1672
1673    drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels,         \
1674      objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, legloc,\
1675      True)
1676
1677    objdomf.close()
1678
1679    return
1680
1681def movievalslice(origslice, dimmovien, framenum):
1682    """ Function to provide variable slice according to a geneation of a movie
1683    movievals(origslice, dimmovien, framenum)
1684      [origslice]= slice original as [dimname1]|[val1],[...,[dimnameN]|[valN]]
1685        ([val] = -1, full length)
1686      [dimmovien]= name of the dimension to produce the movie
1687      [framenum]= value of the frame to substitue in [origslice] as
1688        [dimmovien]|[framenum]
1689    >>> movievalslice('East_West|-1,North_South|-1,Time|2','Time',0)
1690    East_West|-1,North_South|-1,Time|0
1691    """
1692
1693    fname = 'movievalslice'
1694
1695    if origslice == 'h':
1696        print fname + '_____________________________________________________________'
1697        print movievalslice.__doc__
1698        quit()
1699   
1700    dims = origslice.split(',')
1701
1702    movieslice = ''
1703    idim = 0
1704
1705    for dimn in dims:
1706        dn = dimn.split('|')[0]
1707        if dn == dimmovien:
1708            movieslice = movieslice + dn + '|' + str(framenum)
1709        else:
1710            movieslice = movieslice + dimn
1711        if idim < len(dims)-1: movieslice = movieslice + ','
1712
1713        idim = idim + 1
1714
1715    return movieslice
1716
1717class Capturing(list):
1718    """ Class to capture function output as a list
1719    from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call
1720    """
1721#    from cStringIO import StringIO
1722
1723    def __enter__(self):
1724        self._stdout = sys.stdout
1725        sys.stdout = self._stringio = StringIO()
1726        return self
1727    def __exit__(self, *args):
1728        self.extend(self._stringio.getvalue().splitlines())
1729        sys.stdout = self._stdout
1730
1731def create_movie(netcdfile, values, variable):
1732    """ Function to create a movie assuming ImageMagick installed!
1733      values= [graph]#[movie_dimension]#[graph_values]
1734        [graph]: which graphic
1735        [movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval]
1736          [dimnmovie]; name of the dimension from which make the movie
1737          [dimvmovie]; name of the variable with the values of the dimension
1738          [moviedelay]; delay between frames
1739          [interval]; [beg]@[end]@[freq] or -1 (all)
1740        [graph_values]: values to generate the graphic
1741      netcdfile= netCDF file
1742      variable= variable to use (when applicable)
1743    """ 
1744    fname = 'create_movie'
1745
1746    if values == 'h':
1747        print fname + '_____________________________________________________________'
1748        print create_movie.__doc__
1749        quit()
1750
1751    graph = values.split('#')[0]
1752    movie_dim = values.split('#')[1]
1753    graph_vals = values.split('#')[2]
1754
1755    ncobj = NetCDFFile(netcdfile, 'r')
1756
1757# Movie dimension
1758##
1759    dimnmovie = movie_dim.split('@')[0]
1760    dimvmovie = movie_dim.split('@')[1]
1761    moviedelay = movie_dim.split('@')[2]
1762    moviebeg = int(movie_dim.split('@')[3])
1763
1764    if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie):
1765        print errormsg
1766        print '  ' + fname + ": file '" + netcdfile + "' has not dimension '" +      \
1767          dimnmovie + "' !!!"
1768        quit(-1)
1769
1770    objdmovie = ncobj.dimensions[dimnmovie]
1771    dmovie = len(objdmovie)
1772    if moviebeg != -1:
1773        moviend = int(movie_dim.split('@')[4])
1774        moviefreq = int(movie_dim.split('@')[5])
1775    else:
1776        moviebeg = 0
1777        moviend = dmovie
1778        moviefreq = 1
1779
1780    if dimvmovie == 'WRFTimes':
1781        objvdmovie = ncobj.variables['Times']
1782        vdmovieunits = ''
1783        valsdmovie = []
1784        for it in range(objvdmovie.shape[0]):
1785            valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:],           \
1786              'WRFdatetime', 'Y/m/d H-M-S'))
1787    elif dimvmovie == 'CFtime':
1788        objvdmovie = ncobj.variables['time']
1789        vdmovieunits = ''
1790        print objvdmovie.units
1791        valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard',  \
1792          objvdmovie[:])
1793        valsdmovie = []
1794        for it in range(objvdmovie.shape[0]):
1795            valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:],          \
1796              'matYmdHMS', 'Y/m/d H-M-S'))
1797    else:
1798        if  not drw.searchInlist(ncobj.variables.keys(),dimvmovie):
1799            print errormsg
1800            print '  ' + fname + ": file '" + netcdfile + "' has not variable '" +   \
1801              dimvmovie + "' !!!"
1802            quit(-1)
1803        vdmovieunits = objvdmovie.getncattr('units')
1804        objvdmovie = ncobj.variables[dimvmovie]
1805        if len(objvdmovie.shape) == 1:
1806            vasldmovie = objvdmovie[:]
1807        else:
1808            print errormsg
1809            print '  ' + fname + ': shape', objvdmovie.shape, 'of variable with ' +  \
1810              'dimension movie values not ready!!!'
1811            quit(-1)
1812
1813    ncobj.close()
1814    os.system('rm frame_*.png > /dev/null')
1815
1816# graphic
1817##
1818    if graph == 'draw_2D_shad':
1819        graphvals = graph_vals.split(':')
1820
1821        for iframe in range(moviebeg,moviend,moviefreq):
1822            iframeS = str(iframe).zfill(4)
1823
1824            drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \
1825              5, 'frames')
1826            titgraph = dimnmovie + '|=|' + str(valsdmovie[iframe]) + '|' +           \
1827              vdmovieunits
1828
1829            graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe)
1830            graphvals[6] = titgraph
1831            graphvals[7] = 'png'
1832
1833            graphv = drw.numVector_String(graphvals, ":")
1834
1835            with Capturing() as output:
1836                draw_2D_shad(netcdfile, graphv, variable)
1837
1838            os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png')
1839    else:
1840        print errormsg
1841        print '  ' + fname + ": graphic '" +  graph + "' not defined !!!"
1842        quit(-1)
1843
1844    os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif')
1845
1846    print "Succesfuly creation of movie file 'create_movie.gif' !!!"
1847
1848    return
1849
1850def draw_lines(ncfilens, values, varname):
1851    """ Function to draw different lines at the same time from different files
1852    draw_lines(ncfilens, values, varname):
1853      ncfilens= [filen] ',' separated list of netCDF files
1854      values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[legvals]:[colns]:[lines]
1855       [points]:[lwdths]:[psizes]:[freqv]:[figname]:[graphk]:[close]
1856        [dimvname]: ',' list of names of the variable with he values of the common dimension
1857        [valuesaxis]: which axis will be used for the values ('x', or 'y')
1858        [dimtit]: title for the common dimension
1859        [leglabels]: ',' separated list of names for the legend
1860        [vartit]: name of the variable in the graph
1861        [title]: title of the plot ('|' for spaces)
1862        [legvals]=[locleg]|[fontsize]:
1863          [locleg]: location of the legend (0, autmoatic)
1864            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1865            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1866            9: 'upper center', 10: 'center'
1867          [fontsize]: font size for the legend (auto for 12)
1868        [colns]= ',' list of color names ('None' for automatic, single value for all the same)
1869        [lines]= ',' list of style of lines ('None' for automatic, single value for all the same)
1870        [points]= '@' list of style of points ('None' for automatic, single value for all the same)
1871        [lwdths]= ',' list of withs of lines ('None' for automatic, single value for all the same)
1872        [psizes]= ',' list of size of points ('None' for automatic, single value for all the same)
1873        [freqv]= frequency of values ('all' for all values)
1874        [figname]= name of the figure
1875        [graphk]: kind of the graphic
1876        [close]: should figure be closed (finished)
1877      varname= variable to plot
1878      values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1879    """
1880
1881    fname = 'draw_lines'
1882
1883    if values == 'h':
1884        print fname + '_____________________________________________________________'
1885        print draw_lines.__doc__
1886        quit()
1887
1888    expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:' +    \
1889      '[legvals]:[colns]:[lines]:[points]:[lwdths]:[psizes]:[freqv]:[figname]:' +    \
1890      '[graphk]:[close]'
1891    drw.check_arguments(fname,values,expectargs,':')
1892
1893    ncfiles = ncfilens.split(',')
1894    dimvnames = values.split(':')[0]
1895    valuesaxis = values.split(':')[1]
1896    dimtit = values.split(':')[2]
1897    leglabels = gen.latex_text(values.split(':')[3])
1898    vartit = values.split(':')[4]
1899    title = values.split(':')[5].replace('|',' ')
1900    legvals = values.split(':')[6]
1901    colns = gen.str_list(values.split(':')[7], ',')
1902    lines = gen.str_list(values.split(':')[8], ',')
1903    points = gen.str_list(values.split(':')[9], '@')
1904    lwdths = gen.str_list_k(values.split(':')[10], ',', 'R')
1905    psizes = gen.str_list_k(values.split(':')[11], ',', 'R')
1906    freqv0 = values.split(':')[12]
1907    figname = values.split(':')[13]
1908    graphk = values.split(':')[14]
1909    close = gen.Str_Bool(values.split(':')[15])
1910
1911    Nfiles = len(ncfiles)
1912
1913# Getting trajectotries
1914##
1915
1916    varvalues = []
1917    dimvalues = []
1918
1919    print '  ' + fname
1920    ifn = 0
1921    for ifile in ncfiles:
1922        filen = ifile.split('@')[0]
1923
1924        print '    filen:',filen
1925
1926        if not os.path.isfile(filen):
1927            print errormsg
1928            print '  ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1929            quit(-1)
1930
1931        objfile = NetCDFFile(filen, 'r')
1932
1933        if dimvnames.find(',') != -1:
1934            dimvname = dimvnames.split(',')
1935        else:
1936            dimvname = [dimvnames]
1937   
1938        found = False
1939        for dvn in dimvname:
1940            if objfile.variables.has_key(dvn):
1941                found = True
1942                break
1943   
1944        if not found:
1945            print errormsg
1946            print '  ' + fname + ": netCDF file '" + filen +                         \
1947              "' does not have variables '" + dimvnames + "' !!"
1948            quit(-1)
1949
1950        if not objfile.variables.has_key(varname):
1951            print errormsg
1952            print '  ' + fname + ": netCDF file '" + filen +                         \
1953              "' does not have variable '" + varname + "' !!"
1954            quit(-1)
1955
1956        vvobj = objfile.variables[varname]
1957        if len(vvobj.shape) != 1:
1958            print errormsg
1959            print '  ' + fname + ': wrong shape:',vvobj.shape," of variable '" +     \
1960              varname +  "' !!"
1961            quit(-1)
1962
1963        for dimvn in dimvname:
1964            if drw.searchInlist(objfile.variables, dimvn):
1965                vdobj = objfile.variables[dimvn]
1966                if len(vdobj.shape) != 1:
1967                    print errormsg
1968                    print '  ' + fname + ': wrong shape:',vdobj.shape,               \
1969                      " of variable '" + dimvn +  "' !!"
1970                    quit(-1)
1971                break
1972
1973        varvalues.append(vvobj[:])
1974        dimvalues.append(vdobj[:])
1975
1976        if ifn == 0:
1977            varunits = vvobj.units
1978
1979        objfile.close()
1980
1981        ifn = ifn + 1
1982
1983    if freqv0 == 'all':
1984        freqv = None
1985    else:
1986        freqv = int(freqv0)
1987
1988    # Legend
1989    locleg, legfontsize = drw.legend_values(legvals,'|')
1990
1991    drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','),   \
1992      vartit, varunits, title, locleg, legfontsize, colns, lines, points, lwdths,    \
1993      psizes, freqv, figname, graphk, close)
1994
1995    return
1996
1997def draw_lines_time(ncfilens, values, varname0):
1998    """ Function to draw different lines at the same time from different files with times
1999    draw_lines_time(ncfilens, values, varname):
2000      ncfilens= [filen] ',' separated list of netCDF files
2001      values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[rangevals];[timevals];
2002        [legvals];[graphk];[collines];[points];[linewidths];[pointsizes];[pointfreq];[period];[close]
2003        [dimvname]: ',' list of names of the variables with he values of the common dimension
2004        [valuesaxis]: which axis will be used for the values ('x', or 'y')
2005        [dimtit]: title for the common dimension ('|' for spaces)
2006        [leglabels]: ',' separated list of names for the legend ('None', no legend '!' for spaces)
2007        [vartit]: name of the variable in the graph
2008        [title]: title of the plot ('|' for spaces)
2009        [rangevals]: Range of the axis with the values ('None' for 'auto','auto')
2010          [vmin],[vmax]: minimum and maximum values where [vmNN] can also be:
2011            'auto': the computed minimumm or maximum of the values 
2012        [timevals]: [timen]|[units]|[kind]|[tfmt] time labels characteristics
2013           [timen]; name of the time variable
2014           [units]; units string according to CF conventions ([tunits] since
2015             [YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
2016           [kind]; kind of output
2017             'Nval': according to a given number of values as 'Nval',[Nval]
2018             'exct': according to an exact time unit as 'exct',[tunit];
2019               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2020                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2021                'l': milisecond
2022           [tfmt]; desired format
2023        [legvals]=[locleg]|[fontsize]:
2024          [locleg]: location of the legend (0, autmoatic)
2025            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2026            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2027            9: 'upper center', 10: 'center'
2028          [fontsize]: font size for the legend (auto for 12)
2029        [graphk]: kind of the graphic
2030        [lines]: ',' list of type of lines, None for automatic, single value all the same
2031        [collines]: ',' list of colors for the lines, None for automatic, single
2032          value all the same
2033        [points]: '@' list of type of points for the lines, None for automatic, single
2034          value all the same
2035        [linewidths]: ',' list of widths for the lines, None for automatic, single
2036          value all the same
2037        [pointsizes]: ',' list of widths for the lines, None for automatic, single
2038          value all the same
2039        [pointfreq]: frequency of point plotting, 'all' for all time steps
2040        [period]: which period to plot
2041          '-1': all period
2042          [beg],[end]: beginning and end of the period in reference time-units of first file
2043        [close]: Whether figure should be finished or not
2044      varname0= ',' list of variable names to plot (assuming only 1 variable per file)
2045      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'
2046    """
2047
2048    fname = 'draw_lines_time'
2049
2050    if values == 'h':
2051        print fname + '_____________________________________________________________'
2052        print draw_lines_time.__doc__
2053        quit()
2054
2055    expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
2056    expectargs = expectargs + '[rangevals];[timevals];[legvals];[graphk];[lines];'
2057    expectargs = expectargs + '[collines];[points];[linewidths];[pointsizes];'
2058    expectargs = expectargs + '[pointfreq];[period];[close]'
2059    drw.check_arguments(fname,values,expectargs,';')
2060
2061    ncfiles = ncfilens.split(',')
2062    dimvname0 = values.split(';')[0]
2063    valuesaxis = values.split(';')[1]
2064    dimtit = values.split(';')[2].replace('|',' ')
2065    leglabels = gen.latex_text(values.split(';')[3]).replace('!',' ')
2066    vartit = values.split(';')[4]
2067    title = values.split(';')[5].replace('|',' ')
2068    rangevals = values.split(';')[6]
2069    timevals = values.split(';')[7]
2070    legvals = values.split(';')[8]
2071    graphk = values.split(';')[9]
2072    lines0 = values.split(';')[10]
2073    collines0 = values.split(';')[11]
2074    points0 = values.split(';')[12]
2075    linewidths0 = values.split(';')[13]
2076    pointsizes0 = values.split(';')[14]
2077    pointfreq0 = values.split(';')[15]
2078    period = values.split(';')[16]
2079    close = gen.Str_Bool(values.split(':')[17])
2080
2081    Nfiles = len(ncfiles)
2082
2083# Multiple variable-dimension names?
2084    if dimvname0.find(',') != -1:
2085        dimvname = dimvname0.split(',')
2086    else:
2087        dimvname = [dimvname0]
2088
2089# Multiple variables?
2090    if varname0.find(',') != -1:
2091        varname = varname0.split(',')
2092    else:
2093        varname = [varname0]
2094
2095# Multiple lines types?
2096    if lines0.find(',') != -1:
2097        lines = lines0.split(',')
2098    elif lines0 == 'None':
2099        lines = None
2100    else:
2101        lines = []
2102        for il in range(Nfiles):
2103            lines.append(lines0)
2104
2105# Multiple color names?
2106    if collines0.find(',') != -1:
2107        collines = collines0.split(',')
2108    elif collines0 == 'None':
2109        collines = None
2110    else:
2111        collines = []
2112        for ip in range(Nfiles):
2113            collines.append(collines0)
2114
2115# Multiple point types?
2116    if points0.find(',') != -1:
2117        if len(points0) == 1:
2118            points = []
2119            for ip in range(Nfiles):
2120                points.append(points0)
2121        else:
2122            points = points0.split('@')
2123    elif points0 == 'None':
2124        points = None
2125    else:
2126        points = []
2127        for ip in range(Nfiles):
2128            points.append(points0)
2129
2130# Multiple line sizes?
2131    if linewidths0.find(',') != -1:
2132        linewidths = []
2133        Nlines = len(linewidths0.split(','))
2134        for il in range(Nlines):
2135          linewidths.append(np.float(linewidths0.split(',')[il]))
2136    elif linewidths0 == 'None':
2137        linewidths = None
2138    else:
2139        linewidths = [np.float(linewidths0)]
2140
2141# Multiple point sizes?
2142    if pointsizes0.find(',') != -1:
2143        pointsizes = []
2144        Npts = len(pointsizes0.split(','))
2145        for ip in Npts:
2146          pointsizes.append(np.float(pointsizes0.split(',')[ip]))
2147    elif pointsizes0 == 'None':
2148        pointsizes = None
2149    else:
2150        pointsizes = [np.float(pointsizes0)]
2151
2152    timename = timevals.split('|')[0]
2153    timeunit = timevals.split('|')[1].replace('!',' ')
2154    timekind = timevals.split('|')[2]
2155    timefmt = timevals.split('|')[3]
2156
2157    if rangevals == 'None':
2158        valmin = 'auto'
2159        valmax = 'auto'
2160    else:
2161        valmin = rangevals.split(',')[0]
2162        valmax = rangevals.split(',')[1]
2163        if valmin != 'auto': valmin = np.float(valmin)
2164        if valmax != 'auto': valmax = np.float(valmax)
2165
2166    # Legend
2167    locleg, legfontsize = drw.legend_values(legvals,'|')
2168
2169# Getting values
2170##
2171    varvalues = []
2172    dimvalues = []
2173    timvalues = []
2174    timvals0 = timvalues
2175
2176    print '  ' + fname
2177    ifn = 0
2178    mintval = 1.e20
2179    maxtval = -1.e20
2180
2181    for ifile in ncfiles:
2182        filen = ifile.split('@')[0]
2183
2184        print '    filen:',filen
2185
2186        if not os.path.isfile(filen):
2187            print errormsg
2188            print '  ' + fname + ": netCDF file '" + filen + "' does not exist !!"
2189            quit(-1)
2190
2191        objfile = NetCDFFile(filen, 'r')
2192
2193        founddvar = False
2194        for dvar in dimvname:
2195            if objfile.variables.has_key(dvar):
2196                founddvar = True
2197                vdobj = objfile.variables[dvar]
2198                if len(vdobj.shape) != 1:
2199                    print errormsg
2200                    print '  ' + fname + ': wrong shape:',vdobj.shape," of " +       \
2201                      "variable '" + dvar +  "' !!"
2202                    quit(-1)
2203                break
2204        if not founddvar:
2205            print errormsg
2206            print '  ' + fname + ": netCDF file '" + filen +                         \
2207            "' has any variable '", dimvname, "' !!"
2208            quit(-1)
2209
2210        foundvar = False
2211        for var in varname:
2212            if objfile.variables.has_key(var):
2213                foundvar = True
2214                vvobj = objfile.variables[var]
2215                if len(vvobj.shape) != 1:
2216                    print errormsg
2217                    print '  ' + fname + ': wrong shape:',vvobj.shape," of " +       \
2218                      "variable '" + var +  "' !!"
2219                    quit(-1)
2220
2221                break
2222        if not foundvar:
2223            print errormsg
2224            print '  ' + fname + ": netCDF file '" + filen +                         \
2225              "' has any variable '", varname, "' !!"
2226            quit(-1)
2227        if vdobj.units.find('month') != -1:
2228            print warnmsg
2229            print '  ' + fname + ": transforming time units from 'months' to 'days'!!"
2230            timevals0, tunits0 = gen.CFmonthU_daysU(vdobj[:], vdobj.units) 
2231        else:
2232            timevals0 = vdobj[:]
2233            tunits0 = str(vdobj.units) 
2234
2235# Getting period
2236        if ifn > 0: 
2237# Referring all times to the same reference time!
2238            reftvals = gen.coincident_CFtimes(timevals0, timeunit, tunits0)
2239        else:
2240            reftvals = timevals0
2241
2242        dimt = len(vdobj[:])
2243
2244        if period == '-1':
2245            varvalues.append(vvobj[:])
2246            dimvalues.append(reftvals)
2247            mindvals = np.min(reftvals)
2248            maxdvals = np.max(reftvals)
2249        else:
2250            ibeg=-1
2251            iend=-1
2252            tbeg = np.float(period.split(',')[0])
2253            tend = np.float(period.split(',')[1])
2254
2255            for it in range(dimt-1):
2256                if reftvals[it] <= tbeg and reftvals[it+1] > tbeg: ibeg = it
2257                if reftvals[it] <= tend and reftvals[it+1] > tend: iend = it + 1
2258                if ibeg != -1 and iend != -1: break
2259
2260            if ibeg == -1 and iend == -1:
2261                print warnmsg
2262                print '  ' + fname + ': Period:',tbeg,',',tend,'not found!!'
2263                print '    ibeg:',ibeg,'iend:',iend
2264                print '    period in file:',np.min(reftvals), np.max(reftvals)
2265                print '    getting all the period in file !!!'
2266                ibeg = 0
2267                iend = dimt
2268            elif iend == -1:
2269                iend = dimt
2270                print warnmsg
2271                print '  ' + fname + ': end of Period:',tbeg,',',tend,'not found!!'
2272                print '    getting last available time instead'
2273                print '    ibeg:',ibeg,'iend:',iend
2274                print '    period in file:',np.min(reftvals), np.max(reftvals)
2275            elif ibeg == -1:
2276                ibeg = 0
2277                print warnmsg
2278                print '  ' + fname + ': beginning of Period:',tbeg,',',tend,         \
2279                  'not found!!'
2280                print '    getting first available time instead'
2281                print '    ibeg:',ibeg,'iend:',iend
2282                print '    period in file:',np.min(reftvals), np.max(reftvals)
2283
2284            varvalues.append(vvobj[ibeg:iend])
2285            dimvalues.append(reftvals[ibeg:iend])
2286            mindvals = np.min(reftvals[ibeg:iend])
2287            maxdvals = np.max(reftvals[ibeg:iend])
2288
2289            dimt = iend - ibeg
2290
2291        if mindvals < mintval: mintval = mindvals
2292        if maxdvals > maxtval: maxtval = maxdvals
2293        print '  ' + fname + ": file '" + filen + "' period:", mindvals, '->', maxdvals
2294
2295        if ifn == 0:
2296            varunits = drw.units_lunits(vvobj.units)
2297
2298        objfile.close()
2299
2300        ifn = ifn + 1
2301
2302# Times
2303
2304    dtvals = (maxtval - mintval)/dimt
2305#    dti = mintval-dtvals/2.
2306#    dte = maxtval+dtvals/2.
2307    dti = mintval
2308    dte = maxtval
2309    tvals = np.arange(dti, dte, dtvals)
2310
2311    dtiS = drw.datetimeStr_conversion(str(dti) + ',' + timeunit, 'cfTime',           \
2312      'Y/m/d H-M-S')
2313    dteS = drw.datetimeStr_conversion(str(dte) + ',' + timeunit, 'cfTime',           \
2314      'Y/m/d H-M-S')
2315
2316    print '  ' + fname + ': plotting from: ' + dtiS + ' to ' + dteS
2317
2318    timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
2319
2320#    print 'Lluis min/max tval +/- dtval/2:', mintval-dtvals/2., maxtval+dtvals/2.,'dt:', len(tvals)
2321#    for it in range(len(timepos)):
2322#        print timepos[it], timelabels[it]
2323
2324    if leglabels != 'None':
2325        legvals = leglabels.split(',')
2326    else:
2327        legvals = None
2328
2329    if pointfreq0 == 'all':
2330        pointfreq = None
2331    else:
2332        pointfreq = int(pointfreq0)
2333
2334    drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, legvals, vartit,   \
2335      varunits, timepos, timelabels, title, locleg, legfontsize, graphk, valmin,     \
2336      valmax, lines, collines, points, linewidths, pointsizes, pointfreq, close)
2337
2338    return
2339
2340def draw_Neighbourghood_evol(filen, values, variable):
2341    """ Function to draw the temporal evolution of a neighbourghood around a point
2342    draw_Neighbourghood_evol(filen, values, variable)
2343      filen= netCDF file name
2344      values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
2345       [timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
2346        [dimsval]: [dimn1]|[val1]|[dimv1],...,[dimnN]|[valN]|[dimvN] dimension names, values to get
2347          (-1, for all; no name/value pair given full length) and variable with values of the dimension
2348          NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
2349        [neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
2350        [Nneig]: Number of grid points of the full side of the box (odd value)
2351        [Ncol]: Number of columns ('auto': square final plot)
2352        [gvarname]: name of the variable to appear in the graph
2353        [timetits]: [titX],[titY] titles of the axes ('|' for spaces)
2354        [tkinds]: [tkindX]|[tkindY] kinds of time to appear in the graph
2355          'Nval': according to a given number of values as 'Nval',[Nval]
2356          'exct': according to an exact time unit as 'exct',[tunit];
2357            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2358              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2359              'l': milisecond
2360        [timefmts]: [tfmtX],[tfmtY] format of the time labels
2361        [gtitle]: title of the graphic ('|' for spaces)
2362        [shadxtrms]: Extremes for the shading
2363        [cbar]: colorbar to use
2364        [gkind]: kind of graphical output
2365        [ofile]: True/False whether the netcdf with data should be created or not
2366      variable= name of the variable
2367      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'
2368    """ 
2369
2370    fname = 'draw_Neighbourghood_evol'
2371
2372    if values == 'h':
2373        print fname + '_____________________________________________________________'
2374        print draw_Neighbourghood_evol.__doc__
2375        quit()
2376
2377    expectargs = '[gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:' +                 \
2378      '[timetits]:[tkinds]:[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]'
2379 
2380    drw.check_arguments(fname,values,expectargs,':')
2381
2382    gvarname = values.split(':')[0]
2383    dimsval = values.split(':')[1].split(',')
2384    neigdims = values.split(':')[2].split(',')
2385    Nneig = int(values.split(':')[3])
2386    Ncol0 = values.split(':')[4]
2387    timetits = values.split(':')[5].split(',')
2388    timekinds = values.split(':')[6].split('|')
2389    timefmts = values.split(':')[7].split(',')
2390    gtitle = values.split(':')[8].replace('|',' ')
2391    shadxtrms = values.split(':')[9].split(',')
2392    cbar = values.split(':')[10]
2393    gkind = values.split(':')[11]
2394    ofile = values.split(':')[12]
2395
2396    if Ncol0 != 'auto': 
2397        Ncol = int(Ncol0)
2398    else:
2399        Ncol = Ncol0
2400
2401    timetits[0] = timetits[0].replace('|',' ')
2402    timetits[1] = timetits[1].replace('|',' ')
2403
2404    if np.mod(Nneig,2) == 0:
2405        print errormsg
2406        print '  ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
2407        quit(-1)
2408
2409    Nneig2 = int(Nneig/2)
2410
2411# Values to slice the variable
2412    dimvslice = {}
2413    dimvvalues = {}
2414    for dimvs in dimsval:
2415        dimn = dimvs.split('|')[0]
2416        dimv = int(dimvs.split('|')[1])
2417        dimnv = dimvs.split('|')[2]
2418
2419        dimvvalues[dimn] = dimnv
2420        dimvslice[dimn] = dimv
2421
2422    ncobj = NetCDFFile(filen, 'r')
2423
2424    varobj = ncobj.variables[variable]
2425
2426    slicevar = []
2427    newdimn = []
2428    newdimsvar = {}
2429
2430    for dimn in varobj.dimensions:
2431        if not drw.searchInlist(dimvslice.keys(), dimn):
2432            dimsize = len(ncobj.dimensions[dimn])
2433            slicevar.append(slice(0, dimsize+1))
2434            newdimn.append(dimn)
2435            newdimsvar[dimn] = dimseize
2436
2437        for dimslicen in dimvslice.keys():
2438            if dimn == dimslicen:
2439                if dimvslice[dimn] != -1:
2440                    if drw.searchInlist(neigdims, dimn):
2441                        slicevar.append(slice(dimvslice[dimn]-Nneig2,                \
2442                          dimvslice[dimn]+Nneig2+1))
2443                        newdimn.append(dimn)
2444                        newdimsvar[dimn] = Nneig
2445                        break
2446                    else:
2447                        slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
2448                        break
2449                else:
2450                    dimsize = len(ncobj.dimensions[dimn])
2451                    slicevar.append(slice(0, dimsize+1))
2452                    newdimn.append(dimn)
2453                    newdimsvar[dimn] = dimsize
2454                    break
2455 
2456    varv = varobj[tuple(slicevar)]
2457
2458    if len(newdimn) != 3:
2459        print errormsg
2460        print '  ' + fname + ': sliced variable with shape=', varv.shape,            \
2461          ' must have three dimensions',len(varv.shape),'given !!'
2462        quit(-1)
2463
2464    newdims = []
2465    for nwdims in newdimn:
2466        newdims.append(newdimsvar[nwdims])
2467
2468# The dimension which is not in the neighbourhood dimensions must be time!
2469    for dim1 in newdimn:
2470        if not drw.searchInlist(neigdims, dim1):
2471            dimt = newdimsvar[dim1]
2472            dimtime = dim1
2473
2474    if Ncol == 'auto':
2475        dimtsqx = int(np.sqrt(dimt)) + 1
2476        dimtsqy = int(np.sqrt(dimt)) + 1
2477    else:
2478        dimtsqx = int(Ncol)
2479        dimtsqy = dimt/dimtsqx + 1
2480
2481    neighbourghood = np.ones((dimtsqy*Nneig,dimtsqx*Nneig), dtype=np.float)*fillValue
2482
2483    for it in range(dimt):
2484        ity = int(it/dimtsqx)
2485        itx = it-ity*dimtsqx
2486
2487        itty = (dimtsqy - ity - 1)*Nneig + Nneig2
2488        ittx = itx*Nneig + Nneig2
2489
2490        neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]=         \
2491          varv[it,::-1,:]
2492
2493    variablevals = drw.variables_values(variable)
2494    if drw.searchInlist(varobj.ncattrs(), 'units'):
2495        vunits = varobj.units
2496    else:
2497        vunits = variablevals[5]
2498
2499# Time values at the X/Y axes
2500    if ncobj.variables[dimvvalues[dimtime]].dtype == '|S1':
2501        print '    ' + fname + ': WRF time variable!'
2502        refdate = '19491201000000'
2503        tunitsval = 'hours'
2504        dimtvalues = np.zeros((dimt), dtype=np.float)
2505        tvals = ncobj.variables[dimvvalues[dimtime]]
2506        yrref=refdate[0:4]
2507        monref=refdate[4:6]
2508        dayref=refdate[6:8]
2509        horref=refdate[8:10]
2510        minref=refdate[10:12]
2511        secref=refdate[12:14]
2512
2513        refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' +        \
2514          minref + ':' + secref
2515        tunits = tunitsval + ' since ' + refdateS
2516        for it in range(dimt):
2517            wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
2518            dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
2519    else:
2520        dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
2521        tunits = ncobj.variables[newdimsvar[dimtime]].units
2522
2523    dimxv = dimtvalues[0:dimtsqx]
2524    dimyv = dimtvalues[0:dimt:dimtsqx]
2525
2526    dimn = ['time','time']
2527
2528    if ofile == 'True':
2529        ofilen = 'Neighbourghood_evol.nc'
2530        newnc = NetCDFFile(ofilen, 'w')
2531# Dimensions
2532        newdim = newnc.createDimension('time',None)
2533        newdim = newnc.createDimension('y',dimtsqy*Nneig)
2534        newdim = newnc.createDimension('x',dimtsqx*Nneig)
2535# Dimension values
2536        newvar = newnc.createVariable('time','f8',('time'))
2537        newvar[:] = np.arange(dimt)
2538        newattr = drw.basicvardef(newvar, 'time','time',tunits)
2539# Neighbourhghood variable
2540        newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'),        \
2541          fill_value=fillValue)
2542        newvar[:] = neighbourghood
2543
2544        newnc.sync()
2545        newnc.close()
2546        print fname + ": Successfull generation of file '" + ofilen + "' !!"
2547
2548# Time ticks
2549    timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
2550    timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
2551
2552    timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
2553    timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
2554
2555    for i in range(2):
2556        if shadxtrms[i][0:1] != 'S':
2557            shadxtrms[i] = np.float(shadxtrms[i])
2558
2559    drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits,   \
2560      timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
2561
2562def draw_points(filen, values):
2563    """ Function to plot a series of points
2564      [values]= [ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:[pointlabels]:
2565        [legvals]:[figureko]:[figuren]
2566        [ptasciifile]:[file],[comchar],[collon],[collat],[lab]
2567          [file]: column ASCII file with the location of the points
2568          [comchar]: '|' list of characters for commentaries
2569          [collon]: number of column with the longitude of the points
2570          [collat]: number of column with the latitude of the points
2571          [collab]: number of column with the labels of the points ('None', and will get
2572            the values from the [pointlabels] variable
2573        [gtit]: title of the figure ('|' for spaces)
2574        [mapvalues]: drawing coastaline ([proj],[res]) or None
2575          [proj]: projection
2576             * 'cyl', cilindric
2577             * 'lcc', lambert conformal
2578          [res]: resolution:
2579             * 'c', crude
2580             * 'l', low
2581             * 'i', intermediate
2582             * 'h', high
2583             * 'f', full
2584        [kindfigure]: kind of figure
2585          'legend': only points in the map with the legend with the names
2586          'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2587        [pointcolor]: color for the points ('auto' for "red")
2588        [pointlabels]: ',' of labels [only used if [collab]='None'] ('None' for no labels)
2589        [legvals]=[locleg]|[fontsize]:
2590          [locleg]: location of the legend (0, autmoatic)
2591            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2592            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2593            9: 'upper center', 10: 'center'
2594          [fontsize]: font size for the legend (auto for 12)
2595        [figureko]: kind of the output file (pdf, png, ...)
2596        [figuren]: name of the figure
2597      [filen]= [ncfile],[lonvarn],[latvarn][,[varn],[dimvals],[vargn],[min],[max],[cbar],[varu]]
2598        [ncfile]: netCDF to use to geolocalize the points
2599        [lonvarn]: name of the variable with the longitudes
2600        [latvarn]: name of the variable with the latitudes
2601        [varn]: optional variable to add staff into the graph
2602        [dimval]: '@' list of [dimn]|[dimval] to get the values for [varn]
2603          [dimn]: name of the dimension
2604          [dimval]: value of the dimension (no value all range)
2605        [vargn]: name of the variable in the graph
2606        [min]: minimum value for the extra variable
2607        [max]: maximum value for the extra variable
2608        [cbar]: color bar
2609        [varu]: units of the variable
2610    """
2611    fname = 'draw_points'
2612
2613    if values == 'h':
2614        print fname + '_____________________________________________________________'
2615        print draw_points.__doc__
2616        quit()
2617
2618    expectargs = '[ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:' +     \
2619      '[pointlabels]:[locleg]:[figurek]:[figuren]'
2620 
2621    drw.check_arguments(fname,values,expectargs,':')
2622
2623    ptasciifile = values.split(':')[0]
2624    gtit = values.split(':')[1]
2625    mapvalues = values.split(':')[2]
2626    kindfigure = values.split(':')[3]
2627    pointcolor = values.split(':')[4]
2628    pointlabels = values.split(':')[5]
2629    locleg = int(values.split(':')[6])
2630    figureko = values.split(':')[7]
2631    figuren = values.split(':')[8]
2632
2633# Getting station locations
2634##
2635    filev = ptasciifile.split(',')[0]
2636    comchar = ptasciifile.split(',')[1].split('|')
2637    collon = int(ptasciifile.split(',')[2])
2638    collat = int(ptasciifile.split(',')[3])
2639    collab = ptasciifile.split(',')[4]
2640
2641    if not os.path.isfile(filev):
2642        print errormsg
2643        print '  ' + fname + ": file '" + filev + "' does not exist!!"
2644        quit(-1)
2645
2646# Getting points position and labels
2647    oascii = open(filev, 'r')
2648    xptval = []
2649    yptval = []
2650    if collab != 'None':
2651        ptlabels = []
2652        for line in oascii:
2653            if not drw.searchInlist(comchar, line[0:1]):
2654                linevals = drw.reduce_spaces(line)
2655                xptval.append(np.float(linevals[collon].replace('\n','')))
2656                yptval.append(np.float(linevals[collat].replace('\n','')))
2657                ptlabels.append(linevals[int(collab)].replace('\n',''))
2658    else:
2659        ptlabels = None
2660        for line in oascii:
2661            if  not drw.searchInlist(comchar, line[0:1]):
2662                linevals = drw.reduce_spaces(line)
2663                xptval.append(np.float(linevals[collon].replace('\n','')))
2664                yptval.append(np.float(linevals[collat].replace('\n','')))
2665
2666    oascii.close()
2667
2668    if pointlabels != 'None' and collab == 'None':
2669        ptlabels = pointlabels.split(',')
2670
2671# Getting localization of the points
2672##
2673    filev = filen.split(',')
2674    Nvals = len(filev)
2675
2676    ncfile = filev[0]
2677    lonvarn = filev[1]
2678    latvarn = filev[2]
2679    varn = None
2680    varextrav = None
2681    if Nvals == 10:
2682        varn = filev[3]
2683        dimvals = filev[4]
2684        varextrav = [filev[5], np.float(filev[6]), np.float(filev[7]), filev[8],     \
2685          filev[9]]
2686   
2687    if not os.path.isfile(ncfile):
2688        print errormsg
2689        print '  ' + fname + ": file '" + ncfile + "' does not exist!!"
2690        quit(-1)
2691
2692    onc = NetCDFFile(ncfile, 'r')
2693   
2694    lonv, latv = drw.lonlat2D(onc.variables[lonvarn], onc.variables[latvarn])
2695
2696    if varn is not None:
2697        objextra = onc.variables[varn]
2698        vard = objextra.dimensions
2699        dd = {}
2700        for dn in dimvals.split('@'):
2701            ddn = dn.split('|')[0]
2702            ddv = dn.split('|')[1]
2703            dd[ddn] = ddv
2704        slicevar = []
2705        for dv in vard:
2706            found= False
2707            for dn in dd.keys():
2708                if dn == dv:
2709                    slicevar.append(int(dd[dn]))
2710                    found = True
2711                    break
2712            if not found:
2713                slicevar.append(slice(0,len(onc.dimensions[dv])))
2714
2715        varextra = np.squeeze(objextra[tuple(slicevar)])
2716
2717    if mapvalues == 'None':
2718        mapV = None
2719    else:
2720        mapV = mapvalues
2721
2722    drw.plot_points(xptval, yptval, lonv, latv, varextra, varextrav, gtit, mapV,     \
2723      kindfigure, pointcolor, ptlabels, locleg, figureko, figuren)
2724
2725    onc.close()
2726
2727    return
2728
2729def draw_points_lonlat(filen, values):
2730    """ Function to plot a series of lon/lat points
2731      filen= name of the file
2732      values= [lonvarname]:[latvarname]:[gkind]:[gtit]:[ptcolor]:[pttype]:[ptsize]:[labels]:[legvals]:[figureK]
2733        [lonvarname]: name of the variable longitude
2734        [latvarname]: name of the variable latitude
2735        [gkind]: kind of graphical output
2736        [gtit]: graphic title '!' for spaces
2737        [ptcolor]: color of the points ('auto', for "red")
2738        [pttype]: type of point
2739        [ptsize]: size of point
2740        [labels]: ',' list of labels to use
2741        [legvals]=[locleg]|[fontsize]:
2742          [locleg]: location of the legend (0, autmoatic)
2743            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2744            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2745            9: 'upper center', 10: 'center'
2746          [fontsize]: font size for the legend (auto for 12)
2747        [figureK]= kind of figure
2748          'legend': only points in the map with the legend with the names
2749          'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2750    """
2751    fname = 'draw_points_lonlat'
2752
2753    if values == 'h':
2754        print fname + '_____________________________________________________________'
2755        print draw_points_lonlat.__doc__
2756        quit()
2757
2758    expectargs = '[lonvarname]:[latvarname]:[gkind]:[gtit]:[ptcolor]:[pttype]:' +    \
2759      '[ptsize]:[labels]:[locleg]:[figureK]'
2760 
2761    drw.check_arguments(fname,values,expectargs,':')
2762
2763    lonname = values.split(':')[0]
2764    latname = values.split(':')[1]
2765    kindfigure = values.split(':')[2]
2766    gtit = values.split(':')[3].replace('!',' ')
2767    pointcolor = values.split(':')[4]
2768    pointtype = values.split(':')[5]
2769    pointsize = np.float(values.split(':')[6])
2770    labelsv = values.split(':')[7]
2771    loclegend = values.split(':')[8]
2772    figureK = values.split(':')[9]
2773
2774    fname = 'points_lonlat'
2775
2776    onc = NetCDFFile(filen, 'r')
2777    if not onc.variables.has_key(lonname):
2778        print errormsg
2779        print fname + ": file '" + filen + "' does not have longitudes '" + lonname +\
2780          "' !!"
2781        quit(-1)
2782    if not onc.variables.has_key(lonname):
2783        print errormsg
2784        print fname + ": file '" + filen + "' does not have longitudes '" + lonname +\
2785          "' !!"
2786        quit(-1)
2787   
2788    olon = onc.variables[lonname]
2789    olat = onc.variables[latname]
2790
2791    Ndimlon = len(olon.shape)
2792    if Ndimlon == 1:
2793        dx = olon.shape[0]
2794        dy = olat.shape[0]
2795        if dx == dy: 
2796            lonvals = olon[:]
2797            latvals = olat[:]
2798        else: 
2799            lonvals0 = np.zeros((dy,dx), dtype=np.float)
2800            latvals0 = np.zeros((dy,dx), dtype=np.float)
2801            for iL in range(dy):
2802                lonvals0[iL,:] = olon[:]
2803            for il in range(dx):
2804                latvals0[:,il] = olat[:]
2805            lonvals = lonvals0.flatten()
2806            latvals = latvals0.flatten()
2807
2808    elif Ndimlon == 2:
2809        lonvals = olon[:].flatten()
2810        latvals = olat[:].flatten()
2811    elif Ndimlon == 3:
2812        lonvals = olon[1,:,:].flatten()
2813        latvals = olat[1,:,:].flatten()
2814# Playing for Anna
2815#        lonvals = olon[:].flatten()
2816#        latvals = olat[:].flatten()
2817    elif Ndimlon == 4:
2818        lonvals = olon[1,0,:,:].flatten()
2819        latvals = olat[1,0,:,:].flatten()
2820    else:
2821        print errormsg
2822        print '  ' + fname + ': longitude size:',len(olon),' not ready!!'
2823        quit(-1)
2824
2825    if labelsv == 'None':
2826        labels = None
2827    else:
2828        labels = labelsv.split(',')
2829
2830    drw.plot_list_points(lonvals, latvals, lonname, latname, gtit, figureK, pointcolor, pointtype,    \
2831      pointsize, labels, loclegend, kindfigure, fname)
2832
2833    onc.close()
2834
2835    return
2836
2837def draw_timeSeries(filen, values, variables):
2838    """ Function to draw a time-series
2839      draw_timeSeries(filen, values, variable):
2840        filen= name of the file
2841        values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[legvals]:[gkind]:[colorlines]:[pointtype]:[pointfreq]
2842        [gvarname]: name of the variable to appear in the graph
2843        [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
2844        [tkind]: kind of time to appear in the graph (assumed x-axis)
2845          'Nval': according to a given number of values as 'Nval',[Nval]
2846          'exct': according to an exact time unit as 'exct',[tunit];
2847            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2848              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2849              'l': milisecond
2850        [timefmt]: format of the time labels
2851        [title]: title of the graphic ('|' for spaces)
2852        [legvals]=[locleg]|[fontsize]:
2853          [locleg]: location of the legend (0, autmoatic)
2854            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2855            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2856            9: 'upper center', 10: 'center'
2857          [fontsize]: font size for the legend (auto for 12)
2858        [gkind]: kind of graphical output
2859        [colorlines]: ',' list of colors for the lines, None for automatic, single
2860          value all the same
2861        [pointtype]: ',' list of type of points for the lines, None for automatic, single
2862          value all the same
2863        [pointfreq]: frequency of point plotting, 'all' for all time steps
2864      variables= [varname],[timename] names of variable and variable with times
2865      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')
2866    """
2867
2868    fname = 'draw_timeSeries'
2869
2870    if values == 'h':
2871        print fname + '_____________________________________________________________'
2872        print draw_timeSeries.__doc__
2873        quit()
2874
2875    expectargs = '[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:' +                 \
2876      '[locleg]:[gkind]:[colorlines]:[pointtype]:[pointfreq]'
2877 
2878    drw.check_arguments(fname,values,expectargs,':')
2879
2880    gvarname = values.split(':')[0]
2881    timetit = values.split(':')[1].replace('|',' ')
2882    tkind = values.split(':')[2]
2883    timefmt = values.split(':')[3]
2884    title = values.split(':')[4].replace('|',' ')
2885    locleg = int(values.split(':')[5])
2886    gkind = values.split(':')[6]
2887    colorlines = values.split(':')[7]
2888    pointtype = values.split(':')[8]
2889    pointfreq0 = values.split(':')[9]
2890   
2891    ncobj = NetCDFFile(filen, 'r')
2892
2893    variable = variables.split(',')[0]
2894    timevar = variables.split(',')[1]
2895
2896    if not ncobj.variables.has_key(variable):
2897        print errormsg
2898        print '  ' + fname + ": file '" +  filen + "' does not have variable '" +    \
2899          variable + "' !!"
2900        quit(-1)
2901
2902    if not ncobj.variables.has_key(timevar):
2903        print errormsg
2904        print '  ' + fname + ": file '" +  filen + "' does not have variable time '" \
2905          + timevar + "' !!"
2906        quit(-1)
2907
2908    varobj = ncobj.variables[variable]
2909    timeobj = ncobj.variables[timevar]
2910
2911    dimt = len(timeobj[:])
2912    varvals = np.zeros((2,dimt), dtype=np.float)
2913
2914    gunits = varobj.getncattr('units')
2915    tunits = timeobj.getncattr('units')
2916
2917    varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2918    varvals[1,:] = timeobj[:]
2919
2920    tseriesvals = []
2921    tseriesvals.append(varvals)
2922
2923    if colorlines == 'None': 
2924        collines = None
2925    else:
2926        collines = colorlines.split(',')
2927    if pointtype == 'None': 
2928        pttype = None
2929    else:
2930        pttype = pointtype.split(',')
2931
2932    if pointfreq0 == 'all':
2933        pointfreq = None
2934    else:
2935        pointfreq = int(pointfreq0)
2936
2937    drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits,        \
2938      'TimeSeries', gvarname, timetit, tkind, timefmt, title,                        \
2939      gvarname.replace('_','\_'), locleg, gkind, collines, pttype, pointfreq)
2940
2941    return
2942
2943#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:None:None', 'LDQCON,time')
2944
2945def draw_trajectories(trjfilens, values, observations):
2946    """ Function to draw different trajectories at the same time
2947    draw_trajectories(trjfilens, values, observations):
2948      trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2949         time intervals and reference maps (first one will be used to plot)
2950        [filen]: name of the file to use (lines with '#', not readed) as:
2951          [t-step] [x] [y]
2952        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2953        [map]: [file]#[lonname]#[latname]
2954          [file]; with the [lon],[lat] matrices
2955          [lonname],[latname]; names of the longitudes and latitudes variables
2956      values=[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]
2957        [leglabels]: ',' separated list of names for the legend
2958        [lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2959        [title]: title of the plot ('!' for spaces)
2960        [graphk]: kind of the graphic
2961        [mapkind]: drawing coastaline ([proj],[res]) or None
2962          [proj]: projection
2963             * 'cyl', cilindric
2964             * 'lcc', lambert conformal
2965          [res]: resolution:
2966             * 'c', crude
2967             * 'l', low
2968             * 'i', intermediate
2969             * 'h', high
2970             * 'f', full
2971      obsevations= [obsfile],[obsname],[Tint],[null]
2972        [obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2973        [obsname]: name of the observations in the graph
2974        [Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2975        [null]: null value for the observed trajectory
2976    """
2977
2978    fname = 'draw_trajectories'
2979
2980    if values == 'h':
2981        print fname + '_____________________________________________________________'
2982        print draw_trajectories.__doc__
2983        quit()
2984
2985    expectargs = '[leglabels]|[lonlatlims]|[title]|[graphk]|[mapkind]'
2986 
2987    drw.check_arguments(fname,values,expectargs,'|')
2988
2989    trjfiles = trjfilens.split(',')
2990    leglabels = values.split('|')[0]
2991    lonlatlims = values.split('|')[1]
2992    title = values.split('|')[2].replace('!',' ')
2993    graphk = values.split('|')[3]
2994    mapkind = values.split('|')[4]
2995
2996    Nfiles = len(trjfiles)
2997
2998# Getting trajectotries
2999##
3000
3001    lontrjvalues = []
3002    lattrjvalues = []
3003
3004    print '  ' + fname
3005    ifn = 0
3006    for ifile in trjfiles:
3007        filen = ifile.split('@')[0]
3008        Tint = ifile.split('@')[1]
3009
3010        print '    trajectory:',filen
3011
3012        if Tint != '-1':
3013            Tbeg = Tint
3014            Tend = ifile.split('@')[2]
3015            mapv = ifile.split('@')[3]
3016        else:
3017            mapv = ifile.split('@')[2]
3018
3019        if not os.path.isfile(filen):
3020            print errormsg
3021            print '  ' + fname + ": trajectory file '" + filen + "' does not exist !!"
3022            quit(-1)
3023
3024# Charging longitude and latitude values
3025##
3026        lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
3027          mapv.split('#')[2])
3028
3029        if ifn == 0: mapref = mapv
3030        ifn = ifn + 1
3031
3032        objfile = open(filen, 'r')
3033        trjtimev = []
3034        trjxv = []
3035        trjyv = []
3036
3037        for line in objfile:
3038            if line[0:1] != '#':
3039                trjtimev.append(int(line.split(' ')[0]))
3040                trjxv.append(int(line.split(' ')[1]))
3041                trjyv.append(int(line.split(' ')[2]))
3042
3043        objfile.close()
3044
3045        if Tint != '-1':
3046            lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
3047            lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
3048        else:
3049            lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
3050            lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
3051
3052# lonlatlimits
3053##
3054
3055    if lonlatlims == 'None':
3056        lonlatlimsv = None
3057    else:
3058        lonlatlimsv = np.zeros((4), dtype=np.float)
3059        lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
3060        lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
3061        lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
3062        lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
3063
3064# lon/lat objects
3065##
3066    objnc = NetCDFFile(mapref.split('#')[0])
3067    lonobj = objnc.variables[mapref.split('#')[1]]
3068    latobj = objnc.variables[mapref.split('#')[2]]
3069
3070# map
3071##
3072    if mapkind == 'None':
3073        mapkindv = None
3074    else:
3075        mapkindv = mapkind
3076
3077    if observations is None:
3078        obsname = None
3079    else:
3080        obsfile = observations.split(',')[0]
3081        obsname = observations.split(',')[1]
3082        Tint = observations.split(',')[2]
3083        null = np.float(observations.split(',')[3])
3084        print '    observational trajectory:',obsfile
3085
3086        if not os.path.isfile(obsfile):
3087            print errormsg
3088            print '  ' + fname + ": observations trajectory file '" + obsfile +      \
3089              "' does not exist !!"
3090            quit(-1)
3091
3092        objfile = open(obsfile, 'r')
3093        obstrjtimev = []
3094        obstrjxv = []
3095        obstrjyv = []
3096
3097        for line in objfile:
3098            if line[0:1] != '#':
3099                lon = np.float(line.split(' ')[2])
3100                lat = np.float(line.split(' ')[1])
3101                if not lon == null and not lat == null:
3102                    obstrjtimev.append(int(line.split(' ')[0]))
3103                    obstrjxv.append(lon)
3104                    obstrjyv.append(lat)
3105                else:
3106                    obstrjtimev.append(int(line.split(' ')[0]))
3107                    obstrjxv.append(None)
3108                    obstrjyv.append(None)
3109
3110        objfile.close()
3111
3112        if Tint != '-1':
3113            Tint = int(observations.split(',')[2].split('@')[0])
3114            Tend = int(observations.split(',')[2].split('@')[1])
3115            lontrjvalues.append(obstrjxv[Tint:Tend+1])
3116            lattrjvalues.append(obstrjyv[Tint:Tend+1])
3117        else:
3118            lontrjvalues.append(obstrjxv[:])
3119            lattrjvalues.append(obstrjyv[:])
3120
3121    drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','),          \
3122      lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
3123
3124    objnc.close()
3125
3126    return
3127
3128def draw_vals_trajectories(ncfile, values, variable):
3129    """ Function to draw values from the outputs from 'compute_tevolboxtraj'
3130      draw_vals_trajectories(ncfile, values, variable)
3131      ncfile= [ncfile] ',' list of files to use
3132      values= [statisticskind]:[Tint]:[labels]@[legvals]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
3133        [statisticskind]=[statistics][kind]
3134          [statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
3135            'mean2', 'stdev'
3136          [kind]: 'box', 'circle' statistics taking the values from a box or a circle
3137            'trj': value following the trajectory
3138        [Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
3139        [labels]: ',' separated list of labels for the legend
3140        [legvals]=[locleg]|[fontsize]:
3141          [locleg]: location of the legend (0, autmoatic)
3142            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
3143            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
3144            9: 'upper center', 10: 'center'
3145          [fontsize]: font size for the legend (auto for 12)
3146        [gvarname]: name of the variable to appear in the graph
3147        [timetit]: title of the time axis (assumed x-axis, '|' for spaces)
3148        [tkind]: kind of time to appear in the graph (assumed x-axis)
3149          'Nval': according to a given number of values as 'Nval',[Nval]
3150          'exct': according to an exact time unit as 'exct',[tunit];
3151            tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
3152              'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
3153                'l': milisecond
3154        [timefmt]: format of the time labels
3155        [title]: title of the graphic ('|' for spaces)
3156        [gkind]: kind of graphical output
3157      variable= variable to use
3158    """
3159    fname = 'draw_vals_trajectories'
3160
3161    if values == 'h':
3162        print fname + '_____________________________________________________________'
3163        print draw_vals_trajectories.__doc__
3164        quit()
3165
3166    sims = ncfile.split(',')
3167
3168    if len(values.split(':')) != 9:
3169        print errormsg
3170        print '  ' + fname  + ': wrong number of values!', len(values.split(':')),   \
3171          'given 9 needed!!'
3172        print '    ',values.split(':')
3173        quit(-1)
3174
3175    statisticskind = values.split(':')[0]
3176    Tint = values.split(':')[1]
3177    labels = values.split(':')[2]
3178    gvarname = values.split(':')[3]
3179    timetit = values.split(':')[4].replace('|',' ')
3180    tkind = values.split(':')[5]
3181    timefmt = values.split(':')[6]
3182    title = values.split(':')[7].replace('|',' ')
3183    gkind = values.split(':')[8]
3184
3185    leglabels = labels.split('@')[0].split(',')
3186    locleg = int(labels.split('@')[1])
3187
3188    Nsims = len(sims)
3189
3190    if Tint != '-1':
3191        tini = np.float(Tint.split('@')[0])
3192        tend = np.float(Tint.split('@')[1])
3193    else:
3194        tini = -1.
3195        tend = -1.
3196
3197    vartimetrjv = []
3198
3199    print '  ' + fname
3200    for trjfile in sims:
3201        print '    ' + trjfile + ' ...'
3202        if not os.path.isfile(trjfile):
3203            print errormsg
3204            print '  ' + fname + ": trajectory file: '" + trjfile +                  \
3205              "' does not exist !!"
3206            quit(-1)
3207
3208        trjobj = NetCDFFile(trjfile, 'r')
3209        otim = trjobj.variables['time']
3210        if not trjobj.variables.has_key(statisticskind + '_' + variable):
3211            print errormsg
3212            print '  ' + fname + ": file '" + trjfile + "' does not have variable '"+\
3213              statisticskind + '_' + variable + "' !!"
3214            quit(-1)
3215        ovar = trjobj.variables[statisticskind + '_' + variable]
3216        dimt = otim.shape[0]
3217
3218        if trjfile == sims[0]:
3219            gunits = ovar.getncattr('units')
3220            lname = ovar.getncattr('long_name')
3221            tunits = otim.getncattr('units')
3222
3223        if tini != -1:
3224            tiniid = -1
3225            tendid = -1       
3226            for itv in range(dimt):
3227                if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
3228                if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
3229
3230            if tiniid == -1 or tendid == -1:
3231                print errormsg
3232                print '  ' + main + ' time interval ', tini,',',tend,' not found: ',     \
3233                  tendid, ',', tiniid, ' !!'
3234                print '    data interval [',otim[0], otim[dimt-1],']'
3235                quit(-1)
3236            dimt = tendid - tiniid + 1
3237
3238        else:
3239            dimt = otim.shape[0]
3240
3241        valsv = np.zeros((2,dimt), dtype=np.float)
3242# Checking for time consistency
3243        if otim.getncattr('units') != tunits:
3244            print warnmsg
3245            print '  ' + fname + ': different time units in the plot!!'
3246            newtimes = gen.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
3247        else:
3248            newtimes = otim[:]
3249
3250        if tini == -1:
3251            valsv[1,:] = newtimes[:]
3252            valsv[0,:] = ovar[:]
3253        else:
3254            valsv[1,:] = newtimes[tiniid:tendid+1]
3255            valsv[0,:] = ovar[tiniid:tendid+1]
3256
3257        vartimetrjv.append(valsv)
3258        trjobj.close()
3259
3260    drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits,               \
3261      'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
3262      leglabels, locleg, gkind)
3263
3264def variable_values(values):
3265    """ Function to give back values for a given variable
3266      values= [varname] name of the variable
3267    """
3268
3269    fname = 'variable_values'
3270
3271    values = drw.variables_values(values)
3272
3273    print fname,'values:',values
3274    print fname,'variable_name:',values[0]
3275    print fname,'standard_name:',values[1]
3276    print fname,'min,max:',str(values[2]) + ',' + str(values[3])
3277    print fname,'long_name:',values[4]
3278    print fname,'units:',values[5]
3279    print fname,'shad_colors:',values[6]
3280    print fname,'all_values:',drw.numVector_String(values,',')
3281
3282    return
3283
3284def draw_ptZvals(ncfile, values, variable):
3285    """ Function to plot a given list of points and values
3286      ncfile= netCDF file to use
3287      values= [fvname]:[XYvar]:[pointype]:[pointsize]:[graphlimits]:[nxtype]:
3288        [figuretitle]:[colorbar]:[mapvalue]:[kindfig]
3289        fvname: name of the variable in the graph
3290        XYvar: [lon],[lat] variable names
3291        ptype: type of the point
3292        ptsize: size of the point
3293        graphlimits: minLON,minLAT,maxLON,maxLAT limits of the graph 'None' for the full size
3294        nxtype: minimum and maximum type
3295          'auto': values taken from the extrems of the data
3296          [min],[max]: given minimum and maximum values
3297        figtitle: title of the figure
3298        cbar: color bar
3299        mapv: map characteristics: [proj],[res]
3300          see full documentation: http://matplotlib.org/basemap/
3301          [proj]: projection
3302            * 'cyl', cilindric
3303            * 'lcc', lambert-conformal
3304          [res]: resolution:
3305            * 'c', crude
3306            * 'l', low
3307            * 'i', intermediate
3308            * 'h', high
3309            * 'f', full
3310        kfig: kind of figure
3311      variable= name of the variable to plot
3312    """
3313    fname = 'draw_ptZvals'
3314    import numpy.ma as ma
3315
3316    if values == 'h':
3317        print fname + '_____________________________________________________________'
3318        print draw_ptZvals.__doc__
3319        quit()
3320
3321    expectargs = '[fvname]:[XYvar]:[pointype]:[pointsize]:[graphlmits]:[nxtype]:' +  \
3322      '[figuretit]:[colorbar]:[mapvalue]:[kindfig]'
3323 
3324    drw.check_arguments(fname,values,expectargs,':')
3325
3326    fvname = values.split(':')[0]
3327    XYvar = values.split(':')[1]
3328    pointype = values.split(':')[2]
3329    pointsize = int(values.split(':')[3])
3330    graphlimits = values.split(':')[4]
3331    nxtype = values.split(':')[5]
3332    figuretitle = values.split(':')[6].replace('!',' ')
3333    colorbar = values.split(':')[7]
3334    mapvalue = values.split(':')[8]
3335    kindfig = values.split(':')[9]
3336
3337    onc = NetCDFFile(ncfile, 'r')
3338   
3339    if not onc.variables.has_key(variable):
3340        print errormsg
3341        print '  ' + fname + ": file '" + ncfile + "' does not have variable '" +    \
3342          variable + "' !!"
3343        quit(-1)
3344
3345# points
3346    lonvarn = XYvar.split(',')[0]
3347    latvarn = XYvar.split(',')[1]
3348
3349    if not onc.variables.has_key(lonvarn):
3350        print errormsg
3351        print '  ' + fname + ": file '" + ncfile + "' does not have longitude " +    \
3352          "variable '" + lonvarn + "' !!"
3353        quit(-1)
3354   
3355    if not onc.variables.has_key(latvarn):
3356        print errormsg
3357        print '  ' + fname + ": file '" + ncfile + "' does not have latitude " +    \
3358          "variable '" + latvarn + "' !!"
3359        quit(-1)
3360
3361    olonvar = onc.variables[lonvarn]
3362    olatvar = onc.variables[latvarn]
3363    ovarvar = onc.variables[variable]
3364
3365    Lpts = len(olonvar[:].flatten())
3366
3367    pointvalues = ma.masked_array(np.zeros((Lpts,3), dtype=np.float))
3368    pointvalues[:,0] = olonvar[:].flatten()
3369    pointvalues[:,1] = olatvar[:].flatten()
3370    pointvalues[:,2] = ovarvar[:].flatten()
3371
3372    varattrs = ovarvar.ncattrs()
3373    if drw.searchInlist(varattrs, 'units'):
3374        fvunits = ovarvar.getncattr('units')
3375    else:
3376        fvunits = drw.variables_values(variable)[5]
3377
3378# map value
3379    if mapvalue == 'None': mapvalue = None
3380
3381# Graph limits
3382    if graphlimits != 'None':
3383        graphlts = np.zeros((4), dtype=np.float)
3384        for il in range(4): graphlts[il] = np.float(graphlimits.split(',')[il])
3385        pointvalues[:,0] = ma.masked_outside(pointvalues[:,0], graphlts[0],          \
3386          graphlts[2])
3387        pointvalues[:,1] = ma.masked_outside(pointvalues[:,1], graphlts[3],          \
3388          graphlts[2])
3389
3390#        for ip in range(Lpts): 
3391#            if pointvalues[ip,0] < graphlts[0] or pointvalues[ip,0] > graphlts[2]    \
3392#              or pointvalues[ip,1] < graphlts[1] or pointvalues[ip,1] > graphlts[3]:
3393#                print ip,pointvalues[ip,0:2], graphlts
3394#                pointvalues[ip,2] = None
3395    else:
3396        graphlts = None
3397
3398    drw.plot_ptZvals(fvname,fvunits,pointvalues,pointype,pointsize,graphlts, nxtype, \
3399      figuretitle,colorbar,mapvalue,kindfig)
3400
3401    return
3402
3403#draw_ptZvals('OBSnetcdf.nc', 'pracc:lon,lat:o:80:2,42,7,47,:values!of!values:Blues:cyl,l:pdf', 'pr')
3404
3405def draw_vectors(ncfile, values, varns):
3406    """ Function to plot wind vectors
3407      values= [dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:[mapvalues]:
3408        [gtit]:[kindfig]:[figuren]
3409        'X/Y/Z/T'|[dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
3410          [dimname]: name of the dimension in the file
3411          [vardimname]: name of the variable with the values for the dimension in the file
3412          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3413          No value takes all the range of the dimension
3414        [vecvals]= [frequency],[color],[length]
3415          [frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
3416            'auto', computed automatically to have 20 vectors along each axis)
3417          [color]: color of the vectors
3418            'singlecol'@[colorn]: all the vectors same color ('auto': for 'red')
3419            'wind'@[colorbar]: color of the vectors according to wind speed sqrt(u^2+v^2) and given [colorbar]
3420              all vectors the same length
3421            '3rdvar'@[colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
3422              all vectors the same length
3423          [length]: length of the wind vectors ('auto', for 9)
3424        [windlabs]= [windname],[windunits]
3425          [windname]: name of the wind variable in the graph
3426          [windunits]: units of the wind variable in the graph ('None', for the value in the file)
3427        [mapvalues]= map characteristics: [proj],[res]
3428          see full documentation: http://matplotlib.org/basemap/
3429          [proj]: projection
3430            * 'cyl', cilindric
3431            * 'lcc', lambert conformal
3432          [res]: resolution:
3433            * 'c', crude
3434            * 'l', low
3435            * 'i', intermediate
3436            * 'h', high
3437            * 'f', full
3438        gtit= title of the graph ('|', for spaces)
3439        kindfig= kind of figure
3440        figuren= name of the figure
3441      ncfile= file to use
3442      varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
3443    """
3444    fname = 'draw_vectors'
3445
3446    if values == 'h':
3447        print fname + '_____________________________________________________________'
3448        print draw_vectors.__doc__
3449        quit()
3450
3451    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:' +  \
3452      '[mapvalues]:[gtit]:[kindfig]:[figuren]'
3453 
3454    drw.check_arguments(fname,values,expectargs,':')
3455
3456    dimvals = values.split(':')[0]
3457    vecvals = values.split(':')[1]
3458    windlabels = values.split(':')[2]
3459    mapvalues = values.split(':')[3]
3460    gtit = values.split(':')[4]
3461    kindfig = values.split(':')[5]
3462    figuren = values.split(':')[6]
3463
3464    of = NetCDFFile(ncfile,'r')
3465
3466    dims = {}
3467    for dimv in dimvals.split(','):
3468        dns = dimv.split('|')
3469        dims[dns[0]] = [dns[1], dns[2], dns[3]]
3470
3471    varNs = []
3472    for dn in dims.keys():
3473        if dn == 'X':
3474            varNs.append(dims[dn][1])
3475            dimx = len(of.dimensions[dims[dn][0]])
3476        elif dn == 'Y':
3477            varNs.append(dims[dn][1])
3478            dimy = len(of.dimensions[dims[dn][0]])
3479
3480    ivar = 0
3481    for wvar in varns.split(','):
3482        if not drw.searchInlist(of.variables.keys(), wvar):
3483            print errormsg
3484            print '  ' + fname + ": file does not have variable '" + wvar + "' !!"
3485            quit(-1)
3486        if ivar == 0:
3487            varNs.append(wvar)
3488        else:
3489            varNs.append(wvar)
3490
3491    ivar = 0
3492    for varN in varNs:
3493        varslice = []
3494
3495        ovarN = of.variables[varN]
3496        vard = ovarN.dimensions
3497        for vdn in vard:
3498            found = False
3499            for dd in dims.keys():
3500                if dims[dd][0] == vdn:
3501                    if dims[dd][2].find('@') != -1:
3502                        rvals = dims[dd][2].split('@')
3503                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
3504                    elif dims[dd][2] == '-1':
3505                        varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3506                    else:
3507                        varslice.append(int(dims[dd][2]))
3508
3509                    found = True
3510                    break
3511            if not found:
3512                varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3513
3514        if varN == dims['X'][1]:
3515            lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3516        elif varN == dims['Y'][1]:
3517            latvals0 = np.squeeze(ovarN[tuple(varslice)])
3518        elif ivar == 2:
3519            uwvals = np.squeeze(np.array(ovarN[tuple(varslice)]))
3520        elif ivar == 3:
3521            vwvals = np.squeeze(ovarN[tuple(varslice)])
3522
3523        ivar = ivar + 1
3524
3525#    print 'Final shapes:',lonvals0.shape,':',latvals0.shape,':',uwvals.shape,':',
3526#      vwvals.shape
3527
3528    if len(uwvals.shape) != 2 or len(vwvals.shape) != 2:
3529        print errormsg
3530        print '  ' + fname + ': wrong size of the wind fields! they must be ' +      \
3531          '2-dimensional!'
3532        print '    u-winds shape:',uwvals.shape,'dims:',of.variables[varNs[2]]
3533        print '    v-winds shape:',vwvals.shape,'dims:',of.variables[varNs[3]]
3534        print '      provide more values for their dimensions!!'
3535        quit(-1)
3536
3537    if len(lonvals0.shape) == 1:
3538        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3539    else:
3540        lonvals = lonvals0
3541        latvals = latvals0
3542
3543# Vector values
3544    if vecvals.split(',')[0] == 'None':
3545        freqv = None
3546    else:
3547        freqv = vecvals.split(',')[0] 
3548    colorvals = vecvals.split(',')[1]
3549    coln = colorvals.split('@')[0]
3550    colv = colorvals.split('@')[1]
3551    if coln == 'singlecol':
3552        colorv = colv
3553    elif coln == 'wind':
3554        colorv = np.sqrt(uwvals**2 + vwvals**2)
3555    elif coln == '3rdvar':
3556        if len(varn.split(',')) != 3:
3557            print errormsg
3558            print '  ' + fname + ": color of vectors should be according to '" +     \
3559              coln + "' but a third varibale is not provided !!"
3560            quit(-1)
3561        ocolvec = of.variables[varNs[4]]
3562        colorv = ocolvec[:]
3563        stdvn, lonvn, unitsvn = drw.var_3desc(ocolvec)
3564        colorvals = colorvals + '@' + stdvn + '@' + unitsvn
3565    else:
3566        print errormsg
3567        print '  ' + fname + ": color type '" + coln + "' not ready !!"
3568        quit(-1)
3569
3570    lengthv = vecvals.split(',')[2]
3571
3572# Vector labels
3573    windname = windlabels.split(',')[0]
3574    windunits = windlabels.split(',')[1]
3575
3576    drw.plot_vector(lonvals, latvals, uwvals, vwvals, freqv, colorvals, colorv,      \
3577      lengthv, windname, windunits, mapvalues, gtit, kindfig, figuren)
3578
3579    of.close()
3580
3581    return
3582
3583def draw_basins(ncfile, values, varns):
3584    """ Function to plot river basins with their discharge vector and basins id (from 'routing.nc')
3585      values= [lonlatbox]:[mapres]:[cbarname]:[xtrmbasin]:[mapdraw]:[veclength]:[freq]:
3586        [ifreq]:[plotcountry]:[basinidn]:[gtit]:[kindfig]:[figuren]
3587        [lonlatbox]= [lonSW],[lonNE],[latSW],[latNE] coordinates of the lon/lat box
3588        [mapres]= resolution of the mapping information
3589        [cbarname]= colorbar name for the colors
3590        [xtrmbasin]= [minbasin],[maxbasin] minimum and maximum basin numbers
3591        [mapdraw]= whether to draw the map (and project the data) or not ('True/False')
3592        [veclength]= length of the vectors of discharge at each grid cell
3593        [freq]= frequency of values allong each axis (None, all grid points;
3594          'auto', computed automatically to have 20 vectors along each axis)
3595        [plotcountry]= whether country lines should be plotted or not ('True/False')
3596        [plotbasinid]= whether id of the basins should be plotted or not ('True/False')
3597        [gtit]= title of the graph ('|', for spaces)
3598        [kindfig]= kind of figure
3599        [figuren]= name of the figure
3600      ncfile= file to use
3601    """
3602    fname = 'draw_basins'
3603
3604    if values == 'h':
3605        print fname + '_____________________________________________________________'
3606        print draw_vectors.__doc__
3607        quit()
3608
3609    expectargs = '[lonlatbox]:[mapres]:[cbarname]:[xtrmbasin]:[mapdraw]:' +          \
3610      '[veclength]:[freq]:[plotcountry]:[basinidn]:[gtit]:[kindfig]:[figuren]'
3611 
3612    drw.check_arguments(fname,values,expectargs,':')
3613
3614    varn='basins'
3615    lonname = 'nav_lon'
3616    latname = 'nav_lat'
3617    flowname = 'trip'
3618
3619    lonlims =[]
3620    latlims =[]
3621
3622    lonlims.append(np.float(values.split(':')[0].split(',')[0]))
3623    lonlims.append(np.float(values.split(':')[0].split(',')[1]))
3624    latlims.append(np.float(values.split(':')[0].split(',')[2]))
3625    latlims.append(np.float(values.split(':')[0].split(',')[3]))
3626    map_res = values.split(':')[1]
3627    cbarname = values.split(':')[2]
3628    vtit = 'basins'
3629    minbasin = np.int(values.split(':')[3].split(',')[0])
3630    maxbasin = np.int(values.split(':')[3].split(',')[1])
3631    mapdraw = gen.Str_Bool(values.split(':')[4])
3632    veclength = np.float(values.split(':')[5])
3633    freq0 = values.split(':')[6]
3634    plotcountry = gen.Str_Bool(values.split(':')[7])
3635    plotbasinid = gen.Str_Bool(values.split(':')[8])
3636    gtit = values.split(':')[9].replace('|',' ')
3637    kindfig = values.split(':')[10]
3638    figuren = values.split(':')[11]
3639
3640    if freq0 == 'None': freq = None
3641
3642    ofile = NetCDFFile(ncfile, 'r')
3643
3644    obasins = ofile.variables[varn]
3645    olon = ofile.variables[lonname]
3646    olat = ofile.variables[latname]
3647    oflow = ofile.variables[flowname]
3648
3649    lons = olon[:]
3650    lats = olat[:]
3651
3652    lon, lat = drw.lonlat2D(lons, lats)
3653
3654    nlon = lonlims[0]
3655    xlon = lonlims[1]
3656    nlat = latlims[0]
3657    xlat = latlims[1]
3658
3659    imin, imax, jmin, jmax = gen.ijlonlat(lon, lat, nlon, xlon, nlat, xlat)
3660
3661    drw.plot_basins(lon[jmin:jmax,imin:imax], lat[jmin:jmax,imin:imax],              \
3662      oflow[jmin:jmax,imin:imax], freq, cbarname+'@basin@-',                         \
3663      obasins[jmin:jmax,imin:imax], veclength, minbasin, maxbasin, 'outflow', '-',   \
3664      'cyl,'+map_res, plotcountry, plotbasinid, gtit, kindfig, figuren)
3665
3666    ofile.close()
3667
3668    return
3669
3670def draw_basinsold(ncfile, values, varns):
3671    """ Function to plot wind basins
3672      values= [dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:[mapvalues]:
3673        [gtit]:[kindfig]:[figuren]
3674        'X/Y/Z/T'|[dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
3675          [dimname]: name of the dimension in the file
3676          [vardimname]: name of the variable with the values for the dimension in the file
3677          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3678          No value takes all the range of the dimension
3679        [vecvals]= [frequency],[color],[length]
3680          [frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
3681            'auto', computed automatically to have 20 vectors along each axis)
3682          [color]: [colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
3683              all vectors the same length
3684          [length]: length of the wind vectors ('auto', for 9)
3685        [windlabs]= [windname],[windunits]
3686          [windname]: name of the wind variable in the graph
3687          [windunits]: units of the wind variable in the graph ('None', for the value in the file)
3688        [mapvalues]= map characteristics: [proj],[res]
3689          see full documentation: http://matplotlib.org/basemap/
3690          [proj]: projection
3691            * 'cyl', cilindric
3692            * 'lcc', lambert conformal
3693          [res]: resolution:
3694            * 'c', crude
3695            * 'l', low
3696            * 'i', intermediate
3697            * 'h', high
3698            * 'f', full
3699        gtit= title of the graph ('|', for spaces)
3700        kindfig= kind of figure
3701        figuren= name of the figure
3702      ncfile= file to use
3703      varns= [lon],[lat],[outflow],[basinID] ',' list of the name of the variables with the lon,lat, the outflow and the basin ID
3704    """
3705    fname = 'draw_basins'
3706
3707    if values == 'h':
3708        print fname + '_____________________________________________________________'
3709        print draw_vectors.__doc__
3710        quit()
3711
3712    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[vecvals]:[windlabs]:' +  \
3713      '[mapvalues]:[gtit]:[kindfig]:[figuren]'
3714 
3715    drw.check_arguments(fname,values,expectargs,':')
3716
3717    dimvals = values.split(':')[0]
3718    vecvals = values.split(':')[1]
3719    windlabels = values.split(':')[2]
3720    mapvalues = values.split(':')[3]
3721    gtit = values.split(':')[4]
3722    kindfig = values.split(':')[5]
3723    figuren = values.split(':')[6]
3724
3725    of = NetCDFFile(ncfile,'r')
3726
3727    dims = {}
3728    for dimv in dimvals.split(','):
3729        dns = dimv.split('|')
3730        dims[dns[0]] = [dns[1], dns[2], dns[3]]
3731
3732    varNs = []
3733    for dn in dims.keys():
3734        if dn == 'X':
3735            varNs.append(dims[dn][1])
3736            dimx = len(of.dimensions[dims[dn][0]])
3737        elif dn == 'Y':
3738            varNs.append(dims[dn][1])
3739            dimy = len(of.dimensions[dims[dn][0]])
3740
3741    ivar = 0
3742    for wvar in varns.split(','):
3743        if not drw.searchInlist(of.variables.keys(), wvar):
3744            print errormsg
3745            print '  ' + fname + ": file does not have variable '" + wvar + "' !!"
3746            quit(-1)
3747        if ivar == 0:
3748            varNs.append(wvar)
3749        else:
3750            varNs.append(wvar)
3751
3752    ivar = 0
3753    for varN in varNs:
3754        varslice = []
3755
3756        ovarN = of.variables[varN]
3757        vard = ovarN.dimensions
3758        for vdn in vard:
3759            found = False
3760            for dd in dims.keys():
3761                if dims[dd][0] == vdn:
3762                    if dims[dd][2].find('@') != -1:
3763                        rvals = dims[dd][2].split('@')
3764                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
3765                    elif dims[dd][2] == '-1':
3766                        varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3767                    else:
3768                        varslice.append(int(dims[dd][2]))
3769
3770                    found = True
3771                    break
3772            if not found:
3773                varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3774
3775        if varN == dims['X'][1]:
3776            lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3777        elif varN == dims['Y'][1]:
3778            latvals0 = np.squeeze(ovarN[tuple(varslice)])
3779
3780        ivar = ivar + 1
3781
3782    if len(lonvals0.shape) == 1:
3783        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3784    else:
3785        lonvals = lonvals0
3786        latvals = latvals0
3787
3788# Vector values
3789    if vecvals.split(',')[0] == 'None':
3790        freqv = None
3791    else:
3792        freqv = vecvals.split(',')[0] 
3793
3794    colorvals = vecvals.split(',')[1]
3795    if len(varn.split(',')) != 3:
3796        print errormsg
3797        print '  ' + fname + ": color of vectors should be according to '" +         \
3798          coln + "' but a third varibale is not provided !!"
3799        quit(-1)
3800
3801    ocolvec = of.variables[varNs[3]]
3802    colorv = ocolvec[:]
3803    stdvn, lonvn, unitsvn = drw.var_3desc(ocolvec)
3804    colorvals = colorvals + '@' + stdvn + '@' + unitsvn
3805
3806    lengthv = vecvals.split(',')[2]
3807
3808# Vector labels
3809    windname = windlabels.split(',')[0]
3810    windunits = windlabels.split(',')[1]
3811
3812# Vector angles
3813    oflow = ofile.variables[varNs[2]]
3814    angle = (oflow[:] - 1)*np.pi/4
3815    xflow = np.where(oflow[:] < 9, np.float(lengthv)*np.sin(angle), 0.)
3816    yflow = np.where(oflow[:] < 9, np.float(lengthv)*np.cos(angle), 0.)
3817
3818    drw.plot_basins(lonvals, latvals, xflow, yflow, freqv, colorvals, colorv,      \
3819      lengthv, windname, windunits, mapvalues, gtit, kindfig, figuren)
3820
3821    of.close()
3822
3823    return
3824
3825def draw_river_desc(ncfile, values, riverns):
3826    """ Function to plot rivers' description from ORCHIDEE's routing scheme
3827      values= [dimname]|[vardimname]|[value]:[basinvals]:[upstreamvals]:[mapvalues]:
3828        [gtit]:[kindfig]:[legvals]:[figuren]
3829        'X/Y'|[dimname]|[vardimname]|[value]: ',', list for each basic dimension '|' separated of:
3830          [dimname]: name of the dimension in the file
3831          [vardimname]: name of the variable with the values for the dimension in the file
3832          [value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3833          No value takes all the range of the dimension
3834        [basinsvals]= [colorline]
3835          [basincolor]: ',' list of colors of the line to use to mark the basins contours (single value also possible)
3836        [upstreamvals]= [upstreamvarn],[colorbar]
3837          [upstreamcolor]: colorbar to use to plot the basins upstream values
3838        [mapvalues]= map characteristics: [proj],[res]
3839          see full documentation: http://matplotlib.org/basemap/
3840          [proj]: projection
3841            * 'cyl', cilindric
3842            * 'lcc', lambert conformal
3843          [res]: resolution:
3844            * 'c', crude
3845            * 'l', low
3846            * 'i', intermediate
3847            * 'h', high
3848            * 'f', full
3849        gtit= title of the graph ('|', for spaces)
3850        kindfig= kind of figure
3851        [legvals]=[locleg]|[fontsize]:
3852          [locleg]: location of the legend (0, autmoatic)
3853            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
3854            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
3855            9: 'upper center', 10: 'center'
3856          [fontsize]: font size for the legend (auto for 12)
3857        figuren= name of the figure
3858      ncfile= file to use
3859      riverns= ',' list of the name of the rivers to plot
3860    """
3861    import numpy.ma as ma
3862    fname = 'draw_river_desc'
3863
3864    if values == 'h':
3865        print fname + '_____________________________________________________________'
3866        print draw_river_desc.__doc__
3867        quit()
3868
3869    expectargs = '[X/Y/Z/T]|[dimname]|[vardimname]|[value]:[basinvals]:' +           \
3870      '[upstreamvals]:[mapvalues]:[gtit]:[kindfig]:[legloc]:[figuren]'
3871 
3872    drw.check_arguments(fname,values,expectargs,':')
3873
3874    dimvals = values.split(':')[0]
3875    basinvals = values.split(':')[1]
3876    upstreamvals = values.split(':')[2]
3877    mapvals = values.split(':')[3]
3878    gtit = values.split(':')[4]
3879    kindfig = values.split(':')[5]
3880    legloc = int(values.split(':')[6])
3881    figuren = values.split(':')[7]
3882
3883    basincol = basinvals
3884    if basincol.find(',') != 1:
3885        basincolor = basincol.split(',')
3886    else:
3887        basincolor = [basincol]
3888
3889    upstreamcolor = upstreamvals
3890
3891    of = NetCDFFile(ncfile,'r')
3892
3893    dims = {}
3894    for dimv in dimvals.split(','):
3895        dns = dimv.split('|')
3896        dims[dns[0]] = [dns[1], dns[2], dns[3]]
3897
3898    varNs = []
3899    for dn in dims.keys():
3900        if dn == 'X':
3901            varNs.append(dims[dn][1])
3902            dimx = len(of.dimensions[dims[dn][0]])
3903        elif dn == 'Y':
3904            varNs.append(dims[dn][1])
3905            dimy = len(of.dimensions[dims[dn][0]])
3906
3907    if riverns.find(',') != -1:
3908        riverns = riverns.split(',')
3909    else:
3910        riverns = [riverns]
3911
3912    rivers = []
3913    riversubbasins = {}
3914    riversupstream = {}
3915    riversoutflow = {}
3916    for rivern in riverns:
3917        print rivern
3918
3919# subBasins
3920        basinvar = rivern + '_coding'
3921        if not drw.searchInlist(of.variables.keys(), basinvar):
3922            print errormsg
3923            print '  ' + fname + ": file does not have variable '" + basinvar + "' !!"
3924            quit(-1)
3925        rivers.append(rivern)
3926        obasin = of.variables[basinvar]
3927        riversubbasins[rivern] = obasin[:]
3928        if rivern == riverns[0]:
3929            finalmask = obasin[:].mask
3930        else:
3931            finalmask = finalmask * obasin[:].mask
3932
3933# upstream
3934        upstreamvar = rivern + '_upstream'
3935        if not drw.searchInlist(of.variables.keys(), upstreamvar):
3936            print errormsg
3937            print '  ' + fname + ": file does not have variable '" + upstreamvar + "' !!"
3938            quit(-1)
3939        oupstream = of.variables[upstreamvar]
3940        riversupstream[rivern] = oupstream[:]
3941        if rivern == riverns[0]:
3942            uunits = oupstream.getncattr('units')
3943
3944# River metadata
3945        fracvar = rivern + '_frac'
3946        if not drw.searchInlist(of.variables.keys(), fracvar):
3947            print errormsg
3948            print '  ' + fname + ": file does not have variable '" + fracvar + "' !!"
3949            quit(-1)
3950        ofrac = of.variables[fracvar]
3951        riversoutflow[rivern] = [ofrac.getncattr('Longitude_of_outflow_point'),      \
3952          ofrac.getncattr('Latitude_of_outflow_point')]
3953
3954    ivar = 0
3955    for varN in varNs:
3956        varslice = []
3957
3958        ovarN = of.variables[varN]
3959        vard = ovarN.dimensions
3960        for vdn in vard:
3961            found = False
3962            for dd in dims.keys():
3963                if dims[dd][0] == vdn:
3964                    if dims[dd][2].find('@') != -1:
3965                        rvals = dims[dd][2].split('@')
3966                        varslice.append(slice(int(rvals[0]), int(rvals[1])))
3967                    elif dims[dd][2] == '-1':
3968                        varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3969                    else:
3970                        varslice.append(int(dims[dd][2]))
3971
3972                    found = True
3973                    break
3974            if not found:
3975                varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3976
3977        if varN == dims['X'][1]:
3978            lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3979        elif varN == dims['Y'][1]:
3980            latvals0 = np.squeeze(ovarN[tuple(varslice)])
3981
3982        ivar = ivar + 1
3983
3984    if len(lonvals0.shape) == 1:
3985        lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3986    else:
3987        lonvals = lonvals0
3988        latvals = latvals0
3989
3990# Masking only the lon,lat with rivers
3991    malonvals = ma.masked_array(lonvals, mask=finalmask)
3992    malatvals = ma.masked_array(latvals, mask=finalmask)
3993
3994    if mapvals == 'None':
3995        mapvalues = None
3996    else:
3997        mapvalues = mapvals
3998
3999    drw.plot_river_desc(malonvals, malatvals, rivers, riversubbasins, riversupstream, riversoutflow,  \
4000      basincolor, upstreamcolor, uunits, mapvalues, gtit, kindfig, legloc, figuren)
4001
4002    of.close()
4003
4004def draw_vertical_levels(ncfile, values, varn):
4005    """ plotting vertical levels distribution
4006    draw_topo_geogrid(ncfile, values, varn)
4007      ncfile= file to use
4008      values= [zlogs]:[plogs]:[title]:[graphic_kind]:[legvals]
4009        zlogs= zlog,dzlog
4010        zlog: to use logarithmic scale on the height axis ('true/false')
4011        dzlog: to use logarithmic scale on the difference of height between levels axis ('true/false')
4012        plogs= plog,dplog
4013        plog: to use logarithmic scale on the height axis ('true/false')
4014        dplog: to use logarithmic scale on the difference of height between levels axis ('true/false')
4015        title: title of the graph ('!' for spaces)
4016        graphic_kind: kind of figure (jpg, pdf, png)
4017        [legvals]=[locleg]|[fontsize]:
4018          [locleg]: location of the legend (0, autmoatic)
4019            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4020            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4021            9: 'upper center', 10: 'center'
4022          [fontsize]: font size for the legend (auto for 12)
4023      varn= [varnheight],[varnpres]
4024        varnheight: name of the variable with the height of the vertical levels
4025          'WRFz': for WRF z-levels (computed as (PH + PHB)/g, from a PHB(0,i,j) = 0)
4026        varnpres: name of the variable with the pressure of the vertical levels ('None', for no pressure plot)
4027          'WRFp': for WRF p-levels (computed as P + PB, from a PHB(0,i,j) = 0)
4028    """
4029    fname = 'draw_vertical_levels'
4030
4031    if values == 'h':
4032        print fname + '_____________________________________________________________'
4033        print draw_vertical_levels.__doc__
4034        quit()
4035
4036    expectargs = '[zlogs]:[plogs]:[title]:[graphic_kind]:[legloc]'
4037 
4038    drw.check_arguments(fname,values,expectargs,':')
4039
4040    zlog = values.split(':')[0].split(',')[0]
4041    dzlog = values.split(':')[0].split(',')[1]
4042    plog = values.split(':')[1].split(',')[0]
4043    dplog = values.split(':')[1].split(',')[1]
4044    title = values.split(':')[2].replace('!',' ')
4045    kindfig = values.split(':')[3]
4046    legloc = int(values.split(':')[4])
4047
4048    if varn.find(',') == -1:
4049        varnheight = varn
4050        varnpres = None
4051        pvals = None
4052        print warnmsg
4053        print '  ' + fname + ': assuming no pressure variable!!'
4054    else:
4055        varnheight = varn.split(',')[0]
4056        varnpres = varn.split(',')[1]
4057        if varnpres == 'None': 
4058            varnpres = None
4059            pvals = None
4060
4061    if not os.path.isfile(ncfile):
4062        print errormsg
4063        print '  ' + fname + ': file "' + ncfile + '" does not exist !!'
4064        quit(-1)   
4065
4066    objf = NetCDFFile(ncfile, 'r')
4067
4068    if varnheight == 'WRFz': 
4069        if not gen.searchInlist(objf.variables,'PH'):
4070            print errormsg
4071            print '  ' + fname + ": WRF file '" + ncfile + "' does not have " +      \
4072              "variable 'PH' !!"
4073            quit(-1)
4074        if not gen.searchInlist(objf.variables,'PHB'):
4075            print errormsg
4076            print '  ' + fname + ": WRF file '" + ncfile + "' does not have " +      \
4077              "variable 'PHB' !!"
4078            quit(-1)
4079
4080        objph = objf.variables['PH']
4081        objphb = objf.variables['PHB']
4082        geop = objph[:] + objphb[:]
4083
4084        ijz0 = gen.index_mat(geop[0,], 0.)
4085        zvals = geop[0, :, ijz0[0], ijz0[1]] / 9.8
4086    else:
4087        if not gen.searchInlist(objf.variables, varnheight):
4088            print errormsg
4089            print '  ' + fname + ": file '" + ncfile + "' does not have height " +   \
4090              " variable '" + varnheight + "' !!"
4091            quit(-1)
4092        objvar = objf.variables[varn]
4093        if len(objvar.shape) == 4:
4094            print warnmsg
4095            print '  ' + fname + ": assuming that height variable '" + varnheight +  \
4096              "' with shape: dt,dz,dy,dx. Tacking first time-step"
4097
4098            ijz0 = gen.index_mat(objvar[0,0,], 0.)
4099            zvals = objvar[0, :, ijz0[0], ijz0[1]]
4100        elif len(objvar.shape) == 3:
4101            print warnmsg
4102            print '  ' + fname + ": assuming that height variable '" + varnheight +  \
4103              "' with shape: dz,dy,dx"
4104
4105            ijz0 = gen.index_mat(objvar[0,], 0.)
4106            zvals = objvar[:, ijz0[0], ijz0[1]]
4107       
4108        elif len(objvar.shape) == 2:
4109            print warnmsg
4110            print '  ' + fname + ": assuming that height variable '" + varnheight +  \
4111              "' with shape: dz,dyx"
4112
4113            ijz0 = gen.index_mat(objvar[0,], 0.)
4114            zvals = objvar[:, ijz0[0]]
4115        else:
4116            zvals = objvar[:]
4117
4118# Pressure
4119    if varnpres is not None:
4120        if varnpres == 'WRFp': 
4121            if not gen.searchInlist(objf.variables,'P'):
4122                print errormsg
4123                print '  ' + fname + ": WRF file '" + ncfile + "' does not have " +      \
4124                  "variable 'P' !!"
4125                quit(-1)
4126            if not gen.searchInlist(objf.variables,'PB'):
4127                print errormsg
4128                print '  ' + fname + ": WRF file '" + ncfile + "' does not have " +      \
4129                  "variable 'PB' !!"
4130                quit(-1)
4131
4132            objph = objf.variables['P']
4133            objphb = objf.variables['PB']
4134            pres = objph[:] + objphb[:]
4135
4136            pvals = pres[0, :, ijz0[0], ijz0[1]]
4137        else:
4138            if not gen.searchInlist(objf.variables, varnpres):
4139                print errormsg
4140                print '  ' + fname + ": file '" + ncfile + "' does not have pressure " + \
4141                  " variable '" + varnpres + "' !!"
4142                quit(-1)
4143            objvar = objf.variables[varnpres]
4144            if len(objvar.shape) == 4:
4145                print warnmsg
4146                print '  ' + fname + ": assuming that pressure variable '" + varnpres +  \
4147                  "' with shape: dt,dz,dy,dx. Tacking first time-step"
4148   
4149                pvals = objvar[0, :, ijz0[0], ijz0[1]]
4150            elif len(objvar.shape) == 3:
4151                print warnmsg
4152                print '  ' + fname + ": assuming that pressure variable '" + varnpres +  \
4153                  "' with shape: dz,dy,dx"
4154   
4155                pvals = objvar[:, ijz0[0], ijz0[1]]
4156           
4157            elif len(objvar.shape) == 2:
4158                print warnmsg
4159                print '  ' + fname + ": assuming that pressure variable '" + varnpres +  \
4160                  "' with shape: dz,dyx"
4161   
4162                pvals = objvar[:, ijz0[0]]
4163            else:
4164                pvals = objvar[:]
4165
4166# Logarithmic axes
4167    if zlog == 'true':
4168        zlogv = True
4169    elif zlog == 'false':
4170        zlogv = False
4171    else:
4172        print errormsg
4173        print '  ' + fname + ": wrong value for zlog: '" + zlog + "' !!"
4174        print "    must be either: 'true' or 'false'"
4175        quit(-1)
4176
4177    if dzlog == 'true':
4178        dzlogv = True
4179    elif dzlog == 'false':
4180        dzlogv = False
4181    else:
4182        print errormsg
4183        print '  ' + fname + ": wrong value for dzlog: '" + dzlog + "' !!"
4184        print "    must be either: 'true' or 'false'"
4185        quit(-1)
4186
4187    if pvals is not None:
4188        if plog == 'true':
4189            plogv = True
4190        elif plog == 'false':
4191            plogv = False
4192        else:
4193            print errormsg
4194            print '  ' + fname + ": wrong value for plog: '" + plog + "' !!"
4195            print "    must be either: 'true' or 'false'"
4196            quit(-1)
4197        if dplog == 'true':
4198            dplogv = True
4199        elif dplog == 'false':
4200            dplogv = False
4201        else:
4202            print errormsg
4203            print '  ' + fname + ": wrong value for dplog: '" + dplog + "' !!"
4204            print "    must be either: 'true' or 'false'"
4205            quit(-1)
4206
4207    drw.plot_vertical_lev(zvals, pvals, zlogv, dzlogv, plogv, dplogv, title, kindfig, legloc)
4208
4209    objf.close()
4210
4211    return
4212
4213def draw_subbasin(ncfile, values):
4214    """ Function to plot subbasin from 'routnig.nc' ORCDHIEE
4215      ncfile= file to use produced with nc_var.py#subbasin function
4216      values= [subasiname]:[rangecolors]:[mapv]:[basinlinewidth]:[drawsubid]:[gtit]:[figkind]:[legvals]:[figurename]
4217        [subasiname]= name of the subbasin ('!' for spaces)
4218        [rcolor]= '@', list of 'r|g|b' 1-based colors (as much as first level sub-flow). 'None' for automatic
4219        [mapv]= map characteristics: [proj],[res]
4220          see full documentation: http://matplotlib.org/basemap/
4221          [proj]: projection
4222            * 'cyl', cilindric
4223            * 'lcc', lambert conformal
4224          [res]: resolution:
4225            * 'c', crude
4226            * 'l', low
4227            * 'i', intermediate
4228            * 'h', high
4229            * 'f', full
4230        [basinlinewidth]= with of the line to draw the basin
4231        [drawsubid]= wehther sub-flow ids should be plot or not
4232        [graphtit]= title of the graph ('|', for spaces)
4233        [legvals]=[locleg]|[fontsize]:
4234          [locleg]: location of the legend (0, autmoatic)
4235            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4236            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4237            9: 'upper center', 10: 'center'
4238          [fontsize]: font size for the legend (auto for 12)
4239        [figname]= name of the figure
4240    """
4241    fname = 'draw_subbasin'
4242
4243    if values == 'h':
4244        print fname + '_____________________________________________________________'
4245        print draw_subbasin.__doc__
4246        quit()
4247
4248    expectargs = '[subasiname]:[rangecolors]:[mapv]:[basinlinewidth]:[drawsubid]:' + \
4249      '[gtit]:[figkind]:[legloc]:[figurename]'
4250 
4251    drw.check_arguments(fname,values,expectargs,':')
4252
4253    subbasiname = values.split(':')[0].replace('!',' ')
4254    rangecolors = values.split(':')[1]
4255    mapv = values.split(':')[2]
4256    basinlinewidth = np.float(values.split(':')[3])
4257    drawsubid = gen.Str_Bool(values.split(':')[4])
4258    gtit = values.split(':')[5].replace('!',' ')
4259    figkind = values.split(':')[6]
4260    legloc = int(values.split(':')[7])
4261    figurename = values.split(':')[8]
4262
4263    if not os.path.isfile(ncfile):
4264        print errormsg
4265        print '  ' + fname + ': file "' + ncfile + '" does not exist !!'
4266        quit(-1)   
4267
4268    objf = NetCDFFile(ncfile, 'r')
4269
4270    searchvars = ['lon', 'lat', 'lonsubflow', 'latsubflow', 'outsubflow']
4271    for searchvar in searchvars: 
4272        if not gen.searchInlist(objf.variables,searchvar):
4273            print errormsg
4274            print '  ' + fname + ": WRF file '" + ncfile + "' does not have " +      \
4275              "variable '" + searchvar + "' !!"
4276            quit(-1)
4277   
4278# lon,lat
4279    olon = objf.variables['lon']
4280    olat = objf.variables['lat']
4281    lon = olon[:]
4282    lat = olat[:]
4283
4284# sub-flow names
4285    osubnames = objf.variables['subflow']
4286    subnames = drw.get_str_nc(osubnames, osubnames.shape[1])
4287
4288# sub-flow lat, lon
4289    latlonsub = {}
4290    outflowsub = {}
4291    osublon = objf.variables['lonsubflow']
4292    osublat = objf.variables['latsubflow']
4293    oNsubflow = objf.variables['Nsubflow']
4294    ooutsubflow = objf.variables['outsubflow']
4295    Nsubflow = oNsubflow[:]
4296    isub = 0
4297    for Ssub in subnames:
4298        sublatlon = []
4299        suboutflow = []
4300        for igrid in range(Nsubflow[isub]):
4301            sublatlon.append([osublat[isub,igrid], osublon[isub,igrid]])
4302            suboutflow.append(ooutsubflow[isub,igrid])
4303        latlonsub[Ssub] = sublatlon
4304        outflowsub[Ssub] = suboutflow
4305        isub = isub + 1
4306
4307# colors
4308    if rangecolors == 'None':
4309        rangecols = None
4310    else:
4311        cols = rangecolors.split('@')
4312        Ncols = len(cols)
4313        rangecols = []
4314        for icol in range(Ncols):
4315            cval = cols[icol].split('|')
4316            rangecols.append([np.float(cval[0]),np.float(cval[1]),np.float(cval[2])])
4317
4318    drw.plot_subbasin(subbasiname, lon, lat, subnames, latlonsub, outflowsub,        \
4319      rangecols, mapv, basinlinewidth, drawsubid, gtit, figkind, legloc, figurename)
4320
4321    objf.close()
4322
4323    return
4324
4325def draw_2lines(ncfiles, values, varnames):
4326    """ Fucntion to plot two lines in different axes (x/x2 or y/y2)
4327      values= [commonvardim]:[varangeA]:[varangeB]:[varangeaxis]:[axisvals]:[figvarns]:[colors]:
4328       [widths]:[styles]:[sizemarks]:[marks]:[graphtitle]:[labelaxis]:[legvals]:[figname]:[figkind]:[close]
4329        [commonvardim]: name of the common variable-dimension
4330        [varangeA]: ',' separated list of range (min,max) for A values ('None', automatic; 'Extrs' from values extremes)
4331        [varangeB]: ',' separated list of range (min,max) for B values ('None', automatic; 'Extrs' from values extremes)
4332        [varangeaxis]: ',' separated list of range (min,max) for common axis values ('None', automatic; 'Extrs' from
4333          values extremes)
4334        [axisvals]: which is the axis to plot the values ('x' or 'y')
4335        [figvarns]: ',' separated list of names of the variables in the  plot
4336        [colors]: ',' list with color names of the lines for the variables ('None', automatic)
4337        [widths]: ',' list with widths of the lines for the variables ('None', automatic)
4338        [styles]: ',' list with the styles of the lines ('None', automatic)
4339        [sizemarks]: ',' list with the size of the markers of the lines ('None', automatic)
4340        [marks]: ',' list with the markers of the lines ('None', automatic)
4341        [graphtitle]: title of the figure ('!' for spaces)
4342        [labelaxis]: label in the figure of the common axis ('!' for spaces)
4343        [legvals]=[locleg]|[fontsize]:
4344          [locleg]: location of the legend (0, autmoatic)
4345            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4346            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4347            9: 'upper center', 10: 'center'
4348          [fontsize]: font size for the legend (auto for 12)
4349        [figname]: name of the figure
4350        [figkind]: kind of figure
4351        [close]: Whether figure should be finished or not
4352      ncfiles= ',' separated list of files to use
4353      varnames=  ',' separated list of variables names in the files to plot
4354    """
4355    fname = 'draw_2lines'
4356
4357    if values == 'h':
4358        print fname + '_____________________________________________________________'
4359        print draw_2lines.__doc__
4360        quit()
4361
4362    expectargs = '[commonvardim]:[varangeA]:[varangeB]:' +           \
4363     '[varangeaxis]:[axisvals]:[figvarns]:[colors]:[widths]:[styles]:[sizemarks]:' + \
4364     '[marks]:[graphtitle]:[labelaxis]:[lloc]:[figname]:[figkind]:[close]'
4365 
4366    drw.check_arguments(fname,values,expectargs,':')
4367
4368    commonvardim = values.split(':')[0]
4369    varangeA0 = values.split(':')[1]
4370    varangeB0 = values.split(':')[2]
4371    varangeaxis0 = values.split(':')[3]
4372    axisvals = values.split(':')[4]
4373    figvarns = values.split(':')[5].split(',')
4374    colors = gen.str_list(values.split(':')[6],',')
4375    widths = gen.str_list_k(values.split(':')[7],',','np.float')
4376    styles = gen.str_list(values.split(':')[8],',')
4377    sizemarks = gen.str_list_k(values.split(':')[9],',','np.float')
4378    marks = gen.str_list(values.split(':')[10],',')
4379    graphtitle = values.split(':')[11].replace('!',' ')
4380    labelaxis = values.split(':')[12].replace('!',' ')
4381    legvals = values.split(':')[13]
4382    figname = values.split(':')[14]
4383    figkind = values.split(':')[15]
4384    close = gen.Str_Bool(values.split(':')[16])
4385
4386    files = ncfiles.split(',')
4387    invarns = varnames.split(',')
4388
4389    varunits = []
4390
4391    # Values line A
4392    if not os.path.isfile(files[0]):
4393        print errormsg
4394        print '  ' + fname + ": file '" + files[0] + "' does not exist !!"
4395        quit(-1)
4396
4397    oncA = NetCDFFile(files[0], 'r')
4398
4399    if not gen.searchInlist(oncA.variables.keys(), invarns[0]):
4400        print errormsg
4401        print '  ' + fname + ": A file '" + files[0] + "' does not have variable '" +\
4402          invarns[0] + "' !!"
4403        quit(-1)
4404
4405    objvA = oncA.variables[invarns[0]]
4406    varvalsA = objvA[:]
4407    varangeA = np.zeros((2),dtype=np.float)
4408
4409    if gen.searchInlist(objvA.ncattrs(), 'units'):
4410        varunits.append(drw.units_lunits(objvA.getncattr('units')))
4411    else:
4412        valsA = gen.variables_values(invarns[0])
4413        varunits.append(drw.units_lunits(valsA[5]))
4414    if varangeA0 == 'None':
4415        varangeA = [valsA[2], valsA[3]]
4416    elif varangeA0 == 'Extrs':
4417        varangeA = [np.min(varvalsA), np.max(varvalsA)]
4418    else:
4419        for iv in range(2): varangeA[iv] = np.float(varangeA0.split(',')[iv])
4420
4421    if not gen.searchInlist(oncA.variables.keys(), commonvardim):
4422        print errormsg
4423        print '  ' + fname + ": A file '" + files[0] + "' does not have common " +   \
4424          "dimvar '" + commonvardim + "' !!"
4425        quit(-1)
4426    objvd = oncA.variables[commonvardim]
4427    varvalsaxisA = objvd[:]   
4428
4429    oncA.close()
4430
4431    # Values line B
4432    if not os.path.isfile(files[1]):
4433        print errormsg
4434        print '  ' + fname + ": file '" + files[1] + "' does not exist !!"
4435        quit(-1)
4436
4437    oncB = NetCDFFile(files[1], 'r')
4438
4439    if not gen.searchInlist(oncB.variables.keys(), invarns[1]):
4440        print errormsg
4441        print '  ' + fname + ": B file '" + files[1] + "' does not have variable '" +\
4442          invarns[1] + "' !!"
4443        quit(-1)
4444
4445    objvB = oncB.variables[invarns[1]]
4446    varvalsB = objvB[:]
4447    varangeB = np.zeros((2),dtype=np.float)
4448
4449    if gen.searchInlist(objvB.ncattrs(), 'units'):
4450        varunits.append(drw.units_lunits(objvB.getncattr('units')))
4451    else:
4452        valsB = gen.variables_values(invarns[1])
4453        varunits.append(drw.units_lunits(valsB[5]))
4454    if varangeB0 == 'None':
4455        varangeB = [valsB[2], valsB[3]]
4456    elif varangeB0 == 'Extrs':
4457        varangeA = [np.min(varvalsA), np.max(varvalsA)]
4458    else:
4459        for iv in range(2): varangeB[iv] = np.float(varangeB0.split(',')[iv])
4460
4461    # Common vardim
4462    if not gen.searchInlist(oncB.variables.keys(), commonvardim):
4463        print errormsg
4464        print '  ' + fname + ": B file '" + files[1] + "' does not have common " +   \
4465          "dimvar '" + commonvardim + "' !!"
4466        quit(-1)
4467    objvd = oncB.variables[commonvardim]
4468    varvalsaxisB = objvd[:]
4469
4470    # Range of the axis
4471    varangeaxis = np.zeros((2),dtype=np.float)
4472
4473    valsVD = gen.variables_values(commonvardim)
4474    if gen.searchInlist(objvd.ncattrs(), 'units'):
4475        dimvarunits = drw.units_lunits(objvd.getncattr('units'))
4476    else:
4477        dimvarunits = drw.units_lunits(valsVD[5])
4478    if varangeaxis0 == 'None':
4479        varangeaxis = [valsVD[2], valsVD[3]]
4480    elif varangeaxis0 == 'Extrs':
4481        varangeaxis[0] = np.min([np.min(varvalsaxisA), np.min(varvalsaxisB)])
4482        varangeaxis[1] = np.max([np.max(varvalsaxisA), np.max(varvalsaxisB)])
4483    else:
4484        for iv in range(2): varangeaxis[iv] = np.float(varangeaxis0.split(',')[iv])
4485   
4486    oncB.close()
4487
4488    labelaxis = valsVD[0] + ' (' + dimvarunits + ')'
4489
4490    # Lines characteristics
4491    colvalues, linekinds, pointkinds, lwidths, psizes = drw.ColorsLinesPointsStyles( \
4492      2, colors, styles, marks, widths, sizemarks, 'None')
4493
4494    # legend
4495    lloc, lsize = drw.legend_values(legvals,'|')
4496
4497    drw.plot_2lines(varvalsA, varvalsB, varvalsaxisA, varvalsaxisB, varangeA,        \
4498      varangeB, varangeaxis, axisvals, figvarns, varunits, colvalues, lwidths,       \
4499      linekinds, psizes, pointkinds, graphtitle, labelaxis, lloc, lsize, figname,    \
4500      figkind, close)
4501
4502def draw_2lines_time(ncfiles, values, varnames):
4503    """ Function to plot two time-lines in different axes (x/x2 or y/y2)
4504      values= [timevardim]:[varangeA]:[varangeB]:[timeaxisfmt]:[timeaxis]:[figvarns]:[colors]:
4505       [widths]:[styles]:[sizemarks]:[marks]:[graphtitle]:[labelaxis]:[legvals]:[figname]:[figkind]:[close]
4506        [timevardim]: name of the common variable-dimension time
4507        [varangeA]: ',' separated list of range (min,max) for A values ('None', automatic; 'Extrs' from values extremes)
4508        [varangeB]: ',' separated list of range (min,max) for B values ('None', automatic; 'Extrs' from values extremes)
4509        [timeaxisfmt]=[tkind];[tfmt]: format of the ticks for the time axis:
4510           [kind]: kind of time to appear in the graph
4511             'Nval': according to a given number of values as 'Nval',[Nval]
4512             'exct': according to an exact time unit as 'exct',[tunit];
4513               tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
4514                'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
4515                'l': milisecond
4516           [tfmt]; desired format
4517        [timeaxis]: which is the time axis in the plot ('x' or 'y')
4518        [figvarns]: ',' separated list of names of the variables in the  plot
4519        [colors]: ',' list with color names of the lines for the variables ('None', automatic)
4520        [widths]: ',' list with widths of the lines for the variables ('None', automatic)
4521        [styles]: ',' list with the styles of the lines ('None', automatic)
4522        [sizemarks]: ',' list with the size of the markers of the lines ('None', automatic)
4523        [marks]: ',' list with the markers of the lines ('None', automatic)
4524        [graphtitle]: title of the figure ('!' for spaces)
4525        [labelaxis]: label in the figure of the common axis ('!' for spaces)
4526        [legvals]=[locleg]|[fontsize]:
4527          [locleg]: location of the legend (0, autmoatic)
4528            1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4529            5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4530            9: 'upper center', 10: 'center'
4531          [fontsize]: font size for the legend (auto for 12)
4532        [figname]: name of the figure
4533        [figkind]: kind of figure
4534        [close]: Whether figure should be finished or not
4535      ncfiles= ',' separated list of files to use
4536      varnames=  ',' separated list of variables names in the files to plot
4537    """
4538    fname = 'draw_2lines_time'
4539
4540    if values == 'h':
4541        print fname + '_____________________________________________________________'
4542        print draw_2lines_time.__doc__
4543        quit()
4544
4545    expectargs = '[timevardim]:[varangeA]:[varangeB]:[timeaxisfmt]:[timeaxis]:' +    \
4546     '[figvarns]:[colors]:[widths]:[styles]:[sizemarks]:[marks]:[graphtitle]:' +     \
4547     '[labelaxis]:[lloc]:[figname]:[figkind]:[close]'
4548 
4549    drw.check_arguments(fname,values,expectargs,':')
4550
4551    timevardim = values.split(':')[0]
4552    varangeA0 = values.split(':')[1]
4553    varangeB0 = values.split(':')[2]
4554    timeaxisfmt = values.split(':')[3]
4555    timeaxis = values.split(':')[4]
4556    figvarns = values.split(':')[5].split(',')
4557    colors = gen.str_list(values.split(':')[6],',')
4558    widths = gen.str_list_k(values.split(':')[7],',','np.float')
4559    styles = gen.str_list(values.split(':')[8],',')
4560    sizemarks = gen.str_list_k(values.split(':')[9],',','np.float')
4561    marks = gen.str_list(values.split(':')[10],',')
4562    graphtitle = values.split(':')[11].replace('!',' ')
4563    labelaxis = values.split(':')[12].replace('!',' ')
4564    legvals = values.split(':')[13]
4565    figname = values.split(':')[14]
4566    figkind = values.split(':')[15]
4567    close = gen.Str_Bool(values.split(':')[16])
4568
4569    files = ncfiles.split(',')
4570    invarns = varnames.split(',')
4571
4572    varunits = []
4573
4574    # Values line A
4575    if not os.path.isfile(files[0]):
4576        print errormsg
4577        print '  ' + fname + ": file '" + files[0] + "' does not exist !!"
4578        quit(-1)
4579
4580    oncA = NetCDFFile(files[0], 'r')
4581
4582    if not gen.searchInlist(oncA.variables.keys(), invarns[0]):
4583        print errormsg
4584        print '  ' + fname + ": A file '" + files[0] + "' does not have variable '" +\
4585          invarns[0] + "' !!"
4586        quit(-1)
4587    if not gen.searchInlist(oncA.variables.keys(), timevardim):
4588        print errormsg
4589        print '  ' + fname + ": A file '" + files[0] + "' does not have time " +     \
4590          "variable '" + timevardim + "' !!"
4591        quit(-1)
4592
4593    objvA = oncA.variables[invarns[0]]
4594    varvalsA = objvA[:]
4595    varangeA = np.zeros((2),dtype=np.float)
4596    objtA = oncA.variables[timevardim]
4597    timevalsA = objtA[:]
4598    trangeA = [np.min(timevalsA), np.max(timevalsA)]
4599    tunitsA = objtA.getncattr('units')
4600
4601    if len(varvalsA.shape) != 1:
4602        print errormsg
4603        print '  ' + fname + ": variable '" + invarns[0] + "' has wrong shape:",     \
4604          varvalsA.shape, 'it must be 1D !!'
4605        quit(-1)
4606
4607    if gen.searchInlist(objvA.ncattrs(), 'units'):
4608        varunits.append(drw.units_lunits(objvA.getncattr('units')))
4609    else:
4610        valsA = gen.variables_values(invarns[0])
4611        varunits.append(drw.units_lunits(valsA[5]))
4612    if varangeA0 == 'None':
4613        varangeA = [valsA[2], valsA[3]]
4614    elif varangeA0 == 'Extrs':
4615        varangeA = [np.min(varvalsA), np.max(varvalsA)]
4616    else:
4617        for iv in range(2): varangeA[iv] = np.float(varangeA0.split(',')[iv])
4618
4619    oncA.close()
4620
4621    # Values line B
4622    if not os.path.isfile(files[1]):
4623        print errormsg
4624        print '  ' + fname + ": file '" + files[1] + "' does not exist !!"
4625        quit(-1)
4626
4627    oncB = NetCDFFile(files[1], 'r')
4628
4629    if not gen.searchInlist(oncB.variables.keys(), invarns[1]):
4630        print errormsg
4631        print '  ' + fname + ": B file '" + files[1] + "' does not have variable '" +\
4632          invarns[1] + "' !!"
4633        quit(-1)
4634    if not gen.searchInlist(oncB.variables.keys(), timevardim):
4635        print errormsg
4636        print '  ' + fname + ": B file '" + files[1] + "' does not have time " +     \
4637          "variable '" + timevardim + "' !!"
4638        quit(-1)
4639
4640    objvB = oncB.variables[invarns[1]]
4641    varvalsB = objvB[:]
4642    varangeB = np.zeros((2),dtype=np.float)
4643    objtB = oncB.variables[timevardim]
4644    timevalsB = objtB[:]
4645    tunitsB = objtB.getncattr('units')
4646
4647    if gen.searchInlist(objvB.ncattrs(), 'units'):
4648        varunits.append(drw.units_lunits(objvB.getncattr('units')))
4649    else:
4650        valsB = gen.variables_values(invarns[1])
4651        varunits.append(drw.units_lunits(valsB[5]))
4652    if varangeB0 == 'None':
4653        varangeB = [valsB[2], valsB[3]]
4654    elif varangeB0 == 'Extrs':
4655        varangeA = [np.min(varvalsA), np.max(varvalsA)]
4656    else:
4657        for iv in range(2): varangeB[iv] = np.float(varangeB0.split(',')[iv])
4658   
4659    oncB.close()
4660
4661    # Time axis taking time units in line A as reference
4662    varvalsaxisB = gen.coincident_CFtimes(timevalsB, tunitsA, tunitsB)
4663    trangeB = [np.min(varvalsaxisB), np.max(varvalsaxisB)]
4664
4665    varangeaxis = [np.min([trangeA[0],trangeB[0]]), np.max([trangeA[1],trangeB[1]])]
4666
4667    timevals = np.arange(varangeaxis[0],varangeaxis[1])
4668    tkind = timeaxisfmt.split(';')[0]
4669    tformat = timeaxisfmt.split(';')[1]
4670    tpos, tlabels = drw.CFtimes_plot(timevals, tunitsA, tkind, tformat)
4671
4672    # Lines characteristics
4673    colvalues, linekinds, pointkinds, lwidths, psizes = drw.ColorsLinesPointsStyles( \
4674      2, colors, styles, marks, widths, sizemarks, 'None')
4675
4676    # legend
4677    lloc, lsize = drw.legend_values(legvals,'|')
4678
4679    drw.plot_2lines_time(varvalsA, varvalsB, timevalsA, varvalsaxisB, varangeA,      \
4680      varangeB, tpos, tlabels, timeaxis, figvarns, varunits, colvalues, lwidths,     \
4681      linekinds, psizes, pointkinds, graphtitle, labelaxis, lloc, lsize, figname,    \
4682      figkind, close)
4683
4684#quit()
4685
4686####### ###### ##### #### ### ## #
4687
4688ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
4689 
4690### Options
4691##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
4692string_operation="""operation to make:
4693  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]
4694  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]
4695    [ckind]:
4696      'cmap': as it gets from colorbar
4697      'fixc,[colname]': fixed color [colname], all stright lines
4698      'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed  line
4699"""
4700####### ###### ##### #### ### ## #
4701
4702# Not checking file operation
4703Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time',                    \
4704  'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_2lines', 'draw_2lines_time',  \
4705  'draw_lines',                                                                      \
4706  'draw_lines_time', 'draw_points', 'draw_topo_geogrid_boxes', 'draw_trajectories',  \
4707  'draw_vals_trajectories', 'variable_values']
4708
4709errormsg='ERROR -- error -- ERROR -- error'
4710
4711# From: http://stackoverflow.com/questions/4041238/why-use-def-main
4712def main():
4713#######    #######
4714## MAIN
4715    #######
4716
4717    parser = OptionParser()
4718    parser.add_option("-f", "--netCDF_file", dest="ncfile", 
4719     help="file to use", metavar="FILE")
4720    parser.add_option("-o", "--operation", type='choice', dest="operation", 
4721      choices=namegraphics, help="operation to make: " + ngraphics, metavar="OPER")
4722    parser.add_option("-S", "--valueS", dest="values", 
4723      help="[WHEN APPLICABLE] values to use according to the operation", 
4724      metavar="VALUES")
4725    parser.add_option("-v", "--variable", dest="varname",
4726      help="[WHEN APPLICABLE] variable to check", metavar="VAR")
4727
4728    (opts, args) = parser.parse_args()
4729
4730    varn=opts.varname
4731    oper=opts.operation
4732
4733    if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and                   \
4734      not drw.searchInlist(Notcheckingfile, oper):
4735        print errormsg
4736        print '  ' + mainn + ': File ' + opts.ncfile + ' does not exist !!'
4737        quit(-1)
4738
4739    if oper == 'create_movie':
4740        create_movie(opts.ncfile, opts.values, opts.varname)
4741    elif oper == 'draw_2D_shad':
4742        draw_2D_shad(opts.ncfile, opts.values, opts.varname)
4743    elif oper == 'draw_2D_shad_time':
4744        draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
4745    elif oper == 'draw_2D_shad_cont':
4746        draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
4747    elif oper == 'draw_2D_shad_cont_time':
4748        draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
4749    elif oper == 'draw_2D_shad_line':
4750        draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
4751    elif oper == 'draw_2D_shad_line_time':
4752        draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
4753    elif oper == 'draw_barbs':
4754        draw_barbs(opts.ncfile, opts.values, opts.varname)
4755    elif oper == 'draw_basins':
4756        draw_basins(opts.ncfile, opts.values, opts.varname)
4757    elif oper == 'draw_Neighbourghood_evol':
4758        draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
4759    elif oper == 'draw_2lines':
4760        draw_2lines(opts.ncfile, opts.values, opts.varname)
4761    elif oper == 'draw_2lines_time':
4762        draw_2lines_time(opts.ncfile, opts.values, opts.varname)
4763    elif oper == 'draw_lines':
4764        draw_lines(opts.ncfile, opts.values, opts.varname)
4765    elif oper == 'draw_lines_time':
4766        draw_lines_time(opts.ncfile, opts.values, opts.varname)
4767    elif oper == 'draw_points':
4768        draw_points(opts.ncfile, opts.values)
4769    elif oper == 'draw_points_lonlat':
4770        draw_points_lonlat(opts.ncfile, opts.values)
4771    elif oper == 'draw_ptZvals':
4772        draw_ptZvals(opts.ncfile, opts.values, opts.varname)
4773    elif oper == 'draw_river_desc':
4774        draw_river_desc(opts.ncfile, opts.values, opts.varname)
4775    elif oper == 'draw_subbasin':
4776        draw_subbasin(opts.ncfile, opts.values)
4777    elif oper == 'draw_timeSeries':
4778        draw_timeSeries(opts.ncfile, opts.values, opts.varname)
4779    elif oper == 'draw_topo_geogrid':
4780        draw_topo_geogrid(opts.ncfile, opts.values)
4781    elif oper == 'draw_topo_geogrid_boxes':
4782        draw_topo_geogrid_boxes(opts.ncfile, opts.values)
4783    elif oper == 'draw_trajectories':
4784        draw_trajectories(opts.ncfile, opts.values, opts.varname)
4785    elif oper == 'draw_vals_trajectories':
4786        draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
4787    elif oper == 'draw_vectors':
4788        draw_vectors(opts.ncfile, opts.values, opts.varname)
4789    elif oper == 'draw_vertical_levels':
4790        draw_vertical_levels(opts.ncfile, opts.values, opts.varname)
4791    elif oper == 'list_graphics':
4792    # From: http://www.diveintopython.net/power_of_introspection/all_together.html
4793        import drawing as myself
4794        object = myself
4795        for opern in namegraphics:
4796            if  opern != 'list_graphics': 
4797                print opern + '_______ ______ _____ ____ ___ __ _'
4798                print getattr(object, opern).__doc__
4799    elif oper == 'variable_values':
4800        variable_values(opts.values)
4801    else:
4802        print errormsg
4803        print '  ' + main + ": the graphic '" + oper + "' is not ready !!"
4804        print errormsg
4805        quit()
4806
4807if __name__ == '__main__':
4808    main()
4809
Note: See TracBrowser for help on using the repository browser.