Context Navigation

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

Last change on this file since 232 was 232, checked in by aslmd, 14 years ago
MESOSCALE: corrections petits bugs graphiques.
Property svn:executable set to ``*
File size: 19.6 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,\
	25	zoom=None):
[184]	26
[192]	27	####################################################################################################################
	28	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
	29
[184]	30	#################################
	31	### Load librairies and functions
[180]	32	from netCDF4 import Dataset
[202]	33	from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
[232]	34	fmtvar,definecolorvec,getwinds,defcolorb,getprefix,putpoints
[199]	35	from mymath import deg,max,min,mean
[202]	36	from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor
[192]	37	from matplotlib.cm import get_cmap
[182]	38	import numpy as np
[202]	39	from numpy.core.defchararray import find
[184]	40
[202]	41	###
	42	#rcParams['text.usetex'] = True
	43	#rcParams['cairo.format'] = 'svg'
	44
[184]	45	######################
	46	### Load NETCDF object
[180]	47	nc = Dataset(namefile)
[184]	48
	49	###################################
	50	### Recognize predefined file types
	51	if 'controle' in nc.variables: typefile = 'gcm'
[191]	52	elif 'vert' in nc.variables: typefile = 'mesoapi'
[184]	53	elif 'U' in nc.variables: typefile = 'meso'
[225]	54	elif 'HGT_M' in nc.variables: typefile = 'geo'
[184]	55	else:
	56	print "typefile not supported."
[191]	57	print nc.variables
[182]	58	exit()
[184]	59
	60	##############################################################
	61	### Try to guess the projection from wrfout if not set by user
[225]	62	if typefile in ['mesoapi','meso','geo']:
[184]	63	if proj == None: proj = getproj(nc)
	64	### (il faudrait passer CEN_LON dans la projection ?)
	65	elif typefile in ['gcm']:
	66	if proj == None: proj = "cyl"
	67	## pb avec les autres (de trace derriere la sphere ?)
	68
[188]	69	############################################
	70	#### Choose underlying topography by default
	71	if not back:
	72	if not var: back = "mola" ## if no var: draw mola
[225]	73	elif typefile in ['mesoapi','meso','geo'] \
	74	and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
[188]	75	else: pass ## else: draw None
[184]	76
	77	####################################################
	78	### Get geographical coordinates and plot boundaries
	79	if typefile in ['mesoapi','meso']:
	80	[lon2d,lat2d] = getcoord2d(nc)
	81	lon2d = dumpbdy(lon2d)
	82	lat2d = dumpbdy(lat2d)
	83	elif typefile in ['gcm']:
	84	[lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
[225]	85	elif typefile in ['geo']:
	86	[lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
[207]	87	if proj == "npstere": [wlon,wlat] = latinterv("North_Pole")
	88	elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)
	89	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
[225]	90	print wlon
	91	print wlat
[204]	92	if zoom:
	93	dlon = abs(wlon[1]-wlon[0])/2.
	94	dlat = abs(wlat[1]-wlat[0])/2.
	95	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
	96	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
	97	print "zoom %",zoom,wlon,wlat
[184]	98
	99	##############################################################################
	100	### Get winds and know if those are meteorological winds (ie. zon, mer) or not
[191]	101	if winds:
	102	if typefile is 'mesoapi':
	103	[u,v] = getwinds(nc)
	104	metwind = True ## meteorological (zon/mer)
	105	elif typefile is 'gcm':
	106	[u,v] = getwinds(nc,charu='u',charv='v')
	107	metwind = True ## meteorological (zon/mer)
	108	elif typefile is 'meso':
	109	[u,v] = getwinds(nc,charu='U',charv='V')
	110	metwind = False ## geometrical (wrt grid)
	111	print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
[225]	112	elif typefile is 'geo':
	113	winds = None
[184]	114
[191]	115	#####################################################
[184]	116	### Load the chosen variables, whether it is 2D or 3D
	117	if var:
	118	if var not in nc.variables:
	119	print "not found in file:",var
	120	exit()
	121	else:
	122	dimension = np.array(nc.variables[var]).ndim
	123	if dimension == 2: field = nc.variables[var][:,:]
	124	elif dimension == 3: field = nc.variables[var][:,:,:]
	125	elif dimension == 4: field = nc.variables[var][:,nvert,:,:]
[232]	126	fieldcalc = field[ field < 9e+35 ]
	127	dev = np.std(fieldcalc)*2.0
	128	if vmin is None: zevmin = mean(fieldcalc) - dev
[199]	129	else: zevmin = vmin
[232]	130	if vmax is None: zevmax = mean(fieldcalc) + dev
[199]	131	else: zevmax = vmax
[225]	132	if vmin == vmax:
[232]	133	zevmin = mean(fieldcalc) - dev ### for continuity
	134	zevmax = mean(fieldcalc) + dev ### for continuity
	135	print "field ", min(fieldcalc), max(fieldcalc)
[199]	136	print "bounds ", zevmin, zevmax
[202]	137	### some already defined colormaps
	138	if colorb is True: colorb = defcolorb(var)
[207]	139	else:
	140	dimension = 0
[201]	141
[191]	142	###########################
	143	### Get length of time axis
	144	if winds: nt = len(u[:,0,0,0])
	145	elif var:
	146	if dimension == 2: nt = 1
	147	else : nt = len(field[:,0,0])
	148
[186]	149	#########################################
	150	### Name for title and graphics save file
[225]	151	if var and winds: basename = var + '_UV'
	152	elif var: basename = var
	153	elif winds: basename = 'UV'
	154	else: exit()
[202]	155	###
[225]	156	if dimension == 4 or winds:
[207]	157	if typefile is 'meso': stralt = "_lvl" + str(nvert)
	158	elif typefile is 'mesoapi':
	159	zelevel = int(nc.variables['vert'][nvert])
[232]	160	if abs(zelevel) < 10000.: strheight=str(zelevel)+"m"
	161	else: strheight=str(int(zelevel/1000.))+"km"
	162	if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR"
	163	elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS"
	164	elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight
[207]	165	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
	166	else: stralt = ""
	167	else:
	168	stralt = ""
[202]	169	###
[188]	170	basename = basename + stralt
[186]	171
[184]	172	##################################
	173	### Open a figure and set subplots
[180]	174	fig = figure()
[225]	175	rcParams['font.size'] = 12.
	176	if typefile in ['geo']: numplot = 1
[186]	177	if numplot > 0:
	178	if numplot == 4:
	179	sub = 221
[201]	180	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
[186]	181	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
	182	elif numplot == 2:
	183	sub = 121
[232]	184	fig.subplots_adjust(wspace = 0.35)
[186]	185	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
	186	elif numplot == 3:
	187	sub = 131
[202]	188	fig.subplots_adjust(wspace = 0.5)
	189	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[186]	190	elif numplot == 6:
	191	sub = 231
	192	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
[202]	193	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
[186]	194	elif numplot == 8:
	195	sub = 331 #241
[202]	196	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
[186]	197	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	198	elif numplot == 9:
[207]	199	sub = 331
[202]	200	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
[186]	201	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	202	elif numplot == 1:
[189]	203	sub = 99999
[186]	204	else:
[207]	205	print "supported: 1,2,3,4,6,8,9"
[186]	206	exit()
	207	### Prepare time loop
	208	if nt <= numplot or numplot == 1:
	209	tabrange = [0]
	210	numplot = 1
	211	else:
	212	tabrange = range(0,nt,int(nt/numplot)) #nt-1
	213	tabrange = tabrange[0:numplot]
	214	else:
	215	tabrange = range(0,nt,1)
	216	sub = 99999
	217	print tabrange
[184]	218
	219	#################################
	220	### Time loop for plotting device
[186]	221	found_lct = False
[184]	222	for i in tabrange:
	223
[197]	224	### Which local time ?
	225	ltst = ( interv[0] + 0.5(wlon[0]+wlon[1])/15.) + iinterv[1]
	226	ltst = int (ltst * 10) / 10.
	227	ltst = ltst % 24
	228
[184]	229	### General plot settings
[186]	230	if numplot > 1:
	231	subplot(sub)
	232	found_lct = True
	233	elif numplot == 1:
	234	found_lct = True
[188]	235	### If only one local time is requested (numplot < 0)
[186]	236	elif numplot <= 0:
[201]	237	if int(ltst) + numplot != 0: continue
	238	else: found_lct = True
[184]	239
	240	### Map projection
[180]	241	m = define_proj(proj,wlon,wlat,back=back)
	242	x, y = m(lon2d, lat2d)
[184]	243
	244	#### Contour plot
	245	if var:
	246	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:])
[225]	247	elif typefile in ['geo']: what_I_plot = field[0,:,:]
[184]	248	elif typefile in ['gcm']:
	249	if dimension == 2: what_I_plot = field[:,:]
	250	elif dimension == 3: what_I_plot = field[i,:,:]
[204]	251	palette = get_cmap(name=colorb)
	252	#palette.set_over('b', 1.0)
[225]	253	#print np.array(x).shape
	254	#print np.array(y).shape
	255	#print np.array(what_I_plot).shape
[201]	256	if not tile:
	257	zelevels = np.linspace(zevmin,zevmax)
[232]	258	hole = True
	259	if not hole:
	260	what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7)
	261	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
	262	what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
[204]	263	contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels )
[201]	264	else:
[204]	265	pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
[232]	266	#putpoints(m,fig)
[225]	267	if var in ['HGT']: pass
[204]	268	elif colorb:
[225]	269	ndiv = 10
	270	colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\
	271	ticks=np.linspace(zevmin,zevmax,ndiv+1),\
	272	extend='max',spacing='proportional')
	273	# both min max neither
[184]	274	### Vector plot
[191]	275	if winds:
	276	if typefile in ['mesoapi','meso']:
	277	[vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
	278	key = True
	279	elif typefile in ['gcm']:
	280	[vecx,vecy] = [ u[i,nvert,:,:] , v[i,nvert,:,:] ]
	281	key = False
	282	if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
[202]	283	if var == None: colorvec = definecolorvec(back)
	284	else: colorvec = definecolorvec(colorb)
[191]	285	vectorfield(vecx, vecy,\
[225]	286	x, y, stride=stride, csmooth=2,\
[188]	287	scale=15., factor=300., color=colorvec, key=key)
[225]	288	#200. ## or csmooth=stride
[184]	289
	290	### Next subplot
[202]	291	plottitle = basename
[225]	292	if typefile in ['mesoapi','meso']:
	293	if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst)
	294	else: plottitle = plottitle + "_LT"+str(ltst)
[195]	295	ptitle( plottitle )
[180]	296	sub += 1
[184]	297
	298	##########################################################################
	299	### Save the figure in a file in the data folder or an user-defined folder
[202]	300	if typefile in ['meso','mesoapi']: prefix = getprefix(nc)
	301	elif typefile in ['gcm']: prefix = 'LMD_GCM_'
	302	else: prefix = ''
	303	###
	304	zeplot = prefix + basename
	305	if addchar: zeplot = zeplot + addchar
	306	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
	307	###
	308	if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
	309	else: zeplot = target + "/" + zeplot
	310	###
[186]	311	if found_lct: makeplotpng(zeplot,pad_inches_value=0.35)
	312	else: print "Local time not found"
[202]	313
	314	###############
	315	### Now the end
[193]	316	return zeplot
[180]	317
[184]	318	###########################################################################################
	319	###########################################################################################
	320	### What is below relate to running the file as a command line executable (very convenient)
[180]	321	if __name__ == "__main__":
	322	import sys
[184]	323	from optparse import OptionParser ### to be replaced by argparse
[186]	324	from api_wrapper import api_onelevel
[195]	325	from netCDF4 import Dataset
	326	from myplot import getlschar
[225]	327	from numpy import ones
[193]	328
	329	#############################
	330	### Get options and variables
[182]	331	parser = OptionParser()
[225]	332	parser.add_option('-f', action='append', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)')
[188]	333	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)')
[184]	334	parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection')
[182]	335	parser.add_option('-b', action='store', dest='back', type="string", default=None, help='background')
	336	parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder')
[186]	337	parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)')
[225]	338	parser.add_option('-v', action='append', dest='var', type="string", default=None, help='variable contoured (append)')
	339	parser.add_option('-n', action='store', dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -numplot)')
[188]	340	parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)')
[202]	341	parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap')
[207]	342	parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors')
[225]	343	parser.add_option('-m', action='append', dest='vmin', type="float", default=None, help='bounding minimum value (append)')
	344	parser.add_option('-M', action='append', dest='vmax', type="float", default=None, help='bounding maximum value (append)')
[204]	345	parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)')
	346	parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %')
[207]	347	parser.add_option('-N', action='store_true', dest='nocall', default=False, help='do not recreate api file')
[201]	348	#parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables')
[182]	349	(opt,args) = parser.parse_args()
	350	if opt.namefile is None:
	351	print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
	352	exit()
	353	print "Options:", opt
[225]	354
	355	if opt.var is None:
	356	pass
	357	else:
	358	listvar = ''
	359	zelen = len(opt.var)
	360	if zelen == 1: listvar = opt.var[0] + ','
	361	else :
	362	for jjj in range(zelen): listvar += opt.var[jjj] + ','
	363	listvar = listvar[0:len(listvar)-1]
	364	if opt.vmin:
	365	if zelen != len(opt.vmin):
	366	print "not enough or too much vmin values... setting same values all variables"
	367	vmintab = ones(zelen) * opt.vmin[0]
	368	else:
	369	vmintab = opt.vmin
	370	else:
	371	vmintab = None
	372	if opt.vmax:
	373	if zelen != len(opt.vmax):
	374	print "not enough or too much vmax values... setting same values all variables"
	375	vmaxtab = ones(zelen) * opt.vmax[0]
	376	else:
	377	vmaxtab = opt.vmax
	378	else:
	379	vmaxtab = None
[193]	380
[225]	381	for i in range(len(opt.namefile)):
[184]	382
[225]	383	zefile = opt.namefile[i]
	384	print zefile
	385	zelevel = opt.nvert
	386	stralt = None
	387	[lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file)
	388
	389	#####################################################
	390	### Call Fortran routines for vertical interpolations
	391	if opt.interp is not None:
	392	if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010
	393	### winds or no winds
	394	if opt.winds : zefields = 'uvmet'
	395	else : zefields = ''
	396	### var or no var
	397	if opt.var is None : pass
	398	elif zefields == '' : zefields = listvar
	399	else : zefields = zefields + "," + listvar
	400	print zefields
	401	zefile = api_onelevel ( path_to_input = '', \
	402	input_name = zefile, \
	403	#path_to_output = opt.target, \
	404	fields = zefields, \
	405	interp_method = opt.interp, \
	406	onelevel = zelevel, \
	407	nocall = opt.nocall )
	408	print zefile
	409	zelevel = 0 ## so that zelevel could play again the role of nvert
[193]	410
[225]	411	if opt.var is None:
	412	#############
	413	### Main call
	414	name = winds (zefile,int(zelevel),\
	415	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
	416	addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom)
	417	print 'Done: '+name
	418	else:
	419	for jjj in range(len(opt.var)):
	420	if vmintab: argvmin = vmintab[jjj]
	421	else: argvmin = None
	422	if vmaxtab: argvmax = vmaxtab[jjj]
	423	else: argvmax = None
	424	#############
	425	### Main call
	426	name = winds (zefile,int(zelevel),\
	427	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var[jjj],numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\
	428	addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom)
	429	print 'Done: '+name
	430
	431	#########################################################
	432	### Generate a .sh file with the used command saved in it
	433	command = ""
	434	for arg in sys.argv: command = command + arg + ' '
	435	f = open(name+'.sh', 'w')
	436	f.write(command)

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

Download in other formats: