%load_ext load_style
%load_style talk.css
%matplotlib inline
windspharm is a Python library developed by Andrew Dawson which provides an pythonic interface to the pyspharm module, which is basically a bindings to the [spherepack] Fortran library
Installation
Download and unpack pyspharm
Download spherepack from http://www2.cisl.ucar.edu/resources/legacy/spherepack and unpack
copy all the Fortran files in [path_to_spherepack]/spherepack3.2/src
to [path_to_pyspharm]/pyspharm-1.0.8/src
install pyspharm:
pyspharm-1.0.8 ᐅ python setup.py build
pyspharm-1.0.8 ᐅ python setup.py install
ᐅ pip install windspharm
windspharm has 3 different interfaces
:
standard
: expects numpy arrays as inputscdms
: cdms2 objects (cdms2 is a class for opening netcdf files [amongst other formats]) part of the cdat-lite package or UV-CDAT distribution)iris
: iris cubesWe are going to use xray here, and thus use the standard interface, passing the underlying numpy arrays
from windspharm.standard import VectorWind
from windspharm.tools import prep_data, recover_data, order_latdim
import os, sys
import pandas as pd
import numpy as np
from numpy import ma
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap as bm
dpath = os.path.join(os.environ.get('HOME'), 'data/NCEP1')
def plot_field(X, lat, lon, vmin, vmax, step, cmap=plt.get_cmap('jet'), ax=False, title=False, grid=False):
if not ax:
f, ax = plt.subplots(figsize=(10, (X.shape[0] / float(X.shape[1])) * 10))
m.ax = ax
im = m.contourf(lons, lats, X, np.arange(vmin, vmax+step, step), latlon=True, cmap=cmap, extend='both', ax=ax)
m.drawcoastlines()
if grid:
m.drawmeridians(np.arange(0, 360, 60), labels=[0,0,0,1])
m.drawparallels([-40, 0, 40], labels=[1,0,0,0])
m.colorbar(im)
if title:
ax.set_title(title)
import xray; print(xray.__version__)
0.4.0rc1
dset_u = xray.open_dataset(os.path.join(dpath,'uwnd.2014.nc'))
dset_v = xray.open_dataset(os.path.join(dpath,'vwnd.2014.nc'))
dset_u = dset_u.sel(level=200)
dset_v = dset_v.sel(level=200)
dset_u = dset_u.mean('time')
dset_v = dset_v.mean('time')
lats = dset_u['lat'].values
lons = dset_u['lon'].values
uwnd = dset_u['uwnd'].values
vwnd = dset_v['vwnd'].values
uwnd, uwnd_info = prep_data(uwnd, 'yx')
vwnd, vwnd_info = prep_data(vwnd, 'yx')
# It is also required that the latitude dimension is north-to-south. Again the
# bundled tools make this easy.
lats, uwnd, vwnd = order_latdim(lats, uwnd, vwnd)
lons, lats = np.meshgrid(lons, lats)
w = VectorWind(uwnd, vwnd)
sf, vp = w.sfvp()
vp = vp * 10e-6
sf = sf * 10e-6
m = bm(projection='cyl',llcrnrlat=-90,urcrnrlat=90,\
llcrnrlon=0,urcrnrlon=360,\
lat_ts=0,resolution='c')
plot_field(sf.squeeze(), lats, lons, -1500, 1500, 100, cmap=plt.get_cmap('RdBu_r'), \
title="Streamfunction at 200 hPa ($10^6$m$^2$s$^{-1}$)", grid=True)
plot_field(vp.squeeze(), lats, lons, -100, 100, 10, cmap=plt.get_cmap('RdBu_r'), \
title="Velocity Potential at 200 hPa ($10^6$m$^2$s$^{-1}$)", grid=True)