Plotting NOAA ERSSTs in Python with Basemap

gistfile1.py

#! /usr/bin/env python

from mpl_toolkits.basemap import Basemap, shiftgrid, cm
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset

# The netCDF file is from ftp://ftp.cdc.noaa.gov/Datasets/noaa.ersst/sst.mnmean.nc
d = Dataset("sst.mon.anom.nc")
sst = d.variables["sst"][:][0, :, :]  # Should take the first time slice.
lon = d.variables["lon"][:]
lat = d.variables["lat"][:]
sst, lon = shiftgrid(185, sst, lon, start = False)
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
# For north Polar Stereographic projection
#m = Basemap(projection='npstere',boundinglat=10,lon_0=270,resolution='l')
# For classic Mercator... in case you hate geography.
m = Basemap(projection='merc', llcrnrlat=-80, urcrnrlat=80,
            llcrnrlon=-180, urcrnrlon=180, lat_ts=20, resolution='c')
# 555660 is because this is a 5.0 x 5.0 degree grid; 111132 meters are in about 111 km appart. This is fairly rough.
nx = int((m.xmax - m.xmin)/555660); ny = int((m.ymax - m.ymin)/555660)
sst = m.transform_scalar(sst, lon, lat, nx, ny)
im = m.imshow(sst, interpolation = "none")

m.drawcoastlines()
parallels = np.arange(-90, 90, 30)
meridians = np.arange(-180, 180, 60)
m.drawparallels(parallels, labels = [1, 0, 0, 1])
m.drawmeridians(meridians, labels = [1, 0, 0, 1])

cb = m.colorbar(im, "right", size = "5%", pad = "2%")
ax.set_title("Kaplan ERSST January 1895")
plt.show()