High-resolution deterministic weather forecasts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# High-resolution deterministic forecasts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates using [netcdf4-python](http://netcdf4-python.googlecode.com) to read some [NOMADS](http://nomads.ncep.noaa.gov/txt_descriptions/NOMADS_doc.shtml) [openDAP](http://opendap.org) GFS forecast data, use [numpy](http://numpy.org) to do some rudimentary analysis, and have [matplotlib](http://matplotlib.org) and [Basemap](http://matplotlib.org/basemap) plot the results. \n",
    "\n",
    "We start by loading the required modules. First up is the Matplotlib Basemap package and related matplotlib modules (pyplot, gridspec, cm and axes_grid1), then Numpy, netCDF4-python, some datetime routines and IPython widgets to make this notebook interactive.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from mpl_toolkits.basemap import Basemap  # import Basemap matplotlib toolkit\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "from mpl_toolkits.basemap import cm as basemapcm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from netCDF4 import Dataset, num2date # netcdf4-python module\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "# Import widgets for interactive notebook\n",
    "from ipywidgets import interact, fixed, Dropdown"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A small function to set the appropriate forecast time. There's usually about a 4-5 hour delay between the initial time and the data becoming available on the NOMADS system. This only picks up the 00Z or 12Z runs - a bit more logic would enable the function to return the 06Z or 18Z time stamps as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_forecast_time(delay=6):\n",
    "    now = datetime.utcnow()\n",
    "    print \"Current time: {0}\".format(now)\n",
    "    fcast_time = datetime.utcnow()\n",
    "    if now.hour < delay:\n",
    "        fcast_time = fcast_time - timedelta(1)\n",
    "        fcast_time = fcast_time.replace(hour=delay+6, minute=0)\n",
    "    elif (now.hour >=delay and now.hour < (delay+12)):\n",
    "        fcast_time = fcast_time.replace(hour=0)\n",
    "    elif now.hour >= (delay+12):\n",
    "        fcast_time = fcast_time.replace(hour=12)\n",
    "    print \"Forecast time: {0}\".format(fcast_time)\n",
    "    return fcast_time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we construct the URL for the latest available dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcast_time = get_forecast_time()\n",
    "\n",
    "fcast_date_string = fcast_time.strftime(\"gfs%Y%m%d/gfs_0p25_%Hz\")\n",
    "data_url = 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p25/{0}'.format(fcast_date_string)\n",
    "print \"Dataset URL: \" +  data_url"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfs_fcst = Dataset(data_url, mode='r')\n",
    "print \"\"\n",
    "print \"Dataset description:\"\n",
    "print gfs_fcst # get some summary information about the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = gfs_fcst.variables['time']\n",
    "print time\n",
    "valid_dates = num2date(time[:], time.units).tolist()\n",
    "\n",
    "timelist = [d.strftime('%Y-%m-%d %H:%M') for d in valid_dates]\n",
    "print \"Valid times:\"\n",
    "print timelist\n",
    "levels = gfs_fcst.variables['lev']\n",
    "print \"Levels:\"\n",
    "print [level for level in levels[:]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_rain(timestamp):\n",
    "    idx = timelist.index(timestamp)\n",
    "    apcp = gfs_fcst.variables['apcpsfc'][idx,:,:]\n",
    "    prmsl = gfs_fcst.variables['prmslmsl'][idx,:,:]\n",
    "    hgt1000 = gfs_fcst.variables['hgtprs'][idx,0,:,:]\n",
    "    hgt500 = gfs_fcst.variables['hgtprs'][idx,12,:,:]\n",
    "    thk = (hgt500 - hgt1000)/10.\n",
    "\n",
    "    lats = gfs_fcst.variables['lat'][:]; lons = gfs_fcst.variables['lon'][:]\n",
    "    lons, lats = np.meshgrid(lons, lats)\n",
    "    fig = plt.figure(figsize=(18,10))\n",
    "\n",
    "    m = Basemap(projection='mill', llcrnrlon=130., llcrnrlat=-45, \n",
    "                urcrnrlon=160., urcrnrlat=-25.,resolution='i')\n",
    "    \n",
    "    m.drawcoastlines()\n",
    "    m.drawstates()\n",
    "    m.drawcountries()\n",
    "    x,y = m(lons, lats) # convert lats/lons to map projection coordinates\n",
    "    clevs = [0,.1,.25,.5,.75,1.0,1.5,2.0,3.0,4.0,5.0,7.0,10.0,\n",
    "             15.0,20.0,25.0,30.0,40.0,50.0,60.0,75.0]\n",
    "    \n",
    "\n",
    "    cs = m.contourf(x,y, apcp, clevs, cmap=basemapcm.s3pcpn, extend='max')  # color-filled contours\n",
    "    ct = m.contour(x, y, thk, np.arange(500, 600, 4), colors='0.75', linestyle='--')\n",
    "    cp = m.contour(x, y, prmsl/100., np.arange(900, 1040, 2), colors='k')\n",
    "    cb = m.colorbar(cs, extend='max')  # draw colorbar\n",
    "    parallels = m.drawparallels(np.arange(-50, 0, 10), labels=[1,0,0,0])  # draw parallels, label on left\n",
    "    meridians = m.drawmeridians(np.arange(80, 190, 10), labels=[0,0,0,1]) # label meridians on bottom\n",
    "    \n",
    "    fig.suptitle('Forecast apcp (mm) from %s for %s' % (valid_dates[0],valid_dates[idx]),fontweight='bold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_tsindex(timestamp):\n",
    "    idx = timelist.index(timestamp)\n",
    "    apcp = gfs_fcst.variables['apcpsfc'][idx,:,:]\n",
    "    prmsl = gfs_fcst.variables['prmslmsl'][idx,:,:]\n",
    "    no4lftx = gfs_fcst.variables['no4lftxsfc'][idx,:,:]\n",
    "    cape = gfs_fcst.variables['capesfc'][idx,:,:]\n",
    "    #thk = (hgt500 - hgt1000)/10.\n",
    "\n",
    "    lats = gfs_fcst.variables['lat'][:]; lons = gfs_fcst.variables['lon'][:]\n",
    "    lons, lats = np.meshgrid(lons, lats)\n",
    "    fig = plt.figure(figsize=(18,10))\n",
    "\n",
    "    m = Basemap(projection='mill', llcrnrlon=130., llcrnrlat=-45, \n",
    "                urcrnrlon=160., urcrnrlat=-25.,resolution='i')\n",
    "    \n",
    "    m.drawcoastlines()\n",
    "    m.drawstates()\n",
    "    m.drawcountries()\n",
    "    x,y = m(lons, lats) # convert lats/lons to map projection coordinates\n",
    "    clevs = [0, 10, 25, 50, 75, 100, 150, 200, 300, 400, 500, 750, 1000,\n",
    "             1500, 2000, 2500, 3000, 4000, 5000]\n",
    "\n",
    "    cs = m.contourf(x,y, cape, clevs, cmap=basemapcm.s3pcpn, extend='max')  # color-filled contours\n",
    "    clm = m.contour(x, y, no4lftx, np.arange(-20, 0, 2), colors='0.75', linestyle='--')\n",
    "    clp = m.contour(x, y, no4lftx, np.arange(0, 20, 2), colors='0.75', linestyle='-')\n",
    "    cp = m.contour(x, y, prmsl/100., np.arange(900, 1040, 2), colors='k')\n",
    "    #cc = m.contour(x, y, cape, np.arange(0, 5000, 100), colors='b', linewidth=2)\n",
    "    cb = m.colorbar(cs, extend='max')  # draw colorbar\n",
    "    parallels = m.drawparallels(np.arange(-50, 0, 10), labels=[1,0,0,0])  # draw parallels, label on left\n",
    "    meridians = m.drawmeridians(np.arange(80, 190, 10), labels=[0,0,0,1]) # label meridians on bottom\n",
    "    \n",
    "    fig.suptitle('Forecast CAPE (J/kg), LI (K) for %s' % (valid_dates[idx]),fontweight='bold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_wind(timestamp):\n",
    "    idx = timelist.index(timestamp)\n",
    "    uu = gfs_fcst.variables['ugrdsig995'][idx,:,:]\n",
    "    vv = gfs_fcst.variables['vgrdsig995'][idx,:,:]\n",
    "    wspd = np.sqrt(uu*uu + vv*vv)\n",
    "    pmsl = gfs_fcst.variables['prmslmsl'][idx,:,:]\n",
    "    hgt1000 = gfs_fcst.variables['hgtprs'][idx,0,:,:]\n",
    "    hgt500 = gfs_fcst.variables['hgtprs'][idx,12,:,:]\n",
    "    thk = (hgt500 - hgt1000)/10.\n",
    "\n",
    "\n",
    "    lats = gfs_fcst.variables['lat'][:]; lons = gfs_fcst.variables['lon'][:]\n",
    "    lons, lats = np.meshgrid(lons, lats)\n",
    "    fig = plt.figure(figsize=(18,10))\n",
    "\n",
    "    m = Basemap(projection='mill',llcrnrlon=130., llcrnrlat=-45, \n",
    "                urcrnrlon=160., urcrnrlat=-25.,resolution='i')\n",
    "    m.drawcoastlines()\n",
    "    m.drawstates()\n",
    "    m.drawcountries()\n",
    "    x,y = m(lons, lats) # convert lats/lons to map projection coordinates\n",
    "    clevs = np.arange(5.,41.)\n",
    "\n",
    "    cs = m.contourf(x,y,wspd,clevs,cmap=basemapcm.GMT_haxby_r, extend='both')  # color-filled contours\n",
    "    ct = m.contour(x,y,thk, np.arange(500, 600, 4), colors='0.75', linestyle='--')\n",
    "    cp = m.contour(x,y,pmsl/100.,np.arange(900, 1040, 2), colors='k')\n",
    "    cb = m.colorbar(cs,extend='both')  # draw colorbar\n",
    "    parallels = m.drawparallels(np.arange(-50,0,10),labels=[1,0,0,0])  # draw parallels, label on left\n",
    "    meridians = m.drawmeridians(np.arange(80,190,10),labels=[0,0,0,1]) # label meridians on bottom\n",
    "    \n",
    "    fig.suptitle('Forecast wind speed (m/s) from %s for %s' % (valid_dates[0],valid_dates[idx]),fontweight='bold')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(plot_rain, timestamp=Dropdown(options=timelist, value=timelist[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(plot_wind, timestamp=Dropdown(options=timelist, value=timelist[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(plot_tsindex, timestamp=Dropdown(options=timelist, value=timelist[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}