Index: trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py =================================================================== --- trunk/MESOSCALE/PLOT/PYTHON/scripts/domain.py (revision 199) +++ trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py (revision 207) @@ -1,19 +1,32 @@ #!/usr/bin/env python -### A. Spiga -- LMD -- 30/05/2011 +### A. Spiga -- LMD -- 30/06/2011 -- slight modif early 07/2011 -def domain (namefile,proj="ortho",back="molabw",target=None,var='HGT'): +def domain (namefile,proj="ortho",back="vishires",target=None,var='HGT'): from netCDF4 import Dataset - from myplot import getcoord2d,define_proj,makeplotpng,simplinterv - from matplotlib.pyplot import contourf + from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy + from matplotlib.pyplot import contourf,rcParams + from numpy.core.defchararray import find + ### nc = Dataset(namefile) + ### [lon2d,lat2d] = getcoord2d(nc) + lon2d = dumpbdy(lon2d) + lat2d = dumpbdy(lat2d) [wlon,wlat] = simplinterv(lon2d,lat2d) + ### m = define_proj(proj,wlon,wlat,back=back) x, y = m(lon2d, lat2d) - contourf(x, y, nc.variables[var][0,:,:], 40) - if not target: zeplot = namefile+".domain" - else: zeplot = target+"/domain" + ### + what_I_plot = dumpbdy( nc.variables[var][0,:,:] ) + contourf(x, y, what_I_plot, 40) + ### + zeplot = getprefix(nc) + "domain" + if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot + else: zeplot = target + "/" + zeplot + ### makeplotpng(zeplot,pad_inches_value=0.35) + #rcParams['savefig.facecolor'] = 'black' + #makeplotpng(zeplot+"b",pad_inches_value=0.35) if __name__ == "__main__": @@ -22,9 +35,9 @@ from optparse import OptionParser parser = OptionParser() - parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='name of WRF file [NEEDED]') - parser.add_option('-p', action='store', dest='proj', type="string", default="ortho", help='projection') - parser.add_option('-b', action='store', dest='back', type="string", default="molabw", help='background') - parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') - parser.add_option('-v', action='store', dest='var', type="string", default='HGT', help='variable contoured') + parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='name of WRF file [NEEDED]') + parser.add_option('-p', action='store', dest='proj', type="string", default="ortho", help='projection') + parser.add_option('-b', action='store', dest='back', type="string", default="vishires", help='background') + parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') + parser.add_option('-v', action='store', dest='var', type="string", default='HGT', help='variable contoured') (opt,args) = parser.parse_args() if opt.namefile is None: Index: trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py =================================================================== --- trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py (revision 199) +++ trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (revision 207) @@ -1,5 +1,5 @@ #!/usr/bin/env python -### A. Spiga -- LMD -- 30/06/2011 to 04/07/2011 +### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011 ### Thanks to A. Colaitis for the parser trick @@ -16,10 +16,12 @@ numplot=4,\ var=None,\ - colorb=None,\ + colorb=True,\ winds=True,\ addchar=None,\ interv=[0,1],\ vmin=None,\ - vmax=None): + vmax=None,\ + tile=False,\ + zoom=None): #################################################################################################################### @@ -29,9 +31,15 @@ ### Load librairies and functions from netCDF4 import Dataset - from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy + from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\ + fmtvar,definecolorvec,getwinds,defcolorb,getprefix from mymath import deg,max,min,mean - from matplotlib.pyplot import contourf, subplot, figure, rcParams, savefig, colorbar, pcolor + from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor from matplotlib.cm import get_cmap import numpy as np + from numpy.core.defchararray import find + + ### + #rcParams['text.usetex'] = True + #rcParams['cairo.format'] = 'svg' ###################### @@ -63,5 +71,5 @@ if not var: back = "mola" ## if no var: draw mola elif typefile in ['mesoapi','meso'] \ - and proj not in ['merc','lcc']: back = "molabw" ## if var but meso: draw molabw + and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw else: pass ## else: draw None @@ -74,15 +82,13 @@ elif typefile in ['gcm']: [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) - if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") - elif proj == "lcc": [wlon,wlat] = wrfinterv(lon2d,lat2d) - else: [wlon,wlat] = simplinterv(lon2d,lat2d) - - ################### - #### Get local time - #if typefile in ['mesoapi','meso']: - # ltst = interv[0] - # #ltst = int(getattr(nc, 'GMT') + 0.5*(wlon[0]+wlon[1])/15.) - #elif typefile in ['gcm']: - # ltst = 0 + if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") + elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) + else: [wlon,wlat] = simplinterv(lon2d,lat2d) + if zoom: + dlon = abs(wlon[1]-wlon[0])/2. + dlat = abs(wlat[1]-wlat[0])/2. + [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\ + [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ] + print "zoom %",zoom,wlon,wlat ############################################################################## @@ -111,5 +117,5 @@ elif dimension == 3: field = nc.variables[var][:,:,:] elif dimension == 4: field = nc.variables[var][:,nvert,:,:] - dev = np.std(field)*2.5 + dev = np.std(field)*3.0 if vmin is None: zevmin = mean(field) - dev else: zevmin = vmin @@ -117,5 +123,9 @@ else: zevmax = vmax print "bounds ", zevmin, zevmax - + ### some already defined colormaps + if colorb is True: colorb = defcolorb(var) + else: + dimension = 0 + ########################### ### Get length of time axis @@ -127,15 +137,20 @@ ######################################### ### Name for title and graphics save file - if winds: basename = "WINDS_" + if winds: basename = "UV_" else: basename = "" - if var: - basename = basename + var - if typefile is 'meso': stralt = "_lvl" + str(nvert) - elif typefile is 'mesoapi': - zelevel = int(nc.variables['vert'][nvert]) - if 'altitude_abg' in nc.dimensions: stralt = "_"+str(zelevel)+"m-ALS" - elif 'bottom_top' in nc.dimensions: stralt = "_"+str(zelevel)+"m-AMR" - elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" - else: stralt = "" + if var: basename = basename + var + ### + if dimension == 4: + if typefile is 'meso': stralt = "_lvl" + str(nvert) + elif typefile is 'mesoapi': + zelevel = int(nc.variables['vert'][nvert]) + if 'altitude' in nc.dimensions: stralt = "_"+str(zelevel)+"m-AMR" + elif 'altitude_abg' in nc.dimensions: stralt = "_"+str(zelevel)+"m-ALS" + elif 'bottom_top' in nc.dimensions: stralt = "_"+str(zelevel)+"m" + elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" + else: stralt = "" + else: + stralt = "" + ### basename = basename + stralt @@ -146,5 +161,5 @@ if numplot == 4: sub = 221 - fig.subplots_adjust(wspace = 0.1, hspace = 0.3) + fig.subplots_adjust(wspace = 0.3, hspace = 0.3) rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) elif numplot == 2: @@ -154,22 +169,22 @@ elif numplot == 3: sub = 131 - fig.subplots_adjust(wspace = 0.3) - rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) + fig.subplots_adjust(wspace = 0.5) + rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) elif numplot == 6: sub = 231 fig.subplots_adjust(wspace = 0.4, hspace = 0.0) - rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) + rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) elif numplot == 8: sub = 331 #241 - fig.subplots_adjust(wspace = 0.1, hspace = 0.3) + fig.subplots_adjust(wspace = 0.3, hspace = 0.3) rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) elif numplot == 9: - sub = 331 - fig.subplots_adjust(wspace = 0.1, hspace = 0.3) + sub = 331 + fig.subplots_adjust(wspace = 0.3, hspace = 0.3) rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) elif numplot == 1: sub = 99999 else: - print "supported: 1,2,3,4,6,8" + print "supported: 1,2,3,4,6,8,9" exit() ### Prepare time loop @@ -203,6 +218,6 @@ ### If only one local time is requested (numplot < 0) elif numplot <= 0: - if ltst + numplot != 0: continue - else: found_lct = True + if int(ltst) + numplot != 0: continue + else: found_lct = True ### Map projection @@ -216,8 +231,18 @@ if dimension == 2: what_I_plot = field[:,:] elif dimension == 3: what_I_plot = field[i,:,:] - #contourf( x, y, what_I_plot, 30, cmap = get_cmap(name=colorb), vmin=zevmin, vmax=zevmax ) - pcolor( x, y, what_I_plot, cmap = get_cmap(name=colorb), vmin=zevmin, vmax=zevmax ) - if colorb: colorbar(fraction=0.05,pad=0.1,format='%.1e') - + palette = get_cmap(name=colorb) + #palette.set_over('b', 1.0) + if not tile: + zelevels = np.linspace(zevmin,zevmax) + contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels ) + else: + pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax ) + if var in ['HGT']: pass + elif colorb: + ndiv = 10 + colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\ + ticks=np.linspace(zevmin,zevmax,ndiv+1),\ + extend='max',spacing='proportional') + # both min max neither ### Vector plot if winds: @@ -229,6 +254,6 @@ key = False if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d) - if var == None and back == "vishires": colorvec = 'w' - else: colorvec = 'k' + if var == None: colorvec = definecolorvec(back) + else: colorvec = definecolorvec(colorb) vectorfield(vecx, vecy,\ x, y, stride=stride, csmooth=stride,\ @@ -237,6 +262,7 @@ ### Next subplot - plottitle = basename + "_LT"+str(ltst) - if addchar: plottitle = plottitle + addchar + plottitle = basename + if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) + else: plottitle = plottitle + "_LT"+str(ltst) ptitle( plottitle ) sub += 1 @@ -244,25 +270,21 @@ ########################################################################## ### Save the figure in a file in the data folder or an user-defined folder - if not target: zeplot = namefile+"_"+basename - else: zeplot = target+"/"+basename - if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot)) + if typefile in ['meso','mesoapi']: prefix = getprefix(nc) + elif typefile in ['gcm']: prefix = 'LMD_GCM_' + else: prefix = '' + ### + zeplot = prefix + basename + if addchar: zeplot = zeplot + addchar + if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot)) + ### + if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot + else: zeplot = target + "/" + zeplot + ### if found_lct: makeplotpng(zeplot,pad_inches_value=0.35) else: print "Local time not found" + + ############### + ### Now the end return zeplot - - -#################################################### -#################################################### -### A simple program to get wind vectors' components -def getwinds (nc,charu='Um',charv='Vm'): - import numpy as np - u = nc.variables[charu] - v = nc.variables[charv] - if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1] - if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :] - ### ou alors prendre les coordonnees speciales - return u,v - - ########################################################################################### @@ -288,8 +310,12 @@ parser.add_option('-n', action='store', dest='numplot', type="int", default=4, 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=None, help='change colormap (and draw a colorbar)') - parser.add_option('-x', action='store_false', dest='winds', default=True, help='flag: no wind vectors') + parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap') + parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors') parser.add_option('-m', action='store', dest='vmin', type="float", default=None, help='bounding minimum value for color plot') parser.add_option('-M', action='store', dest='vmax', type="float", default=None, help='bounding maximum value for color plot') + parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)') + parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %') + parser.add_option('-N', action='store_true', dest='nocall', default=False, help='do not recreate api file') + #parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables') (opt,args) = parser.parse_args() if opt.namefile is None: @@ -319,5 +345,7 @@ fields = zefields, \ interp_method = opt.interp, \ - onelevel = zelevel ) + onelevel = zelevel, \ + nocall = opt.nocall ) + print zefile zelevel = 0 ## so that zelevel could play again the role of nvert @@ -326,5 +354,6 @@ name = winds (zefile,int(zelevel),\ proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ - addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax) + addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) + print 'Done: '+name #########################################################