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.

```
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!).

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

```
from windspharm.tools import prep_data
```

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

```
info
```

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

```
from windspharm.tools import order_latdim
```

```
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).

```
from windspharm.standard import VectorWind
```

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

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

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

Below is just a simple plotting function.

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

```
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Ã

```
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}$')
```

```
HTML(html)
```