Matplotlib/Basemap and pygrib example¶
netCDF is easier to deal with, but most operational forecast centers provide data in GRIB format. GRIB is a record format, where every record is a 2D field. In this example we use pygrib to read some ECMWF ensemble forecast data, then use matplotlib and Basemap plot forecast maps. Pygrib uses the ECMWF GRIB_API C library under the hood.¶
In [3]:
# Make the output of plotting commands be displayed inline within the notebook,
%matplotlib inline
from mpl_toolkits.basemap import Basemap # import Basemap matplotlib toolkit
import numpy as np
import matplotlib.pyplot as plt
import pygrib # import pygrib interface to grib_api
Open the grib file (downloaded from the NCAR TIGGE data portal, available here). The pygrib.open object behaves much like a python file object.¶
In [4]:
grbs = pygrib.open('ecmwf_ensemble.grb')
The pygrib.open object is a python iterator.¶
In [5]:
for grb in grbs[:4]:
print grb
1:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000 2:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000:hi res cntl fcst 3:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000:pos ens pert 1 4:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000:pos ens pert 2
The iterator yields a grib message object. Each grib message object has a set of attributes, or 'keys', which can be used for searching.¶
In [6]:
print grb.keys()
[u'parametersVersion', u'hundred', u'globalDomain', u'GRIBEditionNumber', u'grib2divider', u'missingValue', u'ieeeFloats', u'section0Length', u'identifier', u'discipline', u'editionNumber', u'totalLength', u'sectionNumber', u'section1Length', u'numberOfSection', u'centre', u'centreDescription', u'subCentre', u'tablesVersion', u'masterDir', u'localTablesVersion', u'significanceOfReferenceTime', u'year', u'month', u'day', u'hour', u'minute', u'second', u'dataDate', u'julianDay', u'dataTime', u'productionStatusOfProcessedData', u'typeOfProcessedData', u'selectStepTemplateInterval', u'selectStepTemplateInstant', u'stepType', u'sectionNumber', u'grib2LocalSectionPresent', u'sectionNumber', u'gridDescriptionSectionPresent', u'section3Length', u'numberOfSection', u'sourceOfGridDefinition', u'numberOfDataPoints', u'numberOfOctectsForNumberOfPoints', u'interpretationOfNumberOfPoints', u'PLPresent', u'gridDefinitionTemplateNumber', u'shapeOfTheEarth', u'scaleFactorOfRadiusOfSphericalEarth', u'scaledValueOfRadiusOfSphericalEarth', u'scaleFactorOfEarthMajorAxis', u'scaledValueOfEarthMajorAxis', u'scaleFactorOfEarthMinorAxis', u'scaledValueOfEarthMinorAxis', u'radius', u'Ni', u'Nj', u'basicAngleOfTheInitialProductionDomain', u'mBasicAngle', u'angleMultiplier', u'mAngleMultiplier', u'subdivisionsOfBasicAngle', u'angleDivisor', u'latitudeOfFirstGridPoint', u'longitudeOfFirstGridPoint', u'resolutionAndComponentFlags', u'resolutionAndComponentFlags1', u'resolutionAndComponentFlags2', u'iDirectionIncrementGiven', u'jDirectionIncrementGiven', u'uvRelativeToGrid', u'resolutionAndComponentFlags6', u'resolutionAndComponentFlags7', u'resolutionAndComponentFlags8', u'ijDirectionIncrementGiven', u'latitudeOfLastGridPoint', u'longitudeOfLastGridPoint', u'iDirectionIncrement', u'jDirectionIncrement', u'scanningMode', u'iScansNegatively', u'jScansPositively', u'jPointsAreConsecutive', u'alternativeRowScanning', u'iScansPositively', u'scanningMode5', u'scanningMode6', u'scanningMode7', u'scanningMode8', u'g2grid', u'latitudeOfFirstGridPointInDegrees', u'longitudeOfFirstGridPointInDegrees', u'latitudeOfLastGridPointInDegrees', u'longitudeOfLastGridPointInDegrees', u'iDirectionIncrementInDegrees', u'jDirectionIncrementInDegrees', u'latLonValues', u'latitudes', u'longitudes', u'distinctLatitudes', u'distinctLongitudes', u'gridType', u'sectionNumber', u'section4Length', u'numberOfSection', u'NV', u'neitherPresent', u'productDefinitionTemplateNumber', u'parameterCategory', u'parameterNumber', u'parameterUnits', u'parameterName', u'typeOfGeneratingProcess', u'backgroundProcess', u'generatingProcessIdentifier', u'hoursAfterDataCutoff', u'minutesAfterDataCutoff', u'indicatorOfUnitOfTimeRange', u'stepUnits', u'forecastTime', u'startStep', u'endStep', u'stepRange', u'stepTypeInternal', u'validityDate', u'validityTime', u'typeOfFirstFixedSurface', u'unitsOfFirstFixedSurface', u'nameOfFirstFixedSurface', u'scaleFactorOfFirstFixedSurface', u'scaledValueOfFirstFixedSurface', u'typeOfSecondFixedSurface', u'unitsOfSecondFixedSurface', u'nameOfSecondFixedSurface', u'scaleFactorOfSecondFixedSurface', u'scaledValueOfSecondFixedSurface', u'pressureUnits', u'typeOfLevel', u'level', u'bottomLevel', u'topLevel', u'typeOfEnsembleForecast', u'perturbationNumber', u'numberOfForecastsInEnsemble', u'paramIdECMF', u'paramId', u'shortNameECMF', u'shortName', u'unitsECMF', u'units', u'nameECMF', u'name', u'cfNameECMF', u'cfName', u'cfVarNameECMF', u'cfVarName', u'marsExpver', u'marsClass', u'marsModel', u'marsType', u'marsStream', u'ifsParam', u'genVertHeightCoords', u'PVPresent', u'sectionNumber', u'section5Length', u'numberOfSection', u'numberOfValues', u'dataRepresentationTemplateNumber', u'packingType', u'referenceValue', u'referenceValueError', u'binaryScaleFactor', u'decimalScaleFactor', u'bitsPerValue', u'typeOfOriginalFieldValues', u'typeOfCompressionUsed', u'targetCompressionRatio', u'lengthOfHeaders', u'sectionNumber', u'section6Length', u'numberOfSection', u'bitMapIndicator', u'bitmapPresent', u'sectionNumber', u'section7Length', u'numberOfSection', u'codedValues', u'values', u'maximum', u'minimum', u'average', u'numberOfMissing', u'standardDeviation', u'skewness', u'kurtosis', u'isConstant', u'changeDecimalPrecision', u'decimalPrecision', u'setBitsPerValue', u'getNumberOfValues', u'scaleValuesBy', u'offsetValuesBy', u'productType', u'section8Length']
By looping over the iterator and checking the attributes, you can select the grib messages you want. The 'values' key contains the data, and the 'latlons' method returns the latitudes and longitudes of the grid. Here we find the grib messages for all the 2-meter temp ensemble members valid at Super Bowl time and put the data into a numpy array.¶
In [11]:
grbs.rewind() # rewind the iterator
from datetime import datetime
date_valid = datetime(2014,2,3,0)
t2mens = []
for grb in grbs:
if grb.validDate == date_valid and grb.parameterName == 'Temperature' and grb.level == 2:
t2mens.append(grb.values)
t2mens = np.array(t2mens)
print t2mens.shape, t2mens.min(), t2mens.max()
lats, lons = grb.latlons() # get the lats and lons for the grid.
print 'min/max lat and lon',lats.min(), lats.max(), lons.min(), lons.max()
(51, 181, 360) 217.003723145 316.734893799 min/max lat and lon -90.0 90.0 0.0 359.0
Plot all the ensemble members on a lambert conformal projection centered on NYC.¶
In [8]:
fig = plt.figure(figsize=(16,35))
m = Basemap(projection='lcc',lon_0=-74,lat_0=41,width=4.e6,height=4.e6)
x,y = m(lons,lats)
for nens in range(1,51):
ax = plt.subplot(10,5,nens)
m.drawcoastlines()
cs = m.contourf(x,y,t2mens[nens],np.linspace(230,300,41),cmap=plt.cm.jet,extend='both')
t = plt.title('ens member %s' % nens)
What's the probability it will be below freezing during the game? Note the ECMWF ensemble is forecasting a lower probability of freezing temps are gametime than the NCEP ensemble.¶
In [9]:
freezeprob = 100.*((t2mens < 273).sum(axis=0))/t2mens.shape[0]
print freezeprob.min(), freezeprob.max(), freezeprob.shape
m.drawcoastlines(color='r')
from mpl_toolkits.basemap import cm as basemapcm
cs = m.contourf(x,y,freezeprob,np.linspace(0,100.,41),cmap=basemapcm.GMT_haxby_r)
cb = m.colorbar()
xpt, ypt = m(-74,41)
nyc = m.plot([xpt],[ypt],'ro') # plot red dot near the Meadowlands.
t = plt.title('ens prob t2m < 273 (go Broncos!)')
0.0 100.0 (181, 360)
In [9]: