Notebook

Figures Module: Imports¶

In [1]:

from __future__ import division

from cStringIO import StringIO
import datetime
import glob
import os

import arrow
from dateutil import tz
import matplotlib
from matplotlib.backends import backend_agg as backend
import matplotlib.cm as cm
import matplotlib.dates as mdates
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
import pandas as pd
import requests
from scipy import interpolate as interp

from salishsea_tools import (
    nc_tools,
    viz_tools,
    stormtools,
    tidetools,
)

Figures Notebook: Imports¶

In [2]:

from __future__ import division

import datetime
from glob import glob
import os

from IPython.core.display import HTML
import netCDF4 as nc

from salishsea_tools.nowcast import figures

%matplotlib inline

Figures Notebook: Results¶

In [3]:

def results_dataset(period, grid, results_dir):
    """Return the results dataset for period (e.g. 1h or 1d)
    and grid (e.g. grid_T, grid_U) from results_dir.
    """
    filename_pattern = 'SalishSea_{period}_*_{grid}.nc'
    filepaths = glob(os.path.join(results_dir, filename_pattern.format(period=period, grid=grid)))
    return nc.Dataset(filepaths[0])

In [4]:

def winds_dataset(run_date, model_path):
    """Return the winds dataset for the current date from model_path."""
    year, month, day = run_date.strftime('%G'), run_date.strftime('%m'), run_date.strftime('%d')
    filename_pattern = 'ops_y{yyyy}m{mm}d{dd}.nc'
    filepaths = glob(os.path.join(model_path, filename_pattern.format(yyyy=year,mm=month,dd=day)))
    return nc.Dataset(filepaths[0])

In [5]:

#run_date = datetime.datetime(2015,1,8)
run_date = datetime.date.today()

# Results dataset location
results_home = '/data/dlatorne/MEOPAR/SalishSea/nowcast/'
results_dir = os.path.join(results_home, run_date.strftime('%d%b%y').lower())
# model winds
model_path = '/ocean/sallen/allen/research/MEOPAR/Operational/'

In [6]:

grid_T_hr = results_dataset('1h', 'grid_T', results_dir)

In [7]:

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

Figures Module: Other stuff¶

In [8]:

model_c = 'MediumBlue'
observations_c = 'DarkGreen'
predictions_c = 'MediumVioletRed'
stations_c = cm.rainbow(np.linspace(0, 1, 7))

In [9]:

time_shift = datetime.timedelta(hours=-8)
hfmt = mdates.DateFormatter('%m/%d %H:%M')

In [10]:

MSL_DATUMS = {
    'Point Atkinson': 3.10, 'Victoria': 1.90,
    'Campbell River': 2.89, 'Patricia Bay': 2.30}

Threshold image¶

In [11]:

reload(figures)

Out[11]:

<module 'salishsea_tools.nowcast.figures' from '/ocean/imachuca/MEOPAR/tools/SalishSeaTools/salishsea_tools/nowcast/figures.pyc'>

In [12]:

fig = figures.plot_threshold_website(bathy, grid_T_hr, model_path, 0.1)

Website thumbnail¶

In [13]:

import scipy.io as sio

In [14]:

coastline = sio.loadmat('/ocean/rich/more/mmapbase/bcgeo/PNW.mat')

def draw_coast(ax, PNW_coastline): coast={} coast['lat'] = PNW_coastline['ncst'][:,1] coast['lon'] = PNW_coastline['ncst'][:,0] ax.plot(coast['lon'],coast['lat'],'-k',rasterized=True,markersize=1)

12 instances of functions which needed "figures" to be referenced

def website_thumbnail(grid_B, grid_T, model_path, PNW_coastline, scale=0.1, PST=1, figsize = (18, 20)): """ Overview image for Salish Sea website. Plots a map of the Salish Sea with markers indicating extreme water at Point Atkinson, Victoria nd Campbell River. Also plots wind vectors averaged over 4 ours before the max ssh at Point Atkinson. Includes text boxes with max water level and timing.

:arg grid_B: Bathymetry dataset for the Salish Sea NEMO model.
:type grid_B: :class:`netCDF4.Dataset`

:arg grid_T: Hourly tracer results dataset from NEMO.
:type grid_T: :class:`netCDF4.Dataset`

:arg model_path: The directory where the model wind files are stored.
:type model_path: string

:arg scale: scale factor or wind arrows
:type scale: float

:arg PST: Specifies if plot should be presented in PST.
1 = plot in PST, 0 = plot in UTC.
:type PST: 0 or 1

:arg figsize: Figure size (width, height) in inches.
:type figsize: 2-tuple

:returns: matplotlib figure object instance (fig).
"""

title_font_thumb = {
'fontname': 'Bitstream Vera Sans', 'size': '40', 'color': 'black',
'weight': 'medium'
}
axis_font_thumb = {'fontname': 'Bitstream Vera Sans', 'size': '30'}

# Stations information
[lats, lons] = station_coords()

# Bathymetry
bathy, X, Y = tidetools.get_bathy_data(grid_B)

# Time range
t_orig, t_final, t = get_model_time_variables(grid_T)

# Wind time
inds = isolate_wind_timing('Point Atkinson', grid_T, grid_B,
                           model_path, t, 4, average=True)

# Set up loop
names = ['Point Atkinson', 'Campbell River', 'Victoria']

# Set up Information
max_sshs = {}; max_times = {}; max_winds = {}

# Figure
fig = plt.figure(figsize=figsize)
gs = gridspec.GridSpec(2, 3, width_ratios=[1, 1, 1], height_ratios=[6, 1])
gs.update(hspace=0.15, wspace=0.05)
ax = fig.add_subplot(gs[0, :])
ax1 = fig.add_subplot(gs[1, 0])
ax2 = fig.add_subplot(gs[1, 1])
ax3 = fig.add_subplot(gs[1, 2])

# Map
viz_tools.set_aspect(ax)
viz_tools.plot_land_mask(ax, grid_B,color='burlywood',coords='map')
viz_tools.plot_coastline(ax,grid_B,coords='map')
draw_coast(ax, PNW_coastline)
ax.set_xlabel('Longitude',**axis_font_thumb)
ax.set_ylabel('Latitude',**axis_font_thumb)
ax.grid()
ax.set_xlim([-126.4,-121.3])
ax.set_ylim([46.8,51.1])

for name in names:
    # Get sea surface height
    [j, i] = tidetools.find_closest_model_point(lons[name], lats[name],
                                        X, Y, bathy, allow_land=False)
    ssh = grid_T.variables['sossheig']
    ssh_loc = ssh[:, j, i]

    # Get tides and ssh
    ttide = get_tides(name)
    sdt = t_orig.replace(minute=0)
    edt = t_final + datetime.timedelta(minutes=30)
    ssh_corr = stormtools.correct_model(ssh_loc, ttide, sdt, edt)

    # Plot thresholds
    plot_threshold_map(ax, ttide, ssh_corr, 'o', 70, 0.3, name)
    # Plot winds
    twind = plot_wind_vector(ax, name, t_orig, t_final, model_path, inds, scale)

    # Information
    res = compute_residual(ssh_loc, ttide, t_orig, t_final)
    [max_ssh, index_ssh, tmax, max_res, max_wind, ind_w] = get_maxes(ssh_corr, t, res, lons[name], lats[name], model_path)
    max_sshs[name] = max_ssh
    max_times[name] = tmax
    max_winds[name] = max_wind

# Add winds for other stations
for name in ['Neah Bay', 'Cherry Point', 'Sandheads', 'Friday Harbor']:
    plot_wind_vector(ax, name, t_orig, t_final, model_path, inds, scale)

# Reference arrow
ax.arrow(-121.5, 50.9, 0.*scale, -5.*scale,
          head_width=0.05, head_length=0.1, width=0.02,
          color='white',fc='DarkMagenta', ec='black')
ax.text(-121.6, 50.95, "Reference: 5 m/s", rotation=90, fontsize=20)

# Location labels
ax.text(-125.7, 47.8, 'Pacific\nOcean', fontsize=30, color='DimGray')
ax.text(-122.8, 50.1, '  British\nColumbia', fontsize=30, color='DimGray')
ax.text(-124.2, 47.3, 'Washington\n    State', fontsize=30, color='DimGray')
ax.text(-122.2, 47.6, ' Puget\nSound', fontsize=20, color='DimGray')
ax.text(-124.6, 48.1, 'Strait of\nJuan de Fuca', fontsize=20, color='DimGray', rotation=-18)
ax.text(-124.5, 49.1, 'Strait of \n Georgia', fontsize=20, color='DimGray', rotation=-18)

# Figure format
t = (twind[0]+PST*time_shift).strftime('%A, %B %d, %Y')
ax.set_title('Marine and Atmospheric Conditions\n {time}'.format(time=t), **title_font_thumb)
fig.patch.set_facecolor('#2B3E50')
axis_colors(ax, 'gray')

# Legend
axs = [ax1, ax2, ax3]
cs=['green','Gold','red']
for ax, name,thresh_c in zip (axs, names,cs):
    plt.setp(ax.spines.values(), visible=False)
    ax.xaxis.set_visible(False); ax.yaxis.set_visible(False)
    axis_colors(ax, 'blue')
    display_time=(max_times[name]+PST*time_shift).strftime('%H:%M')
    ax.set_xlim([0, 1]); ax.set_ylim([0, 1])
    ax.plot(0.2,0.5,marker='o', color=thresh_c,markersize=70,markeredgewidth=2, alpha=0.6)
ax1.text(0.4, 0.2, 'Green:\nNo flooding\nrisk', fontsize=25, color='w')
ax2.text(0.4, 0.2, 'Yellow:\nRisk of\nhigh water', fontsize=25, color='w')
ax3.text(0.4, 0.2, 'Red:\nExtreme risk\nof flooding', fontsize=25, color='w')


return fig

In [15]:

fig = figures.website_thumbnail(bathy, grid_T_hr, model_path, coastline, 0.1)

In [15]: