3_creating_input_datasets.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating input dataset for MetSim\n",
    "\n",
    "Alright, now that you've gotten a feel for how to work with MetSim you may be wondering how to bring in new data.\n",
    "In this portion of the tutorial we will put together a new MetSim setup for Reynolds Creek.\n",
    "This is an experimental watershed in southwestern Idaho.\n",
    "We will be looking at a single location for a year.\n",
    "The data was downloaded in csv format, which we will transform into NetCDF inputs.\n",
    "As usual, we begin with some standard imports."
   ]
  },
  {
   "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",
    "from metsim import MetSim\n",
    "\n",
    "pylab.rcParams['figure.figsize'] = (10, 6)\n",
    "plt.rcParams['figure.dpi'] = 96\n",
    "hv.notebook_extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(gv.tile_sources.StamenTerrainRetina\n",
    " * gv.Points([(-116.51, 43.17)]).opts(style=dict(size=16, color='red'))\n",
    " * gv.Points([(-122.3321, 49.9),(-122.8443, 47.2529), (-122.5, 42.0)]).opts(style=dict(size=0))\n",
    ").opts(width=900, height=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Put together the required meteorological data\n",
    "We've gathered some data form the Reynolds Creek site that we will use as input.\n",
    "There are two sets of CSVs with precipitation and temperature data, respectively.\n",
    "We data for 2009 and 2010, and will be generating MetSim input for 2010.\n",
    "To do this we must first convert it into an `xarray` dataset.\n",
    "Before doing that though, let's just open up one of the temperature files and see what we're working with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "head -n 10 ./reynolds_creek_data/daily/temp_2010.csv\n",
    "echo \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see that we have several pieces of information to weed through.\n",
    "Luckily, pandas makes this quite easy. We can load this into a dataframe and easily select out the columns that we want.\n",
    "As we can see from the header, we are going to be interested in the daily values `TMAX.D-1` and `TMIN.D-1`.\n",
    "Similarly, the column we are interested in, in the precipitation data is `PREC.I-1`.\n",
    "So, let's load things up.\n",
    "\n",
    "## Create the dataset with the relevant dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dates = pd.date_range('1/1/2010', '12/31/2010')\n",
    "shape = (len(dates), 1, 1, )\n",
    "dims = ('time', 'lat', 'lon', )\n",
    "\n",
    "# We are running only one site, at these coordinates\n",
    "lats = [43.17]\n",
    "lons = [-116.51]\n",
    "elev = 1706.90 # meters\n",
    "coords = {'time': dates, 'lat': lats, 'lon': lons}\n",
    "\n",
    "# Create the initial met data input data structure\n",
    "met_data = xr.Dataset(coords=coords)\n",
    "met_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the actual data arrays to put the data into."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for varname in ['prec', 't_min', 't_max']:\n",
    "    met_data[varname] = xr.DataArray(data=np.full(shape, np.nan),\n",
    "                                     coords=coords, dims=dims,\n",
    "                                     name=varname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in the data and put it into the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the precipitation data\n",
    "df = pd.read_csv(\"./reynolds_creek_data/daily/precip_2010.csv\", skiprows=[0,1])\n",
    "# Convert to mm\n",
    "prec_vals = df['PREC.I-1 (in) '].diff().values[1:] * 25.4\n",
    "met_data['prec'].values[:, 0, 0] = prec_vals\n",
    "\n",
    "# And now temp data\n",
    "df = pd.read_csv(\"./reynolds_creek_data/daily/temp_2010.csv\", skiprows=[0,1])\n",
    "tmin_vals = df['TMIN.D-1 (degC) '].values[1:]\n",
    "tmax_vals = df['TMAX.D-1 (degC) '].values[1:]\n",
    "met_data['t_min'].values[:, 0, 0] = tmin_vals\n",
    "met_data['t_max'].values[:, 0, 0] = tmax_vals\n",
    "\n",
    "# This is what we have now\n",
    "met_data.to_netcdf('./input/rc_forcing.nc')\n",
    "met_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Put together the required domain "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We form the domain in a similar fashion\n",
    "# First, by creating the data structure\n",
    "coords = {'lat': lats, 'lon': lons}\n",
    "domain = xr.Dataset(coords=coords)\n",
    "domain['elev'] = xr.DataArray(data=np.full((1,1,), np.nan),\n",
    "                          coords=coords,\n",
    "                          dims=('lat', 'lon', ))\n",
    "domain['mask'] = xr.DataArray(data=np.full((1,1,), np.nan),\n",
    "                          coords=coords,\n",
    "                          dims=('lat', 'lon', ))\n",
    "\n",
    "# Add the data\n",
    "domain['elev'][0, 0] = elev\n",
    "domain['mask'][0, 0] = 1\n",
    "domain.to_netcdf('./input/rc_domain.nc')\n",
    "domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Put together the required state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally, we create the state file - the dates are 90 days prior to \n",
    "# the MetSim run dates - as usual, create an empty data structure to\n",
    "# read the data into\n",
    "dates = pd.date_range('10/3/2009', '12/31/2009')\n",
    "shape = (len(dates), 1, 1, )\n",
    "dims = ('time', 'lat', 'lon', )\n",
    "lats = [43.17]\n",
    "lons = [-116.51]\n",
    "elev = 1706.90 # meters\n",
    "coords = {'time': dates, 'lat': lats, 'lon': lons}\n",
    "state = xr.Dataset(coords=coords)\n",
    "for varname in ['prec', 't_min', 't_max']:\n",
    "    state[varname] = xr.DataArray(data=np.full(shape, np.nan),\n",
    "                               coords=coords, dims=dims,\n",
    "                               name=varname)\n",
    "\n",
    "# Do precip data\n",
    "df = pd.read_csv(\"./reynolds_creek_data/daily/precip_2009.csv\", skiprows=[0,1])\n",
    "prec_vals = df['PREC.I-1 (in) '].diff().values[-90:] * 25.4\n",
    "state['prec'].values[:, 0, 0] = prec_vals\n",
    "\n",
    "# And now temp data\n",
    "df = pd.read_csv(\"./reynolds_creek_data/daily/temp_2009.csv\", skiprows=[0,1])\n",
    "tmin_vals = df['TMIN.D-1 (degC) '].values[-90:]\n",
    "tmax_vals = df['TMAX.D-1 (degC) '].values[-90:]\n",
    "state['t_min'].values[:, 0, 0] = tmin_vals\n",
    "state['t_max'].values[:, 0, 0] = tmax_vals\n",
    "state.to_netcdf('./input/rc_state.nc')\n",
    "state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Registering parameters and building the driver\n",
    "Now that we've built all of the input files we need, we can run MetSim with our new setup.\n",
    "Again, we build a simple configuration and run just as we did in the previous notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dates = pd.date_range('1/1/2010', '12/31/2010')\n",
    "params = {\n",
    "    'time_step'    : \"60\",       \n",
    "    'start'        : dates[0],\n",
    "    'stop'         : dates[-1],\n",
    "    'forcing'      : './input/rc_forcing.nc',     \n",
    "    'domain'       : './input/rc_domain.nc',\n",
    "    'state'        : './input/rc_state.nc',\n",
    "    'forcing_fmt'  : 'netcdf',\n",
    "    'out_dir'      : './output',\n",
    "    'out_prefix': 'reynolds',\n",
    "    'scheduler'    : 'threading',\n",
    "    'chunks'       : \n",
    "        {'lat': 1, 'lon': 1},\n",
    "    'forcing_vars' : \n",
    "        {'prec' : 'prec', 't_max': 't_max', 't_min': 't_min'},\n",
    "    'state_vars'   : \n",
    "        {'prec' : 'prec', 't_max': 't_max', 't_min': 't_min'},\n",
    "    'domain_vars'  : \n",
    "        {'elev': 'elev', 'lat': 'lat', 'lon': 'lon', 'mask': 'mask'}\n",
    "    }               \n",
    "\n",
    "ms = MetSim(params)\n",
    "ms.run()\n",
    "output = ms.open_output().load()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's look at what we've got\n",
    "First, let's just look at our output by plotting the shortwave for the year.\n",
    "We can see a clear annual cycle, as we would hope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output['shortwave'].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finally, a simple comparison against some observations\n",
    "While we're at it, let's compare how MetSim simulated shortwave radiation to some observations.\n",
    "We will just look at some for January. \n",
    "First, let's load in the observations and plot the two timeseries.\n",
    "We see that, generally MetSim has the right timing, but is occasionally off on magnitude.\n",
    "To look at this a little bit further we also show a scatter plot comparing these values.\n",
    "Again, we see that MetSim and the observations are fairly well correlated, although there is a decent amount of spread."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"./reynolds_creek_data/solar_rad_jan_2010.csv\", skiprows=[0,1])\n",
    "df.head()\n",
    "df.index = pd.DatetimeIndex(df['Date'] + \" \" + df['Time'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sliced = output.sel(time=slice('01/01/2010', '01/31/2010')).isel(lat=0, lon=0, drop=True)\n",
    "sliced['shortwave'].plot()\n",
    "df['SRADV.H-1 (watt) '].plot(marker = 'o', ls=':')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "ax.scatter(df['SRADV.H-1 (watt) '], sliced['shortwave'].sel(time=df.index))\n",
    "ax.set_xlabel(r'Observed shortwave $(W/m^2)$')\n",
    "ax.set_xlim([0,500])\n",
    "ax.set_ylabel(r'MetSim shortwave $(W/m^2)$')\n",
    "ax.set_ylim([0,500])\n",
    "ax.plot([0, 500], [0, 500], color='k', linestyle='--')"
   ]
  }
 ],
 "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}