In [2]:
from IPython.core.display import HTML
HTML('<iframe src=http://scitools.org.uk/iris/ width=800 height=350></iframe>')
Out[2]:
In [3]:
import numpy
import matplotlib.pyplot as plt
import datetime as dt
import iris
import iris.plot as iplt
import cartopy.crs as ccrs
In [4]:
def time_near(cube,start):
timevar=cube.coord('time')
itime = timevar.nearest_neighbour_index(timevar.units.date2num(start))
return timevar.points[itime]
In [5]:
def myplot(slice,model=None):
# make the plot
figure(figsize=(12,8))
lat=slice.coord(axis='Y').points
lon=slice.coord(axis='X').points
time=slice.coord('time')[0]
subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
pcolormesh(lon,lat,masked_invalid(slice.data));
colorbar()
grid()
date_str=time.units.num2date(time.points)
plt.title('%s: %s: %s' % (model,slice.long_name,date_str));
In [6]:
# DAP URL: 30 year East Coast wave hindcast (Wave Watch 3 driven by CFSR Winds)
#cubes = iris.load('http://geoport.whoi.edu/thredds/dodsC/fmrc/NCEP/ww3/cfsr/4m/best'); # 4 arc minute resolution
cubes = iris.load('http://geoport.whoi.edu/thredds/dodsC/fmrc/NCEP/ww3/cfsr/10m/best'); # 10 arc minute resolution
In [7]:
print cubes
In [8]:
# use contraints to select geographic subset and nearest time
mytime=dt.datetime(1991,10,31,12) #specified time... Nov 1, 1991 was the "Perfect Storm"
#mytime=dt.datetime.utcnow() # .... or now
In [25]:
cubes2 = cubes.extract(iris.Constraint(time=time_near(cubes[0],mytime),
longitude=lambda cell: -71.5 < cell < -64.0,
latitude=lambda cell: 40.0 < cell < 46.0))
In [26]:
hsig=cubes2[0]
u=cubes2[1]
v=cubes2[2]
In [27]:
print cubes2
In [29]:
print[coord.name() for coord in hsig.coords()]
In [32]:
def myplot(slice,u=None,v=None,model=None):
# make the plot
figure(figsize=(12,8))
lat=slice.coord(axis='Y').points
lon=slice.coord(axis='X').points
time=slice.coord('time')[0]
subplot(111,aspect=(1.0/cos(mean(lat)*pi/180.0)))
pcolormesh(lon,lat,slice.data);
isub = 2
[lon2d,lat2d]=meshgrid(lon,lat)
Q = quiver(lon[::isub,::isub],lat[::isub,::isub],u[::isub,::isub],v[::isub,::isub],scale=20)
maxstr='%3.1f m/s' % maxvel
qk = quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
colorbar()
grid()
date=time.units.num2date(time.points)
date_str=date[0].strftime('%Y-%m-%d %H:%M:%S %Z')
plt.title('%s: %s: %s' % (model,slice.long_name,date_str));
In [43]:
lat=hsig.coord(axis='Y').points
lon=hsig.coord(axis='X').points
hsig = ma.masked_invalid(hsig.data)
u = ma.masked_invalid(u.data)
v = ma.masked_invalid(v.data)
In [45]:
[lon2d,lat2d]=meshgrid(lon,lat)
isub=2
In [46]:
Q = quiver(lon2d[::isub,::isub],lat2d[::isub,::isub],u[::isub,::isub],v[::isub,::isub],scale=20)
maxstr='%3.1f m/s' % maxvel
qk = quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
In [38]:
In [39]:
Q = quiver(lon2d[::isub,::isub],lat2d[::isub,::isub],u[::isub,::isub],v[::isub,::isub],scale=20)
maxstr='%3.1f m/s' % maxvel
qk = quiverkey(Q,0.92,0.08,maxvel,maxstr,labelpos='W')
In [44]:
myplot(cubes[0],u=cubes[1].data,v=cubes[2].data,model='foo')
In [12]:
# first add a Geographic Coord system (required by GRIB2 writer)
# here we add a spherical earth with specified radius
slice.coord(axis='X').coord_system=iris.coord_systems.GeogCS(654321)
slice.coord(axis='Y').coord_system=iris.coord_systems.GeogCS(654321)
In [13]:
# add a forecast_period (required by GRIB2 writer)
slice.add_aux_coord(iris.coords.DimCoord(0, standard_name='forecast_period', units='hours'))
In [14]:
# add a vertical coordinate (required by GRIB2 writer)
slice.add_aux_coord(iris.coords.DimCoord(0, "height", units="m"))
In [15]:
# GRIB2 stores data in long values. Do we have huge values that would create problems?
slice.data.max()
Out[15]:
In [16]:
# ah, we have NaNs. This causes problems for GRIB2 writer
# So convert to masked-array instead.
slice.data=ma.masked_invalid(slice.data)
In [17]:
# Finally... save slice as grib2
iris.save(slice,'hsig.grib2')