Introduction to Jupyter notebooks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 1: Introduction to Jupyter Notebooks\n",
    "### MAQ - 32806, Chiel van Heerwaarden & Imme Benedict, 2016\n",
    "This is the first tutorial. In this tutorial, you will get used to the Jupyter Notebooks, and their work flow. In Jupyter notebooks you can write either text in Markdown format (including mathematics using LaTeX) or code in Python format. In this course, you will work on assignments that you are going to work out in your own notebook.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SOME ESSENTIAL TIPS BEFORE STARTING!\n",
    "* Make a copy of the notebook via the `File` menu. By doing so, you always have the original file to revert to.\n",
    "* Save your work often.\n",
    "* In case of weird behavior, try to restart the notebook via the `Kernel` menu.\n",
    "* Play around, the best way to find out how things work is to try them out!\n",
    "* Use the internet for help on `numpy`, `matplotlib`, `Basemap` or `LaTeX` if you want to try something out.\n",
    "* With `Shift+Enter` you can move down while executing cells.\n",
    "* You can complete commands by pressing `Tab`.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up a Python Basemap and read out 3D fields\n",
    "We start the tutorial by loading the required Python packages and setting the figure properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np              # Numpy is the package for scientific computing in Python.\n",
    "import netCDF4 as nc            # NetCDF is the data format of our meteorological data.\n",
    "import matplotlib.pyplot as pl  # Matplotlib is a scientific plotting package.\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap # Import the map plotting interface.\n",
    "\n",
    "# The statement below enforces the plots to be put inline.\n",
    "%pylab inline\n",
    "pl.rcParams.update({'font.size': 12})          # Set the font size of the plots to 12pt.\n",
    "pl.rcParams.update({'figure.figsize': [13,6]}) # Set the standard figure size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "With the code below, you initialize the world map with the desired coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lon_start = -150.\n",
    "lon_end   =   50.\n",
    "lat_start =   20.\n",
    "lat_end   =   70.\n",
    "m = Basemap(llcrnrlon=lon_start,\n",
    "            llcrnrlat=lat_start,\n",
    "            urcrnrlon=lon_end,\n",
    "            urcrnrlat=lat_end,\n",
    "            projection='mill', resolution='l')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "Now, you load the data and read out the desired variables from a NetCDF file from the ECMWF ERA-Interim data archive. You are loading a file with data from October 1st, 2015."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nc_file = nc.Dataset(\"data/era_data.nc\", \"r\")\n",
    "lat = nc_file.variables[\"latitude\"][:]\n",
    "lon = nc_file.variables[\"longitude\"][:]\n",
    "p = nc_file.variables[\"level\"][:]*100.  # Conversion of pressure to Pa.\n",
    "\n",
    "t = 0 # We only load the first time index from the data.\n",
    "time = nc_file.variables[\"time\"][0]\n",
    "\n",
    "# With the code below, we roll the map, to get the 0 meridian,\n",
    "# rather than the Pacific in the middle.\n",
    "nroll = lon.size//2\n",
    "lon = np.roll(lon, nroll)\n",
    "lon = np.where(lon>=180., lon-360., lon)\n",
    "\n",
    "# We load 3D fields of the two horizontal wind components, the geopotential and the temperature.\n",
    "u     = np.roll(nc_file.variables[\"u\"][t,:,:,:], nroll, -1)\n",
    "v     = np.roll(nc_file.variables[\"v\"][t,:,:,:], nroll, -1)\n",
    "omega = np.roll(nc_file.variables[\"w\"][t,:,:,:], nroll, -1)\n",
    "Phi   = np.roll(nc_file.variables[\"z\"][t,:,:,:], nroll, -1)\n",
    "T     = np.roll(nc_file.variables[\"t\"][t,:,:,:], nroll, -1)\n",
    "\n",
    "nc_file.close() # The file access is closed as no further data is needed now.\n",
    "\n",
    "# We create two-dimensional grids for lons and lats in degrees and radians. These arrays\n",
    "# are used for plotting and computation.\n",
    "lons, lats = m(*np.meshgrid(lon, lat))\n",
    "lonrad, latrad = np.meshgrid(np.pi/180.*lon, np.pi/180.*lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is it. With the actions below, you have done everything to be able to work with ECMWF ERA-Interim weather data.\n",
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with numpy\n",
    "Before plotting the data, we have a quick look at the shape of the loaded data. Our 3D fields have three dimensions: `height`, `lat` and `lon` from which we can select the desired fields. Pressure is the height coordinate of our data. Below we show some simple examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(Phi.shape)     # We have 37 height levels, 241 latitudes and 480 longitudes.\n",
    "print(p[0])        # Our lowest pressure level is 1000 hPa.\n",
    "print(p[6])        # At level 6 we find 850 hPa.\n",
    "Phi_850 = Phi[6,:,:]   # We create a new variable that contains only the geopotential height at 850 hPa.\n",
    "print(Phi_850.shape) # The new data has only two dimensions, lat and lon.\n",
    "\n",
    "# We can also write formatted output to the screen.\n",
    "print(\"The standard deviation is {0}\".format(np.std(Phi_850))) # We can format the output.\n",
    "print(\"The mean, min and max are {0}, {1}, {2}\".format(np.mean(Phi_850), np.min(Phi_850), np.max(Phi_850)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Plotting data\n",
    "We are now going to plot the temperature $T$ and geopotential $\\Phi$ at 500 hPa. Note that our units are SI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "p_plot = 50000.            # Choose the height in pressure level at which you want to plot.\n",
    "n = abs(p-p_plot).argmin() # This is a trick to find the index closest to the desired value.\n",
    "\n",
    "m.pcolormesh(lons, lats, T[n,:,:], cmap=pl.cm.gist_ncar)        # We plot a colormesh using the gist_ncar colormap.\n",
    "m.colorbar()                                                    # We add a colorbar to show the values of temperature.\n",
    "cz = m.contour(lons, lats, Phi[n,:,:], 15, colors='w', linewidths=1.5) # We plot the geopotential in contours.\n",
    "pl.clabel(cz, fmt='%1.0f', fontsize=10.) # We add labels to the contour lines.\n",
    "m.drawcoastlines(color='#222222') # We add coastlines.\n",
    "# Draw parallels.\n",
    "parallels = np.arange(0.,90,30.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# Draw meridians.\n",
    "meridians = np.arange(0.,360.,60.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('T and Phi at p = {0:.0f} hPa'.format(p[n]/100.)); # We add a title to the plot. \n",
    "pl.tight_layout() # With this function we make the figure fit as good as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "It is also possible to make simple line plots. Below we plot the temperature profile in Wageningen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lat_plot = 51.97\n",
    "lon_plot = 5.67\n",
    "i_lon = abs(lon - lon_plot).argmin()\n",
    "j_lat = abs(lat - lat_plot).argmin()\n",
    "\n",
    "pl.subplot(121)\n",
    "pl.plot(T[:,j_lat,i_lon], p, 'bo-')\n",
    "pl.xlabel('T (K)')\n",
    "pl.ylabel('p (Pa)')\n",
    "pl.gca().invert_yaxis()\n",
    "\n",
    "pl.subplot(122)\n",
    "pl.plot(u[:,j_lat,i_lon], p, 'r^:' , label='u')\n",
    "pl.plot(v[:,j_lat,i_lon], p, 'm*--', label='v')\n",
    "pl.xlabel('wind speed (m/s)')\n",
    "pl.ylabel('p (Pa)')\n",
    "pl.legend(loc=0, frameon=False)\n",
    "pl.gca().invert_yaxis()\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Manipulating data\n",
    "In Python, the underlying language of the notebooks, it is very easy to do manipulations on the data. Below, we calculate the 3D field of total wind speed and plot it. Again, we plot the data on the previously defined index of the 500 hPa level `n`. We calculate the total wind speed $U$ as\n",
    "$$ U = \\sqrt{u^2 + v^2} $$\n",
    "\n",
    "The equation above has been rendered using $\\LaTeX$, which is very useful for typesetting equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "U = (u**2 + v**2)**.5\n",
    "\n",
    "m.pcolormesh(lons, lats, U[n,:,:], cmap=pl.cm.gist_ncar)\n",
    "m.colorbar()\n",
    "m.contour(lons, lats, Phi[n,:,:], 15, colors='w', linewidths=1.5)\n",
    "m.drawcoastlines(color='#222222')\n",
    "# draw parallels.\n",
    "parallels = np.arange(0.,90,30.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,60.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('U at p = {0:.0f} hPa'.format(p[n]/100.));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "The `numpy` package provides function for averaging arrays over a chosen axis. Below we calculate and plot the vertical profile of zonal mean temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T_zonal_mean = np.mean(T, axis=2) # Compute the mean over the longitude axis.\n",
    "\n",
    "T_contours = np.arange(180, 320, 10)                       # Define the contour levels.\n",
    "pl.contourf(lat, p, T_zonal_mean, T_contours)              # Plot a colored contour plot\n",
    "pl.colorbar()                                              # Plot a colorbar.\n",
    "pl.contour (lat, p, T_zonal_mean, T_contours, colors='k')  # Plot line contours.\n",
    "pl.gca().invert_yaxis()                                    # Reverse the y-axis.\n",
    "pl.xlabel('lat (deg)')\n",
    "pl.ylabel('p (Pa)')\n",
    "pl.xlim(-90, 90)\n",
    "pl.ylim(1e5, 2e3)\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A WORD OF WARNING\n",
    "The `axis=` specification can be quite confusing. The field `T[:,:,:]` has three dimensions. The first (`axis=0`) is the height in pressure levels, the second (`axis=1`) is the latitude and the third (`axis=2`) is the longitude. The field `T[0,:,:]`, however, has only two dimensions, because the height disappears as soon as we select one level. This means that the axis that corresponds to latitude is `axis=0`, and longitude is `axis=1`. The Python code below illustrates this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('The shape of T[:,:,:] is {0}'.format(T[:,:,:].shape))\n",
    "print('The length of axis=0 is {0}'.format(T[:,:,:].shape[0]))\n",
    "print('The shape of T[0,:,:] is {0}'.format(T[0,:,:].shape))\n",
    "print('The length of axis=0 is {0}'.format(T[0,:,:].shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Computing differential operators\n",
    "With the example below, you can calculate the gradient of the geopotential in the zonal and meridional direction. This is a tricky example that needs to be studied with care. In spherical coordinates, the two horizontal components of the gradient of the geopotential $\\nabla\\Phi$ are\n",
    "\n",
    "$$\n",
    "\\dfrac{\\partial \\Phi}{\\partial x} = \\dfrac{1}{a \\cos \\phi} \\dfrac{\\partial \\Phi}{\\partial \\lambda},\\\\\n",
    "\\dfrac{\\partial \\Phi}{\\partial y} = \\dfrac{1}{a} \\dfrac{\\partial \\Phi}{\\partial \\phi},\n",
    "$$\n",
    "\n",
    "where $a$ is the radius of the earth, $\\phi$ the latitude and $\\lambda$ the longitude.\n",
    "With `numpy`, these terms can be computed using the `np.gradient` function. This function applies centered finite differences and assumes that the distance between nodes is equal to unity. Therefore, the gradients needs to be computed as\n",
    "\n",
    "$$\n",
    "\\dfrac{\\partial \\Phi}{\\partial x} \\approx \\dfrac{1}{a \\cos \\phi} \\dfrac{\\delta \\Phi}{\\delta \\lambda} =  \\dfrac{1}{a \\cos \\phi} \\dfrac{\\dfrac{\\delta \\Phi}{\\delta i}}{\\dfrac{\\delta \\lambda}{\\delta i}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_earth = 6.37e6 # Here, we define the radius of the earth.\n",
    "\n",
    "# Two arrays of zeros are initialized with the same shape as that of Z.\n",
    "dPhidx = np.zeros(Phi.shape)\n",
    "dPhidy = np.zeros(Phi.shape)\n",
    "\n",
    "# We initialize an array with the cosine of the latitude.\n",
    "cos_lat = np.cos(latrad)\n",
    "\n",
    "# Below, we calculate the two partial derivatives defined in the equations above. Note that axis keyword\n",
    "# indicates the direction over which the gradient is computed. We create a loop here over all pressure\n",
    "# levels to have a 3D field.\n",
    "for k in range(p.size):\n",
    "    dPhidx[k,:,:] = 1./(r_earth*cos_lat) * np.gradient(Phi[k,:,:], axis=1) / np.gradient(lonrad, axis=1)\n",
    "\n",
    "for k in range(p.size):\n",
    "    dPhidy[k,:,:] = (1./r_earth) * np.gradient(Phi[k,:,:], axis=0) / np.gradient(latrad, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m.pcolormesh(lons, lats, dPhidy[n,:,:], cmap=pl.cm.gist_ncar)\n",
    "m.colorbar()\n",
    "m.contour(lons, lats, Phi[n,:,:], 15, colors='w', linewidths=1.5)\n",
    "m.drawcoastlines(color='#222222')\n",
    "# draw parallels.\n",
    "parallels = np.arange(0.,90,30.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,60.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('$\\partial \\Phi/\\partial y$ at p = {0:.0f} hPa'.format(p[n]/100));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "In a similar fashion, it is possible to compute vertical gradients. Note that your data is in pressure coordinates. Vertical gradients will therefore expressed in a change per unit of Pascal. In the example below, we calculate the vertical gradient of temperature\n",
    "\n",
    "$$ \\dfrac{\\partial T}{\\partial p} \\approx \\dfrac{\\delta T}{\\delta p} = \\dfrac{\\dfrac{\\delta T}{\\delta k}}{\\dfrac{\\delta p}{\\delta k}}.$$\n",
    "\n",
    "We plot the data again at the previously chosen height. Our pressure levels are only a one-dimensional array, but `numpy` provides a trick to make the array pretend to be three-dimensional like the temperature array. With the `newaxis` keyword, we can add an additional axis that takes automatically the correct values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dTdp = np.gradient(T[:,:,:], axis=0) / np.gradient(p[:, np.newaxis, np.newaxis], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m.pcolormesh(lons, lats, dTdp[n,:,:], cmap=pl.cm.jet)\n",
    "m.colorbar()\n",
    "m.contour(lons, lats, Phi[n,:,:], 15, colors='w', linewidths=1.5)\n",
    "m.drawcoastlines(color='#222222')\n",
    "# draw parallels.\n",
    "parallels = np.arange(0.,90,30.)\n",
    "m.drawparallels(parallels,labels=[1,0,0,0])\n",
    "# draw meridians\n",
    "meridians = np.arange(0.,360.,60.)\n",
    "m.drawmeridians(meridians,labels=[0,0,0,1])\n",
    "pl.title('dT/dp at p = {0:.0f} hPa'.format(p[n]/100.));\n",
    "pl.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "### Exercises\n",
    "You have reached the end of the examples. Now it is time for you to develop your own skills. Below you find exercises that will help you to develop the skills you need in the coming weeks.\n",
    "___\n",
    "#### Exercise 1\n",
    "1. Plot the temperature at 850 hPa and at 300 hPa.\n",
    "2. Calculate and plot the temperature difference between 850 and 300 hPa.\n",
    "3. Plot the temperature at 850 hPa and 300 hPa at a latitude of 45 degrees in an xy-plot that has longitude on the horizontal axis and temperature on the vertical.\n",
    "___\n",
    "#### Exercise 2\n",
    "1. Calculate and plot the zonal gradient in geopotential height at 500 hPa.\n",
    "2. Calculate and plot the meridional gradient in geopotential height at 500 hPa.\n",
    "3. Calculate and plot the magnitude of the horizontal gradient vector.\n",
    "$\\nabla_h \\Phi \\equiv \\left(\\dfrac{\\partial \\Phi}{\\partial x}, \\dfrac{\\partial \\Phi}{\\partial y}, 0 \\right)$\n",
    "___\n",
    "#### Exercise 3\n",
    "1. Calculate the 3D field of potential temperature $\\theta \\equiv T \\left( \\dfrac{p_0}{p} \\right)^{R_d/c_p}$\n",
    "2. Calculate and plot the vertical cross section of zonal mean potential temperature.\n",
    "3. Calculate the zonal mean zonal wind and add it to the previous plot. Do you see a jetstream?\n",
    "___\n",
    "#### Exercise 4\n",
    "1. Calculate all terms of the conservation of mass equation in pressure coordinates $\\dfrac{\\partial u}{\\partial x} + \\dfrac{\\partial v}{\\partial y} + \\dfrac{\\partial \\omega}{\\partial p}$ and check whether they add to zero.\n",
    "2. Calculate the horizontal divergence at 300 hPa.\n",
    "___\n",
    "#### Exercise 5\n",
    "1. Calculate (only in math, either on paper or in the notebook using $\\LaTeX$) the cross product of the nabla operator and the wind vector $\\omega = \\nabla \\times \\boldsymbol{U}$ and explain the meaning of $\\omega$.\n",
    "2. Calculate (only in math, either on paper or in the notebook using $\\LaTeX$) $\\zeta = \\boldsymbol{k} \\cdot \\left( \\nabla \\times \\boldsymbol{U} \\right)$ and explain the meaning of $\\zeta$.\n",
    "3. Calculate and plot $\\zeta$ at the 500 hPa level.\n",
    "\n",
    "#### Exercise 6 (bonus)\n",
    "1. Calculate the derivative in the meridional direction of the geopotential $\\dfrac{\\partial \\Phi}{\\partial y}$ without using a loop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "Work out the exercises below here.\n",
    "___"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}