Glider


In [1]:
import iris

iris.FUTURE.netcdf_promote = True

url = ('http://tds.marine.rutgers.edu:8080/thredds/dodsC/'
       'cool/glider/mab/Gridded/20130911T000000_20130920T000000_gp2013_modena.nc')

glider = iris.load(url)

lon = glider.extract_strict('Longitude').data
lat = glider.extract_strict('Latitude').data
glider = glider.extract_strict('Temperature')
depth = glider.coord('depth').points


/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1301: UserWarning: Ignoring netCDF variable 'salinity' invalid units 'psu'
  warnings.warn(msg.format(msg_name, msg_units))

In [2]:
import numpy as np
import numpy.ma as ma
import seawater as sw
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

from utilities import time_coord

%matplotlib inline


def plot_glider(cube, mask_topo=False, track_inset=False, **kw):
    """Plot glider cube."""
    cmap = kw.pop('cmap', plt.cm.rainbow)
    
    data = ma.masked_invalid(cube.data.squeeze())
    t = time_coord(cube)
    #t = t.units.num2date(t.points.squeeze())
    
    dist, pha = sw.dist(lat, lon, units='km')
    dist = np.r_[0, np.cumsum(dist)]
    
    dist, z = np.broadcast_arrays(dist[..., None], depth)

    try:
        z_range = cube.coord(axis='Z').attributes['actual_range']
    except KeyError:
        z_range = z.min(), z.max()
    try:
        data_range = cube.attributes['actual_range']
    except KeyError:        
        data_range = data.min(), data.max()
    
    condition = np.logical_and(data >= data_range[0], data <= data_range[1])
    data = ma.masked_where(~condition, data)
    
    condition = np.logical_and(z >= z_range[0], z <= z_range[1])
    z = ma.masked_where(~condition, z)

    fig, ax = plt.subplots(figsize=(9, 3.75))
    cs = ax.pcolor(dist, z, data, cmap=cmap, snap=True, **kw)
    if mask_topo:
        h = z.max(axis=1)
        x = dist[:, 0]
        ax.plot(x, h, color='black', linewidth='0.5', zorder=3)
        ax.fill_between(x, h, y2=h.max(), color='0.9', zorder=3)
    #ax.set_title('Glider track from {} to {}'.format(t[0], t[-1]))
    fig.tight_layout()
    
    if track_inset:
        axin = inset_axes(ax, width="25%", height="30%", loc=4)
        axin.plot(lon, lat, 'k.')
        start, end = (lon[0], lat[0]), (lon[-1], lat[-1])
        kw = dict(marker='o', linestyle='none')
        axin.plot(*start, color='g', **kw)
        axin.plot(*end, color='r', **kw)
        axin.axis('off')
    return fig, ax, cs

Models


In [3]:
from utilities import CF_names, quick_load_cubes

models = dict(useast=('http://ecowatch.ncddc.noaa.gov/thredds/dodsC/'
                      'ncom_us_east_agg/US_East_Apr_05_2013_to_Current_best.ncd'),
              hycom=('http://ecowatch.ncddc.noaa.gov/thredds/dodsC/'
                     'hycom/hycom_reg1_agg/HYCOM_Region_1_Aggregation_best.ncd'),
              sabgom=('http://omgsrv1.meas.ncsu.edu:8080/thredds/dodsC/'
                      'fmrc/sabgom/SABGOM_Forecast_Model_Run_Collection_best.ncd'),
              coawst=('http://geoport.whoi.edu/thredds/dodsC/'
                      'coawst_4/use/fmrc/coawst_4_use_best.ncd'))


name_list = CF_names['sea_water_temperature']

coawst = quick_load_cubes(models['coawst'], name_list, strict=True)
useast = quick_load_cubes(models['useast'], name_list, strict=True)
hycom = quick_load_cubes(models['hycom'], name_list, strict=True)


/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/cf.py:1040: UserWarning: Ignoring formula terms variable u'h' referenced by data variable u'v' via variable u's_rho': Dimensions (u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_v', u'xi_v')
  warnings.warn(msg)
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/cf.py:1040: UserWarning: Ignoring formula terms variable u'zeta' referenced by data variable u'v' via variable u's_rho': Dimensions (u'time', u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_v', u'xi_v')
  warnings.warn(msg)
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/cf.py:1040: UserWarning: Ignoring formula terms variable u'h' referenced by data variable u'u' via variable u's_rho': Dimensions (u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_u', u'xi_u')
  warnings.warn(msg)
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/cf.py:1040: UserWarning: Ignoring formula terms variable u'zeta' referenced by data variable u'u' via variable u's_rho': Dimensions (u'time', u'eta_rho', u'xi_rho') do not span (u'time', u's_rho', u'eta_u', u'xi_u')
  warnings.warn(msg)
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1396: UserWarning: Gracefully filling 'time' dimension coordinate masked points
  warnings.warn(msg.format(str(cf_coord_var.cf_name)))
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/netcdf.py:441: UserWarning: Unable to find coordinate for variable u'zeta'
  '{!r}'.format(name))
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/netcdf.py:441: UserWarning: Unable to find coordinate for variable u'h'
  '{!r}'.format(name))
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/netcdf.py:558: UserWarning: Unable to construct Ocean s-coordinate, generic form 1 factory due to insufficient source coordinates.
  warnings.warn('{}'.format(e))
/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/_pyke_rules/compiled_krb/fc_rules_cf_fc.py:1301: UserWarning: Ignoring netCDF variable 'surf_salt_flux' invalid units 'psu-m/s'
  warnings.warn(msg.format(msg_name, msg_units))

In [4]:
from datetime import datetime
from utilities import proc_cube

# Glider info.
start = glider.coord(axis='T').attributes['minimum']
stop = glider.coord(axis='T').attributes['maximum']

start = datetime.strptime(start, '%Y-%m-%d %H:%M:%S') 
stop = datetime.strptime(stop, '%Y-%m-%d %H:%M:%S')

bbox = lon.min(), lat.min(), lon.max(), lat.max()

# Subsetting the cube to the glider limits.
coawst = proc_cube(coawst, bbox=bbox, time=(start, stop), units=glider.units)
useast = proc_cube(useast, bbox=bbox, time=(start, stop), units=glider.units)
hycom = proc_cube(hycom, bbox=bbox, time=(start, stop), units=glider.units)

coawst, useast, hycom


Out[4]:
(<iris 'Cube' of sea_water_potential_temperature / (Celsius) (time: 209; ocean_s_coordinate_g1: 16; -- : 8; -- : 11)>,
 <iris 'Cube' of sea_water_temperature / (degC) (time: 69; depth: 40; latitude: 12; longitude: 14)>,
 <iris 'Cube' of sea_water_temperature / (degC) (time: 67; depth: 40; latitude: 5; longitude: 6)>)

In [5]:
for aux in coawst.aux_factories:
    coawst.remove_aux_factory(aux)

In [6]:
from iris.analysis import trajectory

sample_points = [('latitude', lat),
                 ('longitude', lon),
                 ('time', glider.coord(axis='T').points)]

In [7]:
depth = glider.coord('depth').points
fig, ax, cs = plot_glider(glider, mask_topo=False, track_inset=True)



In [8]:
iuseast = trajectory.interpolate(useast, sample_points)
iuseast.transpose()

depth = -iuseast.coord(axis='Z').points

fig, ax, cs = plot_glider(iuseast, mask_topo=False, track_inset=True)
ax.set_ylim(-120, 0)
t = ax.set_title("USEAST")



In [9]:
ihycom = trajectory.interpolate(hycom, sample_points)
ihycom.transpose()

depth = -ihycom.coord(axis='Z').points

fig, ax, cs = plot_glider(ihycom, mask_topo=False, track_inset=True)
ax.set_ylim(-120, 0)
t = ax.set_title("HYCOM")



In [10]:
icoawst = trajectory.interpolate( coawst, sample_points)

icoawst.transpose()

depth = -icoawst.coord(axis='Z').points

fig, ax, cs = plot_glider(ihycom, mask_topo=False, track_inset=True)
ax.set_ylim(-120, 0)
t = ax.set_title("COAWST")


---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-10-89d02ffa5123> in <module>()
----> 1 icoawst = trajectory.interpolate( coawst, sample_points)
      2 
      3 icoawst.transpose()
      4 
      5 depth = -icoawst.coord(axis='Z').points

/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/analysis/trajectory.py in interpolate(cube, sample_points, method)
    242             column_index = iris.analysis.interpolate._nearest_neighbour_indices_ndcoords(cube, point, cache=cache)
    243             column = cube[column_index]
--> 244             new_cube.data[..., i] = column.data
    245 
    246         # Fill in the empty squashed (non derived) coords.

/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/cube.py in data(self)
   1482         if not isinstance(data, np.ndarray):
   1483             try:
-> 1484                 data = data.masked_array()
   1485             except MemoryError:
   1486                 msg = "Failed to create the cube's data as there was not" \

/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/biggus/__init__.py in masked_array(self)
    824 
    825     def masked_array(self):
--> 826         array = self._apply_keys()
    827         # We want the shape of the result to match the shape of the
    828         # Array, so where we've ended up with an array-scalar,

/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/biggus/__init__.py in _apply_keys(self)
    950     """
    951     def _apply_keys(self):
--> 952         array = self.concrete.__getitem__(self._keys)
    953         return array
    954 

/home/filipe/.virtualenvs/iris/lib/python2.7/site-packages/iris/fileformats/netcdf.py in __getitem__(self, keys)
    269             variable = dataset.variables[self.variable_name]
    270             # Get the NetCDF variable data and slice.
--> 271             data = variable[keys]
    272         finally:
    273             dataset.close()

netCDF4.pyx in netCDF4.Variable.__getitem__ (netCDF4.c:42664)()

netCDF4.pyx in netCDF4.Variable._get (netCDF4.c:51595)()

netCDF4.pyx in netCDF4.Variable.dimensions.__get__ (netCDF4.c:38554)()

netCDF4.pyx in netCDF4.Variable._getdims (netCDF4.c:37622)()

/home/filipe/.virtualenvs/iris/lib64/python2.7/encodings/utf_8.pyc in decode(input, errors)
     13 encode = codecs.utf_8_encode
     14 
---> 15 def decode(input, errors='strict'):
     16     return codecs.utf_8_decode(input, errors, True)
     17 

KeyboardInterrupt: