basemap_pygrib.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Matplotlib/Basemap and pygrib example\n",
      "====================================="
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "netCDF is easier to deal with, but most operational forecast centers provide data in GRIB format. [GRIB](http://en.wikipedia.org/wiki/GRIB) is a record format, where every record is a 2D field. In this example we use [pygrib](http://pygrib.googlecode.com) to read some ECMWF ensemble forecast data, then use matplotlib and Basemap plot forecast maps.  Pygrib uses the ECMWF [GRIB_API](https://software.ecmwf.int/wiki/display/GRIB/Home) C library under the hood."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make the output of plotting commands be displayed inline within the notebook,\n",
      "%matplotlib inline \n",
      "from mpl_toolkits.basemap import Basemap  # import Basemap matplotlib toolkit\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import pygrib # import pygrib interface to grib_api"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Open the grib file (downloaded from the [NCAR TIGGE](http://tigge.ucar.edu/home/home.htm) data portal, available [here](ftp://ftp.cdc.noaa.gov/Public/jwhitaker/ecmwf_ensemble.grb)).  The pygrib.open object behaves much like a python file object."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "grbs = pygrib.open('ecmwf_ensemble.grb')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "The pygrib.open object is a python iterator."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for grb in grbs[:4]:\n",
      "    print grb"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000\n",
        "2:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000:hi res cntl fcst\n",
        "3:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000:pos ens pert 1\n",
        "4:Mean sea level pressure:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201401230000:pos ens pert 2\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "The iterator yields a grib message object.  Each grib message object has a set of attributes, or 'keys', which can be used for searching."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print grb.keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[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']\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "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."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "grbs.rewind() # rewind the iterator\n",
      "from datetime import datetime\n",
      "date_valid = datetime(2014,2,3,0)\n",
      "t2mens = []\n",
      "for grb in grbs:\n",
      "    if grb.validDate == date_valid and grb.parameterName == 'Temperature' and grb.level == 2: \n",
      "        t2mens.append(grb.values)\n",
      "t2mens = np.array(t2mens)\n",
      "print t2mens.shape, t2mens.min(), t2mens.max()\n",
      "lats, lons = grb.latlons()  # get the lats and lons for the grid.\n",
      "print 'min/max lat and lon',lats.min(), lats.max(), lons.min(), lons.max()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(51, 181, 360) 217.003723145 316.734893799\n",
        "min/max lat and lon"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " -90.0 90.0 0.0 359.0\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Plot all the ensemble members on a lambert conformal projection centered on NYC."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(16,35))\n",
      "m = Basemap(projection='lcc',lon_0=-74,lat_0=41,width=4.e6,height=4.e6)\n",
      "x,y = m(lons,lats)\n",
      "for nens in range(1,51):\n",
      "    ax = plt.subplot(10,5,nens)\n",
      "    m.drawcoastlines()\n",
      "    cs = m.contourf(x,y,t2mens[nens],np.linspace(230,300,41),cmap=plt.cm.jet,extend='both')\n",
      "    t = plt.title('ens member %s' % nens)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",