Using Iris to access NCEP CFSR 30-year Wave Hindcast

This demonstrates extracting data from a large aggregated archive via OPeNDAP, taking advantage of CF conventions to enable a plot without specification of lon,lat variables, and save extracted data to GRIB2



In [20]:

    
from IPython.core.display import HTML
HTML('<iframe src=http://scitools.org.uk/iris/ width=800 height=350></iframe>')









    Out[20]:



In [21]:

    
import numpy
import matplotlib.pyplot as plt
import datetime as dt

import iris
import iris.plot as iplt
import iris.quickplot as qplt
import cartopy.crs as ccrs



In [22]:

    
def time_near(cube,start):
    timevar=cube.coord('time')
    itime = timevar.nearest_neighbour_index(timevar.units.date2num(start))
    return timevar.points[itime]



In [23]:

    
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));

Extract some data using OPeNDAP



In [24]:

    
# 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 [25]:

    
print cubes









    



0: Significant height of combined wind waves and swell @ Ground or water surface / (m) (time: 90584; latitude: 481; longitude: 586)
1: u-component of wind @ Ground or water surface / (m/s) (time: 90096; latitude: 481; longitude: 586)
2: v-component of wind @ Ground or water surface / (m/s) (time: 90096; latitude: 481; longitude: 586)
3: Primary wave direction (degree true) @ Ground or water surface / (unknown) (time: 90584; latitude: 481; longitude: 586)
4: Primary wave mean period @ Ground or water surface / (s) (time: 90584; latitude: 481; longitude: 586)



In [26]:

    
hsig=cubes[0]



In [27]:

    
print hsig









    



Significant height of combined wind waves and swell @ Ground or water surface / (m) (time: 90584; latitude: 481; longitude: 586)
     Dimension coordinates:
          time                                                                           x                -               -
          latitude                                                                       -                x               -
          longitude                                                                      -                -               x
     Attributes:
          Analysis_or_forecast_generating_process_identifier_defined_by_originating_centre: Missing
          Conventions: CF-1.6
          GRIB_table_version: 2,1
          Grib2_Generating_Process_Type: Forecast
          Grib2_Level_Type: 1
          Grib2_Parameter: [10  0  3]
          Grib2_Parameter_Category: Waves
          Grib2_Parameter_Discipline: Oceanographic products
          Grib2_Parameter_Name: Significant height of combined wind waves and swell
          Grib_Variable_Id: VAR_10-0-3_L1
          Originating_or_generating_Center: US National Weather Service, National Centres for Environmental Prediction...
          Originating_or_generating_Subcenter: 0
          Type_of_generating_process: Forecast
          _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
          abbreviation: HTSGW
          featureType: GRID
          file_format: NETCDF3_CLASSIC
          history: Read using CDM IOSP Grib2Collection



In [28]:

    
# 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
slice=hsig.extract(iris.Constraint(time=time_near(hsig,mytime),
    longitude=lambda cell: -71.5 < cell < -64.0,
    latitude=lambda cell: 40.0 < cell < 46.0))



In [29]:

    
print slice









    



Significant height of combined wind waves and swell @ Ground or water surface / (m) (latitude: 89; longitude: 112)
     Dimension coordinates:
          latitude                                                                           x              -
          longitude                                                                          -              x
     Scalar coordinates:
          time: 1991-10-31 12:00:00
     Attributes:
          Analysis_or_forecast_generating_process_identifier_defined_by_originating_centre: Missing
          Conventions: CF-1.6
          GRIB_table_version: 2,1
          Grib2_Generating_Process_Type: Forecast
          Grib2_Level_Type: 1
          Grib2_Parameter: [10  0  3]
          Grib2_Parameter_Category: Waves
          Grib2_Parameter_Discipline: Oceanographic products
          Grib2_Parameter_Name: Significant height of combined wind waves and swell
          Grib_Variable_Id: VAR_10-0-3_L1
          Originating_or_generating_Center: US National Weather Service, National Centres for Environmental Prediction...
          Originating_or_generating_Subcenter: 0
          Type_of_generating_process: Forecast
          _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention
          abbreviation: HTSGW
          featureType: GRID
          file_format: NETCDF3_CLASSIC
          history: Read using CDM IOSP Grib2Collection



In [30]:

    
print[coord.name() for coord in slice.coords()]









    



['latitude', 'longitude', u'time']

Plot using Iris Quickplot



In [31]:

    
qplt.contourf(slice)









    Out[31]:





<matplotlib.contour.QuadContourSet instance at 0x537acb0>

Plot using Iris.plot

There is also Iris.quickplot, but I wanted to add my own title here and control the orientation of the colorbar, so I used the more flexible Iris.plot



In [32]:

    
figure(figsize=(12,8))

# set the projection
ax1 = plt.axes(projection=ccrs.Mercator())

# color filled contour plot
h = iplt.contourf(slice,64)

# add coastlines, colorbar and title
plt.gca().coastlines(resolution='10m')
colorbar(h,orientation='vertical');
title('WaveWatch 3 significant wave height driven by CFSR winds');

Save extracted Cube to GRIB2 file



In [33]:

    
# 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 [34]:

    
# add a forecast_period (required by GRIB2 writer)
slice.add_aux_coord(iris.coords.DimCoord(0, standard_name='forecast_period', units='hours'))



In [35]:

    
# add a vertical coordinate (required by GRIB2 writer)
slice.add_aux_coord(iris.coords.DimCoord(0, "height", units="m"))



In [36]:

    
# GRIB2 stores data in long values.  Do we have huge values that would create problems?
slice.data.max()









    Out[36]:





nan



In [37]:

    
# ah, we have NaNs.  This causes problems for GRIB2 writer
# So convert to masked-array instead.
slice.data=ma.masked_invalid(slice.data)



In [38]:

    
# Finally... save slice as grib2
iris.save(slice,'hsig.grib2')