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