NCCOOS Trac Projects: Top | Web | Platforms | Processing | Viz | Sprints | Sandbox | (Wind)

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

Revision 70 (checked in by cbc, 13 years ago)

Add wind distribution script utilizing windrose.

Line 
1 import os,datetime,glob
2 import numpy as np
3 import matplotlib as mpl
4 mpl.use("Agg")
5 from windrose import WindroseAxes
6 from matplotlib import pyplot as plt
7 import pycdf
8
9 ncDir = '/seacoos/data/nccoos/level1/billymitchell/sodar1'
10 ncFileGlob = 'billymitchell_sfas_*.nc'
11  
12 pngDir = '/seacoos/data/nccoos/level1/billymitchell/sodar1/plots/winddist'
13 if not os.path.exists(pngDir):
14     os.mkdir(pngDir)
15
16 pngExt = 'png'
17 htmlExt = 'html'
18  
19 ncFilePattern = os.path.join(ncDir,ncFileGlob)
20 files = glob.glob(ncFilePattern)
21 previous = [files[-1]] + files[:-1]
22 next = files[1:] + [files[0]]
23 files = zip(previous,files,next)
24  
25 first = True
26 for previous,ncFile,next in files:
27     print 'Processing',ncFile
28     ncFileName = os.path.splitext(os.path.basename(ncFile))[0]
29     year = ncFileName[19:23]
30     month = ncFileName[24:26]
31     month = datetime.datetime(int(year),int(month),1).strftime('%B')
32     if previous:
33         previous = os.path.splitext(os.path.basename(previous))[0]
34         previousYear = previous[19:23]
35         previousMonth = previous[24:26]
36         previousMonth = datetime.datetime(int(previousYear),int(previousMonth),1).strftime('%B')
37         previous = previous + os.extsep + htmlExt
38     if next:
39         next = os.path.splitext(os.path.basename(next))[0]
40         nextYear = next[19:23]
41         nextMonth = next[24:26]
42         nextMonth = datetime.datetime(int(nextYear),int(nextMonth),1).strftime('%B')
43         next = next + os.extsep + htmlExt
44  
45     nc = pycdf.CDF(ncFile)
46
47     t = nc.var('time')[:]
48     z = nc.var('z')[:]
49     u = nc.var('u')[:]
50     v = nc.var('v')[:]
51     w = nc.var('w')[:]
52  
53     mu = np.ma.masked_invalid(u)
54     mv = np.ma.masked_invalid(v)
55     mw = np.ma.masked_invalid(w)
56  
57     mumean = mu.mean(axis=0)
58     mvmean = mv.mean(axis=0)
59     mwmean = mw.mean(axis=0)
60  
61     ucount = mu.count(axis=0)
62     vcount = mv.count(axis=0)
63     wcount = mw.count(axis=0)
64
65     mrho = np.sqrt((mu**2) + (mv**2))
66     rho = mrho.T.data
67     mtheta = (180 * np.arctan2(mv, mu)/np.pi)
68     theta = np.piecewise(mtheta,
69                          (mtheta <= 90, mtheta > 90),
70                          (lambda x: 90 - x, lambda x: 450 - x))
71     theta = theta.T
72
73     for x in range(0,rho.shape[0]):
74         if rho[x].any():
75             fig = plt.figure(figsize=(8, 8), dpi=80, facecolor='w', edgecolor='w')
76             rect = [0.1, 0.1, 0.8, 0.8]
77             ax = WindroseAxes(fig, rect, axisbg='w')
78             fig.add_axes(ax)
79             ax.bar(theta[x], rho[x], normed=True, opening=0.8, edgecolor='white')
80             l = ax.legend(axespad=-0.10)
81             plt.setp(l.get_texts(), fontsize=8)
82             if not os.path.exists(os.path.join(pngDir,ncFileName)):
83                 os.mkdir(os.path.join(pngDir,ncFileName))
84             outFile = os.path.join(pngDir,ncFileName,('%dm' % z[x]) + os.extsep + pngExt)
85             print 'Saving', outFile
86             fig.savefig(outFile)
87
Note: See TracBrowser for help on using the browser.