Index: trunk/UTIL/PYTHON/README.SPARSE
===================================================================
--- trunk/UTIL/PYTHON/README.SPARSE	(revision 568)
+++ trunk/UTIL/PYTHON/README.SPARSE	(revision 568)
@@ -0,0 +1,72 @@
+*****************************************
+*****************************************
+*****************************************
+         SPARSE.PY EXAMPLES
+*****************************************
+               Author : TN
+*****************************************
+  DON'T FORGET YOUR BEST FRIEND IS
+   sparse.py -h [or] sparse.py --help
+*****************************************
+*****************************************
+*****************************************
+
+*****************************************
+SIMPLE EXAMPLES on a SAMPLE TES DATA FILE
+*****************************************
+
+sparse.py reads sparse data, i.e. not binned in a lat/lon/vert/time grid and does plots.
+
+
+Goal: I would like to plot surface temperature vs. latitude for a .sav file from IDL.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y surft
+[files in /planeto/lmlmd/TESdata/]
+
+Goal: I would like the same plot, but with local time for points color.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y surft -C loctime
+
+Goal: I would like the same plot, but for local time > 10.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y surft -C loctime -w 'loctime>12'
+
+Goal: I would like the same plot, but for longitude between 170 and 180.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y surft -C loctime -w 'loctime>12,lon>170,lon<180'
+
+Goal: I would like another plot for local time < 12.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y surft -C loctime -w 'loctime>12,lon>170,lon<180' -w 'loctime<12,lon>170,lon<180'
+
+Goal: I would the same plot with the same 'spectral' colorbar.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y surft -C loctime -w 'loctime>12,lon>170,lon<180' -w 'loctime<12,lon>170,lon<180' -c spectral -m 0 -M 24
+
+Goal: I would like to map surface temperature on a polar projection.
+sparse.py -f TES_MY26_90_120.sav -X lon -Y lat -C surft -w 'lat>70,loctime>12' -p npstere --blat 70
+
+Goal: I would like to plot surface temperature vs LS for Northern pole ice deposits for two different files.
+sparse.py -f TES_MY26_90_120.sav,TES_MY26_120_150.sav -X areoc -Y lat -C surft -w 'lat>70,loctime>10,lon>160,lon<180'
+
+Goal: I would like the same plot, but with merged data from the two files.
+sparse.py -f TES_MY26_90_120.sav,TES_MY26_120_150.sav -X areoc -Y lat -C surft -w 'lat>70,loctime>10,lon>160,lon<180' --merge
+
+Goal: I would like to plot surface temperature occurences for day and night.
+sparse.py -f TES_MY26_90_120.sav -X surft -w 'loctime>12,areoc>110' -w  'loctime<12,areoc>110' --div 50
+
+Goal: I would like the plot surface temperature and surface pressure for z-axis.
+sparse.py -f TES_MY26_90_120.sav -X lat -Y lon -C surft -Z surfp -w 'loctime>12,lon>170,lon<180,lat>75'
+
+Goal: I would like to plot data from .txt files, knowing that 1st column is latitude, 2nd is longitude and 3rd is water vapor.
+sparse.py -f media_prum_SPICAM.dat,media_prum_TES.dat,media_prum_OMEGA.dat -X 1 -Y 2 -C 3 -c 'spectral' -w '3>1' -m 0 -M 65
+[files in /planeto/tnalmd/data/water_vapor]
+
+
+sparse.py also supports min and max options (--xmin,--xmax,-m,-M,etc..), background (-b) and 
+transparency (--trans) options from planetoplot.
+
+
+*****************************************
+to do:
+*****************************************
+- read other formats than txt and sav ?
+- variable that depends on 2 variables, such as atmospheric temperatures measurements
+- interface with binned data for comparison with simulations (use of writeascii in planetoplot, etc...)
+
+
+
Index: trunk/UTIL/PYTHON/myplot.py
===================================================================
--- trunk/UTIL/PYTHON/myplot.py	(revision 567)
+++ trunk/UTIL/PYTHON/myplot.py	(revision 568)
@@ -223,6 +223,6 @@
         subh = 99999
     elif numplot == 1: 
-        subv = 99999
-        subh = 99999
+        subv = 1
+        subh = 1
     elif numplot == 2:
         subv = 1 #2
@@ -999,4 +999,21 @@
            
   return zesaxis
+
+## Author: TN
+def readdata(data,datatype,coord1,coord2):
+## Read sparse data
+  from numpy import empty
+  if datatype is 'txt':
+     if len(data[coord1].shape) == 1:
+         return data[coord1][:]
+     elif len(data[coord1].shape) == 2:
+         return data[coord1][:,int(coord2)-1]
+     else:
+         errormess('error in readdata')
+  elif datatype is 'sav':
+     return data[coord1][coord2]
+  else:
+     errormess(datatype+' type is not supported!')
+
 
 ## Author: AS
Index: trunk/UTIL/PYTHON/sparse.py
===================================================================
--- trunk/UTIL/PYTHON/sparse.py	(revision 568)
+++ trunk/UTIL/PYTHON/sparse.py	(revision 568)
@@ -0,0 +1,300 @@
+#!/usr/bin/env python
+
+### T.Navarro
+
+##
+#   This routine reads sparse data, i.e. not binned in a lat/lon/vert/time grid and does plots
+#
+#   There is nothing gnuplot could not be able to do for now, except command line and ease and all python plots possibilities, 
+#   but the goal is to provide a framework to compare with data from simulations on a constant grid.
+###
+
+###########################################################################################
+###########################################################################################
+### What is below relate to running the file as a command line executable (very convenient)
+if __name__ == "__main__":
+   import sys
+   from optparse import OptionParser    ### to be replaced by argparse
+   from os import system,path
+   from scipy.io.idl import readsav
+   import numpy as np
+   from myplot import separatenames, definesubplot, errormess, defcolorb, fmtvar, polarinterv, simplinterv, define_proj, readdata
+   from mymath import min, max, writeascii
+   import matplotlib as mpl
+   from matplotlib.pyplot import subplot, figure, plot, scatter, colorbar, show,  title, close, legend, xlabel, ylabel, axis, hist
+   from matplotlib.cm import get_cmap
+   from mpl_toolkits.mplot3d import Axes3D
+   parser = OptionParser() 
+   
+   #############################
+   ### Options
+   parser.add_option('-f', '--file',    action='store',dest='file',     type="string",  default=None,  help='[NEEDED] filename, comma separated')
+   parser.add_option('--ftype',         action='store',dest='ftype',   default=None, help=' Force data file type [\'sav\',\'txt\']')
+   parser.add_option('-X', '--xvar',    action='store',dest='xvar',     type="string",  default=None,  help='[NEEDED] x variable: a name, a number, etc ...')
+   parser.add_option('-Y', '--yvar',    action='store',dest='yvar',     type="string",  default=None,  help='y variable')
+   parser.add_option('-Z', '--zvar',    action='store',dest='zvar',     type="string",  default=None,  help='z variable')
+   parser.add_option('-C', '--cvar',    action='store',dest='cvar',     type="string",  default=None,  help='c variable for color')
+   parser.add_option('--xmin',          action='store',dest='xmin',     type="float",   default=None,  help='min values for x var, comma separated')
+   parser.add_option('--xmax',          action='store',dest='xmax',     type="float",   default=None,  help='max values for x var, comma separated')
+   parser.add_option('--ymin',          action='store',dest='ymin',     type="float",   default=None,  help='min values for y var, comma separated')
+   parser.add_option('--ymax',          action='store',dest='ymax',     type="float",   default=None,  help='max values for y var, comma separated')
+   parser.add_option('--zmin',          action='store',dest='zmin',     type="float",   default=None,  help='min values for z var, comma separated')
+   parser.add_option('--zmax',          action='store',dest='zmax',     type="float",   default=None,  help='max values for z var, comma separated')
+   parser.add_option('--cmin',          action='store',dest='cmin',     type="string",  default=None,  help='min values for c var, comma separated')
+   parser.add_option('--cmax',          action='store',dest='cmax',     type="string",  default=None,  help='max values for c var, comma separated')
+   parser.add_option('-w','--with',     action='append',dest='cond',    type="string",  default=None,  help='conditions..')
+   parser.add_option('--merge',   action='store_true', dest='merge',                      default=False, help='merge datafiles in a single plot [False]')
+   parser.add_option('-c','--color',  action='store',dest='clb',       type="string",  default="def", help='change colormap (also: nobar,onebar)')
+   parser.add_option('--trans',        action='store',dest='trans',     type="float",   default=1.,    help='shaded plot transparency, 0 to 1 (=opaque) [1]')
+   parser.add_option('-p', '--proj',   action='store',dest='proj',      type="string",  default=None,  help='projection')
+   parser.add_option('--blat',         action='store',dest='blat',      type="int",     default=None,  help='reference lat (or bounding lat for stere) [computed]')
+   parser.add_option('--blon',         action='store',dest='blon',      type="int",     default=None,  help='reference lon [computed]')
+   parser.add_option('-b', '--back',   action='store',dest='back',      type="string",  default=None,  help='background image [None]')
+   parser.add_option('-m', '--min',    action='store',dest='vmin',      type="float",   default=None,  help='bounding minimum value [min]')    
+   parser.add_option('-M', '--max',    action='store',dest='vmax',      type="float",   default=None,  help='bounding maximum value [max]') 
+   parser.add_option('--div',          action='store',dest='ndiv',      type="int",     default=10,    help='number of divisions in histogram [10]')
+   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]')
+
+
+   #############################
+   ### Get options and variables
+   (opt,args) = parser.parse_args()
+
+   #############################
+   ### Load and check data
+
+   #############################
+   ### Load and check data
+
+   if opt.file is None:
+      print "You must specify at least a file to process with -f."
+      exit()
+   if opt.xvar is None:
+      print "You must specify at least a 1st field with -X."
+      exit()
+   if opt.proj is not None and opt.yvar is None:
+      print "Why did you ask a projection with only one variable?"
+      exit()
+      
+   filename=separatenames(opt.file)
+   
+   #print 'conditions', opt.cond
+   
+   print opt.vmax
+   print opt.vmin
+   
+   if opt.cond is None: nslices = 1
+   else:                nslices = len(opt.cond)
+   numplot   = nslices
+   if opt.merge is False: numplot = numplot*len(filename)
+   
+   print ' '
+   if opt.merge is False:
+      print nslices, 'condition(s) and', len(filename), 'file(s) without merging ---->', numplot, 'plot(s)'
+   else:
+      print nslices, 'condition(s) and', len(filename), 'file(s) with merging ---->', numplot, 'plot(s)'
+   print ' '
+
+   all_type = [[]]*len(filename)
+   all_x  = [[]]*len(filename)
+   all_y  = [[]]*len(filename)
+   all_z  = [[]]*len(filename)
+   all_c  = [[]]*len(filename)
+   all_data  = [[]]*len(filename)
+   #index  = [[]]*len(filename)
+
+
+##############################
+##################   READ DATA
+
+   for k in range(len(filename)):
+   
+##### Find file type
+     if filename[k].find('.sav') is not -1:
+        all_type[k] ='sav'
+        print '.sav file detected for', filename[k],'!!'
+     elif filename[k].find('.txt') is not -1:
+        all_type[k] ='txt'
+        #print '.txt file detected for', filename[k],'!!'
+     else:
+        all_type[k] ='txt'
+        print 'no file type detected for', filename[k],'!!, default type is', all_type[k]
+     
+     if all_type[k] != all_type[0]:
+        print 'Two different types were used: ', all_type[k], all_type[0]
+        errormess('Not suported')
+
+
+##### Read file       
+     if all_type[k] is 'sav':
+       print 'reading .sav file', filename[k], '...'
+       data = {}
+       data = readsav(filename[k], idict=data, python_dict=False) #, verbose=True)
+       all_data[k] = data
+     elif all_type[k] is 'txt':
+       print 'reading .text file', filename[k], '...'
+       all_data[k] = np.loadtxt(filename[k])
+       
+     all_x[k] = readdata(all_data,all_type[k],k,opt.xvar)
+     if opt.yvar is not None: all_y[k] = readdata(all_data,all_type[k],k,opt.yvar)
+     else:                    all_y[k] = None
+     if opt.zvar is not None: all_z[k] = readdata(all_data,all_type[k],k,opt.zvar)
+     else:                    all_z[k] = None
+     if opt.cvar is not None: all_c[k] = readdata(all_data,all_type[k],k,opt.cvar)
+     else:                    all_c[k] = None
+       
+     if opt.merge is True and k >=1 :
+          all_x[0] = np.concatenate((all_x[0],all_x[k]))
+          if opt.yvar is not None: all_y[0] = np.concatenate((all_y[0],all_y[k]))
+          if opt.zvar is not None: all_z[0] = np.concatenate((all_z[0],all_z[k]))
+          if opt.cvar is not None: all_c[0] = np.concatenate((all_c[0],all_c[k]))
+       
+   
+   
+##############################
+##################   PLOT DATA
+
+## Open a figure and set subplots
+   fig = figure()
+   subv,subh = definesubplot(numplot,fig)
+   palette = get_cmap(name=opt.clb)
+   
+   
+   for nplot in range(numplot):
+       
+       print ' '
+       
+###### Find which data and which file for plot nplot
+       if opt.merge is False:
+          index_s = ((nplot)//len(filename))%nslices
+          index_f = ((nplot-1)//nslices)%len(filename)
+       else:
+          index_s = ((nplot)//1)%nslices
+          index_f = ((nplot-1)//nslices)%1
+       #print 'nplot,numplot,index_f,index_s', nplot,numplot,index_f,index_s
+       
+       
+###### Select point under conditions defined in -w option
+       index = np.isfinite(all_x[index_f])
+       if opt.cond is not None:
+           zecond = separatenames(opt.cond[index_s])
+           #print 'hello', opt.cond[index_s], zecond
+           for i in range(len(zecond)):
+              zecondi = zecond[i] 
+              if zecondi.find('>') != -1:
+                 variable = zecondi[0:zecondi.find('>')]
+                 value = float(zecondi[zecondi.find('>')+1:len(zecondi)])
+                 if opt.merge is True: # ultra moche de reconcatener a chaque fois, mais bon on s'en fout c'est du python
+                    zedata = []
+                    for k in range(len(filename)):
+                       zedata = np.concatenate((zedata,readdata(all_data,all_type[k],k,variable)))
+                 else : zedata = readdata(all_data,all_type[k],index_f,variable)
+                 #print index.shape,zedata.shape, value
+                 index = index*(zedata>value)
+                 print 'All points such as', variable, '>', value, 'have been selected'
+              elif zecondi.find('<') != -1:
+                 variable = zecondi[0:zecondi.find('<')]
+                 value = float(zecondi[zecondi.find('<')+1:len(zecondi)])
+                 if opt.merge is True:
+                    zedata = []
+                    for k in range(len(filename)):
+                       zedata = np.concatenate((zedata,readdata(all_data,all_type[k],k,variable)))
+                 else : zedata = readdata(all_data,all_type[k],index_f,variable)
+                 #print index.shape,zedata.shape, value
+                 index = index*(zedata<value)
+                 print 'All points such as', variable, '<', value, 'have been selected'
+              else:
+                 print ''
+                 print 'I do not understand that condition :', zecondi
+                 errormess('')
+       
+       
+       if np.sum(index) == 0:
+           print '***********  WARNING  ***********'
+           print '***********  NO DATA  ***********'
+           errormess('')
+       else:
+           print np.sum(index),'points have been selected among', len(all_x[index_f]), \
+           'that is to say %2.2f' %(float(100*np.sum(index))/float(len(all_x[index_f]))), '%.'
+       
+          
+       #print 'numplot', numplot
+       changesubplot = (numplot > 1) #and (len(what_I_plot.shape) != 1)
+       if changesubplot: subplot(subv,subh,nplot+1)
+       #print 'subv,subh,nplot', subv,subh,nplot
+       
+###### Projection
+       if opt.proj is not None: # Nota : NEVER TRY TO DO A MESHGRID ON SPARSE DATA. WAY TOO MUCH POINTS.
+          wlon = [min(all_x[index_f][index]),max(all_x[index_f][index])]
+          wlat = [min(all_y[index_f][index]),max(all_y[index_f][index])]
+          m = define_proj(opt.proj,wlon,wlat,back=opt.back,blat=opt.blat,blon=opt.blon)
+          x, y = m(all_x[index_f][index],all_y[index_f][index])
+       else:
+          x = all_x[index_f][index]
+          if opt.yvar is not None: y = all_y[index_f][index]                              
+       
+###### Plot: 1D histogram
+       if opt.yvar is None:
+          hist(x,bins=opt.ndiv,histtype='step',linewidth=2)
+          if opt.save == 'txt':
+              print 'saving file profile'+str(nplot+1)+'.txt'
+              writeascii(np.transpose(x),'profile'+str(nplot+1)+'.txt')
+###### Plot: 2D cloud
+       elif opt.zvar is None and opt.cvar is None :
+          plot(x,y,'.b')
+          if opt.save == 'txt':
+               print 'saving file profile'+str(nplot+1)+'.txt'
+               writeascii(np.transpose(np.array([x,y])),'profile'+str(nplot+1)+'.txt')
+###### Plot: 2D cloud with color
+       elif opt.zvar is None and opt.cvar is not None :
+          if opt.save == 'txt':
+               print 'saving file profile'+str(nplot+1)+'.txt'
+               writeascii(np.transpose(np.array([x,y,all_c[index_f][index]])),'profile'+str(nplot+1)+'.txt')
+          scatter(x,y,c=all_c[index_f][index],\
+          marker='o', edgecolor='None',cmap = palette, alpha=opt.trans, vmin = opt.vmin,vmax = opt.vmax)
+###### Plot: 3D cloud
+       elif opt.zvar is not None and opt.cvar is None :
+          ax = fig.add_subplot(subv,subh,nplot+1, projection='3d')
+          ax.plot(x,y,all_z[index_f][index],'.b')
+###### Plot: 3D cloud with color
+       else :
+          ax = fig.add_subplot(subv,subh,nplot+1, projection='3d')
+          ax.scatter(x,y,all_z[index_f][index],c=all_c[index_f][index],\
+          marker='o', edgecolor='None', cmap = palette, alpha=opt.trans,vmin = opt.vmin,vmax = opt.vmax)
+   
+###### Colorbar: http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
+       if opt.clb != 'nobar' and all_c[index_f] is not None and all_z[index_f] is None: # pourquoi la colorbar marche pas en 3d?
+          colorbar( fraction=0.05,pad=0.03,format='%0.2f',\
+          extend='neither',spacing='proportional' )
+        
+###### Plot limits (options --xmin, --xmax, etc ..) 
+       if opt.proj is None:
+          xlabel(opt.xvar)
+          ylabel(opt.yvar)
+          if opt.zvar is None :
+             if opt.xmin is not None: mpl.pyplot.xlim(xmin=opt.xmin)
+             else:                    mpl.pyplot.xlim(xmin=min(all_x[index_f][index]))
+             if opt.xmax is not None: mpl.pyplot.xlim(xmax=opt.xmax)
+             else:                    mpl.pyplot.xlim(xmax=max(all_x[index_f][index]))
+             if opt.yvar is not None:
+                if opt.ymin is not None: mpl.pyplot.ylim(ymin=opt.ymin)
+                else:                    mpl.pyplot.ylim(ymin=min(all_y[index_f][index]))
+                if opt.ymax is not None: mpl.pyplot.ylim(ymax=opt.ymax)
+                else:                    mpl.pyplot.ylim(ymax=max(all_y[index_f][index]))
+
+          
+       if opt.cond is not None:
+          title(opt.cond[index_s])
+       else:
+          title('all point selected')
+       
+
+
+   
+   show()
+   print ''
+
+   
+
+
+   
