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 |
files = files[:-1] # sodar was broken last month |
---|
22 |
previous = [files[-1]] + files[:-1] |
---|
23 |
next = files[1:] + [files[0]] |
---|
24 |
files = zip(previous,files,next) |
---|
25 |
|
---|
26 |
firstMonth = True |
---|
27 |
for previous,ncFile,next in files: |
---|
28 |
print 'Processing',ncFile |
---|
29 |
ncFileName = os.path.splitext(os.path.basename(ncFile))[0] |
---|
30 |
year = ncFileName[19:23] |
---|
31 |
month = ncFileName[24:26] |
---|
32 |
month = datetime.datetime(int(year),int(month),1).strftime('%B') |
---|
33 |
if previous: |
---|
34 |
previous = os.path.splitext(os.path.basename(previous))[0] |
---|
35 |
previousYear = previous[19:23] |
---|
36 |
previousMonth = previous[24:26] |
---|
37 |
previousMonth = datetime.datetime(int(previousYear),int(previousMonth),1).strftime('%B') |
---|
38 |
previous = previous + os.extsep + htmlExt |
---|
39 |
if next: |
---|
40 |
next = os.path.splitext(os.path.basename(next))[0] |
---|
41 |
nextYear = next[19:23] |
---|
42 |
nextMonth = next[24:26] |
---|
43 |
nextMonth = datetime.datetime(int(nextYear),int(nextMonth),1).strftime('%B') |
---|
44 |
next = next + os.extsep + htmlExt |
---|
45 |
|
---|
46 |
nc = pycdf.CDF(ncFile) |
---|
47 |
|
---|
48 |
t = nc.var('time')[:] |
---|
49 |
z = nc.var('z')[:] |
---|
50 |
u = nc.var('u')[:] |
---|
51 |
v = nc.var('v')[:] |
---|
52 |
w = nc.var('w')[:] |
---|
53 |
|
---|
54 |
mu = np.ma.masked_invalid(u) |
---|
55 |
mv = np.ma.masked_invalid(v) |
---|
56 |
mw = np.ma.masked_invalid(w) |
---|
57 |
|
---|
58 |
mumean = mu.mean(axis=0) |
---|
59 |
mvmean = mv.mean(axis=0) |
---|
60 |
mwmean = mw.mean(axis=0) |
---|
61 |
|
---|
62 |
ucount = mu.count(axis=0) |
---|
63 |
vcount = mv.count(axis=0) |
---|
64 |
wcount = mw.count(axis=0) |
---|
65 |
|
---|
66 |
mrho = np.sqrt((mu**2) + (mv**2)) |
---|
67 |
rho = mrho.T.data |
---|
68 |
mtheta = (180 * np.arctan2(mv, mu)/np.pi) |
---|
69 |
theta = np.piecewise(mtheta, |
---|
70 |
(mtheta <= 90, mtheta > 90), |
---|
71 |
(lambda x: 90 - x, lambda x: 450 - x)) |
---|
72 |
theta = theta.T |
---|
73 |
|
---|
74 |
firstElevation = True |
---|
75 |
for x in range(0,rho.shape[0]): |
---|
76 |
if rho[x].any(): |
---|
77 |
fig = plt.figure(figsize=(8, 8), dpi=80, facecolor='w', edgecolor='w') |
---|
78 |
spt = fig.suptitle('Billy Mitchell Sodar :: ' + month + " " + year) |
---|
79 |
rect = [0.16, 0.16, 0.68, 0.68] |
---|
80 |
ax = WindroseAxes(fig, rect, axisbg='w') |
---|
81 |
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) |
---|
82 |
fig.add_axes(ax) |
---|
83 |
ax.bar(theta[x], rho[x], normed=True, opening=0.8, edgecolor='white') |
---|
84 |
l = ax.legend(axespad=-0.20) |
---|
85 |
plt.setp(l.get_texts(), fontsize=8) |
---|
86 |
if not os.path.exists(os.path.join(pngDir,ncFileName)): |
---|
87 |
os.mkdir(os.path.join(pngDir,ncFileName)) |
---|
88 |
outFile = os.path.join(pngDir,ncFileName,('%dm' % z[x]) + os.extsep + pngExt) |
---|
89 |
print 'Saving', outFile |
---|
90 |
fig.savefig(outFile) |
---|
91 |
|
---|