# Using iris trajectory class to extract a bathymetric profile

Iris has a built-in Trajectory class that is very easy to use. You have to enter the waypoints and the number of linear spaced points between them. Imagine a np.linspace for lon, lat points. Below I will demonstrate with a simple use to extract bathymetric values from ETOPO2.

In [2]:
import iris

url = "http://opendap.ccst.inpe.br/Misc/etopo2/ETOPO2v2c_f4.nc"

print(cube)
z / (1)                             (y: 5400; x: 10800)
Dimension coordinates:
y                           x        -
x                           -        x
Attributes:
Conventions: COARDS
actual_range: [-10791.   8440.]
node_offset: 1
source:                         -Rd -I2m -ZTLf
title:

The cube has global data, let's reduce it to a smaller area.

In [3]:
lon = iris.Constraint(x=lambda cell: -54 <= cell <= -20)
lat = iris.Constraint(y=lambda cell: -30 <= cell <= -12)

cube = cube.extract(lon & lat)

print(cube)
z / (1)                             (y: 540; x: 1020)
Dimension coordinates:
y                           x       -
x                           -       x
Attributes:
Conventions: COARDS
actual_range: [-10791.   8440.]
node_offset: 1
source:                         -Rd -I2m -ZTLf
title:

Here is how to construct a trajectory.

In [4]:
from iris.analysis import trajectory

waypoints = [{'x': -38.5, 'y': -20.65}, {'x': -28.5, 'y': -20.65}]

#[(-38.594367271025789, -20.40487006995896),
#(-28.695012274806029, -20.531784877602803)]

traj = trajectory.Trajectory(waypoints, sample_count=30)

traj
Out[4]:
Trajectory([{'y': -20.65, 'x': -38.5}, {'y': -20.65, 'x': -28.5}], sample_count=30)
In [5]:
traj.sampled_points[:5]  # First 5 points
Out[5]:
[{'x': -38.5, 'y': -20.65},
{'x': -38.1551724137931, 'y': -20.65},
{'x': -37.810344827586206, 'y': -20.65},
{'x': -37.46551724137931, 'y': -20.65},
{'x': -37.12068965517241, 'y': -20.65}]
In [6]:
traj.length * 111  # Approximate distance in km.
Out[6]:
1110.0

Here is where iris is a little bit annoying. The trajectory.interpolate seems to be the intended consumer of this trajectory object, but interpolate takes a tuple with x and y points, while the Trajectory object is a list of dictionaries. I am sure there is a reason for that! (I just cannot see it ;-)

Anyways, we can quickly transform this into the input expected by interpolate.

In [7]:
lon = [d['x'] for d in traj.sampled_points]
lat = [d['y'] for d in traj.sampled_points]

sampled_points = [('x', lon),
('y', lat)]

section = trajectory.interpolate(cube, sampled_points)

print(section)
z / (1)                             (-- : 30)
Auxiliary coordinates:
x                             x
y                             x
Attributes:
Conventions: COARDS
actual_range: [-10791.   8440.]
node_offset: 1
source:                         -Rd -I2m -ZTLf
title:

Not that the interpolated section data is a cube. Lovely! Now all that is to plot the section we extracted.

In [8]:
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from palettable import colorbrewer
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

land = colorbrewer.get_map('Greens', 'sequential', 9)
ocean = colorbrewer.get_map('Blues', 'sequential', 9, reverse=True)

colors = np.array(ocean.mpl_colors + land.mpl_colors)
levels = [-9000, -6000, -5000, -4000, -3000, -2500, -1000, -700, -400, -145, -10, 0, 10, 145, 400, 800, 1200, 1600]

def make_map(projection=ccrs.PlateCarree()):
fig, ax = plt.subplots(figsize=(8, 8), 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
In [9]:
import seawater as sw
from matplotlib.ticker import MaxNLocator

lons = cube.coord('x').points
lats = cube.coord('y').points

dist = np.r_[0, sw.dist(lat, lon)[0].cumsum()]

fig, ax = make_map(projection=ccrs.PlateCarree())
offset = 1
ax.set_extent([lons.min()+offset, lons.max()-offset,
lats.min()+offset, lats.max()-offset])
cs = ax.contourf(lons, lats, cube.data, levels=levels, colors=colors)
l, = ax.plot(lon, lat, 'r', linewidth=2, alpha=0.7)

inset_ax = fig.add_axes([0.47, 0.4, 0.26, 0.125], frameon=False)
c = inset_ax.fill_between(dist, -section.data, y2=5000, color="none",
hatch="//", edgecolor='k', interpolate=True)
inset_ax.yaxis.tick_right()
inset_ax.axis([0, 1090, -400, 4000])
inset_ax.xaxis.tick_bottom()
inset_ax.yaxis.tick_right()
inset_ax.invert_yaxis()

inset_ax.yaxis.set_major_locator(MaxNLocator(4))
inset_ax.xaxis.set_major_locator(MaxNLocator(5))

If you are curious here is the feature I am using for the example.

Have fun!

In [10]:
HTML(html)
Out[10]:

This post was written as an IPython notebook. It is available for download or as a static html.