An example of plotting the zonal average T anomaly w.r.t. WOA2005 climatology

plot_zonal_average_T_anomaly.py

import netCDF4
import numpy
import matplotlib
import matplotlib.pyplot as plt

[ny,nx] = netCDF4.Dataset('/archive/gold/datasets/OM4_025/mosaic.v20140610.unpacked/ocean_topog.nc').variables['depth'].shape
area = netCDF4.Dataset('/archive/gold/datasets/OM4_025/mosaic.v20140610.unpacked/ocean_hgrid.nc').variables['area'][:].reshape(ny,2,nx,2).sum(axis=-1).sum(axis=1)
lat = netCDF4.Dataset('/archive/gold/datasets/OM4_025/mosaic.v20140610.unpacked/ocean_hgrid.nc').variables['y'][1::2,1::2]

annualFile = '/archive/bls/tikal_201406_mom6_2014.07.01/OM4_SIS_baseline/gfdl.ncrc2-intel-prod/pp/ocean_annual/av/annual_5yr/ocean_annual.1900-1904.ann.nc'
zInterface = netCDF4.Dataset(annualFile).variables['e'][0,:]
T1 = netCDF4.Dataset(annualFile).variables['temp'][0,:]
To = netCDF4.Dataset('/archive/gold/datasets/OM4_025/obs/WOA05_ptemp_salt_annual.nc').variables['temp'][:]

hThickness = zInterface[:-1] - zInterface[1:]
volumes = area * hThickness
delTxave = numpy.sum( volumes * (T1 - To), axis=-1 ) / numpy.sum( volumes, axis=-1)

y = numpy.amax( lat, axis=-1) # Note: zonal averaging is inappropriate on tri-polar grid
z = numpy.amin( zInterface, axis=-1 )

plt.gca().set_axis_bgcolor([.5,.5,.5]) # Makes topography appear as grey
cmap = plt.get_cmap('seismic')
cmap.set_over([0.4,0.1,0.1])
cmap.set_under([0.1,0.1,0.4])
ci = numpy.array([-4, -3, -2, -1, -.5, -.2, .2, .5, 1, 2, 3, 4])
norm = matplotlib.colors.BoundaryNorm( ci, cmap.N)
ch = plt.pcolormesh( y, z, delTxave, cmap=cmap, norm=norm ); cbax=plt.colorbar(ch, extend='both')
cbax.set_ticks(ci) # Force CI labels on all contour intervals
plt.xlim(-80, 90); plt.ylim(-6500, 0)
plt.xlabel(r'Latitude [$\degree$N]')
plt.ylabel('z* [m]')
plt.title(r'Zonal average $\theta$ anomaly [$\degree$C]')
plt.show()