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"
cube = iris.load_cube(url)

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.

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/.

python4oceanographers

Turning ripples into waves

Using iris trajectory class to extract a bathymetric profile

Comments