1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
from __future__ import (absolute_import, division, print_function)
# example using matplotlib.animation to create a movie
# reads data over http - needs an active internet connection.
import numpy as np
import matplotlib.pyplot as plt
import numpy.ma as ma
import datetime, time
from mpl_toolkits.basemap import Basemap, shiftgrid
from netCDF4 import Dataset as NetCDFFile, date2index, num2date
import matplotlib.animation as animation
# times for March 1993 'storm of the century'
date1 = datetime.datetime(1993,3,10,0)
date2 = datetime.datetime(1993,3,17,0)
# set OpenDAP server URL.
URL="http://nomad2.ncep.noaa.gov:9090/dods/reanalyses/reanalysis-2/6hr/pgb/pgb"
try:
data = NetCDFFile(URL)
except:
raise IOError('opendap server not providing the requested data')
# read lats,lons,times.
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')
# 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)
# make orthographic basemap.
m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.)
uin = udata[ntime1:ntime2+1,0,:,:]
vin = vdata[ntime1:ntime2+1,0,:,:]
# 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.
nt = 0; date = dates[nt]
CS1 = m.contour(x,y,slp[nt,:,:],clevs,linewidths=0.5,colors='k')
CS2 = m.contourf(x,y,slp[nt,:,:],clevs,cmap=plt.cm.RdBu_r)
# 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 map.
#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 up.
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 map.
Q = m.quiver(xx,yy,urot,vrot,scale=500,zorder=10)
# 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)
txt = plt.title('SLP and Wind Vectors '+str(date))
# 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(CS2,drawedges=True, cax=cax) # draw colorbar
cax.text(0.0,-0.05,'mb')
plt.axes(ax) # reset current axes
def updatefig(nt):
global CS1,CS2,Q
date = dates[nt]
for c in CS1.collections: c.remove()
CS1 = m.contour(x,y,slp[nt,:,:],clevs,linewidths=0.5,colors='k')
for c in CS2.collections: c.remove()
CS2 = m.contourf(x,y,slp[nt,:,:],clevs,cmap=plt.cm.RdBu_r)
ugrid,newlons = shiftgrid(180.,u[nt,:,:],longitudes,start=False)
vgrid,newlons = shiftgrid(180.,v[nt,:,:],longitudes,start=False)
urot,vrot,xx,yy = m.transform_vector(ugrid,vgrid,newlons,latitudes,51,51,returnxy=True,masked=True)
txt.set_text('SLP and Wind Vectors '+str(date))
Q.set_UVC(urot,vrot)
ani = animation.FuncAnimation(fig, updatefig, frames=len(dates))
#ani.save('movie.mp4')
plt.show()
|
