# make a sequence of png files which may then be displayed # as an animation using a tool like imagemagick animate, or # converted to an animate gif (using imagemagick convert). # reads data over http - needs an active internet connection. import numpy as np import matplotlib.pyplot as plt import matplotlib.mlab as mlab import numpy.ma as ma import datetime, sys, time, subprocess from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile, date2index, num2date # times for March 1993 'storm of the century' date1 = datetime.datetime(1993,3,10,0) date2 = datetime.datetime(1993,3,17,0) print date1, date2 # set OpenDAP server URL. URL="http://nomad1.ncep.noaa.gov:9090/dods/reanalyses/reanalysis-2/6hr/pgb/pgb" print URL try: data = NetCDFFile(URL) except: raise IOError, 'opendap server not providing the requested data' # read lats,lons,times. print data.variables.keys() latitudes = data.variables['lat'][:] longitudes = data.variables['lon'][:].tolist() times = data.variables['time'] ntime1 = date2index(date1,times,calendar='standard') ntime2 = date2index(date2,times,calendar='standard') print 'ntime1,ntime2:',ntime1,ntime2 print num2date(times[ntime1],times.units,calendar='standard'), num2date(times[ntime2],times.units,calendar='standard') # get sea level pressure and 10-m wind data. slpdata = data.variables['presmsl'] udata = data.variables['ugrdprs'] vdata = data.variables['vgrdprs'] # mult slp by 0.01 to put in units of millibars. slpin = 0.01*slpdata[ntime1:ntime2+1,:,:] uin = udata[ntime1:ntime2+1,0,:,:] vin = vdata[ntime1:ntime2+1,0,:,:] dates = num2date(times[ntime1:ntime2+1], times.units, calendar='standard') # add cyclic points slp = np.zeros((slpin.shape[0],slpin.shape[1],slpin.shape[2]+1),np.float64) slp[:,:,0:-1] = slpin; slp[:,:,-1] = slpin[:,:,0] u = np.zeros((uin.shape[0],uin.shape[1],uin.shape[2]+1),np.float64) u[:,:,0:-1] = uin; u[:,:,-1] = uin[:,:,0] v = np.zeros((vin.shape[0],vin.shape[1],vin.shape[2]+1),np.float64) v[:,:,0:-1] = vin; v[:,:,-1] = vin[:,:,0] longitudes.append(360.); longitudes = np.array(longitudes) # make 2-d grid of lons, lats lons, lats = np.meshgrid(longitudes,latitudes) print 'min/max slp,u,v' print slp.min(), slp.max() print uin.min(), uin.max() print vin.min(), vin.max() print 'dates' print dates # make orthographic basemaplt. m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.) plt.ion() # interactive mode on. uin = udata[ntime1:ntime2+1,0,:,:] vin = vdata[ntime1:ntime2+1,0,:,:] # make orthographic basemaplt. m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.) plt.ion() # interactive mode on. # create figure, add axes (leaving room for colorbar on right) fig = plt.figure() ax = fig.add_axes([0.1,0.1,0.7,0.7]) # set desired contour levels. clevs = np.arange(960,1061,5) # compute native x,y coordinates of grid. x, y = m(lons, lats) # define parallels and meridians to draw. parallels = np.arange(-80.,90,20.) meridians = np.arange(0.,360.,20.) # number of repeated frames at beginning and end is n1. nframe = 0; n1 = 10 pos = ax.get_position() l, b, w, h = pos.bounds # loop over times, make contour plots, draw coastlines, # parallels, meridians and title. for nt,date in enumerate(dates): CS = m.contour(x,y,slp[nt,:,:],clevs,linewidths=0.5,colors='k',animated=True) CS = m.contourf(x,y,slp[nt,:,:],clevs,cmap=plt.cm.RdBu_r,animated=True) # plot wind vectors on lat/lon grid. # rotate wind vectors to map projection coordinates. #urot,vrot = m.rotate_vector(u[nt,:,:],v[nt,:,:],lons,lats) # plot wind vectors over maplt. #Q = m.quiver(x,y,urot,vrot,scale=500) # plot wind vectors on projection grid (looks better). # first, shift grid so it goes from -180 to 180 (instead of 0 to 360 # in longitude). Otherwise, interpolation is messed uplt. ugrid,newlons = shiftgrid(180.,u[nt,:,:],longitudes,start=False) vgrid,newlons = shiftgrid(180.,v[nt,:,:],longitudes,start=False) # transform vectors to projection grid. urot,vrot,xx,yy = m.transform_vector(ugrid,vgrid,newlons,latitudes,51,51,returnxy=True,masked=True) # plot wind vectors over maplt. Q = m.quiver(xx,yy,urot,vrot,scale=500) # make quiver key. qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W') # draw coastlines, parallels, meridians, title. m.drawcoastlines(linewidth=1.5) m.drawparallels(parallels) m.drawmeridians(meridians) plt.title('SLP and Wind Vectors '+str(date)) if nt == 0: # plot colorbar on a separate axes (only for first frame) cax = plt.axes([l+w-0.05, b, 0.03, h]) # setup colorbar axes fig.colorbar(CS,drawedges=True, cax=cax) # draw colorbar cax.text(0.0,-0.05,'mb') plt.axes(ax) # reset current axes plt.draw() # draw the plot # save first and last frame n1 times # (so gif animation pauses at beginning and end) if nframe == 0 or nt == slp.shape[0]-1: for n in range(n1): plt.savefig('anim%03i'%nframe+'.png') nframe = nframe + 1 else: plt.savefig('anim%03i'%nframe+'.png') nframe = nframe + 1 ax.clear() # clear the axes for the next plot. print """ Now display animation using imagemagick 'animate -delay 10 anim*png' or, create an animated gif using 'convert -delay 10 anim*png anim.gif' and display in a web browser (or Powerpoint). Will try to run imagemagick animate in 5 seconds ...""" time.sleep(5) subprocess.call('animate -delay 10 anim*png',shell=True)