Context Navigation

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

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

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

Download in other formats: