Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py (revision 240) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py (revision 241) @@ -71,45 +71,20 @@ return output, error -def latinterv (area): - if area == "Europe": - wlat = [20.,80.] - wlon = [-50.,50.] - elif area == "Central_America": - wlat = [-10.,40.] - wlon = [230.,300.] - elif area == "Africa": - wlat = [-20.,50.] - wlon = [-50.,50.] - elif area == "Whole": - wlat = [-90.,90.] - wlon = [-180.,180.] - elif area == "Southern_Hemisphere": - wlat = [-90.,60.] - wlon = [-180.,180.] - elif area == "Northern_Hemisphere": - wlat = [-60.,90.] - wlon = [-180.,180.] - elif area == "Tharsis": - wlat = [-30.,60.] - wlon = [-170.,-10.] - elif area == "Whole_No_High": - wlat = [-60.,60.] - wlon = [-180.,180.] - elif area == "Chryse": - wlat = [-60.,60.] - wlon = [-60.,60.] - elif area == "North_Pole": - wlat = [50.,90.] - wlon = [-180.,180.] - elif area == "Close_North_Pole": - wlat = [75.,90.] - wlon = [-180.,180.] - elif area == "South_Pole": - wlat = [-90.,-50.] - wlon = [-180.,180.] - elif area == "Close_South_Pole": - wlat = [-90.,-75.] - wlon = [-180.,180.] - return wlon,wlat +#def latinterv (area): +# if area == "Europe": [wlat,wlon] = [[ 20., 80.],[- 50., 50.]] +# elif area == "Central_America": [wlat,wlon] = [[-10., 40.],[ 230., 300.]] +# elif area == "Africa": [wlat,wlon] = [[-20., 50.],[- 50., 50.]] +# elif area == "Whole": [wlat,wlon] = [[-90., 90.],[-180.,-180.]] +# elif area == "Southern_Hemisphere": [wlat,wlon] = [[-90., 60.],[-180.,-180.]] +# elif area == "Northern_Hemisphere": [wlat,wlon] = [[-60., 90.],[-180.,-180.]] +# elif area == "Tharsis": [wlat,wlon] = [[-30., 60.],[-170.,- 10.]] +# elif area == "Whole_No_High": [wlat,wlon] = [[-60., 60.],[-180., 180.]] +# elif area == "Chryse": [wlat,wlon] = [[-60., 60.],[- 60., 60.]] +# elif area == "North_Pole": [wlat,wlon] = [[ 50., 90.],[-180., 180.]] +# elif area == "Close_North_Pole": [wlat,wlon] = [[ 75., 90.],[-180., 180.]] +# elif area == "Far_South_Pole": [wlat,wlon] = [[-90.,-40.],[-180., 180.]] +# elif area == "South_Pole": [wlat,wlon] = [[-90.,-50.],[-180., 180.]] +# elif area == "Close_South_Pole": [wlat,wlon] = [[-90.,-75.],[-180., 180.]] +# return wlon,wlat def definesubplot ( numplot, fig ): @@ -142,5 +117,5 @@ elif numplot == 1: sub = 99999 - elif numplot < 0: + elif numplot <= 0: sub = 99999 else: @@ -182,5 +157,5 @@ zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) if zeday < 0: lschar="" ## might have crossed a month... fix soon - else: lschar="_Ls"+str( int( sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) ) + else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. ) ### zetime2 = nc.variables['Times'][1] @@ -472,4 +447,11 @@ return what_I_plot +def nolow(what_I_plot): + from mymath import max,min + lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot))) + print "on vire en dessous de ", lim + what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 + return what_I_plot + def zoomset (wlon,wlat,zoom): dlon = abs(wlon[1]-wlon[0])/2. @@ -492,4 +474,5 @@ "TAU_ICE": "%.2f",\ "anomaly": "%.1f",\ + "W": "%.1f",\ } if whichvar not in fmtvar: @@ -509,7 +492,10 @@ "ICETOT": "YlGnBu",\ "TAU_ICE": "Blues",\ + "W": "jet",\ "anomaly": "RdBu_r",\ } +#W --> spectral ou jet #spectral BrBG RdBu_r + print "predefined colorbars" if whichone not in whichcolorb: whichone = "def" @@ -567,2 +553,23 @@ return whichlink +def latinterv (area="Whole"): + list = { \ + "Europe": [[ 20., 80.],[- 50., 50.]],\ + "Central_America": [[-10., 40.],[ 230., 300.]],\ + "Africa": [[-20., 50.],[- 50., 50.]],\ + "Whole": [[-90., 90.],[-180.,-180.]],\ + "Southern_Hemisphere": [[-90., 60.],[-180.,-180.]],\ + "Northern_Hemisphere": [[-60., 90.],[-180.,-180.]],\ + "Tharsis": [[-30., 60.],[-170.,- 10.]],\ + "Whole_No_High": [[-60., 60.],[-180., 180.]],\ + "Chryse": [[-60., 60.],[- 60., 60.]],\ + "North_Pole": [[ 50., 90.],[-180., 180.]],\ + "Close_North_Pole": [[ 75., 90.],[-180., 180.]],\ + "Far_South_Pole": [[-90.,-40.],[-180., 180.]],\ + "South_Pole": [[-90.,-50.],[-180., 180.]],\ + "Close_South_Pole": [[-90.,-75.],[-180., 180.]],\ + } + if area not in list: area = "Whole" + [olat,olon] = list[area] + return olon,olat + Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py (revision 240) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py (revision 241) @@ -5,7 +5,9 @@ def domain (namefile,proj=None,back="vishires",target=None): from netCDF4 import Dataset - from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv - from matplotlib.pyplot import contourf,rcParams + from myplot import getcoord2d,define_proj,makeplotres,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv + from mymath import max,min + from matplotlib.pyplot import contourf,rcParams,pcolor from numpy.core.defchararray import find + from numpy import arange ### nc = Dataset(namefile) @@ -23,6 +25,6 @@ zeplot = getprefix(nc) + "domain" ### - lon2d = dumpbdy(lon2d) - lat2d = dumpbdy(lat2d) + lon2d = dumpbdy(lon2d,5) + lat2d = dumpbdy(lat2d,5) if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) @@ -32,11 +34,18 @@ x, y = m(lon2d, lat2d) ### - what_I_plot = dumpbdy( nc.variables[var][0,:,:] ) - contourf(x, y, what_I_plot, 40) + what_I_plot = dumpbdy(nc.variables[var][0,:,:], 5) + #levinterv = 250. + #zelevels = arange(min(what_I_plot)-levinterv,max(what_I_plot)+levinterv,levinterv) + zelevels = 30 + contourf(x, y, what_I_plot, zelevels) + #pcolor(x,y,what_I_plot) ## on voit trop les lignes ! ### if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot else: zeplot = target + "/" + zeplot ### - makeplotpng(zeplot,pad_inches_value=0.35) + pad_inches_value = 0.35 + makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature + makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False) + #makeplotpng(zeplot,pad_inches_value=0.35) #rcParams['savefig.facecolor'] = 'black' #makeplotpng(zeplot+"b",pad_inches_value=0.35) Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (revision 240) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (revision 241) @@ -16,5 +16,5 @@ numplot=2,\ var=None,\ - colorb=True,\ + colorb="def",\ winds=True,\ addchar=None,\ @@ -28,5 +28,6 @@ hole=False,\ save="gui",\ - anomaly=False): + anomaly=False,\ + var2=None): #################################################################################################################### @@ -38,5 +39,5 @@ from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\ fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\ - zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth + zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow from mymath import deg,max,min,mean from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show @@ -45,5 +46,5 @@ from numpy.core.defchararray import find - #rcParams['backend'] = 'PS' + #if save == 'eps': rcParams['backend'] = 'PS' ###################### @@ -64,5 +65,5 @@ ### Define plot boundaries if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") - elif proj == "spstere": [wlon,wlat] = latinterv("South_Pole") + elif proj == "spstere": [wlon,wlat] = latinterv("Far_South_Pole") elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) else: [wlon,wlat] = simplinterv(lon2d,lat2d) @@ -77,4 +78,5 @@ print nc.variables errormess("please set at least winds or var") + if anomaly: basename = 'd' + basename basename = basename + getstralt(nc,nvert) ## can be moved elsewhere for a more generic routine @@ -90,5 +92,6 @@ nplot = 1 error = False - if itstep is None: itstep = int(24./numplot) + if itstep is None and numplot > 0: itstep = int(24./numplot) + else: itstep = 1 while error is False: @@ -109,5 +112,5 @@ ### If only one local time is requested (numplot < 0) elif numplot <= 0: - if int(ltst) + numplot != 0: + if int(ltst) + numplot != 0: itime += 1 if found_lct is True: break ## because it means LT was found at previous iteration @@ -121,4 +124,14 @@ #### Contour plot + if var2: + what_I_contour, error = reducefield( getfield(nc,var2), d4=itime, d3=nvert ) + if not error: + if typefile in ['mesoapi','meso']: what_I_contour = dumpbdy(what_I_contour,6) + zevmin, zevmax = calculate_bounds(what_I_contour) + zelevels = np.linspace(zevmin,zevmax,num=20) + if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.) + contour( x, y, what_I_contour, zelevels, colors='k', linewidths = 0.33 ) #colors='w' )# , alpha=0.5) + + #### Shaded plot if var: what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d3=nvert ) @@ -132,15 +145,15 @@ if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) - if colorb is True: colorb = defcolorb(fvar) - palette = get_cmap(name=colorb) + if colorb in ["def","nobar"]: palette = get_cmap(name=defcolorb(fvar)) + else: palette = get_cmap(name=colorb) if not tile: if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax) - zelevels = np.linspace(zevmin,zevmax) + zelevels = np.linspace(zevmin,zevmax) #,num=20) contourf( x, y, what_I_plot, zelevels, cmap = palette ) - else: - pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax ) - #putpoints(m,fig) - if var in ['HGT']: pass - elif colorb: + else: + if hole: what_I_plot = nolow(what_I_plot) + pcolor( x, y, what_I_plot, cmap = palette, \ + vmin=zevmin, vmax=zevmax ) + if colorb != 'nobar' or var == 'HGT': ndiv = 10 colorbar(fraction=0.05,pad=0.1,format=fmtvar(fvar),\ @@ -148,5 +161,5 @@ extend='neither',spacing='proportional') # both min max neither - + ### Vector plot if winds: @@ -192,5 +205,5 @@ pad_inches_value = 0.35 if save == 'png': - makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value,erase=True) ## a miniature + makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False) elif save in ['eps','svg','pdf']: @@ -214,5 +227,4 @@ outtab = ones(zelen) * -999999 else: - print zelen, len(tab) if zelen != len(tab): print "not enough or too much values... setting same values all variables" @@ -242,8 +254,8 @@ parser.add_option('-t', action='store',dest='target', type="string", default=None, help='destination folder') parser.add_option('-s', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') - parser.add_option('-v', action='append',dest='var', type="string", default=None, help='variable contoured (append)') + parser.add_option('-v', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') parser.add_option('-n', action='store',dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)') parser.add_option('-i', action='store',dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)') - parser.add_option('-c', action='store',dest='colorb', type="string", default=True, help='change colormap') + parser.add_option('-c', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') parser.add_option('-x', action='store_false',dest='winds', default=True, help='no wind vectors') parser.add_option('-m', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)') @@ -257,4 +269,5 @@ parser.add_option('-S', action='store',dest='save', type="string", default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)') parser.add_option('-a', action='store_true',dest='anomaly', default=False, help='compute and plot relative anomaly in %') + parser.add_option('-w', action='store',dest='var2', type="string", default=None, help='variable contoured') #parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables') (opt,args) = parser.parse_args() @@ -262,4 +275,7 @@ 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 @@ -296,4 +312,5 @@ 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 = '', \ @@ -327,5 +344,5 @@ tile=opt.tile,zoom=opt.zoom,display=opt.display,\ itstep=opt.itstep,hole=opt.hole,save=opt.save,\ - anomaly=opt.anomaly) + anomaly=opt.anomaly,var2=opt.var2) print 'Done: '+name system("rm -f to_be_erased")