Animated SSH on blue marble

demo_animate_ssh_blue_marble.py

import netCDF4
import matplotlib.pyplot as plt
import numpy
from mpl_toolkits.basemap import Basemap

lon = netCDF4.Dataset('/archive/gold/datasets/MOM6z_SIS_025/siena/mosaic.unpacked/ocean_hgrid.nc').variables['x'][::2,::2]
lat = netCDF4.Dataset('/archive/gold/datasets/MOM6z_SIS_025/siena/mosaic.unpacked/ocean_hgrid.nc').variables['y'][::2,::2]

dailyFiles = '/archive/aja/awg/tikal_201403/CM4_c192L48_am4a1_2000/gfdl.ncrc2-intel-prod-openmp/history/unpack/000[3456]0101.ocean_daily.nc'
ssh = netCDF4.MFDataset(dailyFiles).variables['ssh']
dayNum = netCDF4.MFDataset(dailyFiles).variables['time']

def JulianDate(n):
  """Returns a (year/month/day) date for the given Julian day number."""
  a = int(n) + 32044
  b = (4*a + 3)//146097
  c = a - (146097*b)//4
  d = (4*c + 3)//1461
  e = c - (1461*d)//4
  m = (5*e + 2)//153
  day = e + 1 - (153*m + 2)//5
  month = m + 3 - 12*(m//10)
  year = 100*b + d - 4800 + m/10
  return year, month, day

def makePlot(fig, ssh, n, year, month, day):
  """This function plots ssh[n] to fig."""
  plt.clf()
  ax = fig.add_axes([0.03,0.03,0.94,0.94])
  m = Basemap(projection='kav7',lon_0=-120,resolution=None)
  m.drawmapboundary(fill_color='0.3')
  im0 = m.bluemarble(scale=0.25)
  im1 = m.pcolormesh(numpy.minimum(lon,60.),lat,ssh[n],shading='flat',cmap=plt.cm.spectral,latlon=True)
  plt.clim(-2,1.25)
  m.drawparallels(numpy.arange(-90.,99.,30.))
  m.drawmeridians(numpy.arange(-180.,180.,60.))
  cbaxes = fig.add_axes([0.01, 0.91, 0.16, 0.03]) 
  cb = plt.colorbar(im1, cax = cbaxes, orientation='horizontal', extend='both', ticks=[-2.,-1.,0.,1.])
  cb.set_label('[m]',fontsize=10)
  cb.ax.tick_params(labelsize=10)
  plt.figtext(0.09,.975,'Sea surface height',verticalalignment='top',size='medium',horizontalalignment='center')
  plt.figtext(0.81,.975,'Model date = %4.4i.%2.2i.%2.2i'%(year,month,day),verticalalignment='top')
  plt.draw()

fig = plt.figure(figsize=(12.8,7.2),dpi=100) # 720p
#plt.show(block=False)
for n in range(ssh.shape[0]):
  year, month, day = JulianDate(1721426 + dayNum[n])
  print 'record %i (y/m/d %i/%i/%i)'%(n,year,month,day)
  makePlot(fig, ssh, n, year, month, day)
  plt.savefig('figs/ssh.%4.4i.%2.2i.%2.2i.png'%(year,month,day))