Quick peak at new Neah Bay Sea surface height forcing files.

In [1]:

import netCDF4 as NC
import numpy as np

import matplotlib.pyplot as plt
import glob
import os
import datetime

from salishsea_tools import nc_tools
from StringIO import StringIO
import urllib2
import csv

%matplotlib inline

In [2]:

def plot_winds(sax,eax):

    response = urllib2.urlopen('http://www.ndbc.noaa.gov/data/realtime2/46087.cwind')
    html = response.read()
    text_file = open("Output.txt", "w")
    text_file.write(html)
    text_file.close()

    fil= csv.reader(open('Output.txt','rb'),delimiter=' ')
    fil.next()
    fil.next()
    timeair=[]; wind_speed=[]; wind_dir=[]
    for row in fil:
        aList=[]
        for element in row:
            if element != '':
                aList.append(element)
        timeair.append(datetime.datetime(int(aList[0]),int(aList[1]),int(aList[2]),int(aList[3]),int(aList[4])))
        d = -float(aList[5])+270
        d=d+ 360 * (d<0)
        wind_dir.append(d)
        wind_speed.append(float(aList[6]))
    timeair=np.array(timeair)
    wind_dir=np.array(wind_dir)
    wind_speed=np.array(wind_speed)

    fig,axs=plt.subplots(2,1,figsize=(10,10))
    ax=axs[0]
    ax.plot(timeair,wind_speed);
    ax.set_xlim(sax,eax); ax.set_ylim([0,20]); ax.set_ylabel('wind speed (m/s)')
    ax.set_title('Winds at Neah Bay Buoy - 46087')
    ax=axs[1]
    ax.plot(timeair,wind_dir); ax.set_ylim([0,360])
    ax.set_xlim(sax,eax); ax.set_ylabel('wind direction (degrees CCW of East)')

In [3]:

#forecasts
files = glob.glob('/ocean/nsoontie/MEOPAR/sshNeahBay/fcst/ssh*')
files.sort(key=os.path.basename)

#observations
fileso =glob.glob('/ocean/nsoontie/MEOPAR/sshNeahBay/obs/ssh*')
fileso.sort(key=os.path.basename)

Compare Neah Bay Forecasts to Observations¶

In [4]:

#forecasts
sshs=[]; times=[]
start = datetime.datetime(2014,10,29)
day=0
for f in files:
    data=NC.Dataset(f)
    ssh=data.variables['sossheig']
    sshs.append(ssh[:,0,0])
    for h in range(ssh.shape[0]):
        t=start+datetime.timedelta(days=day,hours=h)
        times.append(t);
    day=day+1;
sshs=np.array(sshs).reshape(len(files)*24)
    
#obs
sshso=[]; timeso=[]
start = datetime.datetime(2014,10,28)
day=0
for f in fileso:
    data=NC.Dataset(f)
    ssh=data.variables['sossheig']
    sshso.append(ssh[:,0,0])
    for h in range(ssh.shape[0]):
        t=start+datetime.timedelta(days=day,hours=h)
        timeso.append(t);
    day=day+1;
sshso=np.array(sshso).reshape(len(fileso)*24)
    

fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(times,sshs,label='forecast')
ax.plot(timeso,sshso,label='observation')
ax.legend(loc=0)
ax.set_ylabel('SSH anomaly (m)')
ax.set_title('Neah Bay storm surge')
ax.set_xlim([datetime.datetime(2014,11,26),datetime.datetime(2014,12,8)])
ax.grid()

Out[4]:

<matplotlib.text.Text at 0x7f31d418e190>

Neah Bay Winds¶

In [5]:

eax =times[-1]; sax = timeso[0];
plot_winds(sax,eax)

In [5]: