| 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 |
html1 = """<html> |
|---|
| 27 |
<head> |
|---|
| 28 |
<title> |
|---|
| 29 |
Billy Mitchell Sodar :: """ |
|---|
| 30 |
|
|---|
| 31 |
|
|---|
| 32 |
html2 = """ :: Wind Distribution at """ |
|---|
| 33 |
|
|---|
| 34 |
html3 = """ Meters Above Sea Level |
|---|
| 35 |
</title> |
|---|
| 36 |
</head> |
|---|
| 37 |
<body> |
|---|
| 38 |
<form> |
|---|
| 39 |
<a href=\"""" |
|---|
| 40 |
|
|---|
| 41 |
html4 = """\"><Previous Month</a> |
|---|
| 42 |
<a href=\"""" |
|---|
| 43 |
|
|---|
| 44 |
html5 = """\">Next Month></a> |
|---|
| 45 |
<a href=\"""" |
|---|
| 46 |
|
|---|
| 47 |
html6 = """\"><Previous Elevation</a> |
|---|
| 48 |
<a href=\"""" |
|---|
| 49 |
|
|---|
| 50 |
html7 = """\">Next Elevation></a> |
|---|
| 51 |
<a href=\"all.html\"><All At Once></a> |
|---|
| 52 |
</form> |
|---|
| 53 |
<br/> |
|---|
| 54 |
<img src=\"""" |
|---|
| 55 |
|
|---|
| 56 |
html8 = """\"/> |
|---|
| 57 |
</body> |
|---|
| 58 |
</html>""" |
|---|
| 59 |
|
|---|
| 60 |
firstMonth = True |
|---|
| 61 |
for previous,ncFile,next in files: |
|---|
| 62 |
print 'Processing',ncFile |
|---|
| 63 |
ncFileName = os.path.splitext(os.path.basename(ncFile))[0] |
|---|
| 64 |
year = ncFileName[19:23] |
|---|
| 65 |
month = ncFileName[24:26] |
|---|
| 66 |
month = datetime.datetime(int(year),int(month),1).strftime('%B') |
|---|
| 67 |
if previous: |
|---|
| 68 |
previous = os.path.splitext(os.path.basename(previous))[0] |
|---|
| 69 |
previousYear = previous[19:23] |
|---|
| 70 |
previousMonth = previous[24:26] |
|---|
| 71 |
previousMonth = datetime.datetime(int(previousYear),int(previousMonth),1).strftime('%B') |
|---|
| 72 |
if next: |
|---|
| 73 |
next = os.path.splitext(os.path.basename(next))[0] |
|---|
| 74 |
nextYear = next[19:23] |
|---|
| 75 |
nextMonth = next[24:26] |
|---|
| 76 |
nextMonth = datetime.datetime(int(nextYear),int(nextMonth),1).strftime('%B') |
|---|
| 77 |
|
|---|
| 78 |
nc = pycdf.CDF(ncFile) |
|---|
| 79 |
|
|---|
| 80 |
t = nc.var('time')[:] |
|---|
| 81 |
z = nc.var('z')[:] |
|---|
| 82 |
u = nc.var('u')[:] |
|---|
| 83 |
v = nc.var('v')[:] |
|---|
| 84 |
w = nc.var('w')[:] |
|---|
| 85 |
|
|---|
| 86 |
mu = np.ma.masked_invalid(u) |
|---|
| 87 |
mv = np.ma.masked_invalid(v) |
|---|
| 88 |
mw = np.ma.masked_invalid(w) |
|---|
| 89 |
|
|---|
| 90 |
mumean = mu.mean(axis=0) |
|---|
| 91 |
mvmean = mv.mean(axis=0) |
|---|
| 92 |
mwmean = mw.mean(axis=0) |
|---|
| 93 |
|
|---|
| 94 |
ucount = mu.count(axis=0) |
|---|
| 95 |
vcount = mv.count(axis=0) |
|---|
| 96 |
wcount = mw.count(axis=0) |
|---|
| 97 |
|
|---|
| 98 |
mrho = np.sqrt((mu**2) + (mv**2)) |
|---|
| 99 |
mrho = mrho.T |
|---|
| 100 |
mtheta = (180 * np.arctan2(mv, mu)/np.pi) |
|---|
| 101 |
theta = np.piecewise(mtheta, |
|---|
| 102 |
(mtheta <= 90, mtheta > 90), |
|---|
| 103 |
(lambda x: 90 - x, lambda x: 450 - x)) |
|---|
| 104 |
mtheta = np.ma.array(theta,mask=mtheta.mask) |
|---|
| 105 |
mtheta = mtheta.T |
|---|
| 106 |
|
|---|
| 107 |
rho = [np.array([elem |
|---|
| 108 |
for elem,mask |
|---|
| 109 |
in zip(row.data,row.mask) |
|---|
| 110 |
if not mask]) |
|---|
| 111 |
for row |
|---|
| 112 |
in mrho] |
|---|
| 113 |
theta = [np.array([elem |
|---|
| 114 |
for elem,mask |
|---|
| 115 |
in zip(row.data,row.mask) |
|---|
| 116 |
if not mask]) |
|---|
| 117 |
for row |
|---|
| 118 |
in mtheta] |
|---|
| 119 |
|
|---|
| 120 |
firstElevation = True |
|---|
| 121 |
indices = range(len(rho)) |
|---|
| 122 |
previousX = [indices[-1]] + indices[:-1] |
|---|
| 123 |
nextX = indices[1:] + [indices[0]] |
|---|
| 124 |
indices = zip(previousX,indices,nextX) |
|---|
| 125 |
for previousX,x,nextX in indices: |
|---|
| 126 |
if rho[x].any(): |
|---|
| 127 |
fig = plt.figure(figsize=(8, 8), dpi=80, facecolor='w', edgecolor='w') |
|---|
| 128 |
spt = fig.suptitle('Billy Mitchell Sodar :: ' + month + " " + year) |
|---|
| 129 |
rect = [0.16, 0.16, 0.68, 0.68] |
|---|
| 130 |
ax = WindroseAxes(fig, rect, axisbg='w') |
|---|
| 131 |
ax.set_title('Wind Distribution at %d Meters Above Sea Level\nPercentage of Wind Magnitude From Direction\n\n' % z[x],fontsize=spt.get_fontsize()*0.8) |
|---|
| 132 |
fig.add_axes(ax) |
|---|
| 133 |
ax.bar(theta[x], rho[x], bins = np.arange(0.0,20.0,2.5), normed=True, opening=0.8, edgecolor='white') |
|---|
| 134 |
yticklabels = ['',] + \ |
|---|
| 135 |
[" %s%%" % plt.getp(textobj,'text') |
|---|
| 136 |
for textobj in plt.getp(ax,'yticklabels')[1:]] |
|---|
| 137 |
plt.setp(ax,yticklabels=yticklabels) |
|---|
| 138 |
l = ax.legend(axespad=-0.20,title="Magnitude (m/s)") |
|---|
| 139 |
plt.setp(l.get_texts(), fontsize=8) |
|---|
| 140 |
if not os.path.exists(os.path.join(pngDir,ncFileName)): |
|---|
| 141 |
os.mkdir(os.path.join(pngDir,ncFileName)) |
|---|
| 142 |
outFile = os.path.join(pngDir,ncFileName,('%dm' % z[x]) + os.extsep + pngExt) |
|---|
| 143 |
print 'Saving', outFile |
|---|
| 144 |
fig.savefig(outFile) |
|---|
| 145 |
|
|---|
| 146 |
# htmlFile = os.path.join(pngDir,ncFileName + "_" + ('%dm' % z[x]) + os.extsep + htmlExt) |
|---|
| 147 |
# html = html1 + month + " " + year |
|---|
| 148 |
# html = html + html2 + ('%dm' % z[x]) |
|---|
| 149 |
# html = html + html3 + previous + "_" + ('%dm' % z[x]) + os.extsep + htmlExt |
|---|
| 150 |
# html = html + html4 + next + "_" + ('%dm' % z[x]) + os.extsep + htmlExt |
|---|
| 151 |
# html = html + html5 + ncFileName + "_" + ('%dm' % z[previousX]) + os.extsep + htmlExt |
|---|
| 152 |
# html = html + html6 + ncFileName + "_" + ('%dm' % z[nextX]) + os.extsep + htmlExt |
|---|
| 153 |
# html = html + html7 + os.path.join(ncFileName,os.path.basename(outFile)) + html8 |
|---|
| 154 |
# handle = open(htmlFile,'w') |
|---|
| 155 |
# handle.write(html) |
|---|
| 156 |
# handle.close() |
|---|
| 157 |
|
|---|
| 158 |
# if firstMonth and firstElevation: |
|---|
| 159 |
# firstMonth = False |
|---|
| 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 |
|
|---|
