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

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

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