LAB-1-Operational-SWE-Products-empty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Operational Snow Water Equivalent (SWE) Products \n",
    "In this lab we will explore two operational SWE products to estimate basin wide SWE. These are standard products, produced in near real time (NRT), that are publically avalible. A word of caution: as you might remember from the lectures, estimating SWE in the mountains is hard. Operational SWE products perform better in large, flat areas with consisten depth and densisty of snowpack across many square kilometeres. Fortunatly, ongoing development of these products means they are continuously improving. While these products might not produce completely accurate basin scale SWE, they can often be used as an index along with data from previous years. \n",
    "\n",
    "## Two types of SWE products: Measurements and Models\n",
    "There are two main catagories of operational SWE products: remotely sensed measurements of SWE and process based model estimates of SWE. We will explore an example of both. For measurement of SWE, we will look at data from the unified Advanced Microwave Scanning Radiometer (AMSR) SWE products which provide daily estimates of global SWE. Note that the product may not be available in some mountainous regions. For process based Model estimates of SWE we will explore the newest near real-time (NRT) ECWMF Re-analysis (ERA-5) SWE estimates. \n",
    "\n",
    "Goals of the lab:\n",
    "- Gain familiarity with SWE plots\n",
    "- Understand the two main sources of SWE estimates: remotely sensed measurements and process based models\n",
    "- Gain familiarty with interpolation and smoothing of data\n",
    "- Data science techniques to visualize the SWE in the context of historical records"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PART 1: AMSR-Unified SWE estimates\n",
    "\n",
    "The first product we are going to look at is the AMSR Unifed SWE product. These data are collected from satellites that measure the emission of microwave radiation from Earth. Snow scatters and absorbs the microwave radiation.  Based on the measured reduction in passive microwave signal, we can estimate the snow water equivalent. This method works best in areas that have shallow snowpack, are devoid of forests, and are flat for at least 25 kilometers.\n",
    "\n",
    "This AMSR-E Unified SWE product is produced using data collected by the Advanced Microwave Scanning Radiometer (AMSR-E) on board the Aqua satellite from June 2002 to October 2011, and from the Advanced Microwave Scanning Radiometer 2 (AMSR2) on board the JAXA GCOM-W1 satellite from July 2012 to the present. The purpose of the AMSR Unified Data Set is to provide the science community with intercalibrated climate products from both the AMSR-E and the AMSR2 instruments. \n",
    "\n",
    "The data sets are archived here:\n",
    "\n",
    "AMSR-E SWE: 2002-2011\n",
    "https://nsidc.org/data/AE_DySno/versions/2\n",
    "\n",
    "\n",
    "AMSR2 SWE: 2012-now \n",
    "https://nsidc.org/data/AU_DySno/versions/1\n",
    "\n",
    "\n",
    "These data are avalible as daily, 5-day maximum, and monthly maximum datasets. We will look at the daily data in this lab exercise. \n",
    "\n",
    "## Set up the Lab workspace\n",
    "To explore these datasets we need to import the required packages for geospatial analysis in python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd \n",
    "import rasterio as rio\n",
    "from rasterio.plot import show\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.image as mpimg\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import glob\n",
    "import contextily as ctx\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choose an Area of Interest\n",
    "We have three different AOIs. Choose one of the following to complete your lab: Langtang catchment in Nepal, Panjsher basin in Afganistan, or wangchu watershed in Bhutan. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list all the files in our directory to see what is avalible and what we are working with in the lab\n",
    "\n",
    "#what about just shapefiles or just tiff files? Can we get a list of all the shapefiles or\n",
    "#just the tiff files we have to work with?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lets look at the Area of Interest (AOI) on a map. Pick one of the following:\n",
    "\n",
    "#choose on of the following areas for your labratory. \n",
    "#Feel free to explore each of them!\n",
    "#panjsher_basin_afghanistan\n",
    "#wangchu_watershed_bhutan\n",
    "#langtang_catchment_nepal\n",
    "############## INSERT AOI WATERSHED NAME HERE ###################\n",
    "watershedName='panjsher_basin_afghanistan'\n",
    "#################################################################\n",
    "aoiFile=watershedName+'_wgs84.shp'\n",
    "basinMaskFile=watershedName+'_basinMask.tif'\n",
    "aoi_AMSR_SWE=watershedName+'_AMSR_SWE.tif'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the Area of Interest on a basemap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#load the shapefile with geopandas and display on a basemap with contextily\n",
    "#The 'EPSG' codes are standerdized information that discribe which projection the datasets \n",
    "#are in. We want to add basemap in same projection as the SWE data and AOI data\n",
    "\n",
    "#load the shapefile into python\n",
    "\n",
    "\n",
    "#what is the projection of this shapefile?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the shapefile over a basemap (in the same projection)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the AMSR-E SWE data\n",
    "\n",
    "First, it is importnat to see how the data are packaged. The datasets are packaged for the entire northern hemisphere as a single image each day in the daily dataset. This is not very useful for looking a small basin over many days or months! Lets see what the full dataset look like for one day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load one day of AMSR data 'fullAMSR_SWE_20100218.png' and print an image to the screen\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "interesting? Do you recognize any land features on earth? Where are we looking at earth from? This is a unique projection. More information is here:\n",
    "https://nsidc.org/ease/ease-grid-projection-gt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore AMSR SWE for our Area of Interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have reprojected, cropped, and assembled data cubes of the full AMSR Unifed SWE records for this lab for each area of interest. \n",
    "\n",
    "A data cube is a stack of data, for Earth science datasets, one common type of datacube is a \"time/space\" cube. This is when we have data for a specific area, for many days, all assembled into a single file for analysis. For this lab we will be working with time/space datacubs of AMSR SWE data. \n",
    "\n",
    "To start, we will load the datacube, and visualize a few days to get a sense of the SWE data that we have to work with. \n",
    "## Load SWE datacube and check size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the basin mask and see how big it is\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# is it the correct size? if not, fix it. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the basin mask\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#load the AMSR SWE dataset\n",
    "\n",
    "\n",
    "#see how big the dataset is - what does each dimension represent? \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have loaded the SWE datacube for our study area. the ***'.shape'*** command lets us see the size of our datacube. From this, can you figure out how many days of data are in the datacube? what about how many measurements per day?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot SWE at one location for entire data record\n",
    "Next, we will load up the dates for these data, which are stored seperatly and plot one of the pixels in the datacube for the whole record of data. We need to find a location to plot the SWE at. Looking at the basemap above, and using your knowledge of the basin, we can take latitude and longitude coordinates and with them figure out which pixel in the image represents the SWE at that location. \n",
    "### Pick a location to map SWE and convert back and forth between lat/lon and pixel coordiantes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hint(70,35.75) is within our AOI\n",
    "\n",
    "#longitude of point of interest\n",
    "\n",
    "#latitude of point of interest\n",
    "\n",
    "#convert longitude and latitude to pixel 'row/column' cooridantes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This transformation can also go in the other direction, we can take a pixel coordinate and convert it to latitude and longitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert backwards from row/colum to latitude/longitude\n",
    "\n",
    "#pixel row\n",
    "\n",
    "#pixel column\n",
    "\n",
    "#convert pixel row and column to latitude and longitude\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#load the dates for the dataset - look at the csv and load it into a dataframe with a column name\n",
    "dateFile=watershedName+'AMSR_sweDates.csv'\n",
    "\n",
    "#here we extract the timeseries of SWE data from the pixel that is at the \n",
    "#latitude and longitude of our point of interest\n",
    "\n",
    "#assign the data from our lat/lon point to a 'SWE_mm' column in the dataframe with measurement dates\n",
    "\n",
    "#set the dates as the index for the dataframe\n",
    "\n",
    "#values above 240 are various flags for no data, we will set to NaN\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the data organized into a DataFrame, lets plot one of the pixels in the datacube for the whole record of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create figure and plot space\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "# Add x-axis and y-axis\n",
    "ax.plot(pointSWEdf.index.values,\n",
    "        pointSWEdf['SWE_mm'],\n",
    "        color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "ax.set(xlabel=\"Date\",\n",
    "       ylabel=\"Snow Water Equivalent (mm)\",\n",
    "       title=\"Daily Snow Water Equivalent \\n\" + watershedName)\n",
    "\n",
    "# Rotate tick marks on x-axis\n",
    "plt.setp(ax.get_xticklabels(), rotation=45)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explore one water year of data\n",
    "Often we would like to see how snow accumulates and melts over a year. Lets look closey at one water year of data, the 2004-2005 water year which begins October 1, 2004 and ends September 30, 2005.\n",
    "### Pick a water year to explore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#pick a water year and pull all the data for our point of interest into a new variable for just that water year\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot SWE for the water year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create figure and plot space\n",
    "\n",
    "# Add data\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "anything look weird? Why does the SWE disapear on lots of days within the data series?\n",
    "\n",
    "\n",
    "This is because the satellite does not take a measurement everyday. Sometimes the pixel is not imaged by the satellite, so there are fill values in the datasets on those days. What can we do about this? One solution is to interpolate the value based on measurements on the preceeding days and following days that are valid measurements. This task is very common with satellite snow observations. Many times there are numerious days where we just do not have a good measurement. This can be becuase of cloud cover, no image was taken on that day, or the image is of poor quality for a variety of reasons. (off nadir satellite view angle, noise, etc....)\n",
    "\n",
    "## Fill and smooth AMSR SWE datacube\n",
    "Now we will smooth the data to remove inconsistencies. This is an example of a simple way to smooth the data. In the lectures, Karl will show you some of our more advanced routines and we will see the quality of the research grade SWE products that include more difficult time/space smoothing. For now, we will explore how smoothing works with some simpler methods. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#calculate the cubic spline estimate of the AMSR data \n",
    "\n",
    "# 7 day  moving average of the interpolated data in the cubic spline estimate\n",
    "\n",
    "#plot the smoothed output.\n",
    "# Create figure and plot space\n",
    "\n",
    "# Add data points and line color\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the maximum SWE value at the point durring the water year?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find the maximum SWE value from the water year\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare original SWE with interpolated and smoothed SWE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define plot space\n",
    "\n",
    "#plot original SWE data as a solid green line\n",
    "\n",
    "\n",
    "#plot interpolated and smoothed SWE data as a dashed black line\n",
    "\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "\n",
    "#set axis limits\n",
    "\n",
    "\n",
    "# Rotate tick marks on x-axis\n",
    "\n",
    "\n",
    "#let python choose where to put the legend\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate volume of SWE in the basin each day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#write a function to calculate the volume of SWE wihtin the basin each day\n",
    "\n",
    "\n",
    "#call the function compute_volume and use the dataCube and basinMask to figure out basin SWE.\n",
    "    \n",
    "    #find all the pixels within the basin\n",
    "    \n",
    "    #initalize output as the right size and datatype\n",
    "  \n",
    "    #calculate basin SWE volume each day and store the answer\n",
    "   \n",
    "        \n",
    "        #remeber that values above 240 are fill in the original dataset, so we dont want to add - assume missing data\n",
    "        daySWE[daySWE>240]=0\n",
    "        \n",
    "        #volume of water in each pixel\n",
    "        #SWE*area=volumne of water - we want volume in cubic kilometers\n",
    "        swe2km3=(31**2)/1000000\n",
    "\n",
    "        #volume of water in the basin\n",
    "        #all data for this lab has been reporjected to 31km pixels equivalent to the ERA-5 grid.\n",
    "        #AMSR-U SWE is in 25km native resolution \n",
    "        #calculate basin volume for the day and save in output\n",
    "        \n",
    "        \n",
    "        #finish the fucntion and return the output\n",
    "\n",
    "        \n",
    "#use the funciton to calculate basin volume and see the size of our answer\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets make a new dataframe with the new basin-wide SWE volume data included."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load AMSR SWE dates and set date/time as index\n",
    "dateFile=watershedName+'AMSR_sweDates.csv'\n",
    "basinSWEdf=pd.read_csv(dateFile,header=None,names=['sweDates'])\n",
    "basinSWEdf['sweDates']=pd.to_datetime(basinSWEdf['sweDates'],format='%Y%m%d')\n",
    "basinSWEdf=basinSWEdf.set_index('sweDates')\n",
    "\n",
    "#Add basin wide SWE to the dataframe\n",
    "\n",
    "#look at the head of the dataset to gut check it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create figure and plot space\n",
    "\n",
    "# Add basinwide SWE data to plot as a purple line\n",
    "\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets look at basin-wide water volume for the 2004-2005 water year again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# create a data frame for the 2004-2005 water year\n",
    "\n",
    "# Create figure and plot space\n",
    "\n",
    "# plot the water year data with a purple line\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we can really see those interpolated values in our data which have replaced the fill values. Lets smooth these also and see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "#calculate the cubic spline estimate of the basin wide AMSR data \n",
    "\n",
    "# 7 day  moving average of the interpolated data\n",
    "\n",
    "#plot the smoothed output\n",
    "# Create figure and plot space\n",
    "\n",
    "# plot the smooth data (after the cubic spline or after the rolling average)\n",
    "\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "# Rotate tick marks on x-axis\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare original basin wide SWE with interpolated and smoothed basin wide SWE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define plot space\n",
    "\n",
    "#plot original basin wide SWE data as green\n",
    "\n",
    "#plot interpolated and smoothed SWE data as a dashed black line\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "#set axis limits\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to remeber:\n",
    "SWE is not streamflow. Many physical processes are occuring in the environment for example, sublimation, evaporation, transpiration, recharging groundwater, and human's consumptive use.\n",
    "These are SWE estimates. Operational products are not highly accurate over mountainous terrain.\n",
    "However, the tools and techniques that we can use to work with the data and get information from it (i.e SWE to water volume) are still useful in understanding the relative amounts of SWE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# PART 2: ERA-5 SWE\n",
    "Lets look at the ERA-5 data now. \n",
    "\n",
    "ERA5 is a process based physical model of the earth. Equations that model the evolution of the earth system are combined with measurements (when avalible) to provide the best possible estimate of the state of the earth at each time step. \n",
    "https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5\n",
    "\n",
    "ERA5 provides hourly estimates of many meteorology and hydrology parameters going back to 1979. \n",
    "The volume of data form ERA5 is huge, so we have ordered the data nessesary for the lab from the online data repository:\n",
    "https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview  \n",
    "The webpage above shows the full list of variables that are avalible from ERA5 and another tab for downloading the data locally.\n",
    "\n",
    "We want to look at the snow water equivalent, which is stored as parameter 141:  \n",
    "https://apps.ecmwf.int/codes/grib/param-db?id=141\n",
    "\n",
    "Now we will perfrom some of the same analysis as above on ERA5 data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load ERA-5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load ERA5 SWE as a data frame\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "data frame of era5 + dates for one pixel (same pixel as the AMSR part of the lab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Load dates \n",
    "\n",
    "#load ERA 5 SWE data for same Lat/Lon as part 1\n",
    "\n",
    "#set index for the data frame as the date each day\n",
    "\n",
    "\n",
    "#display the first 5 values in the data frame\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create figure and plot space\n",
    "\n",
    "# plot all the ERA 5 data\n",
    "\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look Different? ERA5 gives a much longer historical record of the SWE for the location. This is a nice feture of the models and a dangerous one. They can calculate values for times when they do not have much data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trim the dataset to full water years of data only\n",
    "We can see from the above graph that the dataset starts midway through a water year and ends midway through a water year. To run statistical analysis on full water years of data, lets trim the dataset to full water years of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#trim the dataset to only cover full water years\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#replot the data with only complete water years of data\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets look at the 2004-2005 water year like we did with AMSR. Again, the water year starts October 1, 2004 and ends September 30, 2005."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a water year and exract the lat/lon point data from that location for the water year\n",
    "\n",
    "\n",
    "# Create figure and plot space\n",
    "\n",
    "#plot the data\n",
    "\n",
    "# Set title and labels for axes\n",
    "\n",
    "# Rotate tick marks on x-axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see there is no need to interpolate or smooth! An advantage of model outputs relative to satellite observations is that they are continous. The physical equations are used to generate estimates between measurements, this is a form of \"data assimmilation\" where the combination of models and measuremtns enable you to make the best estimate possible, considering all sources of information you have available, at any timestep.\n",
    "\n",
    "With the long history of data available from the model we can contextualize any year within the historical context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find out how many years of data are in the ERA5 dataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Create figure and plot space\n",
    "\n",
    "#plot all the background years as grey lines\n",
    "\n",
    "\n",
    "# Add x-axis and y-axis \n",
    "# hint( x axis is not year this time becuase we have lots of years overlaied )\n",
    "\n",
    "\n",
    "#highlight the 2004-2005 water year as a purple line\n",
    "\n",
    "\n",
    "# Set title and labels for axes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This graph gives some context of the highlighted year we are interested in, in the context of other years. But the graph is still kind of messy. Lets look at another way we can visualize a year of interest in the context of the historical record. What is the driest year? The wettest year? To figure this out, and plot them along with some more historical context, we will use on of the more powerful datascience functions in python the \"GroupBy\" function you can use on pandas dataframes. We can use GroupBy to find the maximum SWE of each calendar year, which we can see from the above graph, for these pixels occurs at peak SWE each spring. [add some more info about group by and some links to good online references]  First, lets use GroupBy to see what the peak SWE was each year for this pixel. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use groupby to split the dataframe into tables for each year, and then find the maximum value in the 'SWE_mm' column\n",
    "\n",
    "#the data starts in october of 1979 so the peak in 1979 is not the peak for the water year\n",
    "#drop the 1979 peak\n",
    "\n",
    "#view results\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a list of the peak SWE each year, we can use this information to identify the water year that had the largest snowpack and the water year that had the smallest snowpack for our contextualization graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#what was the driest year?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#what was the wettest year? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save these two years in a list and print to screen\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#how different are the wettest and the driest year? Plot both of them on the same graph\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know the wettest and driest years, lets use that information to add historical context to our graph of the year of interest. In additon, we can look at the interquartile range of SWE through time. This is the bounds for any given day for which the SWE model is within the range 50% of the time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find the interquartile range (25th percintile and 75th percintile of SWE from ERA-5 for this pixel \n",
    "#for each day of year across the record. (day of water year = DOWY)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the wettest year, the driest year, the interquartile range, and the year of interest on a single graph\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is interesting, dicuss with fellow students how you would interpret this graph. Is the wettest year always the most swe on a given day for a pixel? Is the driest year always the driest? \n",
    "lets now calcuate basin wide storage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PART 3: Challenge Quesiton (optional)\n",
    "Above we made a graph of the ERA-5 data showing the minimum year, maximum year, interquartile range, and current year, from data that spanned a 41 year record. \n",
    "## What is the minimum SWE year and maximum SWE year in the AMSR SWE dataset? \n",
    "Below is some code to get you started, with the AMSR data set you need to answer the question interpolated, smoothed, and plotted for you. To answer the quesiton you will have to modify and add to the codebase below based on what we have used in the lab above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#record of AMSR SWE data for the same point as the ERA 5 data\n",
    "#calculate the cubic spline estimate of the AMSR data \n",
    "amsrSWE_interp_fullRecord=pointSWEdf['SWE_mm'].interpolate(method='cubicspline')\n",
    "\n",
    "# 7 day  moving average of the interpolated data\n",
    "amsrSWE_avg_fullRecord=amsrSWE_interp_fullRecord.rolling(7).mean()\n",
    "\n",
    "#trim to full water years\n",
    "amsrSWE_avg_fullRecord = amsrSWE_avg_fullRecord[~(amsrSWE_avg_fullRecord.index < '2002-10-01')]\n",
    "amsrSWE_avg_fullRecord = amsrSWE_avg_fullRecord[~(amsrSWE_avg_fullRecord.index > '2011-08-30')]\n",
    "\n",
    "#plot the smoothed output\n",
    "# Create figure and plot space\n",
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "# Add x-axis and y-axis\n",
    "ax.plot(amsrSWE_avg_fullRecord.index.values,\n",
    "        amsrSWE_avg_fullRecord,\n",
    "        color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "ax.set(xlabel=\"Date\",\n",
    "       ylabel=\"Snow Water Equivalent (km3)\",\n",
    "       title=\"Daily Snow Water Equivalent \\n\" + watershedName)\n",
    "\n",
    "# Rotate tick marks on x-axis\n",
    "plt.setp(ax.get_xticklabels(), rotation=45)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use the code base from the ERA 5 section of the lab to calculate the wettest and driest year\n",
    "hint: the comments in the empty code below match the comments in the ERA 5 section of the code above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use groupby to split the dataframe into tables for each year,\n",
    "#then find the maximum value in the 'SWE_mm' column\n",
    "\n",
    "#what was the driest year?\n",
    "\n",
    "#what was the wettest year? for this we use .idxmax() instead of .idxmin()\n",
    "\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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}