%matplotlib inline
import time
import contextlib
from datetime import datetime, timedelta
import numpy as np
import numpy.ma as ma
import matplotlib.tri as tri
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
import iris
from iris.unit import Unit
from iris.exceptions import CoordinateNotFoundError
import cartopy.crs as ccrs
from cartopy.feature import NaturalEarthFeature, COLORS
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
LAND = NaturalEarthFeature('physical', 'land', '10m', edgecolor='face',
facecolor=COLORS['land'])
iris.FUTURE.netcdf_promote = True
iris.FUTURE.cell_datetime_objects = True # <- TODO!
def time_coord(cube):
"""Return the variable attached to time axis and rename it to time."""
try:
cube.coord(axis='T').rename('time')
except CoordinateNotFoundError:
pass
timevar = cube.coord('time')
return timevar
def z_coord(cube):
"""Heuristic way to return the dimensionless vertical coordinate."""
try:
z = cube.coord(axis='Z')
except CoordinateNotFoundError:
z = cube.coords(axis='Z')
for coord in cube.coords(axis='Z'):
if coord.ndim == 1:
z = coord
return z
def time_near(cube, datetime):
"""Return the nearest index to a `datetime`."""
timevar = time_coord(cube)
try:
time = timevar.units.date2num(datetime)
idx = timevar.nearest_neighbour_index(time)
except IndexError:
idx = -1
return idx
def time_slice(cube, start, stop=None):
"""TODO: Re-write to use `iris.FUTURE.cell_datetime_objects`."""
istart = time_near(cube, start)
if stop:
istop = time_near(cube, stop)
if istart == istop:
raise ValueError('istart must be different from istop!'
'Got istart {!r} and '
' istop {!r}'.format(istart, istop))
return cube[istart:istop, ...]
else:
return cube[istart, ...]
def plot_surface(cube, model='', unstructure=False, **kw):
projection = kw.pop('projection', ccrs.PlateCarree())
figsize = kw.pop('figsize', (8, 6))
cmap = kw.pop('cmap', plt.cm.rainbow)
fig, ax = plt.subplots(figsize=figsize,
subplot_kw=dict(projection=projection))
ax.set_extent(get_bbox(cube))
ax.add_feature(LAND)
ax.coastlines(resolution='10m')
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
z = z_coord(cube)
if z:
positive = z.attributes.get('positive', None)
if positive == 'up':
idx = np.unique(z.points.argmax(axis=0))[0]
else:
idx = np.unique(z.points.argmin(axis=0))[0]
c = cube[idx, ...].copy()
else:
idx = None
c = cube.copy()
c.data = ma.masked_invalid(c.data)
t = time_coord(cube)
t = t.units.num2date(t.points)[0]
if unstructure:
# The following lines would work if the cube is note bbox-sliced.
# lon = cube.mesh.nodes[:, 0]
# lat = cube.mesh.nodes[:, 1]
# nv = cube.mesh.faces
lon = cube.coord(axis='X').points
lat = cube.coord(axis='Y').points
nv = Delaunay(np.c_[lon, lat]).vertices
triang = tri.Triangulation(lon, lat, triangles=nv)
# http://matplotlib.org/examples/pylab_examples/ tricontour_smooth_delaunay.html
if False: # TODO: Test this.
subdiv = 3
min_circle_ratio = 0.01
mask = tri.TriAnalyzer(triang).get_flat_tri_mask(min_circle_ratio)
triang.set_mask(mask)
refiner = tri.UniformTriRefiner(triang)
tri_ref, data_ref = refiner.refine_field(cube.data, subdiv=subdiv)
cs = ax.tricontourf(triang, c.data, cmap=cmap, **kw)
else:
cs = ax.pcolormesh(c.coord(axis='X').points,
c.coord(axis='Y').points,
c.data, cmap=cmap, **kw)
title = (model, t, c.name(), idx)
ax.set_title('{}: {}\nVariable: {} level: {}'.format(*title))
return fig, ax, cs
def get_bbox(cube):
xmin = cube.coord(axis='X').points.min()
xmax = cube.coord(axis='X').points.max()
ymin = cube.coord(axis='Y').points.min()
ymax = cube.coord(axis='Y').points.max()
return [xmin, xmax, ymin, ymax]
@contextlib.contextmanager
def timeit(log=None):
t = time.time()
yield
elapsed = time.strftime("%H:%M:%S", time.gmtime(time.time()-t))
if log:
log.info(elapsed)
else:
print(elapsed)
model = 'NECOFS_FVCOM'
start = datetime.utcnow() - timedelta(days=7)
bbox = [-70.8, 41.4, -69.9, 42.3]
units = Unit('Kelvin')
with timeit():
url = "http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/"
url += "Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc"
cube = iris.load_cube(url, 'sea_water_potential_temperature')
cube = time_slice(cube, start, None)
cube.convert_units(units)
print(cube)
fig, ax, cs = plot_surface(cube, model, unstructure=True)
cbar = fig.colorbar(cs, extend='both', shrink=0.75)
t = cbar.ax.set_title(cube.units)
/home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1004: UserWarning: Ignoring variable u'siglay' referenced by variable u'ww': Dimensions (u'siglay', u'node') do not span (u'time', u'siglay', u'nele') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1004: UserWarning: Ignoring variable u'siglay' referenced by variable u'u': Dimensions (u'siglay', u'node') do not span (u'time', u'siglay', u'nele') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/cf.py:1004: UserWarning: Ignoring variable u'siglay' referenced by variable u'v': Dimensions (u'siglay', u'node') do not span (u'time', u'siglay', u'nele') warnings.warn(msg) /home/usgs/anaconda/lib/python2.7/site-packages/cartopy/io/__init__.py:261: DownloadWarning: Downloading: http://www.nacis.org/naturalearth/10m/physical/ne_10m_land.zip warnings.warn('Downloading: {}'.format(url), DownloadWarning)
sea_water_potential_temperature / (Kelvin) (-- : 10; -- : 98432) Auxiliary coordinates: ocean_sigma_coordinate x x latitude - x longitude - x sea_floor_depth_below_geoid - x sea_surface_height_above_geoid - x Scalar coordinates: time: 2014-09-16 00:00:00 Attributes: Conventions: CF-1.0 CoordinateProjection: init=nad83:1802 CoordinateSystem: Cartesian GroundWater_Forcing: GROUND WATER FORCING IS OFF! River_Forcing: THERE ARE 18 RIVERS IN THIS MODEL. RIVER INFLOW IS ON THE nodes WHERE TEMPERATURE... Surface_Heat_Forcing: FVCOM variable surface heat forcing file: FILE NAME:wrf_for.nc SOURCE:wrf2fvcom... Surface_PrecipEvap_Forcing: FVCOM periodic surface precip forcing: FILE NAME:wrf_for.nc SOURCE:wrf2fvcom... Surface_Wind_Forcing: FVCOM variable surface Wind forcing: FILE NAME:wrf_for.nc SOURCE:wrf2fvcom... Tidal_Forcing: TIDAL ELEVATION FORCING IS OFF! _DODS_Unlimited_Dimension: time cdm_data_type: any coverage_content_type: modelResult grid: fvcom_grid history: Fri Sep 19 09:31:25 2014: ncrcat -O -v x,y,lat,lon,xc,yc,lonc,latc,siglay,siglev,nv,nbe,aw0,awx,awy,h,temp,salinity,u,v,ww,ua,va,zeta,Times... institution: School for Marine Science and Technology location: node mesh: fvcom_mesh nco_openmp_thread_number: 1 references: http://fvcom.smast.umassd.edu, http://codfish.smast.umassd.edu source: FVCOM_3.0 summary: Latest forecast from the FVCOM Northeast Coastal Ocean Forecast System... title: NECOFS Massachusetts (FVCOM) - Massachusetts Coastal - Latest Forecast type: data 00:00:05
/home/usgs/anaconda/lib/python2.7/site-packages/IPython/core/formatters.py:239: FormatterWarning: Exception in image/png formatter: HTTP Error 404: Not Found FormatterWarning,
<matplotlib.figure.Figure at 0x7f4189009e90>
X
and Y
the subset works.¶with timeit():
url = "http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/"
url += "Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc"
cube = iris.load_cube(url, 'sea_water_potential_temperature')
cube = time_slice(cube, start, None)
cube.convert_units(units)
print(cube.coord(axis='Y'))
print(cube.coord(axis='X'))
print(cube.coord(axis='Z'))
print("\n")
cube = cube.intersection(longitude=(bbox[0], bbox[2]),
latitude=(bbox[1], bbox[3]))
print(cube)
fig, ax, cs = plot_surface(cube, model, unstructure=True)
cbar = fig.colorbar(cs, extend='both', shrink=0.75)
t = cbar.ax.set_title(cube.units)
AuxCoord(array([ 43.30540466, 43.29943466, 43.29204559, ..., 42.37276459, 42.37341309, 42.37305069], dtype=float32), standard_name=u'latitude', units=Unit('degrees'), long_name=u'nodal latitude', var_name='lat') AuxCoord(array([-70.56564331, -70.54589081, -70.52242279, ..., -71.13265991, -71.13302612, -71.13330841], dtype=float32), standard_name=u'longitude', units=Unit('degrees'), long_name=u'nodal longitude', var_name='lon') AuxCoord(array([[-0.05000001, -0.04501463, -0.03064007, ..., -0.05000001, -0.05000001, -0.05000001], [-0.15000001, -0.13504389, -0.09192021, ..., -0.15000001, -0.15000001, -0.15000001], [-0.25000003, -0.22507316, -0.15320036, ..., -0.25000003, -0.25000003, -0.25000003], ..., [-0.75000006, -0.73931712, -0.70851445, ..., -0.75000006, -0.75000006, -0.75000006], [-0.85000002, -0.84359032, -0.82510877, ..., -0.85000002, -0.85000002, -0.85000002], [-0.95000005, -0.94786352, -0.94170296, ..., -0.95000005, -0.95000005, -0.95000005]], dtype=float32), standard_name=u'ocean_sigma_coordinate', units=Unit('unknown'), long_name=u'Sigma Layers', var_name='siglay', attributes={'positive': 'up'}) sea_water_potential_temperature / (Kelvin) (-- : 10; -- : 44575) Auxiliary coordinates: ocean_sigma_coordinate x x latitude - x longitude - x sea_floor_depth_below_geoid - x sea_surface_height_above_geoid - x Scalar coordinates: time: 2014-09-16 00:00:00 Attributes: Conventions: CF-1.0 CoordinateProjection: init=nad83:1802 CoordinateSystem: Cartesian GroundWater_Forcing: GROUND WATER FORCING IS OFF! River_Forcing: THERE ARE 18 RIVERS IN THIS MODEL. RIVER INFLOW IS ON THE nodes WHERE TEMPERATURE... Surface_Heat_Forcing: FVCOM variable surface heat forcing file: FILE NAME:wrf_for.nc SOURCE:wrf2fvcom... Surface_PrecipEvap_Forcing: FVCOM periodic surface precip forcing: FILE NAME:wrf_for.nc SOURCE:wrf2fvcom... Surface_Wind_Forcing: FVCOM variable surface Wind forcing: FILE NAME:wrf_for.nc SOURCE:wrf2fvcom... Tidal_Forcing: TIDAL ELEVATION FORCING IS OFF! _DODS_Unlimited_Dimension: time cdm_data_type: any coverage_content_type: modelResult grid: fvcom_grid history: Fri Sep 19 09:31:25 2014: ncrcat -O -v x,y,lat,lon,xc,yc,lonc,latc,siglay,siglev,nv,nbe,aw0,awx,awy,h,temp,salinity,u,v,ww,ua,va,zeta,Times... institution: School for Marine Science and Technology location: node mesh: fvcom_mesh nco_openmp_thread_number: 1 references: http://fvcom.smast.umassd.edu, http://codfish.smast.umassd.edu source: FVCOM_3.0 summary: Latest forecast from the FVCOM Northeast Coastal Ocean Forecast System... title: NECOFS Massachusetts (FVCOM) - Massachusetts Coastal - Latest Forecast type: data 00:00:05
<matplotlib.figure.Figure at 0x7f4187b31b50>
with timeit():
url = "http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/"
url += "Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc"
cube = iris.load_cube(url, 'sea_water_potential_temperature')
cube = time_slice(cube, start, None)
cube.convert_units(units)
cube = cube.intersection(longitude=(bbox[0], bbox[2]),
latitude=(bbox[1], bbox[3]))
print(cube)
fig, ax, cs = plot_surface(cube, model, unstructure=True)
cbar = fig.colorbar(cs, extend='both', shrink=0.75)
t = cbar.ax.set_title(cube.units)
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) <ipython-input-5-856ee578308c> in <module>() 7 cube.convert_units(units) 8 cube = cube.intersection(longitude=(bbox[0], bbox[2]), ----> 9 latitude=(bbox[1], bbox[3])) 10 print(cube) 11 /home/usgs/anaconda/lib/python2.7/site-packages/iris/cube.pyc in intersection(self, *args, **kwargs) 2051 result = result._intersect(*arg) 2052 for name, value in kwargs.iteritems(): -> 2053 result = result._intersect(name, *value) 2054 return result 2055 /home/usgs/anaconda/lib/python2.7/site-packages/iris/cube.pyc in _intersect(self, name_or_coord, minimum, maximum, min_inclusive, max_inclusive) 2071 minimum, maximum, 2072 min_inclusive, -> 2073 max_inclusive) 2074 2075 # By this point we have either one or two subsets along the relevant /home/usgs/anaconda/lib/python2.7/site-packages/iris/cube.pyc in _intersect_modulus(self, coord, minimum, maximum, min_inclusive, max_inclusive) 2153 values = coord.bounds 2154 else: -> 2155 values = coord.points 2156 if values.max() > values.min() + modulus: 2157 raise ValueError("coordinate's range greater than coordinate's" /home/usgs/anaconda/lib/python2.7/site-packages/iris/coords.pyc in points(self) 1491 points = self._points 1492 if isinstance(points, biggus.Array): -> 1493 points = points.ndarray() 1494 self._points = points 1495 return points.view() /home/usgs/anaconda/lib/python2.7/site-packages/biggus/__init__.pyc in ndarray(self) 780 781 def ndarray(self): --> 782 array = self._apply_keys() 783 # We want the shape of the result to match the shape of the 784 # Array, so where we've ended up with an array-scalar, /home/usgs/anaconda/lib/python2.7/site-packages/biggus/__init__.pyc in _apply_keys(self) 877 """ 878 def _apply_keys(self): --> 879 array = self.concrete.__getitem__(self._keys) 880 return array 881 /home/usgs/anaconda/lib/python2.7/site-packages/iris/fileformats/netcdf.pyc in __getitem__(self, keys) 247 variable = dataset.variables[self.variable_name] 248 # Get the NetCDF variable data and slice. --> 249 data = variable[keys] 250 finally: 251 dataset.close() /home/usgs/anaconda/lib/python2.7/site-packages/netCDF4.so in netCDF4.Variable.__getitem__ (netCDF4.c:34997)() /home/usgs/anaconda/lib/python2.7/site-packages/netCDF4.so in netCDF4.Variable._get (netCDF4.c:41689)() RuntimeError: NetCDF: I/O failure