|
| 1 | +import numpy as np |
| 2 | +import matplotlib.pyplot as plt |
| 3 | +import datetime |
| 4 | +from mpl_toolkits.basemap import Basemap, shiftgrid |
| 5 | +from netCDF4 import Dataset, date2index, num2date |
| 6 | +from mpl_toolkits.axes_grid1 import make_axes_locatable |
| 7 | +# specify date to plot. |
| 8 | +date = datetime.datetime(1993,3,14,0) |
| 9 | +# set OpenDAP server URL. |
| 10 | +URL="http://nomad1.ncep.noaa.gov:9090/dods/reanalyses/reanalysis-2/6hr/pgb/pgb" |
| 11 | +data = Dataset(URL) |
| 12 | +# read lats,lons,times. |
| 13 | +latitudes = data.variables['lat'][:] |
| 14 | +longitudes = data.variables['lon'][:].tolist() |
| 15 | +times = data.variables['time'] |
| 16 | +# get time index for desired date. |
| 17 | +ntime = date2index(date,times,calendar='standard') |
| 18 | +# get sea level pressure and 10-m wind data. |
| 19 | +slpdata = data.variables['presmsl'] |
| 20 | +udata = data.variables['ugrdprs'] |
| 21 | +vdata = data.variables['vgrdprs'] |
| 22 | +# mult slp by 0.01 to put in units of hPa. |
| 23 | +slpin = 0.01*slpdata[ntime,:,:] |
| 24 | +uin = udata[ntime,0,:,:] |
| 25 | +vin = vdata[ntime,0,:,:] |
| 26 | +# add cyclic points manually (could use addcyclic function) |
| 27 | +slp = np.zeros((slpin.shape[0],slpin.shape[1]+1),np.float) |
| 28 | +slp[:,0:-1] = slpin; slp[:,-1] = slpin[:,0] |
| 29 | +u = np.zeros((uin.shape[0],uin.shape[1]+1),np.float64) |
| 30 | +u[:,0:-1] = uin; u[:,-1] = uin[:,0] |
| 31 | +v = np.zeros((vin.shape[0],vin.shape[1]+1),np.float64) |
| 32 | +v[:,0:-1] = vin; v[:,-1] = vin[:,0] |
| 33 | +longitudes.append(360.); longitudes = np.array(longitudes) |
| 34 | +# make 2-d grid of lons, lats |
| 35 | +lons, lats = np.meshgrid(longitudes,latitudes) |
| 36 | +# make orthographic basemap. |
| 37 | +m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.) |
| 38 | +# create figure, add axes |
| 39 | +fig1 = plt.figure(figsize=(8,10)) |
| 40 | +ax = fig1.add_axes([0.1,0.1,0.8,0.8]) |
| 41 | +# set desired contour levels. |
| 42 | +clevs = np.arange(960,1061,5) |
| 43 | +# compute native x,y coordinates of grid. |
| 44 | +x, y = m(lons, lats) |
| 45 | +# define parallels and meridians to draw. |
| 46 | +parallels = np.arange(-80.,90,20.) |
| 47 | +meridians = np.arange(0.,360.,20.) |
| 48 | +# plot SLP contours. |
| 49 | +CS1 = m.contour(x,y,slp,clevs,linewidths=0.5,colors='k',animated=True) |
| 50 | +CS2 = m.contourf(x,y,slp,clevs,cmap=plt.cm.RdBu_r,animated=True) |
| 51 | +# plot wind vectors on projection grid. |
| 52 | +# first, shift grid so it goes from -180 to 180 (instead of 0 to 360 |
| 53 | +# in longitude). Otherwise, interpolation is messed up. |
| 54 | +ugrid,newlons = shiftgrid(180.,u,longitudes,start=False) |
| 55 | +vgrid,newlons = shiftgrid(180.,v,longitudes,start=False) |
| 56 | +# transform vectors to projection grid. |
| 57 | +uproj,vproj,xx,yy = \ |
| 58 | +m.transform_vector(ugrid,vgrid,newlons,latitudes,31,31,returnxy=True,masked=True) |
| 59 | +# now plot. |
| 60 | +Q = m.quiver(xx,yy,uproj,vproj,scale=700) |
| 61 | +# make quiver key. |
| 62 | +qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W') |
| 63 | +# draw coastlines, parallels, meridians. |
| 64 | +m.drawcoastlines(linewidth=1.5) |
| 65 | +m.drawparallels(parallels) |
| 66 | +m.drawmeridians(meridians) |
| 67 | +# use axes_grid toolkit to make colorbar axes. |
| 68 | +divider = make_axes_locatable(ax) |
| 69 | +cax = divider.append_axes("bottom", size="5%", pad=0.1) |
| 70 | +cb = plt.colorbar(CS2,orientation='horizontal',cax=cax) |
| 71 | +cb.set_label('hPa') |
| 72 | +# set plot title |
| 73 | +ax.set_title('SLP and Wind Vectors '+str(date)) |
| 74 | +fig1.savefig('plotwindvec1.png') |
| 75 | + |
| 76 | +# create 2nd figure, add axes |
| 77 | +fig2 = plt.figure(figsize=(8,10)) |
| 78 | +ax = fig2.add_axes([0.1,0.1,0.8,0.8]) |
| 79 | +# plot SLP contours |
| 80 | +CS1 = m.contour(x,y,slp,clevs,linewidths=0.5,colors='k',animated=True) |
| 81 | +CS2 = m.contourf(x,y,slp,clevs,cmap=plt.cm.RdBu_r,animated=True) |
| 82 | +# plot wind barbs over map. |
| 83 | +barbs = m.barbs(xx,yy,uproj,vproj,length=5,barbcolor='k',flagcolor='r',linewidth=0.5) |
| 84 | +# draw coastlines, parallels, meridians. |
| 85 | +m.drawcoastlines(linewidth=1.5) |
| 86 | +m.drawparallels(parallels) |
| 87 | +m.drawmeridians(meridians) |
| 88 | +# use axes_grid toolkit to make colorbar axes. |
| 89 | +divider = make_axes_locatable(ax) |
| 90 | +cax = divider.append_axes("bottom", size="5%", pad=0.1) |
| 91 | +cb = plt.colorbar(CS2,orientation='horizontal',cax=cax) |
| 92 | +cb.set_label('hPa') |
| 93 | +# set plot title. |
| 94 | +ax.set_title('SLP and Wind Barbs '+str(date)) |
| 95 | +fig2.savefig('plotwindvec2.png') |
0 commit comments