Revision Log

winddist.py

Revision 85 (checked in by cbc, 11 years ago)	Fix directionality of plot labels.

Line
1	import sys,os,datetime,glob
2	import multiprocessing
3	import numpy as np
4	import matplotlib as mpl
5	mpl.use("Agg")
6	from windrose import WindroseAxes
7	from matplotlib import pyplot as plt
8	import pycdf
9
10	ncDir = '/seacoos/data/nccoos/level1/billymitchell/sodar1'
11	ncFileGlob = 'billymitchell_sfas_*.nc'
12
13	pngDir = '/seacoos/data/nccoos/level1/billymitchell/sodar1/plots/winddist'
14	if not os.path.exists(pngDir):
15	os.mkdir(pngDir)
16
17	pngExt = 'png'
18	htmlExt = 'html'
19
20	ncFilePattern = os.path.join(ncDir,ncFileGlob)
21	files = glob.glob(ncFilePattern)
22	files.sort()
23	# files=files[-3:]
24	previous = [files[-1]] + files[:-1]
25	next = files[1:] + [files[0]]
26	files = zip(previous,files,next)
27
28	html1 = """<html>
29	<head>
30	<title>
31	Billy Mitchell Sodar :: """
32
33
34	html2 = """ :: Wind Distribution at """
35
36	html3 = """ Meters Above Sea Level
37	</title>
38	</head>
39	<body>
40	<form>
41	<a href=\""""
42
43	html4 = """\"><Previous Month</a>
44	<a href=\""""
45
46	html5 = """\">Next Month></a>
47	<a href=\""""
48
49	html6 = """\"><Previous Elevation</a>
50	<a href=\""""
51
52	html7 = """\">Next Elevation></a>
53	</form>
54	<br/>
55	<img src=\""""
56
57	html8 = """\"/>
58	</body>
59	</html>"""
60
61	def _components(fileSpec):
62	"Figure out year and month from file specification."
63
64	fileName = os.path.splitext(os.path.basename(fileSpec))[0]
65	year = fileName[-7:-3]
66	month = fileName[-2:]
67	monthName = datetime.datetime(int(year),int(month),1).strftime('%B')
68	return (fileName, year, month, monthName)
69
70	def winddist((previous, ncFile, next), pngDir, genHtml=False, lastMonth=False):
71	print 'Processing',ncFile
72	ncFileName, year, month, monthName = _components(ncFile)
73	previous, previousYear, previousMonth, previousMonthName = _components(previous)
74	next, nextYear, nextMonth, nextMonthName = _components(next)
75
76	nc = pycdf.CDF(ncFile)
77
78	t = nc.var('time')[:]
79	z = nc.var('z')[:]
80	u = nc.var('u')[:]
81	v = nc.var('v')[:]
82	w = nc.var('w')[:]
83
84	mu = np.ma.masked_invalid(u)
85	mv = np.ma.masked_invalid(v)
86	mw = np.ma.masked_invalid(w)
87
88	mumean = mu.mean(axis=0)
89	mvmean = mv.mean(axis=0)
90	mwmean = mw.mean(axis=0)
91
92	ucount = mu.count(axis=0)
93	vcount = mv.count(axis=0)
94	wcount = mw.count(axis=0)
95
96	mrho = np.sqrt((mu2) + (mv2))
97	mrho = mrho.T
98	mtheta = (180 * np.arctan2(mv, mu)/np.pi)
99	theta = np.piecewise(mtheta,
100	(mtheta <= 90, mtheta > 90),
101	(lambda x: 90 - x, lambda x: 450 - x))
102	mtheta = np.ma.array(theta,mask=mtheta.mask)
103	mtheta = mtheta.T
104
105	rho = [np.array([elem
106	for elem,mask
107	in zip(row.data,row.mask)
108	if not mask])
109	for row
110	in mrho]
111	theta = [np.array([elem
112	for elem,mask
113	in zip(row.data,row.mask)
114	if not mask])
115	for row
116	in mtheta]
117
118	firstElevation = True
119	indices = range(len(rho))
120	previousX = [indices[-1]] + indices[:-1]
121	nextX = indices[1:] + [indices[0]]
122	indices = zip(previousX,indices,nextX)
123	for previousX,x,nextX in indices:
124	if rho[x].any():
125	fig = plt.figure(figsize=(8, 8), dpi=80, facecolor='w', edgecolor='w')
126	spt = fig.suptitle('Billy Mitchell Sodar :: ' + monthName + " " + year)
127	rect = [0.16, 0.16, 0.68, 0.68]
128	ax = WindroseAxes(fig, rect, axisbg='w')
129	ax.set_title('Wind Distribution at %d Meters Above Sea Level\nPercentage of Wind Magnitude Towards Direction\n\n' % z[x],fontsize=spt.get_fontsize()*0.8)
130	fig.add_axes(ax)
131	ax.bar(theta[x], rho[x], bins = np.arange(0.0,20.0,2.5), normed=True, opening=0.8, edgecolor='white')
132	yticklabels = ['',] + \
133	[" %s%%" % plt.getp(textobj,'text')
134	for textobj in plt.getp(ax,'yticklabels')[1:]]
135	plt.setp(ax,yticklabels=yticklabels)
136	l = ax.legend(axespad=-0.20,title="Magnitude (m/s)")
137	plt.setp(l.get_texts(), fontsize=8)
138	dirName = "%4u_%02u" % (int(year), int(month))
139	if not os.path.exists(os.path.join(pngDir, dirName)):
140	os.mkdir(os.path.join(pngDir, dirName))
141	outFile = os.path.join(pngDir,dirName,('%03um' % z[x]) + os.extsep + pngExt)
142	print 'Saving', outFile
143	fig.savefig(outFile)
144	fig.clear()
145
146	if genHtml:
147	htmlFile = os.path.join(pngDir,ncFileName + "_" + ('%dm' % z[x]) + os.extsep + htmlExt)
148	html = html1 + monthName + " " + year
149	html = html + html2 + ('%dm' % z[x])
150	html = html + html3 + previous + "_" + ('%dm' % z[x]) + os.extsep + htmlExt
151	html = html + html4 + next + "_" + ('%dm' % z[x]) + os.extsep + htmlExt
152	html = html + html5 + ncFileName + "_" + ('%dm' % z[previousX]) + os.extsep + htmlExt
153	html = html + html6 + ncFileName + "_" + ('%dm' % z[nextX]) + os.extsep + htmlExt
154	html = html + html7 + os.path.join(dirName,os.path.basename(outFile)) + html8
155	handle = open(htmlFile,'w')
156	handle.write(html)
157	handle.close()
158
159	if lastMonth and firstElevation:
160	firstElevation = False
161	htmlFile = os.path.join(pngDir,"index" + os.extsep + htmlExt)
162	handle = open(htmlFile,'w')
163	handle.write(html)
164	handle.close()
165
166	if __name__ == "__main__":
167	lastMonth = False
168	for previous,ncFile,next in files:
169	if next <= ncFile:
170	lastMonth = True
171	p = multiprocessing.Process(target=winddist,
172	args=((previous, ncFile, next),
173	pngDir,
174	True,
175	lastMonth))
176	p.start()
177	p.join()
178	if p.exitcode:
179	print "Exitcode %s from processing %s" % (p.exitcode, ncFile)

Note: See TracBrowser for help on using the browser.

root/sodarplot/trunk/sodarplot/scintec/winddist.py

Download in other formats: