Notebook

Nancy's first tries for looking at nowcasts

In [1]:

%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np
from salishsea_tools import viz_tools, stormtools, tidetools, nc_tools
import datetime
import pandas as pd
import pytz
import requests
import arrow
from StringIO import StringIO
import matplotlib.dates as mdates

In [2]:

#load up data
directory='19nov14'
filepart='20141119_20141119'

fil = '/data/dlatorne/MEOPAR/SalishSea/nowcast/{}/SalishSea_1h_{}_grid_T.nc'.format(directory,filepart)
fil_d ='/data/dlatorne/MEOPAR/SalishSea/nowcast/{}/SalishSea_1d_{}_grid_T.nc'.format(directory,filepart)
filU = '/data/dlatorne/MEOPAR/SalishSea/nowcast/{}/SalishSea_1d_{}_grid_U.nc'.format(directory,filepart)
filV = '/data/dlatorne/MEOPAR/SalishSea/nowcast/{}/SalishSea_1d_{}_grid_V.nc'.format(directory,filepart)

gridT = NC.Dataset(fil)
gridT_d = NC.Dataset(fil_d)
gridU = NC.Dataset(filU)
gridV = NC.Dataset(filV)

bathy = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

Sea Surface Height¶

Tidal Cycle¶

In [3]:

def PA_tidal_predictions(gridT,figsize=(10,5)):
    """ Plots the tidal cycle at Point Atkinson during a 4 week period centred around the dsimulation start date.
    Assumes that a tidal prediction file exists in a specific directory.
    
    :arg gridT: Hourly tracer results dataset from NEMO.
    :type gridT: :class:`netCDF4.Dataset`

    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    :returns: Matplotlib figure object instance
    """
    
    #begining and end time of the simulation file.
    t_orig=(nc_tools.timestamp(gridT,0)).datetime
    t_end =((nc_tools.timestamp(gridT,-1)).datetime)
    
    #set axis limits 2 weeks before and after start date
    ax_start = t_orig - datetime.timedelta(weeks=2)
    ax_end = t_orig + datetime.timedelta(weeks=2)
    ylims=[-3,3]
    
    #load the tidal prediciton file. 
    #Note: How can I do this better? Should I run ttide to make sure this file exists? 
    #Where should this file live?
    path='/data/nsoontie/MEOPAR/analysis/Susan/'
    filename='Point Atkinson_t_tide_compare8_31-Dec-{}_02-Jan-{}.csv'.format(t_orig.year-1,t_orig.year+1)
    tfile = path+filename
    ttide,msl= stormtools.load_tidal_predictions(tfile)
    
    #plotting
    fig,ax=plt.subplots(1,1,figsize=figsize)
    fig.autofmt_xdate()
    ax.plot(ttide.time,ttide.pred_all,'-k')
    #line indicating current date
    ax.plot([t_orig,t_orig],ylims,'-r')
    ax.plot([t_end,t_end],ylims,'-r')
    #axis limits and labels
    ax.set_xlim([ax_start,ax_end])
    ax.set_ylim(ylims)
    ax.set_title('Tidal Predictions at Point Atkinson')
    ax.set_ylabel('SSH [m]')
    
    return fig

In [4]:

fig=PA_tidal_predictions(gridT)

Note: Not sure what to do about where to store the predictions file.

Model Comparisons¶

Canadian Stations¶

In [5]:

def compare_tidal_predictions(name, gridT,gridB,figsize=(10,5) ):
    """ Compares modelled water levels to tidal predictions at a station over one day.
    It is assummed that the tidal predictions were calculated ahead of time and stored in a very specific location.
    Tidal predictions were calculated with the eight consituents used in the model.
    
    :arg name: name of station (e.g Point Atkinson).
    :type name: string
    
    :arg gridT: Hourly tracer results dataset from NEMO.
    :type gridT: :class:`netCDF4.Dataset`
    
    :arg gridB: Model bathymetry file.
    :type gridB: :class:`netCDF4.Dataset`
    
    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    :returns: Matplotlib figure object instance
    """
    
    #Locations of interest.
    lons={'Point Atkinson':-123.26}
    lats={'Point Atkinson': 49.33}
    
    #load bathymetry
    bathy = gridB.variables['Bathymetry'][:, :]
    X = gridB.variables['nav_lon'][:, :]
    Y = gridB.variables['nav_lat'][:, :]
    
    #time stamp of simulation
    t_orig=(nc_tools.timestamp(gridT,0)).datetime
    
    #loading the tidal predictions
    path='/data/nsoontie/MEOPAR/analysis/Susan/'
    filename='_t_tide_compare8_31-Dec-{}_02-Jan-{}.csv'.format(t_orig.year-1,t_orig.year+1)
    tfile = path+name+filename
    ttide,msl= stormtools.load_tidal_predictions(tfile)
    
    #identify grid point of location of interest
    [j,i]=tidetools.find_closest_model_point(lons[name],lats[name],X,Y,bathy,allow_land=True)
    
    #load variables for plotting
    ssh = gridT.variables['sossheig'][:,j,i]
    count=gridT.variables['time_counter'][:]
    t = nc_tools.timestamp(gridT,np.arange(count.shape[0]))
    #convert times to datetimes because that is what the plot wants
    for i in range(len(t)):
        t[i]=t[i].datetime
    
    #plotting
    fig,ax =plt.subplots(1,1,figsize=figsize) 
    ax.plot(t,ssh,label='model')
    ax.plot(ttide.time,ttide.pred_8,label='tidal predictions')
    ax.set_xlim(t_orig,t_orig + datetime.timedelta(days=1))
    ax.set_ylim([-3,3])
    ax.legend(loc=0)
    ax.set_title(name)
    ax.set_ylabel('Water levels wrt MSL (m)')
    ax.set_xlabel('time [UTC]')

    
    return fig

Note: I'm hoping this function will change when we have a systematic way of collecting recent water level observations.

In [6]:

fig =compare_tidal_predictions('Point Atkinson', gridT,bathy)

Note: There is no correction for 8 constituents in the model output. But, the tidal predictions were calculated with only 8 constituents

American Stations¶

In [7]:

def get_NOAA_wlevels(station_no, start_date,end_date):
    """ Retrieves recent, 6 minute interval, NOAA water levels relative to mean sea level
    from a station in a given date range.
    
    :arg station_no: the NOAA station number.
    :type station_no: integer
    
    :arg start_date: the start of the date range eg. 01-Jan-2014.
    :type start_date: string
    
    :arg end_date: the end of the date range ef. 02-Jan-2014.
    :type end_date: string
    
    :returns: DataFrame object with time and wlev columns, among others that are irrelevant.
    """
    
    outfile = 'wlev_'+str(station_no)+'_'+str(start_date)+'_'+str(end_date)+'.csv'
    
    st_ar=arrow.Arrow.strptime(start_date, '%d-%b-%Y')
    end_ar=arrow.Arrow.strptime(end_date, '%d-%b-%Y')
    
    base_url = 'http://tidesandcurrents.noaa.gov/api/datagetter?product=water_level&application=NOS.COOPS.TAC.WL'
    params = {
    'begin_date': st_ar.format('YYYYMMDD'),
    'end_date': end_ar.format('YYYYMMDD'),
    'datum':  'MSL',
    'station': str(station_no),
    'time_zone': 'GMT',
    'units':'metric',
    'format': 'csv',
}
    response = requests.get(base_url, params=params)

    fakefile = StringIO(response.content)
    
    obs = pd.read_csv(fakefile)
    obs=obs.rename(columns={'Date Time': 'time', ' Water Level': 'wlev'})
    
    return obs

In [8]:

def compare_water_levels(name,gridT,gridB,figsize=(10,5) ):
    """ Compares modelled water levels to observed water levels at a NOAA station over one day.
    
    :arg name: name of the NOAA station (e.g NeahBay, CherryPoint, FridayHarbor).
    :type name: string
    
    :arg gridT: Hourly tracer results dataset from NEMO.
    :type gridT: :class:`netCDF4.Dataset`
    
    :arg gridB: Model bathymetry file.
    :type gridB: :class:`netCDF4.Dataset`
    
    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    :returns: Matplotlib figure object instance
    """
    
    stations = {'CherryPoint': 9449424,'NeahBay':9443090, 'FridayHarbor': 9449880 }
    lats={'CherryPoint': 48.866667,'NeahBay': 48.4, 'FridayHarbor': 48.55}
    lons={'CherryPoint': -122.766667, 'NeahBay':-124.6, 'FridayHarbor': -123.016667}
    
    bathy = gridB.variables['Bathymetry'][:, :]
    X = gridB.variables['nav_lon'][:, :]
    Y = gridB.variables['nav_lat'][:, :]
    
    t_orig=(nc_tools.timestamp(gridT,0)).datetime
    start_date = t_orig.strftime('%d-%b-%Y')
    end_date = (t_orig + datetime.timedelta(days=1)).strftime('%d-%b-%Y')
    
    obs=get_NOAA_wlevels(stations[name],start_date,end_date)
    
    [j,i]=tidetools.find_closest_model_point(lons[name],lats[name],X,Y,bathy,allow_land=False)
    
    ssh = gridT.variables['sossheig'][:,j,i]
    count=gridT.variables['time_counter'][:]
    t = nc_tools.timestamp(gridT,np.arange(count.shape[0]))
    for i in range(len(t)):
        t[i]=t[i].datetime
    
    fig,ax =plt.subplots(1,1,figsize=figsize) 
    ax.plot(t,ssh,label='model')
    ax.plot(obs.time,obs.wlev,label='observed water levels')
    ax.set_xlim(t_orig,t_orig + datetime.timedelta(days=1))
    ax.set_ylim([-1.5,1.5])
    ax.legend(loc=0)
    ax.set_title(name)
    ax.set_ylabel('Water levels wrt MSL (m)')

    
    return fig

In [9]:

fig=compare_water_levels('NeahBay',gridT,bathy)

In [10]:

fig=compare_water_levels('CherryPoint',gridT,bathy)

In [11]:

fig=compare_water_levels('FridayHarbor',gridT,bathy)

Note: there is no correction for 8 constituents.

Maximum Water Levels¶

In [12]:

def PA_max_ssh(gridT, gridB,figsize=(15,10)):
    """Function that plots the water level at every hour throughout the day and identifies the maximum.
    It also plots the sea surface height throughout the region for the time when the sea surface height
    was at its maximum at Point Atkinson.
    
    :arg gridT: Hourly tracer results dataset from NEMO.
    :type gridT: :class:`netCDF4.Dataset`
    
    :arg gridB: Model bathymetry file.
    :type gridB: :class:`netCDF4.Dataset`
    
    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    :returns: Matplotlib figure object instance
    """

    #defining Point Atkinson
    lon_PA=-123.26
    lat_PA= 49.33
    bathy = gridB.variables['Bathymetry'][:, :]
    X = gridB.variables['nav_lon'][:, :]
    Y = gridB.variables['nav_lat'][:, :]
    [j,i]=tidetools.find_closest_model_point(lon_PA,lat_PA,X,Y,bathy,allow_land=True)
    
    #loading sea surface height
    ssh = gridT.variables['sossheig']
    #loading sea surface height at Point Atkinson
    ssh_loc = ssh[:,j,i]
    #index when sea surface height is at its maximum at Point Atkinson
    m = max(ssh_loc)
    index = [k for k, l in enumerate(ssh_loc) if l == m]
    index = index[0]
    #sea surface height when there is a maximum at Point Atkinson
    ssh_max = np.ma.masked_values(ssh[index], 0)
    
    #time for curve
    count=gridT.variables['time_counter'][:]
    t = nc_tools.timestamp(gridT,np.arange(count.shape[0]))
    for ind in range(len(t)):
        t[ind]=t[ind].datetime
    t_orig=(nc_tools.timestamp(gridT,0)).datetime
    start_date = t_orig.strftime('%d-%b-%Y')
    end_date = (t_orig + datetime.timedelta(days=1)).strftime('%d-%b-%Y')
    
    #timestamp
    timestamp = nc_tools.timestamp(gridT,index)
    
    #figure
    fig,(ax1,ax2) =plt.subplots(1,2,figsize=figsize) 
    
    #curve plot
    ax1.plot(t,ssh_loc,label='Model')
    ax1.plot(t[index],ssh_loc[index],color='Indigo',marker='D',markersize=12,label='Maximum SSH')
    ax1.set_xlim(t_orig,t_orig + datetime.timedelta(days=1))
    ax1.set_ylim([-3,3])
    ax1.set_title('Point Atkinson')
    ax1.set_xlabel(timestamp.strftime('%d-%b-%Y'))
    ax1.set_ylabel('Water levels wrt MSL (m)')
    ax1.legend(loc = 0, numpoints = 1)
    ax1.set_position((0, 0.3, 0.55, 0.4))
    
    #ssh profile
    viz_tools.set_aspect(ax2)
    land_colour = 'burlywood'
    ax2.set_axis_bgcolor(land_colour)
    cmap = plt.get_cmap('RdYlGn_r')
    cs = [-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.4,2.8,3]
    mesh=ax2.contourf(ssh_max,cs,cmap=cmap,extend='both')
    cbar = fig.colorbar(mesh,ax=ax2)
    cbar.set_ticks(cs)
    cbar.set_label('[m]')
    ax2.grid()
    ax2.set_xlabel('x Index')
    ax2.set_ylabel('y Index')
    viz_tools.plot_coastline(ax2,gridB)
    ax2.set_title('Sea Suface Height: ' + timestamp.strftime('%d-%b-%Y, %H:%M'))
    ax2.plot(i,j,marker='D',color='Indigo',ms=8)

    return fig

Above: Smaller chunks?

In [13]:

fig =PA_max_ssh(gridT,bathy)

Winds¶

In [14]:

def plot_Sandheads_winds(gridT,figsize=(10,7)):
    """ Plot the observed winds at Sandheads during the simulation.
    
    :arg gridT: Hourly tracer results dataset from NEMO.
    :type gridT: :class:`netCDF4.Dataset`

    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    :returns: Matplotlib figure object instance
    """
    
    #simulation date range.
    t_orig=(nc_tools.timestamp(gridT,0)).datetime
    t_end=(nc_tools.timestamp(gridT,-1)).datetime
    #strings for timetamps of EC data
    start=t_orig.strftime('%d-%b-%Y')
    end=t_end.strftime('%d-%b-%Y')

    [winds,dirs,temps,time, lat,lon] = stormtools.get_EC_observations('Sandheads',start,end)

    fig,axs=plt.subplots(2,1,figsize=figsize)
    #plotting wind speed
    ax=axs[0]
    ax.set_title('Winds at Sandheads ' + start )
    ax.plot(time,winds)
    ax.set_xlim([time[0],time[-1]])
    ax.set_ylim([0,20])
    ax.set_ylabel('Wind speed (m/s)')
    #plotting wind direction
    ax=axs[1]
    ax.plot(time,dirs)
    ax.set_ylabel('Wind direction \n (degress CCW of East)')
    ax.set_ylim([0,360])
    ax.set_xlim([time[0],time[-1]])
    ax.set_xlabel('Time [UTC]')
    
    #Eventually add comparison with modelled winds?
    
    return fig

In [15]:

fig=plot_Sandheads_winds(gridT)

We'll have to think about a way to load the modelled weather too.

Would quivers be better? Map indicating where Sandheads is?

Daily Averaged Fields¶

Thalweg¶

In [16]:

def thalweg_salinity(gridT_d,figsize=(20,8),cs = [26,27,28,29,30,30.2,30.4,30.6,30.8,31,32,33,34]):
    """ Plot the daily averaged salinity field along the thalweg.
    
    :arg gridT_d: Daily tracer results dataset from NEMO.
    :type gridT_d: :class:`netCDF4.Dataset`

    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    :arg cs: list of salinity contour levels for shading.
    :type cs: list
    
    :returns: Matplotlib figure object instance
    """
    
    lon_d = gridT_d.variables['nav_lon']
    lat_d = gridT_d.variables['nav_lat']
    dep_d = gridT_d.variables['deptht']
    sal_d = gridT_d.variables['vosaline']
    
    lines = np.loadtxt('/data/nsoontie/MEOPAR/tools/analysis_tools/thalweg.txt', delimiter=" ", unpack=False)
    lines = lines.astype(int)

    thalweg_lon = lon_d[lines[:,0],lines[:,1]]
    thalweg_lat = lat_d[lines[:,0],lines[:,1]]

    ds=np.arange(0,lines.shape[0],1);
    XX,ZZ = np.meshgrid(ds,-dep_d[:])
    
    salP=sal_d[0,:,lines[:,0],lines[:,1]]
    salP= np.ma.masked_values(salP,0)
    
    fig,ax=plt.subplots(1,1,figsize=figsize)
    land_colour = 'burlywood'
    ax.set_axis_bgcolor(land_colour)
    mesh=ax.contourf(XX,ZZ,salP,cs,cmap='hsv',extend='both')
    
    cbar=fig.colorbar(mesh,ax=ax)
    cbar.set_ticks(cs)
    cbar.set_label('Practical Salinity [psu]')
    
    timestamp = nc_tools.timestamp(gridT_d,0)
    ax.set_title('Salinity field along thalweg: ' + timestamp.format('DD-MMM-YYYY'))
    ax.set_ylabel('Depth [m]')
    ax.set_xlabel('position along thalweg')
    
    return fig

In [17]:

fig = thalweg_salinity(gridT_d)

Surface Fields¶

In [18]:

def plot_surface(gridT_d, gridU, gridV, gridB, figsize=(18,10)):
    """Function that plots the daily average surface salinity, temperature and currents.
    
    :arg gridT_d: Daily tracer results dataset from NEMO.
    :type gridT_d: :class:`netCDF4.Dataset`
    
    :arg gridU: Daily zonal velocity results dataset from NEMO.
    :type gridU: :class:`netCDF4.Dataset`
    
    :arg gridV: Daily meridional velocity results dataset from NEMO.
    :type gridV: :class:`netCDF4.Dataset`
    
    :arg gridB: Model bathymetry file.
    :type gridB: :class:`netCDF4.Dataset`
    
    :arg figsize: figure size (width, height) in inches.
    :type figsize: 2-tuple
    
    """
    
    #loading lon, lat, depth, salinity, temperature
    lon_d = gridT_d.variables['nav_lon']
    lat_d = gridT_d.variables['nav_lat']
    dep_d = gridT_d.variables['deptht']
    sal_d = gridT_d.variables['vosaline']
    tem_d = gridT_d.variables['votemper']
    
    #tracers ready to plot
    t, z = 0, 0
    sal_d = np.ma.masked_values(sal_d[t, z], 0)
    tem_d = np.ma.masked_values(tem_d[t, z], 0)
    
    
    #for loop
    tracers = [sal_d, tem_d]
    titles = ['Average Salinity: ','Average Temperature: ']
    cmaps = ['ocean','jet']
    units = ['[psu]','[degC]']
    
    fig, (ax1,ax2,ax3) = plt.subplots(1, 3, figsize=figsize)
    axs = [ax1, ax2]
    plots = np.arange(1,3,1)
    
    #***temperature and salinity plots
    
    for ax,tracer,title,cmap,unit,plot in zip(axs,tracers,titles,cmaps,units,plots):
        #general set up
        viz_tools.set_aspect(ax)
        land_colour = 'burlywood'
        ax.set_axis_bgcolor(land_colour)
        cmap = plt.get_cmap(cmap)
        
        #colourmap breakdown
        if plot == 1:
            cs = [0,4,8,12,16,20,24,28,32,36]
        if plot == 2:
            cs = [0,2,4,6,8,10,12,14,16,18,20]

        #plotting tracers
        mesh=ax.contourf(tracer,cs,cmap=cmap,extend='both')
        
        #colour bars
        cbar = fig.colorbar(mesh,ax=ax)
        cbar.set_ticks(cs)
        
        #labels
        ax.grid()
        ax.set_xlabel('x Index')
        ax.set_ylabel('y Index')
        timestamp = nc_tools.timestamp(gridT_d,0)
        ax.set_title(title + timestamp.format('DD-MMM-YYYY'))
        cbar.set_label(unit)
    
    #loading velocity components
    ugrid = gridU.variables['vozocrtx']
    vgrid = gridV.variables['vomecrty']
    zlevels = gridU.variables['depthu']
    timesteps = gridU.variables['time_counter']

    #day's average, surface velocity field
    t, zlevel = 0, 0

    #region
    y_slice = np.arange(0, ugrid.shape[2])
    x_slice = np.arange(0, ugrid.shape[3])

    arrow_step = 25
    y_slice_a = y_slice[::arrow_step]
    x_slice_a = x_slice[::arrow_step]

    #masking arrays
    ugrid_tzyx = np.ma.masked_values(ugrid[t, zlevel, y_slice_a, x_slice_a], 0)
    vgrid_tzyx = np.ma.masked_values(vgrid[t, zlevel, y_slice_a, x_slice_a], 0)
    #unstagger velocity values
    u_tzyx, v_tzyx = viz_tools.unstagger(ugrid_tzyx, vgrid_tzyx)
    #velocity magnitudes
    speeds = np.sqrt(np.square(u_tzyx) + np.square(v_tzyx))

    #***velocity plot
    viz_tools.set_aspect(ax3)
    cs = [0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6]

    quiver = ax3.quiver(x_slice_a[1:], y_slice_a[1:], u_tzyx, v_tzyx, speeds, pivot='mid', cmap='gnuplot_r', width=0.015)
    viz_tools.plot_land_mask(ax3, bathy, xslice=x_slice, yslice=y_slice, color='burlywood')
    
    cbar = fig.colorbar(quiver,ax=ax3)
    cbar.set_ticks(cs)
    cbar.set_label('[m / s]')

    ax3.set_xlim(x_slice[0], x_slice[-1])
    ax3.set_ylim(y_slice[0], y_slice[-1])
    ax3.grid()

    ax3.set_xlabel('x Index')
    ax3.set_ylabel('y Index')
    ax3.set_title('Average Velocity Field: ' + timestamp.format('DD-MMM-YYYY') + 
                  u', depth\u2248{d:.2f} {z.units}'.format(d=zlevels[zlevel], z=zlevels))
    ax3.quiverkey(quiver, 355, 850, 1, '1 m/s', coordinates='data', color='Indigo', labelcolor='black')
    
    return fig

In [19]:

fig = plot_surface(gridT_d, gridU, gridV, bathy)

Above: break up into smaller chunks?

In [19]: