In [ ]:
import os
import sys

root = os.path.abspath(os.path.join(os.getcwd(), os.pardir, os.pardir))
sys.path.append(root)

In [ ]:
from datetime import datetime, timedelta

bbox = [-87.4, 24.25, -74.7, 36.70]  # SECOORA: NC, SC GA, FL

time = datetime.utcnow() - timedelta(hours=6)

HF radar


In [ ]:
from iris.unit import Unit
from utilities import get_cubes


def squeeze(cubes):
    if len(cubes) == 1:
        return cubes[0]
    else:
        msg = "> 1 cubes found.  Expected just one.\n {!r}".format
        raise ValueError(msg(res))


url = 'http://hfrnet.ucsd.edu/thredds/dodsC/HFRNet/USEGC/6km/hourly/RTV'

u_name_list = ['x_sea_water_velocity',
               'sea_water_x_velocity',
               'eastward_sea_water_velocity',
               'surface_eastward_sea_water_velocity']

v_name_list = ['y_sea_water_velocity',
               'sea_water_y_velocity',
               'northward_sea_water_velocity',
               'surface_northward_sea_water_velocity']

kw = dict(time=time, bbox=bbox, units=Unit('m s-1'))
u_cube = squeeze(get_cubes(url, name_list=u_name_list, **kw))
v_cube = squeeze(get_cubes(url, name_list=v_name_list, **kw))

In [ ]:
import numpy as np


def subsample_cubes(u_cube, v_cube, n=6):
    """This is just for plotting.  If using cartopy do not
    subsample.  Use cartopy's resample instead."""
    u = u_cube[::n, ::n]
    v = v_cube[::n, ::n]
    lon = u.coord(axis='X').points
    lat = u.coord(axis='Y').points
    if (lon.ndim == 1) and (lat.ndim == 1):
        lon, lat = np.meshgrid(lon, lat)
    time = u.coord('time')
    time.units.num2date(time.points)
    return dict(lon=lon, lat=lat, u=u.data, v=v.data, time=time)


hf_radar = subsample_cubes(u_cube, v_cube, n=6)

SABGOM model


In [ ]:
from utilities import get_roms


model_url = ('http://omgsrv1.meas.ncsu.edu:8080/thredds/dodsC/fmrc/sabgom/'
             'SABGOM_Forecast_Model_Run_Collection_best.ncd')

model = get_roms(model_url, time, n=6)

In [ ]:
import mplleaflet
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
kw = dict(scale=20, headwidth=2, alpha=0.5)

mask = ~hf_radar['u'].mask  # mplleaflet does not support masked array.

Q = ax.quiver(hf_radar['lon'][mask], hf_radar['lat'][mask],
              hf_radar['u'][mask], hf_radar['v'][mask], color='black', **kw)

mask = ~model['u'].mask

Q = ax.quiver(model['lon'][mask], model['lat'][mask],
              model['u'][mask], model['v'][mask], color='red', **kw)

mplleaflet.display(fig=ax.figure)