SNODAS Wrapper.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SNODAS Wrapper\n",
    "\n",
    "This code was written by Jonah Joughin with modifications by A. Arendt and E. Mayorga\n",
    "\n",
    "This code is decomissioned for now. Holding on to it in case we return to the downloading of all snodas data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import re\n",
    "import os\n",
    "import math\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import dask\n",
    "from datetime import datetime, timedelta\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from validation import SNODAS, Elevation, utils as ut, creds\n",
    "import mtnhubsnow as mh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Downloading Data from SNODAS\n",
    "Data can be fetched from SNODAS using the `snodas_ds(date)` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetch data from SNODAS\n",
    "start_date = datetime(2019,11,1)\n",
    "for date in (start_date + timedelta(n) for n in range(31)):\n",
    "    output_path = date.strftime('data/SNODAS/SNODAS_%Y%m%d.nc')\n",
    "    print(output_path)\n",
    "    if not os.path.exists(output_path):\n",
    "        snodas_ds = SNODAS.snodas_ds(date)\n",
    "        ut.save_netcdf(snodas_ds, output_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open a single dataset\n",
    "ds = xr.open_dataset('data/SNODAS/SNODAS_20191201.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merging files into single NetCDF files\n",
    "Merging together NetCDF files for different dates can be accomplished by reading those files into xarray, and then writing the resulting dataset out to a new NetCDF file. It is important to set the proper coordinates along the time axis of the dataset before writing, in order to record the date of each individual layer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_path = 'data/SNODAS/SNODAS_2019_sample.nc'\n",
    "# Only process if file does not already exist\n",
    "if not os.path.exists(output_path):\n",
    "    # Get list of valid files\n",
    "    pattern = re.compile(\"^SNODAS_\\d{8}.nc$\")\n",
    "    # There's something wrong with SNODAS-20170913 - only use files after this for now\n",
    "    files = sorted([os.path.join(\"data/SNODAS\",f) for f in os.listdir(\"data/SNODAS\") if pattern.match(f)])\n",
    "    # Extract dates from files\n",
    "    dates = [ut.date_from_file(f) for f in files]\n",
    "\n",
    "    ds = xr.open_mfdataset(files, concat_dim='time', combine='nested')\n",
    "    ds.coords['time'] = dates\n",
    "\n",
    "# Serializing is causing my local machine to crash - we should explore exporting to S3?\n",
    "    #ds.to_netcdf(output_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MISSING MASKED FILES\n",
    "\n",
    "These notes are from David Hill's [shell script](https://github.com/dfosterhill/SNODAS/blob/master/get_proc_snodas.sh)\n",
    "\n",
    "### The following dates are missing ALL data (YYYY-MM-DD):\n",
    "\n",
    "* 2004-02-25\n",
    "* 2004-08-31\n",
    "* 2004-09-27\n",
    "* 2005-06-25\n",
    "* 2005-08-01\n",
    "* 2005-08-02\n",
    "* 2006-08-26\n",
    "* 2006-08-27\n",
    "* 2006-09-08\n",
    "* 2006-09-30\n",
    "* 2006-10-01\n",
    "* 2007-02-14\n",
    "* 2007-03-26\n",
    "* 2008-03-13\n",
    "* 2008-06-13\n",
    "* 2008-06-18\n",
    "* 2009-08-20\n",
    "* 2012-12-20\n",
    "\n",
    "### The following dates are missing individual files:\n",
    "\n",
    "2003-10-30 is missing one file:\n",
    "\n",
    "* us_ssmv11034tS__T0001TTNATS2003103005HP001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Merged Data in XArray\n",
    "The merged dataset can then be reimported to xarray with the following function. The merging process is important, as it allows xarray to search the dataset much faster, opening up the possibility of API queries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset('data/SNODAS/SNODAS_2019_example.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting snow depth across region\n",
    "Data can be selected in a particular region by using the `sel` and `isel` functions in conjunction with `slice`, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.Band1.sel(time = '2019-11-2', lat = slice(40, 54), lon = slice(-126, -110)).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selecting time series at location\n",
    "A time series can be constructed as follows. It is important that the method parameter is set to `'nearest'` in order to recieve results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region = {\n",
    "    'latmax' : 44.25,\n",
    "    'latmin' : 43.75,\n",
    "    'lonmax': -121.65,\n",
    "    'lonmin': -122.5,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = mh.snow_data(limit=10000000, start=datetime(2019,11,9), bbox=region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Google Elevation API now costs $$ to use\n",
    "\n",
    "# We should depreciate the Elevation.py script\n",
    "\n",
    "#obs = Elevation.merge_el_data(obs)\n",
    "#plot = obs.plot(x='date', y='snow_depth', style='o')\n",
    "#snodas_series = ds.Band1.sel(lat=lat, lon=-122, method='nearest') / 10\n",
    "#snodas_series.plot(ax=plot)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting Data Against SNODAS\n",
    "\n",
    "This is older code that needs to be updated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict points to the continental US\n",
    "region = {\n",
    "    'ymax' : 50,\n",
    "    'ymin' : 25,\n",
    "    'xmax': -65,\n",
    "    'xmin': -125,\n",
    "}\n",
    "\n",
    "obs = MountainHub.snow_data(limit=1000, start=datetime(2017,9,14), end=datetime(2018,3,28), box=region)\n",
    "def snodas_depth(ts, lon, lat):\n",
    "    height = ds.Band1.sel(time = ts.strftime('%Y-%m-%d'), lon=lon, lat=lat, method='nearest').item()\n",
    "    if not np.isnan(height):\n",
    "        height /= 10\n",
    "    return height\n",
    "\n",
    "obs['snodas_depth'] = obs.apply(lambda x: snodas_depth(x['date'], x['long'], x['lat']), axis=1)\n",
    "\n",
    "sns.set(color_codes=True)\n",
    "sns.lmplot(x='snodas_depth',y='snow_depth',data=obs, fit_reg=True, lowess=True) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:CSO-validation]",
   "language": "python",
   "name": "conda-env-CSO-validation-py"
  },
  "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}