Context Navigation

source: trunk/MESOSCALE/PLOT/PYTHON/scripts/winds.py @ 188

Last change on this file since 188 was 188, checked in by aslmd, 14 years ago
MESOSCALE: python script for plots. changed defaults values and added automatic captions.
Property svn:executable set to ``*
File size: 11.8 KB

Rev	Line
[180]	1	#!/usr/bin/env python
	2
[184]	3	### A. Spiga -- LMD -- 30/05/2011
	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,\
	16	numplot=4,\
[187]	17	var=None,\
[188]	18	colorb=None):
[184]	19
	20	#################################
	21	### Load librairies and functions
[180]	22	from netCDF4 import Dataset
[184]	23	from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy
[187]	24	from matplotlib.pyplot import contourf, subplot, figure, rcParams, savefig, colorbar
[182]	25	import numpy as np
[184]	26
	27	######################
	28	### Load NETCDF object
[180]	29	nc = Dataset(namefile)
[184]	30
	31	###################################
	32	### Recognize predefined file types
	33	if 'controle' in nc.variables: typefile = 'gcm'
	34	elif 'Um' in nc.variables: typefile = 'mesoapi'
	35	elif 'U' in nc.variables: typefile = 'meso'
	36	else:
	37	print "typefile not supported."
[182]	38	exit()
[184]	39
	40	##############################################################
	41	### Try to guess the projection from wrfout if not set by user
	42	if typefile in ['mesoapi','meso']:
	43	if proj == None: proj = getproj(nc)
	44	### (il faudrait passer CEN_LON dans la projection ?)
	45	elif typefile in ['gcm']:
	46	if proj == None: proj = "cyl"
	47	## pb avec les autres (de trace derriere la sphere ?)
	48
[188]	49	############################################
	50	#### Choose underlying topography by default
	51	if not back:
	52	if not var: back = "mola" ## if no var: draw mola
	53	elif typefile in ['mesoapi','meso'] \
	54	and proj not in ['merc','lcc']: back = "molabw" ## if var but meso: draw molabw
	55	else: pass ## else: draw None
[184]	56
	57	####################################################
	58	### Get geographical coordinates and plot boundaries
	59	if typefile in ['mesoapi','meso']:
	60	[lon2d,lat2d] = getcoord2d(nc)
	61	lon2d = dumpbdy(lon2d)
	62	lat2d = dumpbdy(lat2d)
	63	elif typefile in ['gcm']:
	64	[lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
	65	if proj == "npstere": [wlon,wlat] = latinterv("North_Pole")
	66	elif proj == "lcc": [wlon,wlat] = wrfinterv(lon2d,lat2d)
	67	else: [wlon,wlat] = simplinterv(lon2d,lat2d)
	68
	69	##################
	70	### Get local time
	71	if typefile in ['mesoapi','meso']: ltst = int(getattr(nc, 'GMT') + 0.5*(wlon[0]+wlon[1])/15.)
	72	elif typefile in ['gcm']: ltst = 0
	73
	74	##############################################################################
	75	### Get winds and know if those are meteorological winds (ie. zon, mer) or not
	76	if typefile is 'mesoapi':
	77	[u,v] = getwinds(nc)
	78	metwind = True ## meteorological (zon/mer)
	79	elif typefile is 'gcm':
	80	[u,v] = getwinds(nc,charu='u',charv='v')
	81	metwind = True ## meteorological (zon/mer)
	82	elif typefile is 'meso':
	83	[u,v] = getwinds(nc,charu='U',charv='V')
	84	metwind = False ## geometrical (wrt grid)
[186]	85	print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)."
[184]	86
	87	####################################################
	88	### Load the chosen variables, whether it is 2D or 3D
	89	if var:
	90	if var not in nc.variables:
	91	print "not found in file:",var
	92	exit()
	93	else:
	94	dimension = np.array(nc.variables[var]).ndim
	95	if dimension == 2: field = nc.variables[var][:,:]
	96	elif dimension == 3: field = nc.variables[var][:,:,:]
	97	elif dimension == 4: field = nc.variables[var][:,nvert,:,:]
[186]	98	nt = len(u[:,0,0,0])
[184]	99
[186]	100	#########################################
	101	### Name for title and graphics save file
	102	basename = "WINDS"
	103	if var:
	104	basename = basename + "_" + var
[188]	105	if typefile is 'meso': stralt = "_lvl" + str(nvert)
	106	elif typefile is 'mesoapi':
	107	zelevel = int(nc.variables['vert'][nvert])
	108	if 'altitude_abg' in nc.dimensions: stralt = "_"+str(zelevel)+"m-ALS"
	109	elif 'bottom_top' in nc.dimensions: stralt = "_"+str(zelevel)+"m-AMR"
	110	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
	111	else: stralt = ""
	112	basename = basename + stralt
[186]	113
[184]	114	##################################
	115	### Open a figure and set subplots
[180]	116	fig = figure()
[186]	117	if numplot > 0:
	118	if numplot == 4:
	119	sub = 221
	120	fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
	121	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
	122	elif numplot == 2:
	123	sub = 121
	124	fig.subplots_adjust(wspace = 0.3)
	125	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
	126	elif numplot == 3:
	127	sub = 131
	128	fig.subplots_adjust(wspace = 0.3)
	129	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
	130	elif numplot == 6:
	131	sub = 231
	132	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
	133	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
	134	elif numplot == 8:
	135	sub = 331 #241
	136	fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
	137	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	138	elif numplot == 9:
	139	sub = 331
	140	fig.subplots_adjust(wspace = 0.1, hspace = 0.3)
	141	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
	142	elif numplot == 1:
	143	pass
	144	else:
	145	print "supported: 1,2,3,4,6,8"
	146	exit()
	147	### Prepare time loop
	148	if nt <= numplot or numplot == 1:
	149	tabrange = [0]
	150	numplot = 1
	151	else:
	152	tabrange = range(0,nt,int(nt/numplot)) #nt-1
	153	tabrange = tabrange[0:numplot]
	154	else:
	155	tabrange = range(0,nt,1)
	156	sub = 99999
	157	print tabrange
[184]	158
	159	#################################
	160	### Time loop for plotting device
[186]	161	found_lct = False
[184]	162	for i in tabrange:
	163
	164	### General plot settings
[186]	165	if numplot > 1:
	166	subplot(sub)
	167	found_lct = True
	168	elif numplot == 1:
	169	found_lct = True
[188]	170	### If only one local time is requested (numplot < 0)
[186]	171	elif numplot <= 0:
	172	if (ltst+i)%24 + numplot != 0: continue
	173	else: found_lct = True
[184]	174
	175	### Map projection
[180]	176	m = define_proj(proj,wlon,wlat,back=back)
	177	x, y = m(lon2d, lat2d)
[184]	178
	179	#### Contour plot
	180	if var:
	181	if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:])
	182	elif typefile in ['gcm']:
	183	if dimension == 2: what_I_plot = field[:,:]
	184	elif dimension == 3: what_I_plot = field[i,:,:]
	185	contourf(x, y, what_I_plot, 30)
[187]	186	if colorb: colorbar(fraction=0.05,pad=0.1)
[184]	187
	188	### Vector plot
	189	if typefile in ['mesoapi','meso']:
	190	[vecx,vecy] = [dumpbdy(u[i,nvert,:,:]), dumpbdy(v[i,nvert,:,:])]
	191	key = True
	192	elif typefile in ['gcm']:
	193	[vecx,vecy] = [ u[i,nvert,:,:] , v[i,nvert,:,:] ]
	194	key = False
	195	if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
[188]	196	if var == None and back == "vishires": colorvec = 'w'
	197	else: colorvec = 'k'
[184]	198	vectorfield(vecx, vecy,\
[186]	199	x, y, stride=stride, csmooth=stride,\
[188]	200	scale=15., factor=300., color=colorvec, key=key)
	201	#200. ## or csmooth=2
[184]	202
	203	### Next subplot
[186]	204	ptitle( basename + "_LT"+str((ltst+i)%24) )
[180]	205	sub += 1
[184]	206
	207	##########################################################################
	208	### Save the figure in a file in the data folder or an user-defined folder
[186]	209	if not target: zeplot = namefile+"_"+basename
	210	else: zeplot = target+"/"+basename
	211	if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot))
	212	if found_lct: makeplotpng(zeplot,pad_inches_value=0.35)
	213	else: print "Local time not found"
[180]	214
[184]	215
	216
	217	####################################################
	218	####################################################
	219	### A simple program to get wind vectors' components
	220	def getwinds (nc,charu='Um',charv='Vm'):
[182]	221	import numpy as np
[180]	222	u = nc.variables[charu]
	223	v = nc.variables[charv]
[182]	224	if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1]
	225	if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :]
[184]	226	### ou alors prendre les coordonnees speciales
[180]	227	return u,v
	228
[184]	229
	230
	231	###########################################################################################
	232	###########################################################################################
	233	### What is below relate to running the file as a command line executable (very convenient)
[180]	234	if __name__ == "__main__":
	235	import sys
[184]	236	from optparse import OptionParser ### to be replaced by argparse
[186]	237	from api_wrapper import api_onelevel
[182]	238	parser = OptionParser()
[188]	239	parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file')
	240	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]	241	parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection')
[182]	242	parser.add_option('-b', action='store', dest='back', type="string", default=None, help='background')
	243	parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder')
[186]	244	parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)')
[184]	245	parser.add_option('-v', action='store', dest='var', type="string", default=None, help='variable contoured')
[187]	246	parser.add_option('-n', action='store', dest='numplot', type="int", default=4, help='number of plots (def=1)(<0: 1 plot of LT -numplot)')
[188]	247	parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)')
[187]	248	parser.add_option('-c', action='store', dest='colorb', type="string", default=None, help='colorbar')
[182]	249	(opt,args) = parser.parse_args()
	250	if opt.namefile is None:
	251	print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h"
	252	exit()
	253	print "Options:", opt
[184]	254
[186]	255	zefile = opt.namefile
	256	zelevel = opt.nvert
[188]	257	stralt = None
[186]	258	if opt.interp is not None:
[187]	259	if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010
[186]	260	if opt.var is None : zefields = 'uvmet'
	261	else : zefields = 'uvmet,'+opt.var
	262	zefile = api_onelevel ( path_to_input = '', \
	263	input_name = zefile, \
	264	path_to_output = opt.target, \
	265	fields = zefields, \
	266	interp_method = opt.interp, \
	267	onelevel = zelevel )
[188]	268	zelevel = 0 ## so that zelevel could play again the role of nvert
[184]	269
[188]	270	winds (zefile,int(zelevel),\
	271	proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb)

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

Download in other formats: