Context Navigation

source: trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py @ 238

Last change on this file since 238 was 238, checked in by aslmd, 13 years ago
MESOSCALE: small adjustments to graphical plots routines, including additional maps.
Property svn:executable set to ``*
File size: 14.7 KB

Rev	Line
[180]	1	#!/usr/bin/env python
	2
[202]	3	### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011
[184]	4	### Thanks to A. Colaitis for the parser trick
[180]	5
[184]	6
	7	####################################
	8	####################################
	9	### The main program to plot vectors
[182]	10	def winds (namefile,\
	11	nvert,\
[184]	12	proj=None,\
[182]	13	back=None,\
	14	target=None,
[184]	15	stride=3,\
[225]	16	numplot=2,\
[187]	17	var=None,\
[201]	18	colorb=True,\
[195]	19	winds=True,\
[199]	20	addchar=None,\
	21	interv=[0,1],\
	22	vmin=None,\
[201]	23	vmax=None,\
[204]	24	tile=False,\
[233]	25	zoom=None,\
	26	display=True,\
[238]	27	itstep=None,\
	28	hole=False):
[184]	29
[192]	30	####################################################################################################################
	31	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
	32
[184]	33	#################################
	34	### Load librairies and functions
[180]	35	from netCDF4 import Dataset
[237]	36	from myplot import getcoord2d,define_proj,makeplotpng,makeplotpngres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
[233]	37	fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
	38	zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield
[199]	39	from mymath import deg,max,min,mean
[202]	40	from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor
[192]	41	from matplotlib.cm import get_cmap
[182]	42	import numpy as np
[202]	43	from numpy.core.defchararray import find
[184]	44
	45	######################
	46	### Load NETCDF object
[180]	47	nc = Dataset(namefile)
[233]	48
	49	##################################
	50	### Initial checks and definitions
	51	typefile = whatkindfile(nc) ## TYPEFILE
	52	if var not in nc.variables: var = False ## VAR
[238]	53	if winds: ## WINDS
	54	[uchar,vchar,metwind] = getwinddef(nc)
	55	if uchar == 'not found': winds = False
[233]	56	[lon2d,lat2d] = getcoorddef(nc) ## COORDINATES, could be moved below
	57	if proj == None: proj = getproj(nc) ## PROJECTION
[184]	58
[233]	59	##########################
	60	### Define plot boundaries
[207]	61	if proj == "npstere": [wlon,wlat] = latinterv("North_Pole")
	62	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
	63	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
[233]	64	if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom)
[184]	65
[186]	66	#########################################
	67	### Name for title and graphics save file
[225]	68	if var and winds: basename = var + '_UV'
	69	elif var: basename = var
	70	elif winds: basename = 'UV'
[207]	71	else:
[233]	72	print nc.variables
	73	errormess("please set at least winds or var")
	74	basename = basename + getstralt(nc,nvert) ## can be moved elsewhere for a more generic routine
[186]	75
[184]	76	##################################
	77	### Open a figure and set subplots
[180]	78	fig = figure()
[233]	79	sub = definesubplot( numplot, fig )
	80
[184]	81	#################################
	82	### Time loop for plotting device
[186]	83	found_lct = False
[233]	84	itime = 0 ## could be an argument
	85	nplot = 1
	86	error = False
	87	if itstep is None: itstep = int(24./numplot)
	88	while error is False:
[184]	89
[197]	90	### Which local time ?
[237]	91	#print interv[0], interv[1], itime
[233]	92	ltst = ( interv[0] + 0.5(wlon[0]+wlon[1])/15.) + itimeinterv[1]
[197]	93	ltst = int (ltst * 10) / 10.
	94	ltst = ltst % 24
	95
[184]	96	### General plot settings
[233]	97	#print itime, int(ltst), numplot, nplot
	98	if numplot >= 1:
	99	if nplot > numplot: break
	100	if numplot > 1:
	101	if typefile not in ['geo']: subplot(sub+nplot-1)
	102
[186]	103	found_lct = True
[233]	104	### If only one local time is requested (numplot < 0)
[186]	105	elif numplot <= 0:
[233]	106	if int(ltst) + numplot != 0:
	107	itime += 1
	108	if found_lct is True: break ## because it means LT was found at previous iteration
	109	else: continue ## continue to iterate to find the correct LT
	110	else:
	111	found_lct = True
[184]	112
	113	### Map projection
[180]	114	m = define_proj(proj,wlon,wlat,back=back)
	115	x, y = m(lon2d, lat2d)
[184]	116
	117	#### Contour plot
	118	if var:
[233]	119	what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d3=nvert )
	120	if not error:
	121	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot)
	122	zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
	123	if colorb is True: colorb = defcolorb(var)
	124	palette = get_cmap(name=colorb)
	125	if not tile:
[238]	126	if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
[233]	127	zelevels = np.linspace(zevmin,zevmax)
	128	contourf( x, y, what_I_plot, zelevels, cmap = palette )
	129	else:
	130	pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
	131	#putpoints(m,fig)
	132	if var in ['HGT']: pass
	133	elif colorb:
	134	ndiv = 10
	135	colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
	136	ticks=np.linspace(zevmin,zevmax,ndiv+1),\
	137	extend='max',spacing='proportional')
	138	# both min max neither
	139
[184]	140	### Vector plot
[191]	141	if winds:
[233]	142	vecx, error = reducefield( getfield(nc,uchar), d4=itime, d3=nvert )
	143	vecy, error = reducefield( getfield(nc,vchar), d4=itime, d3=nvert )
	144	if not error:
	145	if typefile in ['mesoapi','meso']:
	146	[vecx,vecy] = [dumpbdy(vecx,stag=uchar), dumpbdy(vecy,stag=vchar)]
	147	key = True
	148	elif typefile in ['gcm']:
	149	key = False
	150	if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
	151	if var == False: colorvec = definecolorvec(back)
	152	else: colorvec = definecolorvec(colorb)
	153	vectorfield(vecx, vecy,\
	154	x, y, stride=stride, csmooth=2,\
	155	scale=15., factor=300., color=colorvec, key=key)
	156	#200. ## or csmooth=stride
	157
[184]	158	### Next subplot
[202]	159	plottitle = basename
[225]	160	if typefile in ['mesoapi','meso']:
[233]	161	if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst)
	162	else: plottitle = plottitle + "_LT"+str(ltst)
[195]	163	ptitle( plottitle )
[233]	164	itime += itstep
	165	nplot += 1
[184]	166
	167	##########################################################################
	168	### Save the figure in a file in the data folder or an user-defined folder
[202]	169	if typefile in ['meso','mesoapi']: prefix = getprefix(nc)
	170	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
	171	else: prefix = ''
	172	###
	173	zeplot = prefix + basename
	174	if addchar: zeplot = zeplot + addchar
	175	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
	176	###
	177	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
	178	else: zeplot = target + "/" + zeplot
	179	###
[238]	180	if found_lct: makeplotpng(zeplot,pad_inches_value=0.35,disp=display)
[237]	181	#if found_lct: makeplotpngres(zeplot,200,disp=display)
[186]	182	else: print "Local time not found"
[202]	183
	184	###############
	185	### Now the end
[193]	186	return zeplot
[180]	187
[233]	188	##############################
	189	### A specific stuff for below
	190	def adjust_length (tab, zelen):
	191	from numpy import ones
	192	if tab is None:
	193	outtab = ones(zelen) * -999999
	194	else:
	195	print zelen, len(tab)
	196	if zelen != len(tab):
	197	print "not enough or too much values... setting same values all variables"
	198	outtab = ones(zelen) * tab[0]
	199	else:
	200	outtab = tab
	201	return outtab
	202
[184]	203	###########################################################################################
	204	###########################################################################################
	205	### What is below relate to running the file as a command line executable (very convenient)
[180]	206	if __name__ == "__main__":
	207	import sys
[184]	208	from optparse import OptionParser ### to be replaced by argparse
[186]	209	from api_wrapper import api_onelevel
[195]	210	from netCDF4 import Dataset
	211	from myplot import getlschar
[193]	212
	213	#############################
	214	### Get options and variables
[182]	215	parser = OptionParser()
[233]	216	parser.add_option('-f', action='append',dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)')
	217	parser.add_option('-l', action='store',dest='nvert', type="float", default=0, help='vertical level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)')
	218	parser.add_option('-p', action='store',dest='proj', type="string", default=None, help='projection')
	219	parser.add_option('-b', action='store',dest='back', type="string", default=None, help='background image (def: None)')
	220	parser.add_option('-t', action='store',dest='target', type="string", default=None, help='destination folder')
	221	parser.add_option('-s', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)')
	222	parser.add_option('-v', action='append',dest='var', type="string", default=None, help='variable contoured (append)')
	223	parser.add_option('-n', action='store',dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -numplot)')
	224	parser.add_option('-i', action='store',dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)')
	225	parser.add_option('-c', action='store',dest='colorb', type="string", default=True, help='change colormap')
	226	parser.add_option('-x', action='store_false',dest='winds', default=True, help='no wind vectors')
	227	parser.add_option('-m', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)')
	228	parser.add_option('-M', action='append',dest='vmax', type="float", default=None, help='bounding maximum value (append)')
	229	parser.add_option('-T', action='store_true',dest='tile', default=False, help='draw a tiled plot (no blank zone)')
	230	parser.add_option('-z', action='store',dest='zoom', type="float", default=None, help='zoom factor in %')
	231	parser.add_option('-N', action='store_true',dest='nocall', default=False, help='do not recreate api file')
	232	parser.add_option('-d', action='store_false',dest='display', default=True, help='do not pop up created images')
	233	parser.add_option('-e', action='store',dest='itstep', type="int", default=None, help='stride time (def=4)')
[238]	234	parser.add_option('-H', action='store_true',dest='hole', default=False, help='holes above max and below min')
[201]	235	#parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables')
[182]	236	(opt,args) = parser.parse_args()
	237	if opt.namefile is None:
	238	print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
	239	exit()
	240	print "Options:", opt
[233]	241
	242	listvar = ''
[225]	243	if opt.var is None:
[233]	244	zerange = [-999999]
[225]	245	else:
	246	zelen = len(opt.var)
[233]	247	zerange = range(zelen)
	248	#if zelen == 1: listvar = opt.var[0] + ','
	249	#else :
	250	for jjj in zerange: listvar += opt.var[jjj] + ','
[225]	251	listvar = listvar[0:len(listvar)-1]
[233]	252	vmintab = adjust_length (opt.vmin, zelen)
	253	vmaxtab = adjust_length (opt.vmax, zelen)
[193]	254
[225]	255	for i in range(len(opt.namefile)):
[184]	256
[225]	257	zefile = opt.namefile[i]
	258	print zefile
	259	zelevel = opt.nvert
	260	stralt = None
	261	[lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file)
	262
	263	#####################################################
	264	### Call Fortran routines for vertical interpolations
	265	if opt.interp is not None:
	266	if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010
	267	### winds or no winds
	268	if opt.winds : zefields = 'uvmet'
	269	else : zefields = ''
	270	### var or no var
[233]	271	#if opt.var is None : pass
	272	if zefields == '' : zefields = listvar
[225]	273	else : zefields = zefields + "," + listvar
	274	print zefields
	275	zefile = api_onelevel ( path_to_input = '', \
	276	input_name = zefile, \
	277	#path_to_output = opt.target, \
	278	fields = zefields, \
	279	interp_method = opt.interp, \
	280	onelevel = zelevel, \
	281	nocall = opt.nocall )
	282	print zefile
	283	zelevel = 0 ## so that zelevel could play again the role of nvert
[193]	284
[233]	285	if opt.var is None: zerange = [-999999]
	286	else: zerange = range(zelen)
	287	for jjj in zerange:
	288	if jjj == -999999:
	289	argvar = None
	290	argvmin = None
	291	argvmax = None
	292	else:
	293	argvar = opt.var[jjj]
	294	if vmintab[jjj] != -999999: argvmin = vmintab[jjj]
	295	else: argvmin = None
	296	if vmaxtab[jjj] != -999999: argvmax = vmaxtab[jjj]
	297	else: argvmax = None
[225]	298	#############
	299	### Main call
	300	name = winds (zefile,int(zelevel),\
[233]	301	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\
	302	numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
	303	addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,\
	304	tile=opt.tile,zoom=opt.zoom,display=opt.display,\
[238]	305	itstep=opt.itstep,hole=opt.hole)
[225]	306	print 'Done: '+name
	307
	308	#########################################################
	309	### Generate a .sh file with the used command saved in it
	310	command = ""
	311	for arg in sys.argv: command = command + arg + ' '
	312	f = open(name+'.sh', 'w')
	313	f.write(command)

Note: See TracBrowser for help on using the repository browser.

Download in other formats: