Revision Log

billymitchell_sodar_plot.py

Revision 527 (checked in by cbc, 10 years ago)	Remove pcolor labels for w, sigw, bck.

Line
1	#!/opt/env/haines/dataproc/bin/python
2	# Last modified: Time-stamp: <2011-10-11 12:44:06 cbc>
3
4	"""billymitchell_sodar_plot"""
5
6	#################################################
7	# [1a] search and replace XXXX with platform id
8	# [1b] search and replace YYYY with sensor id
9	# [1c] search and replace ZZZZ plot titles
10	#################################################
11
12	# plot each month
13
14	import os
15	# import sys, glob, re
16	import datetime, time, dateutil, dateutil.tz
17	import pycdf
18	import numpy
19
20	# sys.path.append('/opt/env/haines/dataproc/raw2proc')
21	# del(sys)
22
23	os.environ["MPLCONFIGDIR"]="/home/haines/.matplotlib/"
24
25	from pylab import figure, twinx, savefig, setp, getp, cm, colorbar
26	from matplotlib.dates import DayLocator, HourLocator, MinuteLocator, DateFormatter, date2num, num2date
27	import procutil
28
29	def timeseries(pi, si, yyyy_mm, plot_type='latest'):
30	print 'billymitchell_sodar1_plot ...'
31
32	#
33
34	prev_month, this_month, next_month = procutil.find_months(yyyy_mm)
35	yyyy_mm_str = this_month.strftime('%Y_%m')
36
37	################################
38	# [2a] load primary data file
39	################################
40
41	ncFile1='/seacoos/data/nccoos/level2/billymitchell/sodar1/billymitchell_sfas_'+prev_month.strftime('%Y_%m')+'.nc'
42	ncFile2='/seacoos/data/nccoos/level2/billymitchell/sodar1/billymitchell_sfas_'+this_month.strftime('%Y_%m')+'.nc'
43
44	have_ncFile1 = os.path.exists(ncFile1)
45	have_ncFile2 = os.path.exists(ncFile2)
46
47	print ' ... loading data for graph from ...'
48	print ' ... ... ' + ncFile1 + ' ... ' + str(have_ncFile1)
49	print ' ... ... ' + ncFile2 + ' ... ' + str(have_ncFile2)
50
51	# open netcdf data
52	if have_ncFile1 and have_ncFile2:
53	try:
54	nc = pycdf.CDFMF((ncFile1, ncFile2))
55	except: # files may have different dimensions
56	nc = pycdf.CDFMF((ncFile2,))
57	elif not have_ncFile1 and have_ncFile2:
58	nc = pycdf.CDFMF((ncFile2,))
59	elif have_ncFile1 and not have_ncFile2:
60	nc = pycdf.CDFMF((ncFile1,))
61	else:
62	print ' ... both files do not exist -- NO DATA LOADED'
63	return
64
65	# ncvars = nc.variables()
66	# print ncvars
67	es = nc.var('time')[:]
68	units = nc.var('time').units
69	dt = [procutil.es2dt(e) for e in es]
70	# set timezone info to UTC (since data from level2 should be in UTC!!)
71	dt = [e.replace(tzinfo=dateutil.tz.tzutc()) for e in dt]
72	# return new datetime based on computer local
73	dt_local = [e.astimezone(dateutil.tz.tzlocal()) for e in dt]
74	dn = date2num(dt)
75
76
77	################################
78	# [2b] specify variables
79	################################
80	z = nc.var('z')[:]
81	# convert cm/s to m/s
82	uu = nc.var('u')[:]
83	vv = nc.var('v')[:]
84	ww = nc.var('w')[:]
85	sigw = nc.var('sigw')[:]
86	echo = nc.var('bck')[:]
87
88	nc.close()
89
90	# retain original dt for addnan
91	dto = dt
92
93	# last dt in data for labels
94	dtu = dt[-1]
95	dtl = dt_local[-1]
96
97	diff = abs(dtu - dtl)
98	if diff.days>0:
99	last_dt_str = dtu.strftime("%H:%M %Z on %b %d, %Y") + ' (' + dtl.strftime("%H:%M %Z, %b %d") + ')'
100	else:
101	last_dt_str = dtu.strftime("%H:%M %Z") + ' (' + dtl.strftime("%H:%M %Z") + ')' \
102	+ dtl.strftime(" on %b %d, %Y")
103
104	#######################################
105	# Plot setup
106	#######################################
107
108	fig = figure(figsize=(10, 10))
109	fig.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.1, hspace=0.1)
110
111	ax = fig.add_subplot(5,1,1)
112	axs = [ax]
113
114	# range for horizontal wind plots
115	cmin, cmax = (-20., 20.)
116	# print "%s : %g %g" % ('uv wind', cmin, cmax)
117	# use masked array to hide NaN's on plot
118	(dt, uu) = procutil.addnan(dto, uu)
119	dn = date2num(dt)
120	um = numpy.ma.masked_where(numpy.isnan(uu), uu)
121	pc = ax.pcolor(dn, z, um.T, vmin=cmin, vmax=cmax)
122	pc.set_label('True Eastward Wind (m s-1)')
123	ax.text(0.025, 0.1, pc.get_label(), fontsize="small", transform=ax.transAxes)
124	ax.set_ylabel('Height (m)')
125
126	# setup colorbar axes instance.
127	l,b,w,h = ax.get_position().bounds
128	cax = fig.add_axes([l, b+h+0.04, 0.25*w, 0.03])
129
130	cb = colorbar(pc, cax=cax, orientation='horizontal') # draw colorbar
131	cb.set_label('Wind Velocity (m s-1)')
132	cb.ax.xaxis.set_label_position('top')
133	cb.ax.set_xticks([0.1, 0.3, 0.5, 0.7, 0.9])
134	xtl = numpy.round(numpy.linspace(cmin, cmax, 10), decimals=0)
135	cb.ax.set_xticklabels([xtl[1], xtl[3], xtl[5], xtl[7], xtl[9]])
136
137	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
138	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
139	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
140	ax.set_xticklabels([])
141
142	# this only moves the label not the tick labels
143	ax.xaxis.set_label_position('top')
144	ax.set_xlabel('BILLY MITCHELL SODAR -- ' + yyyy_mm_str)
145
146	# right-hand side scale
147	ax2 = twinx(ax)
148	ax2.yaxis.tick_right()
149	# convert (lhs) meters to (rhs) feet
150	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
151	ax2.set_ylim(feet)
152	ax2.set_ylabel('Height (ft)')
153
154	#######################################
155	#
156	ax = fig.add_subplot(5,1,2)
157	axs.append(ax)
158
159	# use masked array to hide NaN's on plot
160	(dt, vv) = procutil.addnan(dto, vv)
161	dn = date2num(dt)
162	vm = numpy.ma.masked_where(numpy.isnan(vv), vv)
163	pc = ax.pcolor(dn, z, vm.T, vmin=cmin, vmax=cmax)
164	pc.set_label('True Northward Wind (m s-1)')
165	ax.text(0.025, 0.1, pc.get_label(), fontsize="small", transform=ax.transAxes)
166	ax.set_ylabel('Height (m)')
167
168	# ax.set_xlim(date2num(dt[0]), date2num(dt[-1]))
169	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
170	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
171	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
172	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
173	ax.set_xticklabels([])
174
175	# right-hand side scale
176	ax2 = twinx(ax)
177	ax2.yaxis.tick_right()
178	# convert (lhs) meters to (rhs) feet
179	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
180	ax2.set_ylim(feet)
181	ax2.set_ylabel('Height (ft)')
182
183	#######################################
184	#
185	ax = fig.add_subplot(5,1,3)
186	axs.append(ax)
187
188	# range for horizontal wind plots
189	cmin, cmax = (-2., 2.)
190	# print "%s : %g %g" % ('w wind', cmin, cmax)
191
192	# use masked array to hide NaN's on plot
193	(dt, ww) = procutil.addnan(dto, ww)
194	dn = date2num(dt)
195	wm = numpy.ma.masked_where(numpy.isnan(ww), ww)
196	pc = ax.pcolor(dn, z, wm.T, vmin=cmin, vmax=cmax)
197	ax.set_ylabel('Height (m)')
198
199	# setup colorbar axes instance.
200	l,b,w,h = ax.get_position().bounds
201	cax = fig.add_axes([l+0.04, b+h-0.06, 0.25*w, 0.03])
202
203	cb = colorbar(pc, cax=cax, orientation='horizontal') # draw colorbar
204	cb.set_label('Upward Wind Velocity (m s-1)')
205	cb.ax.xaxis.set_label_position('top')
206	cb.ax.set_xticks([0.1, 0.3, 0.5, 0.7, 0.9])
207	xtl = numpy.round(numpy.linspace(cmin, cmax, 10), decimals=0)
208	cb.ax.set_xticklabels([xtl[1], xtl[3], xtl[5], xtl[7], xtl[9]])
209
210	# ax.set_xlim(date2num(dt[0]), date2num(dt[-1]))
211	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
212	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
213	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
214	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
215	ax.set_xticklabels([])
216
217	# right-hand side scale
218	ax2 = twinx(ax)
219	ax2.yaxis.tick_right()
220	# convert (lhs) meters to (rhs) feet
221	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
222	ax2.set_ylim(feet)
223	ax2.set_ylabel('Height (ft)')
224
225	#######################################
226	#
227	ax = fig.add_subplot(5,1,4)
228	axs.append(ax)
229
230	cmin, cmax = (0., 2.)
231	# print "%s : %g %g" % ('sigw', cmin, cmax)
232
233	# use masked array to hide NaN's on plot
234	(dt, sigw) = procutil.addnan(dto, sigw)
235	dn = date2num(dt)
236	sm = numpy.ma.masked_where(numpy.isnan(sigw), sigw)
237	pc = ax.pcolor(dn, z, sm.T, vmin=cmin, vmax=cmax)
238	ax.set_ylabel('Height (m)')
239
240	# setup colorbar axes instance.
241	l,b,w,h = ax.get_position().bounds
242	cax = fig.add_axes([l+0.04, b+h-0.06, 0.25*w, 0.03])
243
244	cb = colorbar(pc, cax=cax, orientation='horizontal') # draw colorbar
245	cb.set_label('Std. Dev. of Vertical Wind (m s-1)')
246	cb.ax.xaxis.set_label_position('top')
247	cb.ax.set_xticks([0.0, 0.2, 0.4, 0.6, 0.8])
248	xtl = numpy.round(numpy.linspace(cmin, cmax, 10), decimals=1)
249	cb.ax.set_xticklabels([xtl[0], xtl[2], xtl[4], xtl[6], xtl[8]])
250
251	# ax.set_xlim(date2num(dt[0]), date2num(dt[-1]))
252	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
253	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
254	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
255	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
256	ax.set_xticklabels([])
257
258	# right-hand side scale
259	ax2 = twinx(ax)
260	ax2.yaxis.tick_right()
261	# convert (lhs) meters to (rhs) feet
262	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
263	ax2.set_ylim(feet)
264	ax2.set_ylabel('Height (ft)')
265
266	#######################################
267	#
268	ax = fig.add_subplot(5,1,5)
269	axs.append(ax)
270
271
272	# cmin, cmax = (0., 1000.)
273	cmin, cmax = (2., 6.)
274	# print "%s : %g %g" % ('echo', cmin, cmax)
275
276	# use masked array to hide NaN's on plot
277	(dt, echo) = procutil.addnan(dto, numpy.log10(echo))
278	dn = date2num(dt)
279	em = numpy.ma.masked_where(numpy.isnan(echo), echo)
280	pc = ax.pcolor(dn, z, em.T, vmin=cmin, vmax=cmax)
281	ax.set_ylabel('Height (m)')
282
283	# setup colorbar axes instance.
284	l,b,w,h = ax.get_position().bounds
285	cax = fig.add_axes([l+0.04, b+h-0.06, 0.25*w, 0.03])
286
287	cb = colorbar(pc, cax=cax, orientation='horizontal') # draw colorbar
288	cb.set_label('log10 Backscatter')
289	cb.ax.xaxis.set_label_position('top')
290	cb.ax.set_xticks([0.0, 1.0])
291	xtl = numpy.round(numpy.linspace(cmin, cmax, 10), decimals=0)
292	cb.ax.set_xticklabels([str(cmin), str(cmax)])
293
294	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
295	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
296	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
297	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
298	ax.set_xticklabels([])
299
300	# right-hand side scale
301	ax2 = twinx(ax)
302	ax2.yaxis.tick_right()
303	# convert (lhs) meters to (rhs) feet
304	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
305	ax2.set_ylim(feet)
306	ax2.set_ylabel('Height (ft)')
307
308	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
309	ax.set_xlabel('BILLY MITCHELL SODAR -- ' + yyyy_mm_str)
310	# save figure
311
312	# save figure for this month
313	ofn = '/home/haines/rayleigh/img/billymitchell/billymitchell_sodar1_'+yyyy_mm_str+'.png'
314	print '... ... write: %s' % (ofn,)
315	savefig(ofn)
316
317
318	#######################################
319	# Last 30 days
320	#######################################
321	if plot_type=='latest':
322	print ' ... Last 30 days'
323	for idx, ax in enumerate(axs):
324	ax.set_xlim(date2num(dt[-1])-30, date2num(dt[-1]))
325	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
326	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
327	ax.set_xticklabels([])
328	if idx==0:
329	ax.set_xlabel('BILLY MITCHELL SODAR -- Last 30 days from ' + last_dt_str)
330	elif idx==len(axs)-1:
331	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
332	ax.set_xlabel('BILLY MITCHELL SODAR -- Last 30 days from ' + last_dt_str)
333
334	savefig('/home/haines/rayleigh/img/billymitchell_sodar1_last30days.png')
335
336
337	#######################################
338	# Last 7 days
339	#######################################
340	if plot_type=='latest':
341	print ' ... Last 7 days'
342	for idx, ax in enumerate(axs):
343	ax.set_xlim(date2num(dt[-1])-7, date2num(dt[-1]))
344	ax.xaxis.set_major_locator( DayLocator(range(0,32,1)) )
345	ax.xaxis.set_minor_locator( HourLocator(range(0,25,6)) )
346	ax.set_xticklabels([])
347	if idx==0:
348	ax.set_xlabel('BILLY MITCHELL SODAR -- Last 7 days from ' + last_dt_str)
349	elif idx==len(axs)-1:
350	ax.xaxis.set_major_formatter( DateFormatter('%m/%d') )
351	ax.set_xlabel('BILLY MITCHELL SODAR -- Last 7 days from ' + last_dt_str)
352
353	savefig('/home/haines/rayleigh/img/billymitchell_sodar1_last07days.png')
354
355
356	#######################################
357	# Last 1 day (24hrs)
358	#######################################
359	if plot_type=='latest':
360	print ' ... Last 1 days'
361
362	for idx, ax in enumerate(axs):
363	ax.set_xlim(date2num(dt[-1])-1, date2num(dt[-1]))
364	ax.xaxis.set_major_locator( HourLocator(range(0,25,1)) )
365	ax.xaxis.set_minor_locator( MinuteLocator(range(0,61,30)) )
366	ax.set_xticklabels([])
367	if idx==0:
368	ax.set_xlabel('BILLY MITCHELL SODAR -- Last 24 hours from ' + last_dt_str)
369	elif idx==len(axs)-1:
370	ax.xaxis.set_major_formatter( DateFormatter('%H') )
371	ax.set_xlabel('BILLY MITCHELL SODAR -- Last 24 hours from ' + last_dt_str)
372
373	savefig('/home/haines/rayleigh/img/billymitchell_sodar1_last01days.png')
374
375
376	def wind_vectors(pi, si, yyyy_mm, plot_type='latest'):
377	print 'billymitchell_arrows_plot ...'
378
379	#
380
381	prev_month, this_month, next_month = procutil.find_months(yyyy_mm)
382	yyyy_mm_str = this_month.strftime('%Y_%m')
383
384	################################
385	# [2a] load primary data file
386	################################
387
388	ncFile1='/seacoos/data/nccoos/level2/billymitchell/sodar1/billymitchell_sfas_'+prev_month.strftime('%Y_%m')+'.nc'
389	ncFile2='/seacoos/data/nccoos/level2/billymitchell/sodar1/billymitchell_sfas_'+this_month.strftime('%Y_%m')+'.nc'
390
391	have_ncFile1 = os.path.exists(ncFile1)
392	have_ncFile2 = os.path.exists(ncFile2)
393
394	print ' ... loading data for graph from ...'
395	print ' ... ... ' + ncFile1 + ' ... ' + str(have_ncFile1)
396	print ' ... ... ' + ncFile2 + ' ... ' + str(have_ncFile2)
397
398	# open netcdf data
399	if have_ncFile1 and have_ncFile2:
400	try:
401	nc = pycdf.CDFMF((ncFile1, ncFile2))
402	except: # files may have different dimensions
403	nc = pycdf.CDFMF((ncFile2,))
404	elif not have_ncFile1 and have_ncFile2:
405	nc = pycdf.CDFMF((ncFile2,))
406	elif have_ncFile1 and not have_ncFile2:
407	nc = pycdf.CDFMF((ncFile1,))
408	else:
409	print ' ... both files do not exist -- NO DATA LOADED'
410	return
411
412	# ncvars = nc.variables()
413	# print ncvars
414	es = nc.var('time')[:]
415	units = nc.var('time').units
416	dt = [procutil.es2dt(e) for e in es]
417	# set timezone info to UTC (since data from level2 should be in UTC!!)
418	dt = [e.replace(tzinfo=dateutil.tz.tzutc()) for e in dt]
419	# return new datetime based on computer local
420	dt_local = [e.astimezone(dateutil.tz.tzlocal()) for e in dt]
421	dn = date2num(dt)
422
423
424	################################
425	# [2b] specify variables
426	################################
427	z = nc.var('z')[:]
428	# convert cm/s to m/s
429	uu = nc.var('u')[:]
430	vv = nc.var('v')[:]
431	ww = nc.var('w')[:]
432	sigw = nc.var('sigw')[:]
433	echo = nc.var('bck')[:]
434
435	nc.close()
436
437	# retain original dt for addnan
438	dto = dt
439
440	# last dt in data for labels
441	dtu = dt[-1]
442	dtl = dt_local[-1]
443
444	diff = abs(dtu - dtl)
445	if diff.days>0:
446	last_dt_str = dtu.strftime("%H:%M %Z on %b %d, %Y") + ' (' + dtl.strftime("%H:%M %Z, %b %d") + ')'
447	else:
448	last_dt_str = dtu.strftime("%H:%M %Z") + ' (' + dtl.strftime("%H:%M %Z") + ')' \
449	+ dtl.strftime(" on %b %d, %Y")
450
451	#######################################
452	# Plot setup
453	#######################################
454
455	fig = figure(figsize=(10, 10))
456	fig.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.1, hspace=0.1)
457
458	ax = fig.add_subplot(1,1,1)
459	axs = [ax]
460
461	# range for horizontal wind plots
462	# cmin, cmax = (0., 20.)
463	# print "%s : %g %g" % ('uv wind', cmin, cmax)
464	# use masked array to hide NaN's on plot
465	(dt, uu) = procutil.addnan(dto, uu)
466	dn = date2num(dt)
467	um = numpy.ma.masked_where(numpy.isnan(uu), uu)
468
469	(dt, vv) = procutil.addnan(dto, vv)
470	vm = numpy.ma.masked_where(numpy.isnan(vv), vv)
471	wspd = numpy.sqrt(umum + vmvm)
472
473	X,Y = numpy.meshgrid(dn, z)
474
475	q1 = ax.quiver(X, Y, um.T, vm.T, wspd.T, units='inches', scale=40)
476	qk = ax.quiverkey(q1, 0.1, 0.8, 20, r'20 m s-1')
477
478	ax.set_ylabel('Height (m)')
479
480	# setup colorbar axes instance.
481	l,b,w,h = ax.get_position().bounds
482	cax = fig.add_axes([l+0.02, b+h-0.06, 0.25*w, 0.03])
483
484	cb = colorbar(q1, cax=cax, orientation='horizontal') # draw colorbar
485	cb.set_label('Wind Velocity (m s-1)')
486
487	# lost control of colorbar ... ???
488	#
489	# cb.ax.xaxis.set_label_position('top')
490	# cb.ax.set_xticks([0.0, 0.5, 1.0])
491	# xtl = numpy.round(numpy.linspace(cmin, cmax, 11), decimals=0)
492	# cb.ax.set_xticklabels([xtl[0], xtl[5], xtl[10]])
493
494	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
495	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
496	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
497	ax.set_xticklabels([])
498
499	# this only moves the label not the tick labels
500	ax.xaxis.set_label_position('bottom')
501	ax.set_xlabel('BILLY MITCHELL Wind Profile -- ' + yyyy_mm_str)
502
503	# right-hand side scale
504	ax2 = twinx(ax)
505	ax2.yaxis.tick_right()
506	# convert (lhs) meters to (rhs) feet
507	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
508	ax2.set_ylim(feet)
509	ax2.set_ylabel('Height (ft)')
510
511
512
513	#######################################
514	# Last 1 day (24hrs)
515	#######################################
516	if plot_type=='latest':
517	print ' ... Last 1 days'
518
519	for idx, ax in enumerate(axs):
520	ax.set_xlim(date2num(dt[-1])-1, date2num(dt[-1]))
521	ax.xaxis.set_major_locator( HourLocator(range(0,25,1)) )
522	ax.xaxis.set_minor_locator( MinuteLocator(range(0,61,30)) )
523	ax.set_xticklabels([])
524	if idx==0:
525	ax.set_xlabel('BILLY MITCHELL Wind Profile -- Last 24 hours from ' + last_dt_str)
526	if idx==len(axs)-1:
527	ax.xaxis.set_major_formatter( DateFormatter('%H') )
528	ax.set_xlabel('BILLY MITCHELL Wind Profile -- Last 24 hours from ' + last_dt_str)
529
530	savefig('/home/haines/rayleigh/img/billymitchell_wind_last01days.png')
531
532
533	def wind_barbs(pi, si, yyyy_mm, plot_type='latest'):
534	print 'billymitchell_arrows_plot ...'
535
536	#
537
538	prev_month, this_month, next_month = procutil.find_months(yyyy_mm)
539	yyyy_mm_str = this_month.strftime('%Y_%m')
540
541	################################
542	# [2a] load primary data file
543	################################
544
545	ncFile1='/seacoos/data/nccoos/level2/billymitchell/sodar1/billymitchell_sfas_'+prev_month.strftime('%Y_%m')+'.nc'
546	ncFile2='/seacoos/data/nccoos/level2/billymitchell/sodar1/billymitchell_sfas_'+this_month.strftime('%Y_%m')+'.nc'
547
548	have_ncFile1 = os.path.exists(ncFile1)
549	have_ncFile2 = os.path.exists(ncFile2)
550
551	print ' ... loading data for graph from ...'
552	print ' ... ... ' + ncFile1 + ' ... ' + str(have_ncFile1)
553	print ' ... ... ' + ncFile2 + ' ... ' + str(have_ncFile2)
554
555	# open netcdf data
556	if have_ncFile1 and have_ncFile2:
557	try:
558	nc = pycdf.CDFMF((ncFile1, ncFile2))
559	except: # files may have different dimensions
560	nc = pycdf.CDFMF((ncFile2,))
561	elif not have_ncFile1 and have_ncFile2:
562	nc = pycdf.CDFMF((ncFile2,))
563	elif have_ncFile1 and not have_ncFile2:
564	nc = pycdf.CDFMF((ncFile1,))
565	else:
566	print ' ... both files do not exist -- NO DATA LOADED'
567	return
568
569	# ncvars = nc.variables()
570	# print ncvars
571	es = nc.var('time')[:]
572	units = nc.var('time').units
573	dt = [procutil.es2dt(e) for e in es]
574	# set timezone info to UTC (since data from level2 should be in UTC!!)
575	dt = [e.replace(tzinfo=dateutil.tz.tzutc()) for e in dt]
576	# return new datetime based on computer local
577	dt_local = [e.astimezone(dateutil.tz.tzlocal()) for e in dt]
578	dn = date2num(dt)
579
580
581	################################
582	# [2b] specify variables
583	################################
584	z = nc.var('z')[:]
585	# convert cm/s to m/s
586	uu = nc.var('u')[:]
587	vv = nc.var('v')[:]
588	ww = nc.var('w')[:]
589	sigw = nc.var('sigw')[:]
590	echo = nc.var('bck')[:]
591
592	nc.close()
593
594	# retain original dt for addnan
595	dto = dt
596
597	# last dt in data for labels
598	dtu = dt[-1]
599	dtl = dt_local[-1]
600
601	diff = abs(dtu - dtl)
602	if diff.days>0:
603	last_dt_str = dtu.strftime("%H:%M %Z on %b %d, %Y") + ' (' + dtl.strftime("%H:%M %Z, %b %d") + ')'
604	else:
605	last_dt_str = dtu.strftime("%H:%M %Z") + ' (' + dtl.strftime("%H:%M %Z") + ')' \
606	+ dtl.strftime(" on %b %d, %Y")
607
608	#######################################
609	# Plot setup
610	#######################################
611
612	fig = figure(figsize=(10, 10))
613	fig.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.1, hspace=0.1)
614
615	ax = fig.add_subplot(1,1,1)
616	axs = [ax]
617
618	# range for horizontal wind plots
619	# cmin, cmax = (0., 20.)
620	# print "%s : %g %g" % ('uv wind', cmin, cmax)
621	# use masked array to hide NaN's on plot
622	(dt, uu) = procutil.addnan(dto, uu)
623	dn = date2num(dt)
624	um = numpy.ma.masked_where(numpy.isnan(uu), uu)
625
626	(dt, vv) = procutil.addnan(dto, vv)
627	vm = numpy.ma.masked_where(numpy.isnan(vv), vv)
628	wspd = numpy.sqrt(umum + vmvm)
629
630	X,Y = numpy.meshgrid(dn, z)
631
632	q1 = ax.barbs(X, Y, um.T, vm.T, wspd.T)
633	# qk = ax.quiverkey(q1, 0.1, 0.8, 20, r'20 m s-1')
634
635	ax.set_ylabel('Height (m)')
636
637	# setup colorbar axes instance.
638	l,b,w,h = ax.get_position().bounds
639	cax = fig.add_axes([l+0.02, b+h-0.06, 0.25*w, 0.03])
640
641	cb = colorbar(q1, cax=cax, orientation='horizontal') # draw colorbar
642	cb.set_label('Wind Velocity (m s-1)')
643
644	# lost control of colorbar ... ???
645	#
646	# cb.ax.xaxis.set_label_position('top')
647	# cb.ax.set_xticks([0.0, 0.5, 1.0])
648	# xtl = numpy.round(numpy.linspace(cmin, cmax, 11), decimals=0)
649	# cb.ax.set_xticklabels([xtl[0], xtl[5], xtl[10]])
650
651	ax.set_xlim(date2num(this_month), date2num(next_month-datetime.timedelta(seconds=1)))
652	ax.xaxis.set_major_locator( DayLocator(range(2,32,2)) )
653	ax.xaxis.set_minor_locator( HourLocator(range(0,25,12)) )
654	ax.set_xticklabels([])
655
656	# this only moves the label not the tick labels
657	ax.xaxis.set_label_position('bottom')
658	ax.set_xlabel('BILLY MITCHELL Wind Profile -- ' + yyyy_mm_str)
659
660	# right-hand side scale
661	ax2 = twinx(ax)
662	ax2.yaxis.tick_right()
663	# convert (lhs) meters to (rhs) feet
664	feet = [procutil.meters2feet(val) for val in ax.get_ylim()]
665	ax2.set_ylim(feet)
666	ax2.set_ylabel('Height (ft)')
667
668
669
670	#######################################
671	# Last 1 day (24hrs)
672	#######################################
673	if plot_type=='latest':
674	print ' ... Last 1 days'
675
676	for idx, ax in enumerate(axs):
677	ax.set_xlim(date2num(dt[-1])-1, date2num(dt[-1]))
678	ax.xaxis.set_major_locator( HourLocator(range(0,25,1)) )
679	ax.xaxis.set_minor_locator( MinuteLocator(range(0,61,30)) )
680	ax.set_xticklabels([])
681	if idx==0:
682	ax.set_xlabel('BILLY MITCHELL Wind Profile -- Last 24 hours from ' + last_dt_str)
683	if idx==len(axs)-1:
684	ax.xaxis.set_major_formatter( DateFormatter('%H') )
685	ax.set_xlabel('BILLY MITCHELL Wind Profile -- Last 24 hours from ' + last_dt_str)
686
687	savefig('/home/haines/rayleigh/img/billymitchell_windbarbs_last01days.png')
688

Note: See TracBrowser for help on using the browser.

root/proc2plot/trunk/proc2plot/billymitchell_sodar_plot.py

Download in other formats: