Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py (revision 345) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py (revision 349) @@ -27,5 +27,5 @@ ### LOAD NETCDF DATA - filename = "psfc.nc" + #filename = "/home/aymeric/Big_Data/psfc.nc" nc = Dataset(filename) psfc = nc.variables["PSFC"] @@ -36,4 +36,5 @@ allsizesx = [] allsizesy = [] + depression = [] stride = 1 #5 for i in range(0,shape[0],stride): @@ -52,7 +53,10 @@ fac = 3.5 #fac = 2.5 + fac = 3. lim = ave - fac*std where = np.where(psfc2d < lim) ############### END CRITERION + + depression = np.append(depression,np.ravel(psfc2d[where])-ave) lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine @@ -115,5 +119,5 @@ allsizesy.sort() - return allsizesx, allsizesy + return allsizesx, allsizesy, depression ######################################################################### @@ -125,13 +129,22 @@ import matplotlib.mlab as mlab import mymath as mym +import myplot as myp + +import plfit +import plplot +import randht +import plpva save = True save = False +pression = False +pression = True if save: - allsizesx, allsizesy = getsize("psfc.nc") + allsizesx, allsizesy, depression = getsize("/home/aymeric/Big_Data/psfc.nc") ### sauvegarde texte pour inspection mym.writeascii(allsizesx,'allsizex.txt') mym.writeascii(allsizesy,'allsizey.txt') + mym.writeascii(depression,'alldepression.txt') ### sauvegarde binaire pour utilisation python myfile = open('allsizex.bin', 'wb') @@ -141,4 +154,7 @@ pickle.dump(allsizesy, myfile) myfile.close() + myfile = open('alldepression.bin', 'wb') + pickle.dump(depression, myfile) + myfile.close() ### load files @@ -147,12 +163,21 @@ myfile = open('allsizey.bin', 'r') allsizesy = pickle.load(myfile) - -### append sizes -plothist = np.append(allsizesx,allsizesy) +myfile = open('alldepression.bin', 'r') +depression = pickle.load(myfile) +depression = np.array(abs(depression))#*1000. + +### sizes +#plothist = np.append(allsizesx,allsizesy) plothist = allsizesx -#plothist = allsizesy +if pression: plothist = depression plothist.sort() - -### make bins +print 'mean ', np.mean(plothist,dtype=np.float64) +print 'std ', np.std(plothist,dtype=np.float64) +print 'max ', np.max(plothist) +print 'min ', np.min(plothist) +print 'len ', len(plothist) + + +### MAKE BINS nbins = 100 zebins = [2.0] @@ -165,24 +190,122 @@ zebins = [2./np.sqrt(2.)] ## ne pas tomber sur une dizaine ronde nbins = 200 -zebins = [0.2/np.sqrt(2.)] ## ne pas tomber sur une dizaine ronde + +if pression: zebins = [0.3] for i in range(0,nbins): zebins.append(zebins[i]*np.sqrt(2)) zebins = np.array(zebins) -### select reasonable bins for DD -zebins = zebins [ zebins > 10. ] -zebins = zebins [ zebins < 1000. ] -#print 'bins ',zebins +#### select reasonable bins for DD +if not pression: + zebins = zebins [ zebins > 15. ] + zebins = zebins [ zebins < 1000. ] +else: + zebins = zebins [ zebins < 10. ] print 'corrected bins ',zebins -print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy) - -#plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.] ## pour calibrer #### HISTOGRAM +plt.figure(1) plt.hist( plothist,\ log=True,\ bins=zebins,\ # cumulative=-1,\ + normed=True,\ ) plt.xscale('log') -plt.show() - +if pression: plt.xlabel('Pressure (Pa)') +else: plt.xlabel('Diameter (m)') +plt.ylabel('Population (normalized)') +if pression: prefix="p" +else: prefix="" +myp.makeplotres(prefix+"histogram",res=200,disp=False) +plt.close(1) + +### COMPARED HISTOGRAMS +### --- FIT WITH POWER LAW +if pression: [alpha, xmin, L] = plfit.plfit(plothist,'xmin',0.3) +else: [alpha, xmin, L] = plfit.plfit(plothist,'limit',20.) +print alpha,xmin + +#a = plpva.plpva(plothist,0.75,'xmin',0.75) +#print a + +#### DEUXIEME ROUTINE +####IL FAUT UTILISER LE DISCRET POUR LA TAILLE !!! +#if pression: myplfit = plfit.plfit(plothist,verbose=True,xmin=0.75) +#else: myplfit = plfit.plfit(plothist,verbose=True,xmin=20.) +#myplfit.plotppf() +#plt.show() +#exit() + + +#plt.figure(1) +#h = plplot.plplot(plothist,xmin,alpha) +#myp.makeplotres(prefix+"fit",res=200,disp=False) +plt.figure(2) +### --- POWER LAW (factor does not really matter) +power = (xmin/2.2)*np.array(randht.randht(10000,'powerlaw',alpha)) +#power = (xmin/2.2)*np.array(randht.randht(10000,'cutoff',alpha,10.)) ##marche pas si trop grand +print 'mean ', np.mean(power,dtype=np.float64) +### --- EXPONENTIAL LAW +expo = randht.randht(10000,'exponential',1./(np.mean(power,dtype=np.float64)*1.00)) +print 'mean ', np.mean(expo,dtype=np.float64) +### --- PLOT +plt.hist( [plothist,power,expo],\ + label=['LES vortices','Power law '+'{:.1f}'.format(alpha),'Exponential law'],\ + log=True,\ + bins=zebins,\ +# cumulative=-1,\ + normed=True,\ + ) +plt.legend() +plt.xscale('log') +if pression: plt.xlabel('Pressure (Pa)') +else: plt.xlabel('Diameter (m)') +plt.ylabel('Population (normalized)') +myp.makeplotres(prefix+"comparison",res=200,disp=False) +plt.close(2) + +######################## +######################## +zebins = [30.,42.,60.,84.,120.,170.,240.,340.] +plothist = [] +plothist = np.append(plothist,30 *np.ones(306)) +plothist = np.append(plothist,42 *np.ones(58) ) +plothist = np.append(plothist,60 *np.ones(66) ) +plothist = np.append(plothist,84 *np.ones(41) ) +plothist = np.append(plothist,120*np.ones(19) ) +plothist = np.append(plothist,170*np.ones(9) ) +plothist = np.append(plothist,240*np.ones(2) ) +plothist = np.append(plothist,340*np.ones(1) ) + +#zebins = [50.,71.,100.,141.,200.,282.,400.] +#plothist = [] +#plothist = np.append(plothist,50. *np.ones(36)) +#plothist = np.append(plothist,71. *np.ones(18) ) +#plothist = np.append(plothist,100. *np.ones(12) ) +#plothist = np.append(plothist,141. *np.ones(6) ) +#plothist = np.append(plothist,200.*np.ones(4) ) +#plothist = np.append(plothist,282.*np.ones(1) ) +#plothist = np.append(plothist,400.*np.ones(2) ) + +plt.figure(3) +[alpha, xmin, L] = plfit.plfit(plothist,'xmin',30)#50.) +print alpha,xmin +#a = plpva.plpva(plothist,30,'xmin',30) +h = plplot.plplot(plothist,xmin,alpha) +plt.loglog(h[0], h[1], 'k--',linewidth=2) +plt.hist( plothist,\ + log=True,\ + bins=zebins,\ +# cumulative=-1,\ + normed=True,\ + ) +plt.xscale('log') +plt.xlabel('Pressure (micro Pa)') +plt.ylabel('Population (normalized)') +myp.makeplotres("data",res=200,disp=False) + +#plt.figure(4) +#[alpha, xmin, L] = plfit.plfit(plothist,'xmin',50.) #,'xmin',30.) +#print alpha,xmin +#h = plplot.plplot(plothist,xmin,alpha) +#myp.makeplotres("datafit",res=200,disp=False) Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/gcm.py (revision 349) @@ -0,0 +1,158 @@ +#!/usr/bin/env python + +### T. Navarro + A. Spiga + +########################################################################################### +########################################################################################### +### 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 api_wrapper import api_onelevel + from netCDF4 import Dataset + from myplot import getlschar, separatenames, readslices, adjust_length + from os import system + from planetoplot import planetoplot + import numpy as np + + ############################# + ### Get options and variables + parser = OptionParser() + parser.add_option('-f', '--file', action='append',dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append). Plot files separated by comas in the same figure') + parser.add_option('-l', '--level', action='store',dest='nvert', type="float", default=0, help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)') + parser.add_option('-p', '--proj', action='store',dest='proj', type="string", default=None, help='projection') + parser.add_option('-b', '--back', action='store',dest='back', type="string", default=None, help='background image (def: None)') + parser.add_option('-t', '--target', action='store',dest='target', type="string", default=None, help='destination folder') + parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') + parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') + parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -*numplot*)') + parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)') + parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') + parser.add_option('-x', '--no-vect',action='store_false',dest='winds', default=True, help='no wind vectors') + parser.add_option('-m', '--min', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)') + parser.add_option('-M', '--max', action='append',dest='vmax', type="float", default=None, help='bounding maximum value (append)') + parser.add_option('-T', '--tiled', action='store_true',dest='tile', default=False, help='draw a tiled plot (no blank zone)') + parser.add_option('-z', '--zoom', action='store',dest='zoom', type="float", default=None, help='zoom factor in %') + parser.add_option('-N', '--no-api', action='store_true',dest='nocall', default=False, help='do not recreate api file') + parser.add_option('-d', '--display',action='store_false',dest='display', default=True, help='do not pop up created images') + parser.add_option('-e', '--itstep', action='store',dest='itstep', type="int", default=None, help='stride time (def=4)') + parser.add_option('-H', '--hole', action='store_true',dest='hole', default=False, help='holes above max and below min') + parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)') + parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly', default=False, help='compute and plot relative anomaly in %') + parser.add_option('-w', '--with', action='store',dest='var2', type="string", default=None, help='variable contoured') + parser.add_option('--div', action='store',dest='ndiv', type="int", default=10, help='number of divisions in colorbar (def: 10)') + parser.add_option('-F', '--first', action='store',dest='first', type="int", default=1, help='first subscript to plot (def: 1)') + parser.add_option('--mult', action='store',dest='mult', type="float", default=1., help='a multiplicative factor to plotted field') + parser.add_option('--title', action='store',dest='zetitle', type="string", default="fill",help='customize the whole title') + #parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables') + + ############# T.N. changes + #parser.add_option('-o','--operation',action='store',dest='operation',type="string", default=None ,help='matrix of operations between files (for now see code, sorry)') + parser.add_option('--lat', action='append',dest='slat',type="string", default=None, help='slices along latitude. One value, or two values separated by comas for averaging') + parser.add_option('--lon', action='append',dest='slon', type="string", default=None, help='slices along longitude. One value, or two values separated by comas for averaging') + parser.add_option('--vert', action='append',dest='svert',type="string", default=None, help='slices along vertical axis. One value, or two values separated by comas for averaging') # quelles coordonnees ? + parser.add_option('--time', action='append',dest='stime',type="string", default=None, help='slices along time. One value, or two values separated by comas for averaging') # quelles coordonnees ? + + + (opt,args) = parser.parse_args() + if opt.namefile is None: + print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h" + exit() + if opt.var is None and opt.anomaly is True: + print "Cannot ask to compute anomaly if no variable is set" + exit() + print "Options:", opt + + #listvar = '' + #if opt.var is None: + # zerange = [-999999] + #else: + # zelen = len(opt.var) + # zerange = range(zelen) + # #if zelen == 1: listvar = opt.var[0] + ',' + # #else : + # for jjj in zerange: listvar += opt.var[jjj] + ',' + # listvar = listvar[0:len(listvar)-1] + # vmintab = adjust_length (opt.vmin, zelen) + # vmaxtab = adjust_length (opt.vmax, zelen) + + print "namefile, length", opt.namefile, len(opt.namefile) + + zeslat = readslices(opt.slat) + zeslon = readslices(opt.slon) + zesvert = readslices(opt.svert) + zestime = readslices(opt.stime) + print "slat,zeslat", opt.slat, zeslat + print "slon,zeslon", opt.slon, zeslon + print "svert,zesvert", opt.svert, zesvert + print "stime,zestime", opt.stime, zestime + + for i in range(len(opt.namefile)): + for j in range(len(opt.var)): + + zenamefiles = separatenames(opt.namefile[i]) + print "zenamefiles", zenamefiles + + if opt.vmin is not None : zevmin = opt.vmin[min(i,len(opt.vmin)-1)] + else: zevmin = None + if opt.vmax is not None : zevmax = opt.vmax[min(i,len(opt.vmax)-1)] + else: zevmax = None + print "vmin, zevmin", opt.vmin, zevmin + print "vmax, zevmax", opt.vmax, zevmax + + zevar = separatenames(opt.var[j]) + zevar = zevar[0] + print "var, zevar", opt.var, zevar + + #checkcoherence(len(zenamefiles),len(opt.slat),len(opt.slon),len(opt.stime)) + + zefile = zenamefiles[0] + + zelevel = opt.nvert + stralt = None + [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) + + ##################################################### + ### Call Fortran routines for vertical interpolations + if opt.interp is not None: + if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 + ### winds or no winds + if opt.winds : zefields = 'uvmet' + else : zefields = '' + ### var or no var + #if opt.var is None : pass + if zefields == '' : zefields = listvar + else : zefields = zefields + "," + listvar + if opt.var2 is not None : zefields = zefields + "," + opt.var2 + print zefields + zefile = api_onelevel ( path_to_input = '', \ + input_name = zefile, \ + fields = zefields, \ + interp_method = opt.interp, \ + onelevel = zelevel, \ + nocall = opt.nocall ) + print zefile + zelevel = 0 ## so that zelevel could play again the role of nvert + + + ############# + ### Main call + name = planetoplot (zenamefiles,int(zelevel),\ + proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\ + numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ + addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\ + tile=opt.tile,zoom=opt.zoom,display=opt.display,\ + itstep=opt.itstep,hole=opt.hole,save=opt.save,\ + anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first,\ + mult=opt.mult,zetitle=opt.zetitle,\ + slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime) + print 'Done: '+name + system("rm -f to_be_erased") + + ######################################################### + ### Generate a .sh file with the used command saved in it + command = "" + for arg in sys.argv: command = command + arg + ' ' + name = 'pycommand' + f = open(name+'.sh', 'w') + f.write(command) Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/little.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/little.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/little.py (revision 349) @@ -0,0 +1,19 @@ +#! /usr/bin/env python + +import matplotlib.pyplot as plt +import numpy as np + +vvv = np.array([ 8.9, 7.3, 3.9, 2.0]) +rrr = np.array([ 2.0, 5.0, 8.0,11.0]) + +logv = np.log(vvv) +logr = np.log(rrr) + +vvv1 = 1./rrr +vvv2 = 1./np.sqrt(rrr) + +logv1 = 3. + np.log(vvv1) +logv2 = 3. + np.log(vvv2) + +plt.plot(logv,logr,'bo',logv1,logr,'r',logv2,logr,'y') +plt.show() Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py (revision 349) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/meso.py (revision 349) @@ -0,0 +1,152 @@ +#!/usr/bin/env python + +### A. Spiga + +########################################################################################### +########################################################################################### +### 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 api_wrapper import api_onelevel + from netCDF4 import Dataset + from myplot import getlschar, separatenames, readslices, adjust_length + from os import system + from planetoplot import planetoplot + import numpy as np + + ############################# + ### Get options and variables + parser = OptionParser() + parser.add_option('-f', '--file', action='append',dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)') + parser.add_option('-l', '--level', action='store',dest='nvert', type="float", default=0, help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)') + parser.add_option('-p', '--proj', action='store',dest='proj', type="string", default=None, help='projection') + parser.add_option('-b', '--back', action='store',dest='back', type="string", default=None, help='background image (def: None)') + parser.add_option('-t', '--target', action='store',dest='target', type="string", default=None, help='destination folder') + parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') + parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') + parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -*numplot*)') + parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)') + parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') + parser.add_option('-x', '--no-vect',action='store_false',dest='winds', default=True, help='no wind vectors') + parser.add_option('-m', '--min', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)') + parser.add_option('-M', '--max', action='append',dest='vmax', type="float", default=None, help='bounding maximum value (append)') + parser.add_option('-T', '--tiled', action='store_true',dest='tile', default=False, help='draw a tiled plot (no blank zone)') + parser.add_option('-z', '--zoom', action='store',dest='zoom', type="float", default=None, help='zoom factor in %') + parser.add_option('-N', '--no-api', action='store_true',dest='nocall', default=False, help='do not recreate api file') + parser.add_option('-d', '--display',action='store_false',dest='display', default=True, help='do not pop up created images') + parser.add_option('-e', '--itstep', action='store',dest='itstep', type="int", default=None, help='stride time (def=4)') + parser.add_option('-H', '--hole', action='store_true',dest='hole', default=False, help='holes above max and below min') + parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)') + parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly', default=False, help='compute and plot relative anomaly in %') + parser.add_option('-w', '--with', action='store',dest='var2', type="string", default=None, help='variable contoured') + parser.add_option('--div', action='store',dest='ndiv', type="int", default=10, help='number of divisions in colorbar (def: 10)') + parser.add_option('-F', '--first', action='store',dest='first', type="int", default=1, help='first subscript to plot (def: 1)') + parser.add_option('--mult', action='store',dest='mult', type="float", default=1., help='a multiplicative factor to plotted field') + parser.add_option('--title', action='store',dest='zetitle', type="string", default="fill",help='customize the whole title') + #parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables') + + ############# T.N. changes + #parser.add_option('-o','--operation',action='store',dest='operation',type="string", default=None ,help='matrix of operations between files (for now see code, sorry)') + parser.add_option('--lat', action='append',dest='slat',type="string", default=None, help='slices along latitude. One value, or two values separated by comas for averaging') + parser.add_option('--lon', action='append',dest='slon', type="string", default=None, help='slices along longitude. One value, or two values separated by comas for averaging') + parser.add_option('--vert', action='append',dest='svert',type="string", default=None, help='slices along vertical axis. One value, or two values separated by comas for averaging') # quelles coordonnees ? + parser.add_option('--time', action='append',dest='stime',type="string", default=None, help='slices along time. One value, or two values separated by comas for averaging') # quelles coordonnees ? + + + (opt,args) = parser.parse_args() + if opt.namefile is None: + print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h" + exit() + if opt.var is None and opt.anomaly is True: + print "Cannot ask to compute anomaly if no variable is set" + exit() + print "Options:", opt + + listvar = '' + if opt.var is None: + zerange = [-999999] + else: + zelen = len(opt.var) + zerange = range(zelen) + #if zelen == 1: listvar = opt.var[0] + ',' + #else : + for jjj in zerange: listvar += opt.var[jjj] + ',' + listvar = listvar[0:len(listvar)-1] + vmintab = adjust_length (opt.vmin, zelen) + vmaxtab = adjust_length (opt.vmax, zelen) + + print "namefile, length", opt.namefile, len(opt.namefile) + + zeslat = readslices(opt.slat) + zeslon = readslices(opt.slon) + zesvert = readslices(opt.svert) + zestime = readslices(opt.stime) + print "slat,zeslat", opt.slat, zeslat + print "slon,zeslon", opt.slon, zeslon + print "svert,zesvert", opt.svert, zesvert + print "stime,zestime", opt.stime, zestime + + for i in range(len(opt.namefile)): + + zefile = opt.namefile[i] + print zefile + zelevel = opt.nvert + stralt = None + [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) + + ##################################################### + ### Call Fortran routines for vertical interpolations + if opt.interp is not None: + if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 + ### winds or no winds + if opt.winds : zefields = 'uvmet' + else : zefields = '' + ### var or no var + #if opt.var is None : pass + if zefields == '' : zefields = listvar + else : zefields = zefields + "," + listvar + if opt.var2 is not None : zefields = zefields + "," + opt.var2 + print zefields + zefile = api_onelevel ( path_to_input = '', \ + input_name = zefile, \ + fields = zefields, \ + interp_method = opt.interp, \ + onelevel = zelevel, \ + nocall = opt.nocall ) + print zefile + zelevel = 0 ## so that zelevel could play again the role of nvert + + if opt.var is None: zerange = [-999999] + else: zerange = range(zelen) + for jjj in zerange: + if jjj == -999999: + argvar = None + zevmin = None + zevmax = None + else: + argvar = opt.var[jjj] + if vmintab[jjj] != -999999: zevmin = vmintab[jjj] + else: zevmin = None + if vmaxtab[jjj] != -999999: zevmax = vmaxtab[jjj] + else: zevmax = None + ############# + ### Main call + name = planetoplot (zefile,int(zelevel),\ + proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\ + numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ + addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\ + tile=opt.tile,zoom=opt.zoom,display=opt.display,\ + itstep=opt.itstep,hole=opt.hole,save=opt.save,\ + anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first,\ + mult=opt.mult,zetitle=opt.zetitle,\ + slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime) + print 'Done: '+name + system("rm -f to_be_erased") + + ######################################################### + ### Generate a .sh file with the used command saved in it + command = "" + for arg in sys.argv: command = command + arg + ' ' + f = open(name+'.sh', 'w') + f.write(command) Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py (revision 345) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py (revision 349) @@ -1,5 +1,6 @@ def min (field,axis=None): import numpy as np - return np.array(field).min(axis=axis) + if field is None: return None + else: return np.array(field).min(axis=axis) def max (field,axis=None): @@ -10,5 +11,6 @@ def mean (field,axis=None): import numpy as np - return np.array(field).mean(axis=axis) + if field is None: return None + else: return np.array(field).mean(axis=axis) def deg (): Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 345) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 349) @@ -5,4 +5,17 @@ exit() return + +## Author: AS +def adjust_length (tab, zelen): + from numpy import ones + if tab is None: + outtab = ones(zelen) * -999999 + else: + if zelen != len(tab): + print "not enough or too much values... setting same values all variables" + outtab = ones(zelen) * tab[0] + else: + outtab = tab + return outtab ## Author: AS @@ -48,7 +61,8 @@ ### it would be actually better to name d4 d3 d2 d1 as t z y x import numpy as np - from mymath import max + from mymath import max,mean dimension = np.array(input).ndim shape = np.array(input).shape + #print 'd1,d2,d3,d4: ',d1,d2,d3,d4 print 'dim,shape: ',dimension,shape output = input @@ -71,10 +85,10 @@ elif max(d1) >= shape[2]: error = True elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,d2,d1] - elif d4 is not None and d2 is not None: output=np.mean(input[d4,:,:],axis=0); output=np.mean(output[d2,:],axis=0) - elif d4 is not None and d1 is not None: output=np.mean(input[d4,:,:],0); output=np.mean(output[:,d1],1) - elif d2 is not None and d1 is not None: output=np.mean(input[:,d2,:],axis=1); output=np.mean(output[:,d1],axis=1) - elif d1 is not None: output = np.mean(input[:,:,d1],axis=2) - elif d2 is not None: output = np.mean(input[:,d2,:],axis=1) - elif d4 is not None: output = np.mean(input[d4,:,:],axis=0) + elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0) + elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1) + elif d2 is not None and d1 is not None: output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1) + elif d1 is not None: output = mean(input[:,:,d1],axis=2) + elif d2 is not None: output = mean(input[:,d2,:],axis=1) + elif d4 is not None: output = mean(input[d4,:,:],axis=0) elif dimension == 4: if max(d4) >= shape[0]: error = True @@ -87,10 +101,10 @@ elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,:,d2,d1] elif d3 is not None and d2 is not None and d1 is not None: output = input[:,d3,d2,d1] - elif d4 is not None and d3 is not None: output = np.mean(input[d4,:,:,:],0); output = np.mean(output[d3,:,:],0) - elif d4 is not None and d2 is not None: output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,d2,:],1) - elif d4 is not None and d1 is not None: output = np.mean(input[d4,:,:,:],0); output = np.mean(output[:,:,d1],2) - elif d3 is not None and d2 is not None: output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,d2,:],1) - elif d3 is not None and d1 is not None: output = np.mean(input[:,d3,:,:],1); output = np.mean(output[:,:,d1],0) - elif d2 is not None and d1 is not None: output = np.mean(input[:,:,d2,:],2); output = np.mean(output[:,:,d1],2) + elif d4 is not None and d3 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[d3,:,:],axis=0) + elif d4 is not None and d2 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,d2,:],axis=1) + elif d4 is not None and d1 is not None: output = mean(input[d4,:,:,:],axis=0); output = mean(output[:,:,d1],axis=2) + elif d3 is not None and d2 is not None: output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,d2,:],axis=1) + elif d3 is not None and d1 is not None: output = mean(input[:,d3,:,:],axis=1); output = mean(output[:,:,d1],axis=0) + elif d2 is not None and d1 is not None: output = mean(input[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2) elif d1 is not None: output = input[:,:,:,d1] elif d2 is not None: output = input[:,:,d2,:] @@ -759,4 +773,6 @@ # output : all indexes to be taken into account for reducing field import numpy as np + if ( np.array(axis).ndim == 2): + axis = axis[:,0] if saxis is None: zeindex = None @@ -820,2 +836,34 @@ #print "define_axis", x, y return what_I_plot,x,y + +# Author: TN + AS +def determineplot(slon, slat, svert, stime): + nlon = 1 # number of longitudinal slices -- 1 is None + nlat = 1 + nvert = 1 + ntime = 1 + nslices = 1 + if slon is not None: + nslices = nslices*len(slon) + nlon = len(slon) + if slat is not None: + nslices = nslices*len(slat) + nlat = len(slat) + if svert is not None: + nslices = nslices*len(svert) + nvert = len(svert) + if stime is not None: + nslices = nslices*len(stime) + ntime = len(stime) + #else: + # nslices = 2 + + mapmode = 0 + if slon is None and slat is None: + mapmode = 1 # in this case we plot a map, with the given projection + #elif proj is not None: + # print "WARNING: you specified a", proj,\ + # "projection but asked for slices", slon,"in longitude and",slat,"in latitude" + print "mapmode: ", mapmode + + return nlon, nlat, nvert, ntime, mapmode, nslices Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py (revision 345) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/planetoplot.py (revision 349) @@ -1,12 +1,9 @@ -#!/usr/bin/env python - -### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011 -### T.Navarro -- LMD -- 10/2011 -### Thanks to A. Colaitis for the parser trick - - -#################################### -#################################### -### The main program to plot vectors +####################### +##### PLANETOPLOT ##### +####################### + +### A. Spiga -- LMD -- 06~09/2011 -- General building and mapping capabilities +### T. Navarro -- LMD -- 10~11/2011 -- Improved use for GCM and added sections + 1Dplot capabilities + def planetoplot (namefiles,\ vertmode,\ @@ -49,5 +46,5 @@ fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\ zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\ - getname,localtime,polarinterv,getsindex,define_axis + getname,localtime,polarinterv,getsindex,define_axis,determineplot from mymath import deg,max,min,mean from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel @@ -56,42 +53,16 @@ from numpy.core.defchararray import find - #whichmode = "" - nlon = 1 # number of longitudinal slices -- 1 is None - nlat = 1 - nvert = 1 - ntime = 1 - nslices = 1 - if slon is not None: - nslices = nslices*len(slon) - nlon = len(slon) - if slat is not None: - nslices = nslices*len(slat) - nlat = len(slat) - if svert is not None: - nslices = nslices*len(svert) - nvert = len(svert) - if stime is not None: - nslices = nslices*len(stime) - ntime = len(stime) - #else: - # nslices = 2 - numplot = len(namefiles)*nslices + ################################ + ### Which plot needs to be done? + nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime) + if mapmode == 0: winds=False + if not isinstance(namefiles, np.ndarray): namefiles = [namefiles] + zelen = len(namefiles) + numplot = zelen*nslices + print "len(namefiles), nslices, numplot: ", zelen, nslices, numplot + all_var = [[]]*zelen + all_var2 = [[]]*zelen + all_title = [[]]*zelen - mapmode = 0 - if slon is None and slat is None: - mapmode = 1 # in this case we plot a map, with the given projection - elif proj is not None: - print "WARNING: you specified a", proj,\ - "projection but asked for slices", slon,"in longitude and",slat,"in latitude" - print "mapmode", mapmode - - - all_var = [[]]*len(namefiles) - all_var2 = [[]]*len(namefiles) - all_title = [[]]*len(namefiles) - - print "len(namefiles), nslices", len(namefiles), nslices - print "numplot", numplot - ######################### ### Loop over the files initially separated by comas to be plotted on the same figure @@ -106,28 +77,36 @@ ################################## ### Initial checks and definitions - typefile = whatkindfile(nc) ## TYPEFILE + ### ... TYPEFILE + typefile = whatkindfile(nc) if firstfile: typefile0 = typefile elif typefile != typefile0: - errormess("Not the same files !", [typefile0, typefile]) - if var not in nc.variables: var = False ## VAR - if winds: ## WINDS + errormess("Not the same kind of files !", [typefile0, typefile]) + ### ... VAR + if var not in nc.variables: var = False + ### ... WINDS + if winds: [uchar,vchar,metwind] = getwinddef(nc) if uchar == 'not found': winds = False if not var and not winds: errormess("please set at least winds or var",printvar=nc.variables) - [lon2d,lat2d] = getcoorddef(nc) ## COORDINATES, could be moved below - if proj == None: proj = getproj(nc) ## PROJECTION - - lat = nc.variables["latitude"][:] - lon = nc.variables["longitude"][:] - time = nc.variables["Time"][:] - vert = nc.variables["altitude"][:] - if firstfile: - lat0 = lat - elif len(lat0) != len(lat): - errormess("Not the same latitude lengths !", [len(lat0), len(lat)]) - elif sum((lat == lat0) == False) != 0: - errormess("Not the same latitudes !", [lat,lat0]) - # Faire d'autre checks sur les compatibilites entre fichiers!! + ### ... COORDINATES, could be moved below + [lon2d,lat2d] = getcoorddef(nc) + ### ... PROJECTION + if proj == None: proj = getproj(nc) + +########################################################## + if typefile == "gcm": + lat = nc.variables["latitude"][:] + lon = nc.variables["longitude"][:] + time = nc.variables["Time"][:] + vert = nc.variables["altitude"][:] + #if firstfile: + # lat0 = lat + #elif len(lat0) != len(lat): + # errormess("Not the same latitude lengths !", [len(lat0), len(lat)]) + #elif sum((lat == lat0) == False) != 0: + # errormess("Not the same latitudes !", [lat,lat0]) + ## Faire d'autre checks sur les compatibilites entre fichiers!! +########################################################## if firstfile: @@ -143,10 +122,7 @@ ### Name for title and graphics save file basename = getname(var=var,winds=winds,anomaly=anomaly) - basename = basename + getstralt(nc,nvert) ## can be moved elsewhere for a more generic routine + basename = basename + getstralt(nc,vertmode) ## can be moved elsewhere for a more generic routine - print "var", var - print "var2", var2 - #print all_var - #print nc + print "var, var2: ", var, var2 if var: all_var[k] = getfield(nc,var) if var2: all_var2[k] = getfield(nc,var2) @@ -173,9 +149,4 @@ elif numplot <= 0: itstep = 1 - ### Map projection - if mapmode == 1: - m = define_proj(proj,wlon,wlat,back=back) - x, y = m(lon2d, lat2d) - #for nplot in range(numplot): while error is False: @@ -201,25 +172,44 @@ else: found_lct = True - - - ## get all indexes to be taken into account for this subplot and then reduce field - ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice - #print "(nplot-1)%nlon", (nplot-1)%nlon - indexlon = getsindex(slon,(nplot-1)%nlon,lon) - #print "indexlon", indexlon - indexlat = getsindex(slat,((nplot-1)//nlon)%nlat,lat) - indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) - indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) - index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles) - #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime) + + ### Map projection + if mapmode == 1: + m = define_proj(proj,wlon,wlat,back=back) + x, y = m(lon2d, lat2d) + +#################################################################### + if typefile in ['mesoapi','meso']: + indexlon = None + indexlat = None + indexvert = vertmode + indextime = itime + nlon = 1 + nlat = 1 + nvert = 1 + ntime = 1 + index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles) + print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime) + elif typefile in ['gcm']: + ## get all indexes to be taken into account for this subplot and then reduce field + ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice + #print "(nplot-1)%nlon", (nplot-1)%nlon + indexlon = getsindex(slon,(nplot-1)%nlon,lon) + print "indexlon", indexlon + indexlat = getsindex(slat,((nplot-1)//nlon)%nlat,lat) + indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) + indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) + index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%len(namefiles) + #print "index lon,lat,vert,time", max(indexlon), max(indexlat), max(indexvert), max(indextime) +#################################################################### #### Contour plot if var2: what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert ) - what_I_plot = what_I_plot*mult + #what_I_plot = what_I_plot*mult if not error: - if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_contour,6) + if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) zevmin, zevmax = calculate_bounds(what_I_plot) - zelevels = np.linspace(zevmin,zevmax,num=ndiv) + zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #20) + if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.) if mapmode == 0: x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ @@ -227,6 +217,6 @@ ### If we plot a 2-D field if len(what_I_plot.shape) is 2: - cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #colors='w' )# , alpha=0.5) - clabel(cs,fmt = '%1.2e') + cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5) + if typefile in ['gcm']: clabel(cs,fmt = '%1.2e') ### If we plot a 1-D field elif len(what_I_plot.shape) is 1: @@ -256,5 +246,5 @@ if len(what_I_plot.shape) is 2: if mapmode == 0: - what_I_plot,x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ + what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\ indextime,what_I_plot, len(all_var[index_f].shape),vertmode) if not tile: @@ -263,5 +253,5 @@ print "what_I_plot.shape", what_I_plot.shape print "x.shape, y.shape", x.shape, y.shape - zelevels = np.linspace(zevmin,zevmax,num=ndiv) + zelevels = np.linspace(zevmin,zevmax,num=ndiv+1) #contourf(what_I_plot, zelevels, cmap = palette ) contourf(x,y,what_I_plot, zelevels, cmap = palette ) @@ -274,5 +264,5 @@ #colorbar() colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar.upper()),\ - ticks=np.linspace(zevmin,zevmax,min([ndiv,20])),\ + ticks=np.linspace(zevmin,zevmax,min([ndiv+1,20])),\ extend='neither',spacing='proportional') # both min max neither @@ -295,11 +285,32 @@ else: errormess("There is an error in reducing field !") - + + ### Vector plot --- a adapter + if winds: + vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert ) + vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert ) + #what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert ) + if not error: + if typefile in ['mesoapi','meso']: + [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar), dumpbdy(vecy,6,stag=vchar)] + key = True + elif typefile in ['gcm']: + key = False + if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d) + if var == False: colorvec = definecolorvec(back) + else: colorvec = definecolorvec(colorb) + vectorfield(vecx, vecy,\ + x, y, stride=stride, csmooth=2,\ + #scale=15., factor=300., color=colorvec, key=key) + scale=20., factor=250., color=colorvec, key=key) + #200. ## or csmooth=stride ### Next subplot - plottitle = basename+' '+namefiles[index_f] if typefile in ['mesoapi','meso']: + plottitle = basename if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) else: plottitle = plottitle + "_LT"+str(ltst) + else: + plottitle = basename+' '+namefiles[index_f] if mult != 1: plottitle = str(mult) + "*" + plottitle if zetitle != "fill": plottitle = zetitle @@ -314,5 +325,5 @@ ptitle( plottitle ) itime += itstep - if nplot == numplot: + if nplot >= numplot: error = True nplot += 1 @@ -356,173 +367,2 @@ ### Now the end return zeplot - -############################## -### A specific stuff for below -def adjust_length (tab, zelen): - from numpy import ones - if tab is None: - outtab = ones(zelen) * -999999 - else: - if zelen != len(tab): - print "not enough or too much values... setting same values all variables" - outtab = ones(zelen) * tab[0] - else: - outtab = tab - return outtab - -########################################################################################### -########################################################################################### -### 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 api_wrapper import api_onelevel - from netCDF4 import Dataset - from myplot import getlschar, separatenames, readslices - from os import system - import numpy as np - - ############################# - ### Get options and variables - parser = OptionParser() - parser.add_option('-f', '--file', action='append',dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append). Plot files separated by comas in the same figure') - parser.add_option('-l', '--level', action='store',dest='nvert', type="float", default=0, help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)') - parser.add_option('-p', '--proj', action='store',dest='proj', type="string", default=None, help='projection') - parser.add_option('-b', '--back', action='store',dest='back', type="string", default=None, help='background image (def: None)') - parser.add_option('-t', '--target', action='store',dest='target', type="string", default=None, help='destination folder') - parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') - parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') - parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -*numplot*)') - parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)') - parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') - parser.add_option('-x', '--no-vect',action='store_false',dest='winds', default=True, help='no wind vectors') - parser.add_option('-m', '--min', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)') - parser.add_option('-M', '--max', action='append',dest='vmax', type="float", default=None, help='bounding maximum value (append)') - parser.add_option('-T', '--tiled', action='store_true',dest='tile', default=False, help='draw a tiled plot (no blank zone)') - parser.add_option('-z', '--zoom', action='store',dest='zoom', type="float", default=None, help='zoom factor in %') - parser.add_option('-N', '--no-api', action='store_true',dest='nocall', default=False, help='do not recreate api file') - parser.add_option('-d', '--display',action='store_false',dest='display', default=True, help='do not pop up created images') - parser.add_option('-e', '--itstep', action='store',dest='itstep', type="int", default=None, help='stride time (def=4)') - parser.add_option('-H', '--hole', action='store_true',dest='hole', default=False, help='holes above max and below min') - parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)') - parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly', default=False, help='compute and plot relative anomaly in %') - parser.add_option('-w', '--with', action='store',dest='var2', type="string", default=None, help='variable contoured') - parser.add_option('--div', action='store',dest='ndiv', type="int", default=10, help='number of divisions in colorbar (def: 10)') - parser.add_option('-F', '--first', action='store',dest='first', type="int", default=1, help='first subscript to plot (def: 1)') - parser.add_option('--mult', action='store',dest='mult', type="float", default=1., help='a multiplicative factor to plotted field') - parser.add_option('--title', action='store',dest='zetitle', type="string", default="fill",help='customize the whole title') - #parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables') - - ############# T.N. changes - #parser.add_option('-o','--operation',action='store',dest='operation',type="string", default=None ,help='matrix of operations between files (for now see code, sorry)') - parser.add_option('--lat', action='append',dest='slat',type="string", default=None, help='slices along latitude. One value, or two values separated by comas for averaging') - parser.add_option('--lon', action='append',dest='slon', type="string", default=None, help='slices along longitude. One value, or two values separated by comas for averaging') - parser.add_option('--vert', action='append',dest='svert',type="string", default=None, help='slices along vertical axis. One value, or two values separated by comas for averaging') # quelles coordonnees ? - parser.add_option('--time', action='append',dest='stime',type="string", default=None, help='slices along time. One value, or two values separated by comas for averaging') # quelles coordonnees ? - - (opt,args) = parser.parse_args() - if opt.namefile is None: - print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h" - exit() - if opt.var is None and opt.anomaly is True: - print "Cannot ask to compute anomaly if no variable is set" - exit() - print "Options:", opt - - #listvar = '' - #if opt.var is None: - # zerange = [-999999] - #else: - # zelen = len(opt.var) - # zerange = range(zelen) - # #if zelen == 1: listvar = opt.var[0] + ',' - # #else : - # for jjj in zerange: listvar += opt.var[jjj] + ',' - # listvar = listvar[0:len(listvar)-1] - # vmintab = adjust_length (opt.vmin, zelen) - # vmaxtab = adjust_length (opt.vmax, zelen) - - print "namefile, length", opt.namefile, len(opt.namefile) - - zeslat = readslices(opt.slat) - zeslon = readslices(opt.slon) - zesvert = readslices(opt.svert) - zestime = readslices(opt.stime) - print "slat,zeslat", opt.slat, zeslat - print "slon,zeslon", opt.slon, zeslon - print "svert,zesvert", opt.svert, zesvert - print "stime,zestime", opt.stime, zestime - - - for i in range(len(opt.namefile)): - for j in range(len(opt.var)): - - zenamefiles = separatenames(opt.namefile[i]) - print "zenamefiles", zenamefiles - - if opt.vmin is not None : zevmin = opt.vmin[min(i,len(opt.vmin)-1)] - else: zevmin = None - if opt.vmax is not None : zevmax = opt.vmax[min(i,len(opt.vmax)-1)] - else: zevmax = None - print "vmin, zevmin", opt.vmin, zevmin - print "vmax, zevmax", opt.vmax, zevmax - - zevar = separatenames(opt.var[j]) - zevar = zevar[0] - print "var, zevar", opt.var, zevar - - #checkcoherence(len(zenamefiles),len(opt.slat),len(opt.slon),len(opt.stime)) - - zefile = zenamefiles[0] - - zelevel = opt.nvert - stralt = None - [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) - #print "lschar ",lschar - #print "zehour ",zehour - #print "zehourin ",zehourin - - ##################################################### - ### Call Fortran routines for vertical interpolations - if opt.interp is not None: - if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 - ### winds or no winds - if opt.winds : zefields = 'uvmet' - else : zefields = '' - ### var or no var - #if opt.var is None : pass - if zefields == '' : zefields = listvar - else : zefields = zefields + "," + listvar - if opt.var2 is not None : zefields = zefields + "," + opt.var2 - print zefields - zefile = api_onelevel ( path_to_input = '', \ - input_name = zefile, \ - fields = zefields, \ - interp_method = opt.interp, \ - onelevel = zelevel, \ - nocall = opt.nocall ) - print zefile - zelevel = 0 ## so that zelevel could play again the role of nvert - - - ############# - ### Main call - name = planetoplot (zenamefiles,int(zelevel),\ - proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=zevar,\ - numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ - addchar=lschar,interv=[zehour,zehourin],vmin=zevmin,vmax=zevmax,\ - tile=opt.tile,zoom=opt.zoom,display=opt.display,\ - itstep=opt.itstep,hole=opt.hole,save=opt.save,\ - anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first,\ - mult=opt.mult,zetitle=opt.zetitle,\ - slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime) - print 'Done: '+name - system("rm -f to_be_erased") - - ######################################################### - ### Generate a .sh file with the used command saved in it - command = "" - for arg in sys.argv: command = command + arg + ' ' - name = 'pycommand' - f = open(name+'.sh', 'w') - f.write(command)