2_scripting_with_metsim.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scripting with MetSim\n",
    "\n",
    "This notebook provides a tutorial for running MetSim in a scripting environment. \n",
    "The scripting environment allows for a more interactive workflow that will make it easier to manipulate the data that MetSim generates.\n",
    "\n",
    "Then, we will begin to expand on this to show how MetSim can be used programmatically to generate multiple sets of meteorologic data.\n",
    "We will run multiple longwave radiation parameterizations and compare their relative qualities.\n",
    "Finally, we will look at how changing some of the various parameters which MetSim exposes affect the output.\n",
    "\n",
    "First, we begin with some standard imports. For this portion of the tutorial we will be making use of several libraries for analysis and plotting:\n",
    "\n",
    " - `xarray` is used under the hood for managing the output from MetSim\n",
    " - `pandas` is used for setting up the dates for our MetSim run\n",
    " - `holoviews` and `geoviews` will be used for some interactive plotting\n",
    " - `geopandas` is used for used for managing the shapefiles, also used for plotting\n",
    " - `cartopy` is also used for managing the shapefiles, also used for plotting\n",
    " - `metsim`, obviously, will be used for data processing\n",
    "\n",
    "We initialize the notebook by importing these, as well as the `pylab` cell magic, which imports `numpy` and `matplotlib`. \n",
    "We also set up `holoviews` to use the `bokeh` interactive plotting library as the backend.\n",
    "\n",
    "Note that there may be some warnings from various Python libraries throughout this tutorial which are not issues as of the presentation date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import cartopy\n",
    "import geoviews as gv\n",
    "import geopandas as gpd\n",
    "import holoviews as hv\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "from metsim import MetSim\n",
    "hv.notebook_extension('bokeh')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Our study location: the Green River basin\n",
    "\n",
    "For this example we will be showing how MetSim can be used on unstructured mesh arrangements as well as the structured latitude-longitude grid that we ran in the command line examples.\n",
    "We will use a setup which consists of twelve \"hydrologic response units\" or HRU.\n",
    "A HRU typically delineates a watershed by topography, soil type, land use, or other defining feature.\n",
    "The HRU delineation we are using comes from the USGS hydrologic unit codes (HUC), with 12 digit codes. \n",
    "Commonly this is referred to as HUC12, which is the finest resolution.\n",
    "Coarser resolutions are available as HUC10 and HUC8.\n",
    "\n",
    "The 12 HRU that we will be running are located along the Green river, just southeast of Seattle. \n",
    "It originates on the western slopes of the Cascade mountains.\n",
    "The upper portion of the Green river provides much of the water for the city of Tacoma, which is south of Seattle.\n",
    "The Green river ultimately ends up as a tributary to the Duwamish river, which flows through Seattle and into Puget Sound.\n",
    "\n",
    "Let's take a look. \n",
    "Because this is an unstructured mesh we will be defining the mesh elements by their respective HUC12 identifiers. \n",
    "We can draw boundaries using a shapefile, which we have obtained from the USGS datasets and included in the tutorial data.\n",
    "There are several tools for plotting maps in Python using shapefiles, most notably Cartopy, GeoPandas, and GeoViews.\n",
    "Other tools for more general GIS functionality include GDAL, Rasterio, Shapely, and Fiona.\n",
    "We will be using Cartopy and GeoPandas to define the geometries and then a combination of GeoViews and HoloViews to actually draw the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapefile = './shapefile/shapefile.shp'\n",
    "shapes = cartopy.io.shapereader.Reader(shapefile)\n",
    "(gv.tile_sources.StamenTerrainRetina\n",
    " * gv.Shape.from_shapefile(shapefile, crs=cartopy.crs.PlateCarree()).opts(\n",
    "     style=dict(fill_color='honeydew', line_color='darkolivegreen', alpha=0.8))\n",
    " * gv.Points([(-122.3321, 47.9),(-122.8443, 47.2529), (-122.5, 47.0)]).opts(style=dict(size=0))\n",
    ").opts(width=900, height=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting set up: Parameters and MetSim instantiation\n",
    "\n",
    "Unlike the command line usage of MetSim, while in the scripting enviornment we will provide the confiugration as a Python dictionary. \n",
    "The information that this dictionary, which we will call `params`, is similar to that which is found in the configuration file.\n",
    "However, there are some notable differences which we will walk through.\n",
    "\n",
    "First, note that we must provide dates as `datetime` objects. \n",
    "This is easily handled via `pandas`. \n",
    "Second, we point out a new parameter, the `scheduler`; we have set it to `threading` for a specific purpose.\n",
    "As mentioned in the command line portion of the tutorial, the scheduler is how MetSim manages jobs under the hood.\n",
    "The default value is the `distributed` scheduler, which does not currently support running in this sort of interactive environment.\n",
    "By setting `scheduler` to `threading` we are able to work in this environment.\n",
    "A notable drawback to this, is that parallelism is significantly hindered. \n",
    "We hope that future versions of `dask` (the library which we use as a job manager) improve this compatibility, and we can make it available in MetSim.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rest of the workflow\n",
    "dates = pd.date_range('1971/01/01', '1971/12/31')\n",
    "params = {\n",
    "    # Run parameters, timestep, start, and stop dates\n",
    "    'time_step'    : \"60\",       \n",
    "    'start'        : dates[0],\n",
    "    'stop'         : dates[-1],\n",
    "    # Where all of our input files are\n",
    "    'forcing'      : './input/green_river_forcing.nc',     \n",
    "    'domain'       : './input/green_river_domain.nc',\n",
    "    'state'        : './input/green_river_state.nc',\n",
    "    # The format of the forcing data\n",
    "    'forcing_fmt'  : 'netcdf',\n",
    "    # Specification of the output\n",
    "    'out_dir'      : './output',\n",
    "    'out_prefix': 'green_river',\n",
    "    # Scheduler to use for job management - see text above for more information\n",
    "    'scheduler'    : 'threading',\n",
    "    # This is analagous to the [chunks] section from the config file\n",
    "    'chunks'       : \n",
    "        {'hru': 3},\n",
    "    # This is analagous to the [forcing_vars] section from the config file\n",
    "    'forcing_vars' : \n",
    "        {'prcp' : 'prec', 'tmax': 't_max', 'tmin': 't_min'},\n",
    "    # This is analagous to the [state_vars] section from the config file\n",
    "    'state_vars'   : \n",
    "        {'prcp' : 'prec', 'tmax': 't_max', 'tmin': 't_min'},\n",
    "    # This is analagous to the [domain_vars] section from the config file\n",
    "    'domain_vars'  : \n",
    "        {'elevation': 'elev', 'latitude': 'lat', 'longitude': 'lon', 'mask': 'mask'}\n",
    "    }               "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Okay, let's run MetSim!\n",
    "\n",
    "## First, let's create our MetSim object\n",
    "For this step, we simply give the parameters as an argument to the MetSim object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms = MetSim(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now, we can run the simulation\n",
    "This is as easy as calling the object's `run` method. Try it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finally, let's open up the output and see what we've got\n",
    "Here we can see that we've got a dataset that looks like we would expect. \n",
    "It's got the twelve HRU each with 8760 timesteps (equal to 365 * 24).\n",
    "We also have 6 variables:\n",
    " - temperature\n",
    " - precipitation\n",
    " - shortwave\n",
    " - longwave\n",
    " - vapor pressure\n",
    " - relative humidity\n",
    " \n",
    "The output is simply an `xarray` dataset, which was read in from the underlying NetCDF file that the `ms.run()` command generated.\n",
    "You can do all of the standard things that were taught in the cyberseminars.\n",
    "Let's start by just printing out what's there like before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = ms.open_output().load()\n",
    "out = out.sortby('hru')\n",
    "gdf = gpd.read_file(shapefile)\n",
    "print(out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interacting with the output\n",
    "\n",
    "Now that we've run MetSim in a notebook we can look at how you might interact with the output in a more meaningful way.\n",
    "Because we are using an unstructured mesh for the domain we will have to use a bit fancier of techniques than the basic `xarray` commands provide.\n",
    "Since this is a workshop and we've got a real live audience it might be nice to build something interactive.\n",
    "We will start by building a map with a time slider that let's us see the values of the output at various times.\n",
    "Before doing so, however, let's take a brief detour to explain _dictionary comprehensions_.\n",
    "\n",
    "## Dictionary comprehensions\n",
    "\n",
    "Dictionary comprehensions are a Python built in functionality which are shorthand for writing for loops.\n",
    "They, along with _list comprehensions_ are very powerful and concise ways to write code.\n",
    "We provide a simple example of how they work below.\n",
    "Following that, we will define a function which will take a variable name and timestep and will produce a map of the variable at that timestep from our MetSim output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphabet = 'abcdefg'\n",
    "\n",
    "# First, a standard for filling up a dictionary:\n",
    "index_of_letter = {}\n",
    "for i in range(len(alphabet)):\n",
    "    index_of_letter[alphabet[i]] = i\n",
    "    \n",
    "# Versus the dictionary comprehension\n",
    "index2 = {alphabet[i]: i for i in range(len(alphabet))}\n",
    "\n",
    "print(index_of_letter)\n",
    "print(index2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A mapping function for our output\n",
    "Below we define the plotting function that we will use to build our interactive plots.\n",
    "It is heavily commented to make things clear, but the basic ideas is we will take a specific variable at a specific time and plot it using our shapefile for the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def var_map(var, timestep):\n",
    "    \"\"\"\n",
    "    Create a map for a variable at a given timestep\n",
    "    Warning: This function is state dependent! \n",
    "    \n",
    "    Params:\n",
    "        var: the variable of interest\n",
    "        timestep: the timestep to plot at\n",
    "    \n",
    "    Returns:\n",
    "        poly: A geoviews \"shape\" plot with the desired data\n",
    "    \"\"\"\n",
    "    # Convert the data from an xarray dataset to a pandas dataframe \n",
    "    # This is necessary for geoviews to be able to plot it\n",
    "    out_df = out[var].sel(time=timestep).to_dataframe()\n",
    "    # Make sure we have some metadata to join with the shapefile\n",
    "    out_df['HUC_12'] = gdf['HUC_12'].values\n",
    "    out_df['hru_name'] = gdf['FIRST_HU_3'].values\n",
    "    # Create the shape plot - some keys:\n",
    "    #  - shapes.records() provides the geometry from the shapefiles\n",
    "    #  - out_df provides the data\n",
    "    #  - value=var chooses which column of the dataframe to use as data for coloring\n",
    "    #  - index=['hru_name'] will provide a name on mouse hover\n",
    "    #  - on='HUC_12' is the key to join the `shapes` with `out_df`\n",
    "    poly = gv.Shape.from_records(shapes.records(), out_df, value=var, index=['hru_name'], on='HUC_12')\n",
    "    # Add some options to make things a bit nicer\n",
    "    #  - width=700, sets the width, as we expect\n",
    "    #  - cmap='plasma' sets the colormap to 'plasma'\n",
    "    #  - tools=['hover'] provides information on mouse hover over\n",
    "    #  - colorbar=True adds a colorbar\n",
    "    #  - alpha=0.7 adds a bit of transparency\n",
    "    poly = poly.opts(width=700, cmap='plasma', tools=['hover'], colorbar=True, alpha=0.7)\n",
    "    return poly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WARNING: Occasionally these plots don't like to render - try rerunning a few times and they should show up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's try it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_map('temp', out.time.values[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Okay, that was nice, but not very explorable...\n",
    "\n",
    "We can do one better by making use of the dictionary comprehensions we just discussed.\n",
    "First we construct a dictionary containing timestampes as keys and maps from the previous plot function as values.\n",
    "HoloViews provides a nice `HoloMap` function which will then take this mapping and let us interact with the whole thing via a slider.\n",
    "Note here that the `Map` in `HoloMap` is `Map` [in the programming sense](https://en.wikipedia.org/wiki/Map_(higher-order_function)), not the cartographic sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapping = {t: var_map('temp', t) for t in out.time.values[0:12]}\n",
    "hv.HoloMap(mapping, kdims=['timestep'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One more layer\n",
    "\n",
    "We can continue to improve this by adding a timeseries to the side, which will show the mean of the entire domain along with the minimum and maximum for a given variable.\n",
    "This will let us quickly explore the output dataset and see what's going on in both space and time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def interactive_plot(var, start, stop, step):\n",
    "    # Just as before, use the `var_map` function to build our map\n",
    "    poly = hv.HoloMap({t: var_map(var, t) for t in out.time.values[start:stop:step]}, kdims=['timestep'])\n",
    "    # Vlines will give us an indicator of which time slice we are looking at\n",
    "    vlines = hv.HoloMap({t: hv.VLine(t) for t in out.time.values[start:stop:step]}, kdims=['timestep'])\n",
    "    # This seems to make the plot appear more often - probably a holo/geoviews bug somewhere \n",
    "    poly\n",
    "    \n",
    "    # Calculate min, max, and mean over the domain for each time\n",
    "    vmin = out[var].isel(time=slice(start, stop)).min(dim='hru').rolling(time=step).min().values\n",
    "    vmax = out[var].isel(time=slice(start, stop)).max(dim='hru').rolling(time=step).max().values\n",
    "    vmean = out[var].isel(time=slice(start, stop)).mean('hru').rolling(center=True, time=step).mean().dropna('time')\n",
    "   \n",
    "    # Build the complete interactive plot. Layers are as follows\n",
    "    #  - gv.tile_sources.EsriTerrain- Add a background with topography\n",
    "    #  - * poly - overlay our map onto the background\n",
    "    #  - + hv.Area... - add the area plot to the right\n",
    "    #  - * hv.Curve... - overlay the mean curve onto the area plot\n",
    "    #  - * vlines... - overlay the timestep indicator\n",
    "    return (gv.tile_sources.EsriTerrain\n",
    "            * poly\n",
    "            + (hv\n",
    "               .Area((out.time[start:stop], vmin, vmax), vdims=['vmin', 'vmax'])\n",
    "               .opts(alpha=0.5, color='gold', line_color=None)\n",
    "               .redim.label(x='Date', vmin=var.capitalize()))\n",
    "            * hv.Curve(vmean).opts(color='purple', alpha=0.8)\n",
    "            * vlines.opts(alpha=0.4, color='red'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interactive_plot('temp', start=0, stop=120, step=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploring the effects of different options\n",
    "\n",
    "Great, so now we've run MetSim in our notebook and have gotten a handle on what's going on in the output.\n",
    "We will now dig in slightly deeper and start playing with some of the options which MetSim provides.\n",
    "First we will look at different longwave radiation parameterizations, then we will follow that with a sensitivity analysis of some of the parameters.\n",
    "\n",
    "## Longwave radiation parameterizations\n",
    "We have implemented 6 different longwave parameterizations in MetSim.\n",
    "By default we use the Prata parameterization, but let's see what the others look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Longwave parameterizations to try\n",
    "lw_types = ['TVA', 'Anderson', 'Brutsaert', 'Idso', 'Satterlund', 'Prata']\n",
    "lw_vals = {}\n",
    "for lw in lw_types:\n",
    "    print(lw)\n",
    "    params['lw_type'] = lw.lower()\n",
    "    ms = MetSim(params)\n",
    "    ms.run()\n",
    "    ds = ms.open_output().load()\n",
    "    lw_vals[lw] = ds['longwave']\n",
    "    ds.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,6))\n",
    "for name, lw in lw_vals.items():\n",
    "    lw.isel(hru=2, time=slice(90,130)).plot(label=name, linewidth=2)\n",
    "plt.legend()\n",
    "plt.ylabel('Longwave Radiation ($W/m^2$)')\n",
    "plt.xlabel('Date')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sensitivity analysis: Cloudiness parameterization\n",
    "\n",
    "The jumps in the longwave plot from above are due to a precipitation event, which causes an estimation of cloud cover to affect the estimated radiation.\n",
    "Let's try varying the `sw_prec_thresh` and `rain_scalar` parameters, which control corrections during periods of precipitation, effectively parameterizing cloud cover.\n",
    "The `sw_prec_thresh` determines the amount of precipitation that must be exceeded to be considered \"cloudy\".\n",
    "By default this is set to 0, so any amount of precipitation scales the shortwave radiation.\n",
    "It is scaled by the `rain_scalar` factor, which is 0.75 by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params['lw_type'] = 'prata'\n",
    "sw_vals = {}\n",
    "spt_vals = [0.0, 12*0.02, 12*0.05]\n",
    "rs_vals = [0.6, 0.75, 0.9]\n",
    "for spt in spt_vals:\n",
    "    for rs in rs_vals:\n",
    "        params['sw_prec_thresh'] = spt\n",
    "        params['rain_scalar'] = rs\n",
    "        ms = MetSim(params)\n",
    "        ms.run()\n",
    "        output = ms.open_output().load()\n",
    "        rs_dict = sw_vals.get(spt, {})\n",
    "        rs_dict[rs] = output['shortwave']\n",
    "        sw_vals[spt] = rs_dict\n",
    "        output.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def param_map(spt, rs, timestep):\n",
    "    \"\"\"\n",
    "    Create a map for a variable at a given timestep\n",
    "    Warning: This function is state dependent! \n",
    "    \n",
    "    Params:\n",
    "        spt: the value for sw_prec_thresh\n",
    "        rs: the value for rain_scalar\n",
    "        timestep: the timestep to plot at\n",
    "    \n",
    "    Returns:\n",
    "        poly: A geoviews \"shape\" plot with the desired data\n",
    "    \"\"\"\n",
    "    out_df = sw_vals[spt][rs].sel(time=timestep).to_dataframe()\n",
    "    out_df['HUC_12'] = gdf['HUC_12'].values\n",
    "    out_df['hru_name'] = gdf['FIRST_HU_3'].values\n",
    "    poly = gv.Shape.from_records(shapes.records(), out_df, value='shortwave', index=['hru_name'], on='HUC_12')\n",
    "    poly = poly.opts(width=700, cmap='plasma', tools=['hover'], colorbar=True, alpha=0.7)\n",
    "    return poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = ms.open_output().load()\n",
    "output['prec'].isel(time=slice(90,130)).plot.line(x='time', add_legend=False, figsize=(12, 6));\n",
    "output.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapping = {(spt, rs, t): param_map(spt, rs, t) for t in out.time.values[90:130] \n",
    "                                               for spt in spt_vals \n",
    "                                               for rs in rs_vals}\n",
    "(gv.tile_sources.EsriTerrain\n",
    " * hv.HoloMap(mapping, kdims=['sw_prec_thresh', 'rain_scalar', 'timestep']))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "metsim_whw",
   "language": "python",
   "name": "metsim_whw"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}