Lab-3-Decision-Support-Tools-SWE-Reconstruction-and-fSCA-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "imposed-bargain",
   "metadata": {},
   "source": [
    "# LAB 3: Snow Cover Analysis and SWE Reconstruction\n",
    "\n",
    "In this lab we will look at high resolution reserach grade snow cover and SWE products for one of the AOI basins. In additon to the Snow water equivalent, a main theme of these trainings, we will also look at the snow covered area. Snow covered area is far easier to accuratly measure in near real time,can give useful clues to water supplies from seasonal snow, and contextualizes the current year. Our view is that alongisde SWE information,  high quality snow cover maps have a large roll to play in decision support tools for water management in snow domination basins.  \n",
    "  \n",
    "For this lab, we will be looking at 7 years of Snow Covered Area and Snow Water Equivalent data at 500m spatial resolution for the Panjsher basin in Afganistan. The snow covered area data is derived from MODIS satellite data and the Snow Water Equivalent data is generated using SWE reconstruciton of the MODIS snow covered area data combined with energy balance estimats. In lecture # 3 we give a detailed account of how these datasets are created. This lab is focsed on using these datasets and understand how they can be used in decision support.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "swedish-funds",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set up the lab workspace\n",
    "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 copy\n",
    "from datetime import date\n",
    "import contextily as ctx\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "determined-wonder",
   "metadata": {},
   "source": [
    "# MODIS datasets\n",
    "This lab uses data from MODIS, the Moderate Resolution Imageing Spectrometer. MODIS images almost the whole earth everyday. A great way to see an example of what MODIS data look like is to browse NASA worldview:  \n",
    "\n",
    "https://worldview.earthdata.nasa.gov/  \n",
    "\n",
    "\n",
    "Browsing Worldview you can see the swaths of the images MODIS aquires onboard the AQUA and TERRA satellites. Browse and look for snow cover in the mountains each day. Some questions to help you think about estimating snow with optical satellites:  \n",
    "\n",
    "- Can you find a mountainous region you are interested in the snowpack of on Worldview?\n",
    "- Can you see all the snow each day?\n",
    "- Are there clouds in the way often?\n",
    "- How does the cloud cover obsurcing the snow vary by time of year? By location in High Mountain Asia?\n",
    "\n",
    "It is always good to look at the \"raw\" data that we generate our snow products from to improve your understanding of how we get to the final estimates of snow properties. This helps develop your understanding of the strengths and limitations of various apporaches to estiamting SWE, as we have mentioned in prior labs. \n",
    "\n",
    "For SWE and snow cover from MODIS, we have to remember that we are relying on an optical satellite for the two crutial snow measuremnts we need for reconsturciton:\n",
    "\n",
    "- Snow covered area\n",
    "- Snow albedo (via grain size and contaminat concentration) \n",
    "\n",
    "As lecture 3 goes into detail on, these are two critical measurements that we need from MODIS to drive the SWE reconstrution. \n",
    "\n",
    "\n",
    "## MODIS data and the sinusodial projection\n",
    "Another important point about MODIS is that data are in the sinusodial- equal area projection from the MODIS satellite. You will see in this lab, that we have left the data in the standard MODIS projection. The basin shape may look different than you expect it to, that is due to the projection of the data. Best practice is to perform all calculations on the data in the projection that it is delivered in, or in the raw satellite swath level format. Everytime the data is reprojected, there is informaton loss and the quality of the estimates decrease. We stringly recommend that users perform all MODIS snow analysis in the sinusodial projection the data is delivered in, and then reproject final maps for qualitativle visualization only in whichever projection end users prefer.\n",
    "\n",
    "Here are some resources on the MODIS sensor that we will be using:\n",
    "\n",
    "First, there are two MODIS sensors in space, on two different satellites. \n",
    "\n",
    "TERRA (morning MODIS images each day): https://terra.nasa.gov/about/terra-instruments/modis\n",
    "\n",
    "AQUA (afternoon MODIS images each day): https://oceancolor.gsfc.nasa.gov/data/aqua/\n",
    "\n",
    "Aqua suffered a loss of band 6 early on. This band is crutial for estimating the albedo of the snow. The loss of band 6 and the fact that aqua aquires images in the afternoon, when convection in the atmosphere tends to make mountainous regions more cloudy, means that we do not use AQUA MODIS data for our SWE analysis. We instead only use the TERRA datasets as they are higher quality and have less cloud cover issues. \n",
    "\n",
    "General information about MODIS: https://modis.gsfc.nasa.gov/about/    \n",
    "\n",
    "Information on the MODIS Grid: https://modis-land.gsfc.nasa.gov/MODLAND_grid.html  \n",
    "\n",
    "Informaton on the Sinusoidal Projection: https://en.wikipedia.org/wiki/Sinusoidal_projection  \n",
    "\n",
    "For this lab, we have prepared snow cover and snow water equivalent data from MODIS for 2003-2011 to explore how these data can be useful for water supply management. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "celtic-collector",
   "metadata": {},
   "source": [
    "### Visualize the basin we will explore snow cover and SWE within for the lab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "criminal-repair",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#lets look at the Area of Interest (AOI) on a map \n",
    "# to remind us where we will be looking at the data from.\n",
    "aoiFile='panjsher_basin_afghanistan_wgs84.shp'\n",
    "\n",
    "#load the shapefile of the basin boundary\n",
    "AOI = gpd.read_file(aoiFile)\n",
    "AOI.crs= \"EPSG:4326\" #({'init': 'epsg:4326'})\n",
    "\n",
    "#plot - set the plot to show a greater extent around the basin than we saw in Lab 1\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "AOI.plot(edgecolor=\"purple\", facecolor=\"None\", ax=ax, linewidth=4)\n",
    "plt.xlim([67, 72]) #longitude\n",
    "plt.ylim([33, 38]) #latitude\n",
    "ctx.add_basemap(ax, zoom=7, crs=\"EPSG:4326\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "iraqi-opening",
   "metadata": {},
   "source": [
    "# Logical Basin Mask - how we calcualte estimates for a specific basin\n",
    "\n",
    "Logical masks are a power tool for analysizing raster data. Instead of a shapefil like we see above, which is a vector line around the bains, a logical mask is a raster that is the same extent as the underlying remote sensing data, and says exactly which pixels are in or out of the basin, we can use this logical mask to color our visualizaing in useful ways, perform calculatons on the SWE and snow cover data, and focus our anaysis on the basin of interest. Below is an image of the mask we will use for this basin. On the backend, when we do analysis of individual regions or what to constrain our analysis to specific subsets of pixels, we use logical masks to perfrom these selections. \n",
    "\n",
    "## In the code below, how do the numbers represent which pixels are within the basin and which pixels are outside of the basin?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "specialized-talent",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the basin mask\n",
    "basinMaskFile='panjsher_basin_afghanistan_basinMask.tif'\n",
    "bMask=rio.open(basinMaskFile)\n",
    "basinMask = bMask.read()\n",
    "basinMask=np.squeeze(basinMask)\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(basinMask)\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "fig.colorbar(im,ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "australian-advice",
   "metadata": {},
   "source": [
    "# Fractional Snow Covered Area: Why fractions/percent cover? \n",
    "\n",
    "\n",
    "Why do we need an estimate of the fractional cover of snow wihtin each pixel? This is very important for remote sening of seasonal snow in the mountains.   \n",
    "  \n",
    "##### Snow varies at a finer spatial resolution than the size of the pixels in the remote sensing images we have of the snow  \n",
    "  \n",
    "Because of this, a binary map of snow for each pixel (snow/no snow) is insufficent. Becuase the snow can cover just a small portion of a pixel, or a large portion of a pixel, and becuase the percent cover of a pixel can change over the course of many weeks as the snow melts, we need to understand exactly what fraction of the pixel is covered in snow for our snow cover and SWE analysis work. \n",
    "\n",
    "Here is a good paper explaining the prevalence of mixed pixels in the mountains and the need to use fractional snow cover estiamtes: https://doi.org/10.3390/rs61212478\n",
    "\n",
    "### Lets look at Figure 6 and Figure 7 in the above paper to see why we need to estimate the percent cover of snow within each MODIS pixel. \n",
    "\n",
    "Now that we have an understanding of the need for fracitona snow cover estimates, lets look at maps of snow cover an then the corresponding maps of snow water equivalent.   \n",
    "\n",
    "# Maps of Snow Covered Area\n",
    "For this section of the lab we will visualize maps of fractional snow covered area (also refered to as snow cover percent) Let's pick a year of data to explore first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "precious-scale",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a calendar year to explor fractional snow covered area from\n",
    "calendarYear='2011'\n",
    "\n",
    "#file with the fsca maps for the calendar year of interest\n",
    "fname='panjsher_basin_afghanistan_fsca_'+calendarYear+'.tif'\n",
    "\n",
    "#load that year of fractional snow cover data with rasterio\n",
    "#open the file with rasterio\n",
    "src=rio.open(fname)\n",
    "\n",
    "#load the image data into our workspace with rasterio\n",
    "FSCA = src.read()\n",
    "\n",
    "#see how big the dataset is\n",
    "numel=FSCA.shape\n",
    "print(numel)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amended-friendly",
   "metadata": {},
   "source": [
    "# Look at one day of Fsca data\n",
    "\n",
    "To start our exploration of these data, lets visualize one day of fsca data and see what it looks like. An imporant part of visualization is choosing good colors for maping the raster data. matplotlib offers great default colormaps. For future work it is always good to spend some time thinking about how you want to present the data. The colors that are choosen can have a signifcant impact on how users respond to the data. \n",
    "\n",
    "\n",
    "# Choosing a colormap\n",
    "\n",
    "Look at the default options in matplot lib: \n",
    "https://matplotlib.org/stable/tutorials/colors/colormaps.html\n",
    "\n",
    "## QUESTION: What are six general catagoires of colormaps for various types of data? \n",
    "\n",
    "\n",
    "# useful resource for picking colormaps\n",
    "\n",
    "The link below is a great tool to help figure out possible colormaps to use for various datasets:\n",
    "\n",
    "https://colorbrewer2.org\n",
    "\n",
    "Things to consider when choosing how to visualize data:\n",
    "- colorblindness\n",
    "- sequential, diverging, or qualitative values\n",
    "- one huse vs multi hue\n",
    "- does the colormap work when printed in black and white?\n",
    "- contextual information colors and inclusion\n",
    "\n",
    "# Visualize the snow cover data with a good colormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "major-victim",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick a day of year to visualize snow cover percent from \n",
    "doy=11"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "graphic-maldives",
   "metadata": {},
   "source": [
    "set the colormap - we will use the \"Blues\" colormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "brief-batman",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap=copy.copy(plt.get_cmap(\"Blues\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "injured-agreement",
   "metadata": {},
   "source": [
    "reverse the standard blue colormap - we want \"white\" colors to represent high fsca values \n",
    "and \"blue\" colors to represent low to no fsca pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subtle-lobby",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap=cmap.reversed()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vertical-youth",
   "metadata": {},
   "source": [
    "nans are set to the maximum integer value in an 8 bit unsigned integer dataset ( a value of 255) \n",
    "we do not want to color them blue, we want them not included in the colormap, we will set them to grey\n",
    "a different color that contrasts with our choosen colormap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "contained-museum",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap.set_over('grey')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "invalid-niagara",
   "metadata": {},
   "source": [
    "Now we can plot our map!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "further-shelf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(src.read(doy),cmap=cmap,vmax=100)\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='snow cover percent',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "clean-brighton",
   "metadata": {},
   "source": [
    "# Plot timeseries for one pixel for one year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ready-liverpool",
   "metadata": {},
   "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(FSCA[:,100,125],\n",
    "        color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "ax.set(xlabel=\"Day of Year\",\n",
    "       ylabel=\"Snow cover fraction of the pixel\",\n",
    "       title=\"Daily Snow Cover Fraction\")\n",
    "\n",
    "# Rotate tick marks on x-axis\n",
    "plt.setp(ax.get_xticklabels(), rotation=45)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "applicable-projection",
   "metadata": {},
   "source": [
    "# Plot basin wide snow cover for each day of the year\n",
    "\n",
    "We now want to look at the total snow cover of the basin by calculating the basin wide fsca each day. below is our function that performs this calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breeding-pacific",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_basin_fsca(dataCube,basinMask):\n",
    "    idx=np.nonzero(basinMask)\n",
    "    numBasinPxls=np.count_nonzero(basinMask)\n",
    "    numValid=np.count_nonzero(basinMask)\n",
    "    numel=dataCube.shape\n",
    "    output=np.zeros(numel[0])   \n",
    "    for i in range(numel[0]):\n",
    "        day_fsca=np.squeeze(dataCube[i,:,:])\n",
    "        \n",
    "        #remeber that values above 100 are fill, so we dont want to add (snow cover can only go to 100% cover of the pixel)\n",
    "        day_fsca[day_fsca>100]=0\n",
    "        \n",
    "    \n",
    "        #fraction of basin covered by snow each day     \n",
    "        output[i]=day_fsca[idx].sum()/numBasinPxls\n",
    "    return output\n",
    "\n",
    "dailyBasinSnowCover=compute_basin_fsca(FSCA,basinMask)\n",
    "dailyBasinSnowCover.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "respected-affair",
   "metadata": {},
   "source": [
    "Now we can take the output and plot it on a graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "removable-miami",
   "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(dailyBasinSnowCover,\n",
    "        color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "ax.set(xlabel=\"Day of Year\",\n",
    "       ylabel=\"Snow cover fraction (%)\",\n",
    "       title=\"Daily Snow Cover Fraction of the Entire Basin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "invisible-rugby",
   "metadata": {},
   "source": [
    "# Look at basin wide snow cover percent each calendar year\n",
    "\n",
    "What if we want to look at all the years at the same time? This can be very useful for contextualizing snow informaton. The plot below is of all years basin wide snow cover overlaid on top of each other. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bright-subsection",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Create figure and plot space\n",
    "plt.figure(figsize=(15,10));\n",
    "\n",
    "#plot all the background years\n",
    "yrs=list(range(2003,2012))\n",
    "\n",
    "for i in range(2003,2012):\n",
    "  \n",
    "    #load the data from the current year\n",
    "    tmpName='panjsher_basin_afghanistan_fsca_'+str(i)+'.tif'\n",
    "    #load the year of fractional snow cover data\n",
    "    tmpSrc=rio.open(tmpName)\n",
    "    tmpFSCA=tmpSrc.read()\n",
    "    \n",
    "    tmp_dailyBasinSnowCover=compute_basin_fsca(tmpFSCA,basinMask)\n",
    "\n",
    "    \n",
    "# Add x-axis and y-axis - x axis is not year this tim ebecuase we have lots of years overlaied\n",
    "    plt.plot(tmp_dailyBasinSnowCover,\n",
    "            color='grey')\n",
    "\n",
    "#highlight the 2004-2005 water year\n",
    "#plt.plot(list(range(len(pointWY_era5))),\n",
    "#         pointWY_era5,\n",
    "#         color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "titleText='Daily Basin Wide Snow Cover \\n' + 'Panjsher Basin Afghanistan'\n",
    "plt.title(titleText);\n",
    "plt.xlabel(\"Calendar Year Days\");\n",
    "plt.ylabel(\"basin wide snow cover (%)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "illegal-seeker",
   "metadata": {},
   "source": [
    "# Look at basin wide snow cover percent each water year\n",
    "\n",
    "calendar years are not the best frame of view for looking at snow cover and swe, lets rearrange the data to look at it by water year instead. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subject-baking",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create figure and plot space\n",
    "plt.figure(figsize=(15,10));\n",
    "\n",
    "# k will step our index forward\n",
    "# the loop runs by water year, \n",
    "# we want the k index to start at zero and go up by one each loop\n",
    "k=0\n",
    "for i in range(2004,2011):\n",
    "  \n",
    "    #pick a water year and a day of the year to create the snow covered days map for\n",
    "    WaterYear=i\n",
    "\n",
    "    # Specific date\n",
    "    day_of_year_wy_start = date(WaterYear-1, 10, 1).timetuple().tm_yday\n",
    "    day_of_year_wy_end = date(WaterYear, 9, 30).timetuple().tm_yday\n",
    "\n",
    "    #create an empty arrary for snow covered days calculation to go into\n",
    "    fname1='panjsher_basin_afghanistan_fsca_'+str(WaterYear-1)+'.tif'\n",
    "    fname2='panjsher_basin_afghanistan_fsca_'+str(WaterYear)+'.tif'\n",
    "\n",
    "    #load data and create a water year array\n",
    "    #load the year of fractional snow cover data\n",
    "    wy1=rio.open(fname1)\n",
    "    wy2=rio.open(fname2)\n",
    "\n",
    "    #image data\n",
    "    fsca1=wy1.read()\n",
    "    fsca2=wy2.read()\n",
    "\n",
    "    #merge the correct portion of each calendar years dataset to create the water year datacube\n",
    "    #days we dont need from calendar year at start of the water year\n",
    "    idx1=list(range(0,day_of_year_wy_start-1))\n",
    "    #days we dont need from calendar year at end of the water year\n",
    "    idx2=list(range(day_of_year_wy_end,365))\n",
    "\n",
    "    #extract the valid dates from the water year from each calendar year array\n",
    "    fsca1=np.delete(fsca1,idx1,axis=0)\n",
    "    fsca2=np.delete(fsca2,idx2,axis=0)\n",
    "    #concatinate the valid dates into a single water year array\n",
    "    wy_fsca=np.concatenate((fsca1,fsca2))\n",
    "      \n",
    "    tmp_dailyBasinSnowCover=compute_basin_fsca(wy_fsca,basinMask)\n",
    "\n",
    "   \n",
    "    #step index forward - the loop runs by year, we want the k index to start at zero and go up by one each loop\n",
    "    k=k+1\n",
    "\n",
    "# Add x-axis and y-axis - x axis is not year this time becuase we have lots of years overlaied\n",
    "    plt.plot(tmp_dailyBasinSnowCover,\n",
    "            color='grey')\n",
    "\n",
    "# Set title and labels for axes\n",
    "titleText='Daily Basin Wide Snow Cover \\n' + 'Panjsher Basin Afghanistan'\n",
    "plt.title(titleText);\n",
    "plt.xlabel(\"Calendar Year Days\");\n",
    "plt.ylabel(\"basin wide snow cover (%)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spiritual-anaheim",
   "metadata": {},
   "source": [
    "# MODIS Fracitonal Snow Cover Decision Support Tools\n",
    "\n",
    "Now that we have the basin wide and per pixel snow cover through time for each year, we can take this data and turn it into useful information for decision making. Below are a few graphs that help turn snow cover data into actionable information."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exposed-memory",
   "metadata": {},
   "source": [
    "# Snow Covered Days\n",
    "\n",
    "We are going to look at how many snow covered days each pixel has in each water year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "naked-instruction",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a water year and a day of the year to create the snow covered days map for\n",
    "WaterYear=2004\n",
    "\n",
    "# Specific date\n",
    "day_of_year_wy_start = date(WaterYear-1, 10, 1).timetuple().tm_yday\n",
    "day_of_year_wy_end = date(WaterYear, 9, 30).timetuple().tm_yday\n",
    "\n",
    "#create an empty arrary for snow covered days calculation to go into\n",
    "fname1='panjsher_basin_afghanistan_fsca_'+str(WaterYear-1)+'.tif'\n",
    "fname2='panjsher_basin_afghanistan_fsca_'+str(WaterYear)+'.tif'\n",
    "\n",
    "#load data and create a water year array\n",
    "#load the year of fractional snow cover data\n",
    "wy1=rio.open(fname1)\n",
    "wy2=rio.open(fname2)\n",
    "\n",
    "#image data\n",
    "fsca1=wy1.read()\n",
    "fsca2=wy2.read()\n",
    "\n",
    "#merge the correct portion of each calendar years dataset to create the water year datacube\n",
    "#days we dont need from calendar year at start of the water year\n",
    "idx1=list(range(0,day_of_year_wy_start-1))\n",
    "#days we dont need from calendar year at end of the water year\n",
    "idx2=list(range(day_of_year_wy_end,365))\n",
    "\n",
    "#extract the valid dates from the water year from each calendar year array\n",
    "fsca1=np.delete(fsca1,idx1,axis=0)\n",
    "fsca2=np.delete(fsca2,idx2,axis=0)\n",
    "#concatinate the valid dates into a single water year array\n",
    "wy_fsca=np.concatenate((fsca1,fsca2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cheap-adult",
   "metadata": {},
   "source": [
    "Now that we have a water year of daily snow covered area assembled, we can do a snow covered days analysis on the data cube. For snow covered days, we are only interested in how many days the pixel was snow covered this water year. So we need to convert the fractional snow cover data into a new variable that has a binary snow mask that identifies each day a pixel has snow cover."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hollywood-wrapping",
   "metadata": {},
   "outputs": [],
   "source": [
    "#snow cover logical mask is any pixel where snow cover is greater than 0%, and not any of the fill pixels (fsca>100)\n",
    "snow_cover=(wy_fsca>0)&(wy_fsca<=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "neither-french",
   "metadata": {},
   "source": [
    "Now that we have this mask, we can look at it. This is similar to what a buinary snow cover product would look like. While it is not as useful as fractional snow cover for a daily view, we can add up the total number of days of snow cover and gain useful informaton about the basin. Another important point is that binary snow maps usually cannot reliabaly identify snow cover for pixels that are less than 50% snow covered so our snow cover map from fractional snow cover data is more extensive than a traditional binary snow cover map becuase the fractional snow cover can be detected when the snow cover is as low at 10% cover for a pixel. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alpha-collapse",
   "metadata": {},
   "outputs": [],
   "source": [
    "#View one days of snow cover\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(snow_cover[10,:,:])\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "turned-router",
   "metadata": {},
   "source": [
    "### Compare the binary map above to the fractional snow cover map earlier in the lab, see how much more information is in the fractional snow cover map.  \n",
    "\n",
    "# Plot Snow Covered Days\n",
    "Now that we see what snow cover looks like, we can add up how many days so far in the water year each pixel has had snow cover. To do this, we need to pick a day of the year that we want to see the snow covered days estimate for. \n",
    "## Convert a date to a day of water year date\n",
    "To do this we will take a date ( month and day of month) and convert it to a dayof water year date. This is the number of the day within the water year that we want to see the snow covered days estimate for. We use water year for this analsyis becuase we want to capture the date the pixel became snow covered - which usually occures prior to January first for many pixels, and is also usually some day after the start of the water year. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "polar-position",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a date we want to see snow covered days for\n",
    "#month\n",
    "m=4\n",
    "#day of month\n",
    "dom=25\n",
    "#number of days from prior calendar year in the water year array\n",
    "numel=fsca1.shape\n",
    "doy=date(WaterYear, m, dom).timetuple().tm_yday\n",
    "#day in the water year array we want to look at per pixel snow covered days\n",
    "analysisDay=numel[0]+doy-1\n",
    "print('the calendar date '+ str(m) + '-' + str(dom) + ' is the water year day: ' + str(analysisDay))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fluid-expense",
   "metadata": {},
   "source": [
    "now that we have picked our analysis day and have the index for this day within the water year we can add up the total number of snow covered days for each pixel so far in the water year leading up to the analysis day and visualize the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indie-isolation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill the array with the snow covered days for each pixel\n",
    "snow_covered_days=np.count_nonzero(snow_cover[0:analysisDay-1,:,:], axis=0, keepdims=True)\n",
    "\n",
    "#View one days of snow cover\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(np.squeeze(snow_covered_days, axis=0))\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='snow covered days',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prime-polish",
   "metadata": {},
   "source": [
    "# Snow Covered Days - compared to the mean snow covered days\n",
    "\n",
    "Now that we have looked at the snow covered days for one year of data, another useful way to analyize the snow covered days is to look at the relationship between the current year snow covered days and the mean number of snow covered days for that pixel for the period of record. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bridal-zealand",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an array that is the number of snow covered days for all years for the date of interest. \n",
    "#pick a date we want to see snow covered days for:\n",
    "#month\n",
    "m=4\n",
    "#day of month\n",
    "dom=25\n",
    "\n",
    "#empty array for 7 years of snow covered day data for the size of our image\n",
    "snow_covered_days_record=np.empty([7,196,169])\n",
    "snow_covered_days_record.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "medieval-dancing",
   "metadata": {},
   "source": [
    "A common way to loop through data and perfrom calculations is to first \"initalize\" the variable that will eventually hold the answers to your calculations. We can see from above that we have initalized an empty array that can hold 7 years of snow covered days analysis for our date of interest and that the empty arrary has the same dimesions as the daily snow cover data. Now we will loop through all 7 water years and calculate the number of snow covered days for each pixel each water year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "weekly-merchant",
   "metadata": {},
   "outputs": [],
   "source": [
    "# k will step our index forward\n",
    "# the loop runs by water year, \n",
    "# we want the k index to start at zero and go up by one each loop\n",
    "k=0\n",
    "for i in range(2004,2011):\n",
    "  \n",
    "    #pick a water year and a day of the year to create the snow covered days map for\n",
    "    WaterYear=i\n",
    "\n",
    "    # Specific date\n",
    "    day_of_year_wy_start = date(WaterYear-1, 10, 1).timetuple().tm_yday\n",
    "    day_of_year_wy_end = date(WaterYear, 9, 30).timetuple().tm_yday\n",
    "\n",
    "    #create an empty arrary for snow covered days calculation to go into\n",
    "    fname1='panjsher_basin_afghanistan_fsca_'+str(WaterYear-1)+'.tif'\n",
    "    fname2='panjsher_basin_afghanistan_fsca_'+str(WaterYear)+'.tif'\n",
    "\n",
    "    #load data and create a water year array\n",
    "    #load the year of fractional snow cover data\n",
    "    wy1=rio.open(fname1)\n",
    "    wy2=rio.open(fname2)\n",
    "\n",
    "    #image data\n",
    "    fsca1=wy1.read()\n",
    "    fsca2=wy2.read()\n",
    "\n",
    "    #merge the correct portion of each calendar years dataset to create the water year datacube\n",
    "    #days we dont need from calendar year at start of the water year\n",
    "    idx1=list(range(0,day_of_year_wy_start-1))\n",
    "    #days we dont need from calendar year at end of the water year\n",
    "    idx2=list(range(day_of_year_wy_end,365))\n",
    "\n",
    "    #extract the valid dates from the water year from each calendar year array\n",
    "    fsca1=np.delete(fsca1,idx1,axis=0)\n",
    "    fsca2=np.delete(fsca2,idx2,axis=0)\n",
    "    #concatinate the valid dates into a single water year array\n",
    "    wy_fsca=np.concatenate((fsca1,fsca2))\n",
    "    \n",
    "    #snow cover logical mask is any pixel where snow cover is greater than 0%, and not any of the fill pixels (fsca>100)\n",
    "    snow_cover=(wy_fsca>0)&(wy_fsca<=100)\n",
    "    \n",
    "    #number of days from prior calendar year in the water year array\n",
    "    doy=date(WaterYear, m, dom).timetuple().tm_yday\n",
    "    #day in the water year array we want to look at per pixel snow covered days\n",
    "    analysisDay=numel[0]+doy-1\n",
    "    \n",
    "    # fill the array with the snow covered days for each pixel\n",
    "    snow_covered_days=np.count_nonzero(snow_cover[0:analysisDay-1,:,:], axis=0, keepdims=True)\n",
    "\n",
    "    #fill our array we initalized before the loop\n",
    "    snow_covered_days_record[k,:,:]=np.squeeze(snow_covered_days, axis=0)\n",
    "    \n",
    "    #step index forward - the loop runs by year, we want the k index to start at zero and go up by one each loop\n",
    "    k=k+1\n",
    "    \n",
    "snow_covered_days_record.shape "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subsequent-replication",
   "metadata": {},
   "source": [
    "After looping through all the data, our array is the same size, but now it is filled with all the data. Now we can view the snow covered days on our day of interest for any water year from our record.  try a few out with the code below.\n",
    "\n",
    "# Question: Do you see any differences in each water year snow covered days analysis?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "colored-connection",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a number 0-6 for which year you want to plot the snow covered days analysis for. \n",
    "year2plot=0\n",
    "\n",
    "#View one days of snow cover\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(snow_covered_days_record[year2plot,:,:])\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='snow covered days',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "colonial-accused",
   "metadata": {},
   "source": [
    "# Snow Covered Days - Comparing to the mean\n",
    "\n",
    "Now we can generate the mean number of snow covered days for each pixel to get the average length of time each pixel is covered by snow from our record. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "integrated-trail",
   "metadata": {},
   "outputs": [],
   "source": [
    "#average number of snow covered days\n",
    "avg_scd=np.mean(snow_covered_days_record,axis=0)\n",
    "\n",
    "#plot the average number of snow covered days\n",
    "#View one days of snow cover\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(avg_scd)\n",
    "\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='mean number of snow covered days',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "computational-dimension",
   "metadata": {},
   "source": [
    "A useful tool is to then compare the current year to the average to understand how snow cover compares to past years. \n",
    "\n",
    "## Question: What type of colormap would be best for this type of comparison?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "outer-camera",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a number 0-6 for which year you want to plot the snow covered days analysis for. \n",
    "year2plot=0\n",
    "\n",
    "# difference between current year and average snow covered days\n",
    "scd_diff=snow_covered_days_record[year2plot,:,:]-avg_scd\n",
    "#set values outside the basin to fill\n",
    "scd_diff[basinMask==0]=999\n",
    "\n",
    "# plot the results\n",
    "\n",
    "#pick a good colormap for visualzing  the data!\n",
    "#set the colormap\n",
    "cmap=copy.copy(plt.get_cmap(\"PiYG\"))\n",
    "#nans are maxint uint8 (255) we do not want to color them blue, we want them not included in the colormap\n",
    "cmap.set_over('grey')\n",
    "\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(scd_diff,cmap=cmap,vmax=40)\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='differnece from mean snow covered days',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "variable-fisher",
   "metadata": {},
   "source": [
    "# Part 2: SWE Reconstruciton \n",
    "\n",
    "We just looked at snow covered area - which we have avalible in real time and helps us understand some of the differences between the current year and past years. Now we will look at SWE reconstruction, which givens us detailed data on past years water stored in the snowpack, and is a useful tool for training machine learning SWE forecasting models.\n",
    "\n",
    "We will learn how to do the following with SWE reconstructions:\n",
    "\n",
    "- one year of swe\n",
    "- plot for one point\n",
    "- visualize on difffernt days\n",
    "- show the melt curve for each year overlaid\n",
    "\n",
    "#### SWE Units: Raw swe data for this lab has been saved in cm to make the files smaller (uint8 data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "printable-internet",
   "metadata": {},
   "source": [
    "# Maps of Snow Water Equivalent\n",
    "For this section of the lab we will visualize maps of Snow Water Equivalent (SWE) Let's pick a year of data to explore first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "manufactured-forestry",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#pick a calendar year to explore snow water equivalent\n",
    "calendarYear='2011'\n",
    "\n",
    "#file with the fsca maps for the calendar year of interest\n",
    "fname='panjsher_basin_afghanistan_ParBal_SWE_'+calendarYear+'.tif'\n",
    "\n",
    "#load that year of fractional snow cover data with rasterio\n",
    "#open the file with rasterio\n",
    "src=rio.open(fname)\n",
    "\n",
    "#load the image data into our workspace with rasterio\n",
    "SWE = src.read()\n",
    "\n",
    "#see how big the dataset is\n",
    "numel=SWE.shape\n",
    "print(numel)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "infinite-fortune",
   "metadata": {},
   "source": [
    "# Look at one day of SWE data\n",
    "\n",
    "To start our exploration of these data, lets visualize one day of SWE data and see what it looks like. An imporant part of visualization is choosing good colors for maping the raster data. Refer back to the first section of the lab on colormaps to think about what a good colormap would be for the SWE dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "forced-malaysia",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick a day of year to visualize snow cover percent from \n",
    "doy=100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "invalid-treatment",
   "metadata": {},
   "source": [
    "set the colormap - we will use the \"viridis\" perceptually uniform sequential colormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eligible-australia",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap=copy.copy(plt.get_cmap(\"viridis\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "listed-universal",
   "metadata": {},
   "source": [
    "nans are set to the maximum integer value in an 8 bit unsigned integer dataset ( a value of 255) \n",
    "we do not want to color them blue, we want them not included in the colormap, we will set them to grey\n",
    "a different color that contrasts with our choosen colormap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "transparent-hamilton",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap.set_under('grey')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "economic-mistake",
   "metadata": {},
   "source": [
    "Now we can plot our map!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "powered-enforcement",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the image to plot and set fill values to zero\n",
    "plotSWE=src.read(doy)\n",
    "plotSWE[plotSWE==255]=0\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(plotSWE,cmap=cmap,vmin=1)\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='snow water equivalent (cm)',size=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "smoking-focus",
   "metadata": {},
   "outputs": [],
   "source": [
    "pointSWE=SWE[:,125,20]\n",
    "pointSWE[pointSWE==255]=0\n",
    "\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(pointSWE,\n",
    "        color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "ax.set(xlabel=\"Day of Year\",\n",
    "       ylabel=\"Snow water equivalent of the pixel (cm)\",\n",
    "       title=\"Daily Snow water equivalent\")\n",
    "\n",
    "# Rotate tick marks on x-axis\n",
    "plt.setp(ax.get_xticklabels(), rotation=45)\n",
    "plt.xlim(80,160)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "greatest-belgium",
   "metadata": {},
   "source": [
    "# Plot basin wide SWE  for each day of the year\n",
    "\n",
    "We are now going to calculate the total water volume in the snow pack each year in the basin. This is the same calculations we did on AMSR and ERA-5 data in Lab 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rental-stockholm",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_volume(dataCube,basinMask):\n",
    "    idx=np.nonzero(basinMask)\n",
    "    numValid=np.count_nonzero(basinMask)\n",
    "    numel=dataCube.shape\n",
    "    output=np.zeros(numel[0])   \n",
    "    for i in range(numel[0]):\n",
    "        daySWE=np.squeeze(dataCube[i,:,:])\n",
    "        \n",
    "        #remeber that values of 255 are fill in the original dataset, so we dont want to add - assume missing data\n",
    "        daySWE[daySWE==255]=0\n",
    "        \n",
    "        #volume of water in each pixel - calculate in cubic meters and convert to cubic kilometers\n",
    "        #SWE*area=volumne of water - we want volume in cubic kilometers\n",
    "        \n",
    "        #calculate volume in each pixel\n",
    "        #convert cm to km and then multiply depth (in km) by area (in km2) to get vol (km3)\n",
    "        pxlVolSWE=daySWE*1e-5*0.463*0.463\n",
    "        \n",
    "        #sum up all pixel volumnes\n",
    "        output[i]=pxlVolSWE[idx].sum()\n",
    "\n",
    "    return output\n",
    "\n",
    "dailyBasinSWE=compute_volume(SWE,basinMask)\n",
    "dailyBasinSWE.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "canadian-penny",
   "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(dailyBasinSWE,\n",
    "        color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "ax.set(xlabel=\"Day of Year\",\n",
    "       ylabel=\"SWE km3\",\n",
    "       title=\"Daily Snow Water Storage Volumne of the Entire Basin\")\n",
    "plt.xlim(50,200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "consecutive-guinea",
   "metadata": {},
   "source": [
    "# Look at basin wide SWE each calendar year\n",
    "\n",
    "Now we will take the melt curve from above and plot each years melt curve on the same graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "numeric-hours",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create figure and plot space\n",
    "plt.figure(figsize=(15,10));\n",
    "\n",
    "#plot all the background years\n",
    "yrs=list(range(2003,2012))\n",
    "\n",
    "for i in range(2003,2012):\n",
    "  \n",
    "    #load the data from the current year\n",
    "    tmpName='panjsher_basin_afghanistan_ParBal_SWE_'+str(i)+'.tif'\n",
    "    #load the year of fractional snow cover data\n",
    "    tmpSrc=rio.open(tmpName)\n",
    "    tmpSWE=tmpSrc.read()\n",
    "    \n",
    "    tmp_dailyBasinSWE=compute_volume(tmpSWE,basinMask)\n",
    "\n",
    "    \n",
    "# Add x-axis and y-axis - x axis is not year this tim ebecuase we have lots of years overlaied\n",
    "    plt.plot(tmp_dailyBasinSWE,\n",
    "            color='grey')\n",
    "\n",
    "#highlight the 2004-2005 water year\n",
    "#plt.plot(list(range(len(pointWY_era5))),\n",
    "#         pointWY_era5,\n",
    "#         color='purple')\n",
    "\n",
    "# Set title and labels for axes\n",
    "titleText='Daily Basin Wide Snow Storage volume (km3) \\n' + 'Panjsher Basin Afghanistan'\n",
    "plt.title(titleText);\n",
    "plt.xlabel(\"Calendar Year Days\");\n",
    "plt.ylabel(\"water volume (km3)\");\n",
    "plt.xlim(1,250)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "standing-maria",
   "metadata": {},
   "source": [
    "# SWE in a given year compared to the mean SWE\n",
    "\n",
    "Now that we have looked at the snow covered days for one year of data, another useful way to analyize the snow covered days is to look at the relationship between the current year snow covered days and the mean number of snow covered days for that pixel for the period of record. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vulnerable-cincinnati",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an array that is the number of snow covered days for all years for the date of interest. \n",
    "#pick a date we want to see snow covered days for:\n",
    "#month\n",
    "m=4\n",
    "#day of month\n",
    "dom=1\n",
    "\n",
    "#get the day of year of the analysis\n",
    "doy=date(2003, m, dom).timetuple().tm_yday\n",
    "\n",
    "#empty array for 7 years of swe data for the size of our image\n",
    "swe_record=np.empty([9,196,169])\n",
    "swe_record.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "divided-hepatitis",
   "metadata": {},
   "outputs": [],
   "source": [
    "#average SWE for this date\n",
    "\n",
    "#loop through and grab SWE from each year for this date\n",
    "k=0\n",
    "for i in range(2003,2012):\n",
    "    print(i)\n",
    "    #load the data from the current year\n",
    "    tmpName='panjsher_basin_afghanistan_ParBal_SWE_'+str(i)+'.tif'\n",
    "    #load the year of SWE data\n",
    "    tmpSrc=rio.open(tmpName)\n",
    "    tmpSWE=tmpSrc.read(doy)\n",
    "    tmpSWE[tmpSWE==255]=0\n",
    "    \n",
    "    #fill our array we initalized before the loop\n",
    "    swe_record[k,:,:]=tmpSWE\n",
    "\n",
    "    k=k+1\n",
    "    \n",
    "swe_record.shape\n",
    "\n",
    "#calculate the average swe for each pixel on this date, across the record.\n",
    "avg_swe=np.mean(swe_record,axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "civilian-blocking",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the mean swe for this date\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(avg_swe)\n",
    "\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='mean swe (cm)',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "improved-intelligence",
   "metadata": {},
   "source": [
    "A useful tool is to then compare the current year to the average to understand how SWE compares to past years. \n",
    "\n",
    "## Question: What type of colormap would be best for this type of comparison?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "incorrect-trail",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# difference between current year and average snow covered days\n",
    "swe_diff=swe_record[1,:,:]-avg_swe\n",
    "\n",
    "#set values outside the basin to fill\n",
    "swe_diff[basinMask==0]=999\n",
    "\n",
    "# plot the results\n",
    "#pick a good colormap for visualzing  the data!\n",
    "#set the colormap\n",
    "cmap=copy.copy(plt.get_cmap(\"cividis\"))\n",
    "#nans are maxint uint8 (255) we do not want to color them blue, we want them not included in the colormap\n",
    "cmap.set_over('white')\n",
    "\n",
    "#set up the plot\n",
    "fig,ax=plt.subplots(figsize=(15,15))\n",
    "#plot the fsca image and the correct colormap with maxcolor value set to fsca=100\n",
    "im=ax.imshow(swe_diff,cmap=cmap,vmax=50)\n",
    "#plt.show()\n",
    "#add the colorbar\n",
    "cbar=fig.colorbar(im,ax=ax)\n",
    "cbar.set_label(label='differnece from mean SWE (cm)',size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fifteen-philadelphia",
   "metadata": {},
   "source": [
    "## We can also look at all the years at once and see the spatial and temporal variabiliy of SWE across the basin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "downtown-module",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    " #pick a good colormap for visualzing  the data!\n",
    "#set the colormap\n",
    "cmap=copy.copy(plt.get_cmap(\"cividis\"))\n",
    "#nans are maxint uint8 (255) we do not want to color them blue, we want them not included in the colormap\n",
    "cmap.set_over('white')\n",
    "\n",
    "#set up the plot\n",
    "fig,axs=plt.subplots(3,3,figsize=(15,15))\n",
    "\n",
    "#list of the water years\n",
    "WY=list(range(2003,2012))\n",
    "    \n",
    "k=0\n",
    "   \n",
    "for col in range(3):\n",
    "    for row in range(3):\n",
    "        \n",
    "        # difference between current year and average snow covered days\n",
    "        swe_diff=swe_record[k,:,:]-avg_swe\n",
    "\n",
    "        #set values outside the basin to fill\n",
    "        swe_diff[basinMask==0]=999\n",
    "        \n",
    "        ax = axs[row, col]\n",
    "        ax.set_title('Water Year: ' + str(WY[k]))\n",
    "        \n",
    "        #vmax and vmin keep the colorbars the same between year, colors are same values each year\n",
    "        pcm = ax.imshow(swe_diff,cmap=cmap,vmax=50,vmin=-40)\n",
    "        fig.colorbar(pcm, ax=ax)\n",
    "        k=k+1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "naughty-frame",
   "metadata": {},
   "source": [
    "We can do this same sort of annual comparison plot wiht the actuall volume of SWE each year on this date, which is also useful:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "diagnostic-question",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pick a good colormap for visualzing  the data!\n",
    "#set the colormap\n",
    "cmap=copy.copy(plt.get_cmap(\"viridis\"))\n",
    "#nans are maxint uint8 (255) we do not want to color them blue, we want them not included in the colormap\n",
    "cmap.set_over('white')\n",
    "\n",
    "#set up the plot\n",
    "fig,axs=plt.subplots(3,3,figsize=(15,15))\n",
    "\n",
    "#list of the water years\n",
    "WY=list(range(2003,2012))\n",
    "    \n",
    "k=0\n",
    "   \n",
    "for col in range(3):\n",
    "    for row in range(3):\n",
    "        \n",
    "        # difference between current year and average snow covered days\n",
    "        yr_swe=swe_record[k,:,:]\n",
    "\n",
    "        #set values outside the basin to fill\n",
    "        yr_swe[basinMask==0]=999\n",
    "        \n",
    "        ax = axs[row, col]\n",
    "        ax.set_title('Water Year: ' + str(WY[k]))\n",
    "        \n",
    "        #vmax and vmin keep the colorbars the same between year, colors are same values each year\n",
    "        pcm = ax.imshow(yr_swe,cmap=cmap,vmax=100)\n",
    "        fig.colorbar(pcm, ax=ax)\n",
    "        k=k+1\n",
    "plt.show()"
   ]
  }
 ],
 "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": 5
}