It has been a while that I wanted to try out windspharm with gridded ocean velocities, but windspharm was available only for Python 2. Well, I got to work and after a few commits the PR for python 3 got merged.

windspharm uses pyspharm, a python wrapper around SPHEREPACK. SPHEREPACK is a collection of functions for modeling geophysical processes.

With windspharm one can calculate:

divergence
vorticity (relative and absolute)
streamfunction
velocity potential
irrotational and non-divergent components (Helmholtz decomposition)
vector gradient of a scalar function
magnitude (speed)

Let's try to calculate the Rossby Wave Source (RWS) term by Sardeshmukh and Hoskins (1988) with windspharm.

In [2]:

from netCDF4 import Dataset

I will use the Mean Dynamic Topographic from a previous post, but here I will reduce the data density a little before the computations and plotting (that reminds me I need a new computer!).

In [3]:

dx = dy = slice(0, -1, 10)

with Dataset('./data/mdt_cnes_cls2009_global_v1.1.nc') as nc:
    h = nc.variables['Grid_0001'][dx, dy]
    u = nc.variables['Grid_0002'][dx, dy]
    v = nc.variables['Grid_0003'][dx, dy]
    lat = nc.variables['NbLatitudes'][dy]
    lon = nc.variables['NbLongitudes'][dx]
    
mask = u.mask

Notice that I stored the continents mask, we will use it later.

windspharm requires latitude and longitude to be the leading dimensions, and that the velocity components must be at least 2D. There are tools in the windspharm module that makes re-shaping the data a lot easier.

In [4]:

from windspharm.tools import prep_data

In [5]:

u, info = prep_data(u, 'xy')
v, _ = prep_data(v, 'xy')

By the way, I loved the declarative way to reshape using txy (time, longitude, latitude). Makes the code much more readable. There are even tools that read the info dict and get your data back into the old shape.

In [6]:

info

Out[6]:

{'intermediate_order': 'yx',
 'intermediate_shape': (72, 144),
 'original_order': 'xy'}

It is also required that the latitude is north-to-south instead of south-to-north. Again windspharm bring tools to make this easy.

In [7]:

from windspharm.tools import order_latdim

In [8]:

lat, u, v = order_latdim(lat, u, v)

The VectorWind object is windspharm "workhorse". There are 3 version, for iris, cdms and standard NumPy arrays. Here I'm using the standard interface because iris is not py3k compatible (yet).

In [9]:

from windspharm.standard import VectorWind

In [10]:

U = VectorWind(u.filled(fill_value=0), v.filled(fill_value=0))

Now we can compute all the components of the RWS: absolute vorticity, divergence, irrotational (divergent) velocity components, and gradients of absolute vorticity.

In [11]:

eta = U.absolutevorticity()
div = U.divergence()
uchi, vchi = U.irrotationalcomponent()
etax, etay = U.gradient(eta)

And a one-liners to combine the components and form the RWS term.

In [12]:

S = -eta * div - (uchi * etax + vchi * etay)

Below is just a simple plotting function.

In [13]:

import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap, addcyclic


def make_map(projection='robin', resolution='c'):
    fig, ax = plt.subplots(figsize=(10, 6))
    m = Basemap(projection='robin', resolution='c',
                lon_0=0, ax=ax)
    m.drawcoastlines()
    m.fillcontinents(color='0.85')
    parallels = np.arange(-60, 90, 30.)
    meridians = np.arange(-360, 360, 60.)
    m.drawparallels(parallels, labels=[1, 0, 0, 0])
    m.drawmeridians(meridians, labels=[0, 0, 0, 1])
    return fig, m

VectorWind does not accept missing values (NaN) nor masked arrays, that is why I had to fill the continents missing data with zeros. It is a dangerous strategy and should be used with care pay attention to:

Never fill missing that as land;
Careful interpretation of the results near the border and,
Always re-apply the same mask.

In [14]:

import numpy as np
import numpy.ma as ma

S = S.squeeze()
S = ma.masked_array(S, np.flipud(mask.T))
S, lon = addcyclic(S, lon)  # Add a cyclic point for plotting.
lon, lat = np.meshgrid(lon, lat)

And Voilà

In [15]:

from palettable import colorbrewer

clevs = [-10, -8, -4, -2, -1, 0, 1, 2, 4, 8, 10]
cmap = colorbrewer.get_map('RdBu', 'diverging',
                           len(clevs), reverse=True).mpl_colormap

fig, m = make_map()
cs = m.contourf(lon, lat, S*1e11, clevs, cmap=cmap,
                extend='both', latlon=True)
cbar = fig.colorbar(cs, orientation='vertical', shrink=0.7, extend='both')
m.ax.set_title('Rossby Wave Source', fontsize=16)
t = cbar.ax.set_title('$10^{-11}$s$^{-1}$')

In [16]:

HTML(html)

Out[16]:

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

Rossby wave source term with windspharm

Comments