After copying-and-pasting so many of my iris code snippets I decided to package everything in a new module: TARDIS.
Tardis is nothing more than a hackish wrapper around iris that expands some of the iris methods to work with 2D-Coords and to facilitate some of its methods calls.
Virtually everything that tardis does could be accomplished using the raw
iris API, but it will be a little bit easier with tardis ;-)
For example, how to load cubes based on a variable list of CF standard names?
from tardis import quick_load_cubes
name_list = ['sea_surface_height',
'sea_surface_elevation',
'sea_surface_height_above_geoid',
'sea_surface_height_above_sea_level',
'water_surface_height_above_reference_datum',
'sea_surface_height_above_reference_ellipsoid']
url = ("http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/his_Best/"
"ESPRESSO_Real-Time_v2_History_Best_Available_best.ncd")
cube = quick_load_cubes(url, name_list, callback=None, strict=True)
Aside: Note all those warnings above? That is one of the reasons why it is complicated to deal with all ROMS variables as a single dataset. Hopefully pysgrid will save us =)
Let's check this cube out.
print(cube)
The module is called tardis because it allows you to "travel" in space, time, and (variables) dimension.
The workhorse function is proc_cube. This function takes a cube and a
series of optional constrains to "conform" a cube into the desired time span,
space bounding box, and/or automatic units conversion.
import iris
import pytz
from datetime import datetime, timedelta
from tardis import proc_cube
# Time.
stop = datetime(2014, 7, 7, 12)
stop = stop.replace(tzinfo=pytz.utc)
start = stop - timedelta(days=1)
# Space.
bbox = [-87.40, 24.25, -74.70, 36.70]
# Units.
units = iris.unit.Unit('meters')
# New cube.
cube = proc_cube(cube, bbox=bbox, time=(start, stop), units=units)
print(cube)
One of the advantages of using iris is the ability to save the data in a CF-compliant netCDF file.
iris.save(cube, "subset.nc")
cube = iris.load_cube("subset.nc")
print(cube)
or,
!ncdump -h subset.nc
The tardis module jewel is the get_nearest_water function. (Which I just
realized should be renamed to get_valid_data!) Together with make_tree
it is possible to easily find any valid data point near to a certain
position.
Note that iris nearest methods only works with 1D-Coords grids, while
tardis will work with virtually any kind of grid: Ugrid, ROMS-Curvilinear
(2D Coords), etc.
Right now make_tree uses Scipy's KDTree, but I am re-writing it to use
rtree instead. In addition the tree
object will be folded into get_nearest_waterget_valid_data and
cached, so we can skip one function call.
from tardis import get_nearest_water, make_tree
tree, lon, lat = make_tree(cube)
obs = dict(lon=-76.67, lat=34.72)
kw = dict(k=10, max_dist=0.04, min_var=0.01)
series, dist, idx = get_nearest_water(cube, tree, obs['lon'], obs['lat'], **kw)
Here is the result.
%matplotlib inline
from matplotlib import style
style.use('ggplot')
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
name='admin_1_states_provinces_shp')
def make_map(projection=ccrs.PlateCarree()):
fig, ax = plt.subplots(figsize=(8, 6),
subplot_kw=dict(projection=projection))
ax.coastlines(resolution='50m')
gl = ax.gridlines(draw_labels=True)
gl.xlabels_top = gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
return fig, ax
import numpy.ma as ma
# Latest time stamp for contouring.
c = cube[-1, ...]
c.data = ma.masked_invalid(c.data)
lon = c.coord(axis='X').points
lat = c.coord(axis='Y').points
extent = (lon.min(), lon.max(),
lat.min(), lat.max())
from oceans import cm
from matplotlib.ticker import MaxNLocator
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
fig, ax = make_map()
ax.set_extent(extent)
cs = ax.pcolormesh(lon, lat, c.data, cmap=cm.rscolmap, alpha=0.5)
ax.plot(obs['lon'], obs['lat'], 'o', color='darkgreen', label='Beaufort, NC', alpha=0.5)
ax.plot(lon[idx], lat[idx], 'o', color='crimson', label='ESPRESSO', alpha=0.5)
ax.set_title(c.attributes.get('title', None))
ax.add_feature(states, edgecolor='gray')
leg = ax.legend(numpoints=1, fancybox=True, framealpha=0, loc="upper left")
axins = zoomed_inset_axes(ax, 20, loc=5,
bbox_to_anchor=(1.1, 0.75),
bbox_transform=ax.figure.transFigure)
axins.pcolormesh(lon, lat, c.data, cmap=cm.rscolmap, alpha=0.5)
axins.plot(obs['lon'], obs['lat'], 'o', color='darkgreen', label='Beaufort, NC', alpha=0.5)
axins.plot(lon[idx], lat[idx], 'o', color='crimson', label='ESPRESSO', alpha=0.5)
axins.axis([-76.7, -76.64,
34.67, 34.73])
mark_inset(ax, axins, loc1=2, loc2=4, ec="0.5")
# Remove ticks
_ = axins.axes.get_xaxis().set_ticks([])
_ = axins.axes.get_yaxis().set_ticks([])
That way we can easily find the nearest ROMS grid point to the Beaufort, NC station.
Let's take a look at the time-series at that point.
import iris.quickplot as qplt
fig, ax = plt.subplots(figsize=(9, 2.75))
l, = qplt.plot(series)
ax.grid(True)
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.