Plotting NCEP/NCAR geopotential height 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
import datetime

datetime_origin = datetime.datetime(1, 1, 1, 0, 0, 0, 0)
# In this case, this is 500 mb. Can check below with d.variables["level"][:]
level = 5
time = 0
# 277830 is because this is a ~2.5 x 2.5 degree grid; 111132 meters are in about 111 km appart. This is fairly rough.
meters_per_grid = 277830

# Download the dataset from: ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.mean.nc
d = Dataset("hgt.mon.mean.nc")
hgt = d.variables["hgt"][:][time, level, :, :]  # Not sure about the dims, here.
lon = d.variables["lon"][:]
# need to reverse direction of lat dimension so it's increasing.
lat = d.variables["lat"][:][::-1]
hgt = hgt[::-1, :]
hgt, lon = shiftgrid(180, hgt, lon, start = False)
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
# Lambert Azimuthal Equal Area
# Larger area:
#m = Basemap(width = 12000000, height = 12000000, resolution = "l", projection = "laea", lat_ts = 50, lat_0 = 50, lon_0 = -107.0)
# Smaller area:
#m = Basemap(width = 12000000, height = 8000000, resolution = "l", projection = "laea", lat_ts = 50, lat_0 = 50, lon_0 = -107.0)
# For north Polar Stereographic projection
m = Basemap(projection='npstere',boundinglat=10,lon_0=270,resolution='l')
nx = int((m.xmax - m.xmin)/meters_per_grid); ny = int((m.ymax - m.ymin)/meters_per_grid)
hgt = m.transform_scalar(hgt, lon, lat, nx, ny)
im = m.imshow(hgt, 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("500 mb Geopotential Height")
plt.show()