source: trunk/UTIL/PYTHON/sparse.py @ 796

Last change on this file since 796 was 585, checked in by acolaitis, 13 years ago

sorry, small typo in sparse.py

  • Property svn:executable set to *
File size: 15.6 KB
Line 
1#!/usr/bin/env python
2
3### T.Navarro
4
5##
6#   This routine reads sparse data, i.e. not binned in a lat/lon/vert/time grid and does plots
7#
8#   There is nothing gnuplot could not be able to do for now, except command line and ease and all python plots possibilities,
9#   but the goal is to provide a framework to compare with data from simulations on a constant grid.
10###
11
12###########################################################################################
13###########################################################################################
14### What is below relate to running the file as a command line executable (very convenient)
15if __name__ == "__main__":
16   import sys
17   from optparse import OptionParser    ### to be replaced by argparse
18   from os import system,path
19   import numpy as np
20   from sparse import sparse
21
22   parser = OptionParser()
23
24   #############################
25   ### Options
26   parser.add_option('-f', '--file',    action='store',dest='file',     type="string",  default=None,  help='[NEEDED] filename, comma separated')
27   parser.add_option('--ftype',         action='store',dest='ftype',   default=None, help=' Force data file type [\'sav\',\'txt\']')
28   parser.add_option('-X', '--xvar',    action='store',dest='xvar',     type="string",  default=None,  help='[NEEDED] x variable: a name, a number, etc ...')
29   parser.add_option('-Y', '--yvar',    action='store',dest='yvar',     type="string",  default=None,  help='y variable')
30   parser.add_option('-Z', '--zvar',    action='store',dest='zvar',     type="string",  default=None,  help='z variable')
31   parser.add_option('-C', '--cvar',    action='store',dest='cvar',     type="string",  default=None,  help='c variable for color')
32   parser.add_option('--xmin',          action='store',dest='xmin',     type="float",   default=None,  help='min values for x var, comma separated')
33   parser.add_option('--xmax',          action='store',dest='xmax',     type="float",   default=None,  help='max values for x var, comma separated')
34   parser.add_option('--ymin',          action='store',dest='ymin',     type="float",   default=None,  help='min values for y var, comma separated')
35   parser.add_option('--ymax',          action='store',dest='ymax',     type="float",   default=None,  help='max values for y var, comma separated')
36   parser.add_option('--zmin',          action='store',dest='zmin',     type="float",   default=None,  help='min values for z var, comma separated')
37   parser.add_option('--zmax',          action='store',dest='zmax',     type="float",   default=None,  help='max values for z var, comma separated')
38   parser.add_option('--cmin',          action='store',dest='cmin',     type="string",  default=None,  help='min values for c var, comma separated')
39   parser.add_option('--cmax',          action='store',dest='cmax',     type="string",  default=None,  help='max values for c var, comma separated')
40   parser.add_option('-w','--with',     action='append',dest='cond',    type="string",  default=None,  help='conditions..')
41   parser.add_option('--merge',   action='store_true', dest='merge',                      default=False, help='merge datafiles in a single plot [False]')
42   parser.add_option('-c','--color',  action='store',dest='clb',       type="string",  default="def", help='change colormap (also: nobar,onebar)')
43   parser.add_option('--trans',        action='store',dest='trans',     type="float",   default=1.,    help='shaded plot transparency, 0 to 1 (=opaque) [1]')
44   parser.add_option('-p', '--proj',   action='store',dest='proj',      type="string",  default=None,  help='projection')
45   parser.add_option('--blat',         action='store',dest='blat',      type="int",     default=None,  help='reference lat (or bounding lat for stere) [computed]')
46   parser.add_option('--blon',         action='store',dest='blon',      type="int",     default=None,  help='reference lon [computed]')
47   parser.add_option('-b', '--back',   action='store',dest='back',      type="string",  default=None,  help='background image [None]')
48   parser.add_option('-m', '--min',    action='store',dest='vmin',      type="float",   default=None,  help='bounding minimum value [min]')
49   parser.add_option('-M', '--max',    action='store',dest='vmax',      type="float",   default=None,  help='bounding maximum value [max]')
50   parser.add_option('--div',          action='store',dest='ndiv',      type="int",     default=10,    help='number of divisions in histogram [10]')
51   parser.add_option('-S', '--save',   action='store',dest='save',      type="string",  default="gui", help='save mode (gui,png,eps,svg,pdf,txt,html,avi) [gui]')
52
53
54   #############################
55   ### Get options and variables
56   (opt,args) = parser.parse_args()
57
58   sparse(oplot=None,file=opt.file,ftype=opt.ftype,xvar=opt.xvar,yvar=opt.yvar,zvar=opt.zvar,cvar=opt.cvar,xmin=opt.xmin,xmax=opt.xmax,ymin=opt.ymin,ymax=opt.ymax,zmin=opt.zmin,zmax=opt.zmax,cmin=opt.cmin,cmax=opt.cmax,cond=opt.cond,merge=opt.merge,clb=opt.clb,trans=opt.trans,proj=opt.proj,blat=opt.blat,blon=opt.blon,back=opt.back,vmin=opt.vmin,vmax=opt.vmax,ndiv=opt.ndiv,save=opt.save)
59
60def sparse(oplot=None,file=None,ftype=None,xvar=None,yvar=None,zvar=None,cvar=None,xmin=None,xmax=None,ymin=None,ymax=None,zmin=None,zmax=None,cmin=None,cmax=None,cond=None,merge=False,clb="def",trans=1.,proj=None,blat=None,blon=None,back=None,vmin=None,vmax=None,ndiv=10,save="gui"):
61
62   import sys
63   from os import system,path
64   from scipy.io.idl import readsav
65   import numpy as np
66   from myplot import separatenames, definesubplot, errormess, defcolorb, fmtvar, polarinterv, simplinterv, define_proj, readdata, makeplotres
67   from mymath import min, max, writeascii
68   import matplotlib as mpl
69   from matplotlib.pyplot import subplot, figure, plot, scatter, colorbar, show,  title, close, legend, xlabel, ylabel, axis, hist
70   from matplotlib.cm import get_cmap
71   from mpl_toolkits.mplot3d import Axes3D
72
73
74   #############################
75   ### Load and check data
76
77   #############################
78   ### Load and check data
79
80   if file is None:
81      print "You must specify at least a file to process with -f."
82      exit()
83   if xvar is None:
84      print "You must specify at least a 1st field with -X."
85      exit()
86   if proj is not None and yvar is None:
87      print "Why did you ask a projection with only one variable?"
88      exit()
89     
90   filename=separatenames(file)
91   
92   #print 'conditions', cond
93   
94   print vmax
95   print vmin
96   
97   if cond is None: nslices = 1
98   else:                nslices = len(cond)
99   numplot   = nslices
100   if merge is False: numplot = numplot*len(filename)
101   
102   print ' '
103   if merge is False:
104      print nslices, 'condition(s) and', len(filename), 'file(s) without merging ---->', numplot, 'plot(s)'
105   else:
106      print nslices, 'condition(s) and', len(filename), 'file(s) with merging ---->', numplot, 'plot(s)'
107   print ' '
108
109   all_type = [[]]*len(filename)
110   all_x  = [[]]*len(filename)
111   all_y  = [[]]*len(filename)
112   all_z  = [[]]*len(filename)
113   all_c  = [[]]*len(filename)
114   all_data  = [[]]*len(filename)
115   #index  = [[]]*len(filename)
116
117
118##############################
119##################   READ DATA
120
121   for k in range(len(filename)):
122   
123##### Find file type
124     if ftype is None:
125      if filename[k].find('.sav') is not -1:
126         all_type[k] ='sav'
127         print '.sav file detected for', filename[k],'!!'
128      elif filename[k].find('.txt') is not -1:
129         all_type[k] ='txt'
130         #print '.txt file detected for', filename[k],'!!'
131      else:
132         all_type[k] ='txt'
133         print 'no file type detected for', filename[k],'!!, default type is', all_type[k]
134     
135      if all_type[k] != all_type[0]:
136         print 'Two different types were used: ', all_type[k], all_type[0]
137         errormess('Not suported')
138     else:
139        all_type[k] = ftype       
140
141##### Read file       
142     if all_type[k] == 'sav':
143       print 'reading .sav file', filename[k], '...'
144       data = {}
145       data = readsav(filename[k], idict=data, python_dict=False) #, verbose=True)
146       all_data[k] = data
147     elif all_type[k] == 'txt':
148       print 'reading .text file', filename[k], '...'
149       all_data[k] = np.loadtxt(filename[k])
150       
151     all_x[k] = readdata(all_data,all_type[k],k,xvar)
152     if yvar is not None: all_y[k] = readdata(all_data,all_type[k],k,yvar)
153     else:                    all_y[k] = None
154     if zvar is not None: all_z[k] = readdata(all_data,all_type[k],k,zvar)
155     else:                    all_z[k] = None
156     if cvar is not None: all_c[k] = readdata(all_data,all_type[k],k,cvar)
157     else:                    all_c[k] = None
158       
159     if merge is True and k >=1 :
160          all_x[0] = np.concatenate((all_x[0],all_x[k]))
161          if yvar is not None: all_y[0] = np.concatenate((all_y[0],all_y[k]))
162          if zvar is not None: all_z[0] = np.concatenate((all_z[0],all_z[k]))
163          if cvar is not None: all_c[0] = np.concatenate((all_c[0],all_c[k]))
164       
165   
166   
167##############################
168##################   PLOT DATA
169
170## Open a figure and set subplots
171   if oplot is not None: fig=oplot
172   else: fig = figure()
173   subv,subh = definesubplot(numplot,fig)
174   palette = get_cmap(name=clb)
175   
176   
177   for nplot in range(numplot):
178       
179       print ' '
180       
181###### Find which data and which file for plot nplot
182       if merge is False:
183          index_s = ((nplot)//len(filename))%nslices
184          index_f = ((nplot-1)//nslices)%len(filename)
185       else:
186          index_s = ((nplot)//1)%nslices
187          index_f = ((nplot-1)//nslices)%1
188       #print 'nplot,numplot,index_f,index_s', nplot,numplot,index_f,index_s
189       
190       
191###### Select point under conditions defined in -w option
192       index = np.isfinite(all_x[index_f])
193       if cond is not None:
194           zecond = separatenames(cond[index_s])
195           #print 'hello', cond[index_s], zecond
196           for i in range(len(zecond)):
197              zecondi = zecond[i] 
198              if zecondi.find('>') != -1:
199                 variable = zecondi[0:zecondi.find('>')]
200                 value = float(zecondi[zecondi.find('>')+1:len(zecondi)])
201                 if merge is True: # ultra moche de reconcatener a chaque fois, mais bon on s'en fout c'est du python
202                    zedata = []
203                    for k in range(len(filename)):
204                       zedata = np.concatenate((zedata,readdata(all_data,all_type[k],k,variable)))
205                 else : zedata = readdata(all_data,all_type[k],index_f,variable)
206                 #print index.shape,zedata.shape, value
207                 index = index*(zedata>value)
208                 print 'All points such as', variable, '>', value, 'have been selected'
209              elif zecondi.find('<') != -1:
210                 variable = zecondi[0:zecondi.find('<')]
211                 value = float(zecondi[zecondi.find('<')+1:len(zecondi)])
212                 if merge is True:
213                    zedata = []
214                    for k in range(len(filename)):
215                       zedata = np.concatenate((zedata,readdata(all_data,all_type[k],k,variable)))
216                 else : zedata = readdata(all_data,all_type[k],index_f,variable)
217                 #print index.shape,zedata.shape, value
218                 index = index*(zedata<value)
219                 print 'All points such as', variable, '<', value, 'have been selected'
220              else:
221                 print ''
222                 print 'I do not understand that condition :', zecondi
223                 errormess('')
224       
225       
226       if np.sum(index) == 0:
227           print '***********  WARNING  ***********'
228           print '***********  NO DATA  ***********'
229           errormess('')
230       else:
231           print np.sum(index),'points have been selected among', len(all_x[index_f]), \
232           'that is to say %2.2f' %(float(100*np.sum(index))/float(len(all_x[index_f]))), '%.'
233       
234         
235       #print 'numplot', numplot
236       changesubplot = (numplot > 1) #and (len(what_I_plot.shape) != 1)
237       if changesubplot: subplot(subv,subh,nplot+1)
238       #print 'subv,subh,nplot', subv,subh,nplot
239       
240###### Projection
241       if proj is not None: # Nota : NEVER TRY TO DO A MESHGRID ON SPARSE DATA. WAY TOO MUCH POINTS.
242          wlon = [min(all_x[index_f][index]),max(all_x[index_f][index])]
243          wlat = [min(all_y[index_f][index]),max(all_y[index_f][index])]
244          m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)
245          x, y = m(all_x[index_f][index],all_y[index_f][index])
246       else:
247          x = all_x[index_f][index]
248          if yvar is not None: y = all_y[index_f][index]                             
249       
250###### Plot: 1D histogram
251       if yvar is None:
252          hist(x,bins=ndiv,histtype='step',linewidth=2)
253          if save == 'txt':
254              print 'saving file profile'+str(nplot+1)+'.txt'
255              writeascii(np.transpose(x),'profile'+str(nplot+1)+'.txt')
256###### Plot: 2D cloud
257       elif zvar is None and cvar is None :
258          plot(x,y,'.b')
259          if save == 'txt':
260               print 'saving file profile'+str(nplot+1)+'.txt'
261               writeascii(np.transpose(np.array([x,y])),'profile'+str(nplot+1)+'.txt')
262###### Plot: 2D cloud with color
263       elif zvar is None and cvar is not None :
264          if save == 'txt':
265               print 'saving file profile'+str(nplot+1)+'.txt'
266               writeascii(np.transpose(np.array([x,y,all_c[index_f][index]])),'profile'+str(nplot+1)+'.txt')
267          scatter(x,y,c=all_c[index_f][index],\
268          marker='o', edgecolor='None',cmap = palette, alpha=trans, vmin = vmin,vmax = vmax)
269###### Plot: 3D cloud
270       elif zvar is not None and cvar is None :
271          ax = fig.add_subplot(subv,subh,nplot+1, projection='3d')
272          ax.plot(x,y,all_z[index_f][index],'.b')
273###### Plot: 3D cloud with color
274       else :
275          ax = fig.add_subplot(subv,subh,nplot+1, projection='3d')
276          ax.scatter(x,y,all_z[index_f][index],c=all_c[index_f][index],\
277          marker='o', edgecolor='None', cmap = palette, alpha=trans,vmin = vmin,vmax = vmax)
278   
279###### Colorbar: http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
280       if clb != 'nobar' and all_c[index_f] is not None and all_z[index_f] is None: # pourquoi la colorbar marche pas en 3d?
281          colorbar( fraction=0.05,pad=0.03,format='%0.2f',\
282          extend='neither',spacing='proportional' )
283       
284###### Plot limits (options --xmin, --xmax, etc ..)
285       if proj is None:
286          xlabel(xvar)
287          ylabel(yvar)
288          if zvar is None :
289             if xmin is not None: mpl.pyplot.xlim(xmin=xmin)
290             else:                    mpl.pyplot.xlim(xmin=min(all_x[index_f][index]))
291             if xmax is not None: mpl.pyplot.xlim(xmax=xmax)
292             else:                    mpl.pyplot.xlim(xmax=max(all_x[index_f][index]))
293             if yvar is not None:
294                if ymin is not None: mpl.pyplot.ylim(ymin=ymin)
295                else:                    mpl.pyplot.ylim(ymin=min(all_y[index_f][index]))
296                if ymax is not None: mpl.pyplot.ylim(ymax=ymax)
297                else:                    mpl.pyplot.ylim(ymax=max(all_y[index_f][index]))
298
299         
300       if cond is not None:
301          title(cond[index_s])
302       else:
303          title('all point selected')
304       
305   zeplot = "output"
306   if save in ['png','eps','svg','pdf']:     makeplotres(zeplot,ext=save)#,res=resolution,pad_inches_value=pad_inches_value,disp=False,ext=save)
307   elif save == 'gui':                       show()
308   elif save == 'return':                    return mpl.pyplot
309   else:                                         print "INFO: save mode not supported. using gui instead." ; show()
310   print ''
311
312   command = ""
313   for arg in sys.argv: command = command + arg + ' '
314   name = zeplot
315   f = open(name+'.sh', 'w')
316   f.write(command)
317
318   
319
320
321   
Note: See TracBrowser for help on using the repository browser.