notebooks/raiderStats/raiderStats_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Peform statistical analyses of GNSS station locations and tropospheric zenith delays\n",
    "\n",
    "**Author**: Simran Sangha, David Bekaert - Jet Propulsion Laboratory\n",
    "\n",
    "This notebook provides an overview of the functionality included in the **`raiderStats.py`** program. Specifically, we outline examples on how to perform basic statistical analyses of GNSS station location and tropospheric zenith delay information over a user defined area of interest, span of time, and seasonal interval. In this notebook, we query GNSS stations spanning California between 2016 and 2018. \n",
    "\n",
    "We will outline the following statistical analysis and filtering options:\n",
    "- Restrict analyses to range of years\n",
    "- Restrict analyses to range of months (i.e. seasonal interval)\n",
    "- Illustrate station distribution and tropospheric zenith delay mean/standard deviation\n",
    "- Illustrate gridded distribution and tropospheric zenith delay mean/standard deviation\n",
    "- Generate variogram plots across specified time periods\n",
    "- Perform basic seasonal amplitude/phase analyses\n",
    "- Examine residuals between weather-models and collocated GNSS stations\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "Both the initial setup (<b>Prep A</b> section) and download of the data (<b>Prep B</b> section) should be run at the start of the notebook. All other sections (examples and applications) do not need to be run in order.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-danger\">\n",
    "<b>Potential Errors/Warnings:</b> \n",
    "    \n",
    "- RAiDER needs to be installed to run this notebook\n",
    "- Most of the examples below specify the use of all available cpus `--cpus all`. If you anticipate this would be an issue on your machine, please adjust to a more conservative number accordingly (e.g. `--cpus 8`).\n",
    "- Just note that downloading data with less than 8 cpus may be very time-consuming and not recommended, while most other steps should run quickly.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    <b>Terminology:</b>\n",
    "    \n",
    "- *GNSS*: Stands for Global Navigation Satellite System. Describes any satellite constellation providing global or regional positioning, navigation, and timing services.\n",
    "- *tropospheric zenith delay*: The precise atmospheric delay satellite signals experience when propagating through the troposphere.\n",
    "- *variogram*: Characterization of the difference between field values at two locations.\n",
    "- *empirical variogram / semivariance*: Provides a description of how the data are correlated with distance. 1/2 * sqrt(obs1 - obs2)\n",
    "- *experimental variogram*: A discrete function calculated using a measure of variability between pairs of points at various distances.\n",
    "- *sill/variance*: The limit of the variogram representing where the variogram appears to level off, tending to infinity lag distances.\n",
    "- *range*: The distance in which the difference of the variogram from the sill becomes negligible, such that the data are no longer autocorrelated.\n",
    "    </div>\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Table of Contents:\n",
    "<a id='example_TOC'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[**Overview of the raiderStats.py program**](#overview)\n",
    "- [1. Basic user input options](#overview_1)\n",
    "- [2. Run parameters](#overview_2)\n",
    "- [3. Optional controls for spatiotemporal subsetting](#overview_3)\n",
    "- [4. Supported types of individual station scatter-plots](#overview_4)\n",
    "- [5. Supported types of gridded station plots](#overview_5)\n",
    "- [6. Supported types of variogram plots](#overview_6)\n",
    "- [7. Optional controls for plotting](#overview_7)\n",
    "\n",
    "[**Examples of the raiderStats.py program**](#examples)\n",
    "- [Example 1. Generate all individual station scatter-plots, as listed under section #4](#example_1)\n",
    "- [Example 2. Generate all basic gridded station plots, as listed under section #5](#example_2)\n",
    "  - [Example 2a. Redo plots efficiently using generated grid raster files](#example_2a)\n",
    "- [Example 3. Generate gridded mean tropospheric zenith delay plot, with stations superimposed](#example_3)\n",
    "- [Example 4. Generate variogram plots](#example_4)\n",
    "- [Example 5. Generate seasonal phase/amplitude plots](#example_5)\n",
    "- [Example 6. Generate weather model/GNSS residual plots](#example_6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Prep A: Initial setup of the notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below we set up the directory structure for this notebook exercise. In addition, we load the required modules into our python environment using the **`import`** command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "## Defining the home and data directories\n",
    "tutorial_home_dir = os.path.abspath(os.getcwd())\n",
    "work_dir = os.path.abspath(os.getcwd())\n",
    "print(\"Tutorial directory: \", tutorial_home_dir)\n",
    "print(\"Work directory: \", work_dir)\n",
    "\n",
    "# Verifying if RAiDER is installed correctly\n",
    "try:\n",
    "    from RAiDER import statsPlot\n",
    "except:\n",
    "    raise Exception('RAiDER is missing from your PYTHONPATH')\n",
    "\n",
    "os.chdir(work_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Prep B: Download prerequisite GNSS station location and tropospheric zenith delay information with the **`raiderDownloadGNSS.py`** program"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Virtually access GNSS station location and zenith delay information for the years '2016-2018', for every 5 days, at a UTC time of day 'HH:MM:SS' of '00:00:00', and across a geographic bounding box '32 42 -125 -114' spanning over California.\n",
    "\n",
    "The footprint of the specified geographic bounding box is again depicted in **Fig. 1**.\n",
    "\n",
    "In addition to querying for multiple years, we will also experiment with using the maximum number of allowed CPUs to save some time! Recall again that the default number of CPUs used for parallelization is 8.\n",
    "\n",
    "Note these features and similar examples are outlined in more detail in the companion notebook **`raiderDownloadGNSS/raiderDownloadGNSS_tutorial.ipynb`**. The outputs of the last example for that notebook will be used for exercises 1-5. So if this example has already been executed, the following prep sets a symlink pointing to these outputs over from the previous notebook to avoid duplicating efforts and wasting time/coffee."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "## Check if outputs from the relevant raiderDownloadGNSS notebook exercise exist\n",
    "## Only run the raiderDownloadGNSS program if these outputs don't exist\n",
    "gnss_tutorial_home_dir = '../raiderDownloadGNSS/products_ex2/UNRcombinedGPS_ztd.csv'\n",
    "run_command = \"raiderDownloadGNSS.py --out products_ex2 --date 20160101 20181231 5 --returntime '00:00:00' \\\n",
    "              --bounding_box '32 42 -125 -114' --cpus all\"\n",
    "if not os.path.exists(gnss_tutorial_home_dir):\n",
    "    from RAiDER import statsPlot\n",
    "    print('Download data needed for this exercise with the following command: {}'.format(run_command))\n",
    "    os.system(run_command)\n",
    "else:\n",
    "    os.symlink(os.path.dirname(gnss_tutorial_home_dir), 'products')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "All of the extracted tropospheric zenith delay information stored under **`GPS_delays`** is concatenated with the GNSS station location information stored under **`gnssStationList_overbbox.csv`** into a primary comprehensive file **`UNRcombinedGPS_ztd.csv`**\n",
    "\n",
    "**`UNRcombinedGPS_ztd.csv`** may in turn be directly used to perform basic statistical analyses using **`raiderStats.py`**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<img src=\"support_docs/bbox_footprint.png\" alt=\"footprint\" width=\"700\">\n",
    "<center><b>Fig. 1</b> Footprint of geopraphic bounding box used in the examples below. </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Overview of the raiderStats.py program\n",
    "<a id='overview'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The **`raiderStats.py`** program provides a suite of convinient statistical analyses of GNSS station locations and tropospheric zenith delays.\n",
    "\n",
    "Running **`raiderStats.py`** with the **`-h`** option will show the parameter options and outline several basic, practical examples. \n",
    "\n",
    "Let us explore these options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 1. Basic user input options\n",
    "<a id='overview_1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Input CSV file (**`--file FNAME`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "**REQUIRED** argument. Provide a valid CSV file as input through **`--file`** which lists the GNSS station IDs (ID), lat/lon coordinates (Lat,Lon), dates in YYYY-MM-DD format (Date), and the desired data field in units of meters.\n",
    "\n",
    "Note that the complementary **`raiderDownloadGNSS.py`** format generates such a primary CSV file named **`UNRcombinedGPS_ztd.csv`** that contains all such fields and is already formatted as expected by **`raiderStats.py`**. Please refer to the accompanying **`raiderDownloadGNSS/raiderDownloadGNSS_tutorial.ipynb `** for more details and practical examples.\n",
    "\n",
    "Additionally supported are rasters containing gridded arrays generated with a previous `raiderStats.py` call with the `--grid_to_raster` option specified. This is designed to be a quick means to replot previously generated arrays while circumventing the overhead involved with the initial gridding/processing of the original CSV file.\n",
    "E.g. if specified with **`--grid_delay_mean`**, then a raster file named `ZTD_grid_delay_mean.tif` containing the gridded mean delay will be generated\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Data column name (**`--column_name COL_NAME`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify name of data column in input CSV file through **`--column_name `** that you wish to perform statistical analyses on. Input assumed to be in units of meters.\n",
    "\n",
    "Default input column name set to **`ZTD`**, the name assigned to tropospheric zenith delays populated under the **`ZTDCombinedGPS_ztd.csv`** file generated through **`raiderDownloadGNSS.py`**\n",
    "\n",
    "The column name is always prepended to output products (e.g. `ZTD_grid_heatmap.tif` and `ZTD_grid_heatmap.png`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Data column unit (**`--unit UNIT`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify unit for output rasters/graphics through **`--unit`**. Again, input assumed to be in units of meters so it will be converted into meters if not already in meters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 2. Run parameters\n",
    "<a id='overview_2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Output directory (**`--workdir WORKDIR`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify directory to deposit all outputs into with **`--workdir`**. Absolute and relative paths are both supported.\n",
    "\n",
    "By default, outputs will be deposited into the current working directory where the program is launched."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Number of CPUs to be used (**`--cpus NUMCPUS`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify number of cpus to be used for multiprocessing with **`--cpus`**. For most cases, multiprocessing is essential in order to access data and perform statistical analyses within a reasonable amount of time.\n",
    "\n",
    "May specify **`--cpus all`** at your own discretion in order to leverage all available CPUs on your system.\n",
    "\n",
    "By default 8 CPUs will be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Verbose mode (**`--verbose`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify **`--verbose`** to print all statements through entire routine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 3. Optional controls for spatiotemporal subsetting\n",
    "<a id='overview_3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Geographic bounding box (**`--bounding_box BOUNDING_BOX`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "An area of interest may be specified as `SNWE` coordinates using the **`--bounding_box`** option. Coordinates should be specified as a space delimited string surrounded by quotes. The common intersection between the user-specified spatial bounds and the spatial bounds computed from the station locations in the input file is then passed. This example below would restrict the analysis to stations over California:\n",
    "**`--bounding_box '32 42 -125 -114'`**\n",
    "\n",
    "If no area of interest is specified, by default the maximum spatial bounds computed based off of the coordinates in the input file is passed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Gridcell spacing (**`--spacing SPACING`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify degree spacing of grid-cells for statistical analyses through **`--spacing`**\n",
    "\n",
    "By default grid-cell spacing is set to 1°. If the specified grid-cell spacing is not a multiple of the spatial bounds of the dataset, the grid-cell spacing again defaults back to  1°."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Subset in time (**`--timeinterval TIMEINTERVAL`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Define temporal bounds with **`--timeinterval TIMEINTERVAL`** by specifying earliest YYYY-MM-DD date followed by latest date YYYY-MM-DD. For example: **`--timeinterval 2018-01-01 2019-01-01`**\n",
    "\n",
    "By default, bounds set to earliest and latest time found in input file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Seasonal interval (**`--seasonalinterval SEASONALINTERVAL`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Observational error limit (**`--obs_errlimit OBS_ERRLIMIT`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify observation error threshold with **`--obs_errlimit OBS_ERRLIMIT`** to discard observations with large uncertainties, where units reflect specified input **`--unit`** (or default `m`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 4. Supported types of individual station scatter-plots\n",
    "<a id='overview_4'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot station distribution (**`--station_distribution`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate each individual station with black markers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot mean tropospheric zenith delay by station (**`--station_delay_mean`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate the tropospheric zenith delay mean for each station with a **`hot`** colorbar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot standard deviation of tropospheric zenith delay by station  (**`--station_delay_stdev`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate the tropospheric zenith delay standard deviation for each station with a **`hot`** colorbar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot phase/amplitude of tropospheric zenith delay by station  (**`--station_seasonal_phase`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate the phase/amplitude of tropospheric zenith delay for each station with a **`hot`** colorbar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 5. Supported types of gridded station plots\n",
    "<a id='overview_5'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot gridded station heatmap (**`--grid_heatmap`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate heatmap of gridded station array with a **`hot`** colorbar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot gridded mean tropospheric zenith delay (**`--grid_delay_mean`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate gridded mean of station-wise means with a **`hot`** colorbar.\n",
    "Alternatively plot absolute gridded delay mean (i.e. measure biased by stations with higher temporal sampling) with the option `--grid_delay_absolute_mean`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot gridded median tropospheric zenith delay (**`--grid_delay_median`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate gridded mean of station-wise medians with a **`hot`** colorbar.\n",
    "Alternatively plot absolute gridded delay median (i.e. measure biased by stations with higher temporal sampling) with the option `--grid_delay_absolute_median`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot gridded standard deviation of tropospheric zenith delay  (**`--grid_delay_stdev`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate gridded mean of station-wise standard deviations with a **`hot`** colorbar.\n",
    "Alternatively plot absolute gridded delay standard deviation (i.e. measure biased by stations with higher temporal sampling) with the option `--grid_delay_absolute_stdev`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot gridded station-wise delay phase/amplitude  (**`--grid_seasonal_phase`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Illustrate gridded mean of station-wise phase/amplitude/period with a **`hot`** colorbar.\n",
    "Alternatively plot absolute gridded delay phase/amplitude/period (i.e. measure biased by stations with higher temporal sampling) with the option `--grid_seasonal_absolute_phase`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 6. Supported types of variogram plots\n",
    "<a id='overview_6'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot variogram (**`--variogramplot`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Passing **`--variogramplot`** toggles plotting of gridded station variogram, where gridded sill and range values for the experimental variogram fits are illustrated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Apply experimental fit to binned variogram (**`--binnedvariogram`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Pass **`--binnedvariogram`** to apply experimental variogram fit to total binned empirical variograms for each time slice. \n",
    "\n",
    "Default is to pass total unbinned empiricial variogram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Save variogram figures per time-slice (**`--variogram_per_timeslice`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify **`--variogram_per_timeslice`** to generate variogram plots per gridded station AND time-slice.\n",
    "\n",
    "If option not toggled, then variogram plots are only generated per gridded station and spanning entire time-span."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Variogram error limit (**`--variogram_errlimit VARIOGRAM_ERRLIMIT`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify variogram RMSE threshold with **`--variogram_errlimit VARIOGRAM_ERRLIMIT`** to discard grid-cells with large RMSE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### 7. Optional controls for plotting\n",
    "<a id='overview_7'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Save gridded array(s) as raster(s) (**`--grid_to_raster`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Save specified gridded array(s) as raster(s). \n",
    "May directly load/plot in successive script call by passing output grid as argument for **`--file`**.\n",
    "E.g. if specified with **`--grid_delay_mean`**, then a raster file named `ZTD_grid_delay_mean.tif` containing the gridded mean delay will be generated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Save debug figures of station-wise seasonal fit (**`--phaseamp_per_station`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Save debug figures of curve-fit vs data per station for seasonal amplitude/phase analaysis options (i.e. **`--grid_seasonal_phase`** and/or **`--station_seasonal_phase`**)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Minimum TS span and minimum fractional observations for seasonal amplitude/phase analyses (**`--min_span`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Minimum TS span (years) and minimum fractional observations in span (fraction) imposed for seasonal amplitude/phase analyses to be performed for a given station.\n",
    "By default set to 2 years and 0.6 respectively (i.e. **`--min_span 2 0.6`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Period limit for seasonal amplitude/phase analyses (**`--period_limit`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Period limit (years) imposed for seasonal amplitude/phase analyses to be performed for a given station."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Variogram density threshold (**`--densitythreshold DENSITYTHRESHOLD`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "For variogram plots, a given grid-cell is only valid if it contains this specified threshold of stations. \n",
    "\n",
    "By default set to 10 stations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Figure DPI (**`--figdpi FIGDPI`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "DPI to use for saving figures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot title (**`--user_title USER_TITLE`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Specify custom title for plots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Plot format (**`--plot_format PLOT_FMT`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "File format for saving plots. Default is PNG."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Colorbar bounds (**`--color_bounds CBOUNDS`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Set lower and upper-bounds for plot colorbars. For example: **`--color_bounds '0 100'`**\n",
    "\n",
    "By default set to the dynamic range of the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Colorbar percentile limits (**`--colorpercentile COLORPERCENTILE COLORPERCENTILE`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Set lower and upper percentile for plot colorbars. For example: **`--colorpercentile 30 100`**\n",
    "\n",
    "By default set to 25% and 95%, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Superimpose individual stations over gridded array (**`--stationsongrids`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In gridded plots, superimpose your gridded array with a scatterplot of station locations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Draw gridlines (**`--drawgridlines`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "In gridded plots, draw gridlines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Generate all supported plots (**`--plotall`**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Generate all supported plots, as outlined under sections #4, #5, and #6 above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "## Examples of the **`raiderStats.py`** program\n",
    "<a id='examples'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Example 1. Generate all individual station scatter-plots, as listed under [section #4](#overview_4) <a id='example_1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Using the file **`UNRcombinedGPS_ztd.csv`** generated by **`raiderDownloadGNSS.py`** as input, produce plots illustrating station distribution, mean tropospheric zenith delay by station, and standard deviation of tropospheric zenith delay by station.\n",
    "\n",
    "Restrict the temporal span of the analyses to all data acquired between 2016-01-01 and 2018-12-31, and restrict the spatial extent to a geographic bounding box '32 42 -125 -114' spanning over California.\n",
    "**NOTE** if not specified, by default the temporal bounds are set to earliest and latest time found in input file, and by default the maximum spatial bounds computed based off of the coordinates in the input file is passed.\n",
    "\n",
    "The footprint of the specified geographic bounding box is depicted in **Fig. 1**.\n",
    "\n",
    "These basic spatiotemporal constraints will be inherited by all successive examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex1 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --station_distribution --station_delay_mean --station_delay_stdev --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex1/figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we visualize the spatial distribution of stations (*ZTD_station_distribution.png*) as blue markers.\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex1/figures/ZTD_station_distribution.png\" alt=\"ZTD_station_distribution\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex1 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --station_distribution\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the mean tropospheric zenith delay by station (*ZTD_station_delay_mean.png*) with a **`hot`** colorbar. \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex1/figures/ZTD_station_delay_mean.png\" alt=\"ZTD_station_delay_mean\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex1 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --station_delay_mean\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the standard deviation of tropospheric zenith delay by station (*ZTD_station_delay_stdev.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex1/figures/ZTD_station_delay_stdev.png\" alt=\"ZTD_station_delay_stdev\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex1 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --station_delay_stdev\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Example 2. Generate all basic gridded station plots, as listed under [section #5](#overview_5)<a id='example_2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Produce plots illustrating gridded station distribution, gridded mean tropospheric zenith delay, and gridded standard deviation of tropospheric zenith delay.\n",
    "\n",
    "Also save gridded arrays as raster files with **`--grid_to_raster`** so as to more conveniently replot in successive script calls (recommended).\n",
    "\n",
    "Finally, use the maximum number of allowed CPUs to save some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex2 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --grid_heatmap --grid_delay_mean --grid_delay_stdev --grid_to_raster --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated rasters (i.e. the TIF files in the specified output directory):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex2/figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the heatmap of gridded station array (*ZTD_grid_heatmap.png*) with a **`hot`** colorbar.\n",
    "\n",
    "Note that the colorbar bounds are somewhat saturated, which demonstrates the utility of plotting options outlined under section #7 such as **`--color_bounds`** and **`--colorpercentile`**\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex2/figures/ZTD_grid_heatmap.png\" alt=\"ZTD_grid_heatmap\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex2 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --grid_heatmap --grid_to_raster\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded mean tropospheric zenith delay (*ZTD_grid_delay_mean.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex2/figures/ZTD_grid_delay_mean.png\" alt=\"ZTD_grid_delay_mean\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex2 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --grid_delay_mean --grid_to_raster\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded standard deviation of tropospheric zenith delay (*ZTD_grid_delay_stdev.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex2/figures/ZTD_grid_delay_stdev.png\" alt=\"ZTD_grid_delay_stdev\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex2 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --grid_delay_stdev --grid_to_raster\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Example 2a. Redo plots efficiently using generated grid raster files <a id='example_2a'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Since we haved the saved gridded arrays as raster files by specifying the **`--grid_to_raster`** option, we may directly replot these graphics by passing a given output raster file as input (e.g. `--file ZTD_grid_heatmap.tif`).\n",
    "\n",
    "This is a practical, efficient means to adjust/replot graphics and save a great deal of time by avoiding the gridding/prep steps involved with processing the initial input CSV file, especially for cases which span continental/global scales.\n",
    "\n",
    "Note though that since the output rasters are static with respect to the original specified spatiotemporal constraints (e.g. `--bounding_box` and `--timeinterval`), you cannot adjust such options with the rasters as input arguments. These rasters must be re-computed for any adjusted spatiotemporal parameters (if necessary) before replotting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "For this replotting command, let us also adjust the colorbar-bounds (using the `--color_bounds` option)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file maps_ex2/ZTD_grid_heatmap.tif --workdir maps_ex2a --color_bounds '10 40'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the replotted heatmap of gridded station array (*ZTD_grid_heatmap.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex2a/figures/ZTD_grid_heatmap.png\" alt=\"ZTD_grid_heatmap\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Example 3. Generate gridded mean tropospheric zenith delay plot, with stations superimposed <a id='example_3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Produce plot illustrating gridded mean tropospheric zenith delay, superimposed with individual station locations (`--stationsongrids`).\n",
    "\n",
    "Additionally, subset data in time for spring. I.e. **`'03-21 06-21'`**\n",
    "\n",
    "Finally, use the maximum number of allowed CPUs to save some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex3 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --seasonalinterval '03-21 06-21' --grid_delay_mean --stationsongrids --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex3/figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded mean tropospheric zenith delay (*ZTD_grid_delay_mean.png*) with a **`hot`** colorbar, with superimposed station locations denoted by blue markers.\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex3/figures/ZTD_grid_delay_mean.png\" alt=\"ZTD_grid_delay_mean\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Example 4. Generate variogram plots, as listed under [section #6](#overview_6)<a id='example_4'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Variogram analysis overview\n",
    "\n",
    "A variogram is the characterization of the difference between field values at two locations across a specified period of time `--timeinterval` and/or annual interval `--seasonalinterval`.\n",
    "\n",
    "Variogram analysis through this program is performed for each grid-cell containing at least a specified threshold of stations `--densitythreshold` (by default 10 stations) and with <=1000 randomly selected pairs (*hardcoded* as `Nsamp=1000`). Before variogram fitting, the data across each grid-cell is detrended using a linear ramp derived from a least-squares solution.\n",
    "\n",
    "The empirical variogram, or semivariance, provides a description of how the data are correlated with distance and is derived as half of the square difference between two observations:\n",
    "1/2 * sqrt(obs1 - obs2)\n",
    "\n",
    "The empirical variogram is then binned into increments of 20, with the top 3rd percentile discarded. These parameters are *hardcoded*.\n",
    "\n",
    "Finally, the experimental variogram is derived using the binned variogram above as input and utilizing an exponential model fit by default (*hardcoded*). The experimental variogram practically is a discrete function calculated using a measure of variability between pairs of points at various distances. **Note** toggling the `--binnedvariogram` option would lead to the an experimental variogram fit to every binned time-slice for a given grid-cell, whereas the default is to bin all data across all time-slices and then fit the experimental variogram for a given grid-cell\n",
    "\n",
    "From this procedure, we derive and plot the output sill/variance and range values. Respectively, the sill/variance represents the limit of the variogram representing where the variogram appears to level off. And the range represents the distance in which the difference of the variogram from the sill becomes negligible, and should not exceed half the size of a given grid-cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Produce plots illustrating empirical/experimental variogram fits per gridded station and time-slice (**`--variogram_per_timeslice`**) and also spanning the entire time-span. Plots of gridded station experimental variogram-derived sill and range values also generated.\n",
    "\n",
    "Additionally, subset data in time for spring. I.e. **`'03-21 06-21'`**\n",
    "\n",
    "Also save gridded arrays as raster files with **`--grid_to_raster`** so as to more conveniently replot in successive script calls (recommended).\n",
    "\n",
    "Finally, use the maximum number of allowed CPUs to save some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex4 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2019-12-31' --seasonalinterval '03-21 06-21' --variogramplot --variogram_per_timeslice --grid_to_raster --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated variograms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex4/variograms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "There are several subdirectories corresponding to each grid-cell that each contain empirical and experimental variograms generated for each time-slice (e.g. **`grid34_timeslice20160321_justEMPvariogram.eps`** and **`grid34_timeslice20160321_justEXPvariogram.eps`**, respectively) and across the entire sampled time period (**`grid34_timeslice20160321–20180619_justEMPvariogram.eps`** and **`grid34_timeslice20160321–20180619_justEXPvariogram.eps`**, respectively).\n",
    "\n",
    "Recall that the former pair of empirical/experimental variograms per time-slice are generated only if the **`---variogram_per_timeslice`** option is toggled. By default only the latter two pair of empirical/experimental variograms spanning the entire time-span are generated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we visualize the total empirical variogram corresponding to the entire sampled time period for grid-cell 34 in the array (*grid34_timeslice20160321–20180619_justEMPvariogram.eps*). \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex4/variograms/grid34/grid34_timeslice20160321–20180619_justEMPvariogram.png\" alt=\"justEMPvariogram\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we visualize the total experimental variogram corresponding to the entire sampled time period for grid-cell 34 in the array (*grid34_timeslice20160321–20180619_justEXPvariogram.eps*). \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex4/variograms/grid34/grid34_timeslice20160321–20180619_justEXPvariogram.png\" alt=\"justEXPvariogram\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The central coordinates for all grid-nodes that satisfy the specified station density threshold (**`--densitythreshold`**, by default 10 stations per grid-cell) for variogram plots are stored in a lookup table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!head maps_ex4/variograms/gridlocation_lookup.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the other generated figures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex4/figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded experimental variogram range (*ZTD_grid_range.png*) with a **`hot`** colorbar.\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex4/figures/ZTD_grid_range.png\" alt=\"ZTD_grid_range\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded experimental variogram variance (*ZTD_grid_variance.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex4/figures/ZTD_grid_variance.png\" alt=\"ZTD_grid_variance\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Examine extreme range of variogram fits\n",
    "\n",
    "Below are a couple case examples of total experimental variograms corresponding to grid-cells which are associated with extreme range values. This may demonstrate the limitations of this approach and/or input dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we visualize the total experimental variogram corresponding to the entire sampled time period for grid-cell 12 in the array (*grid12_timeslice20160321–20180619_justEXPvariogram.eps*). \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex4/variograms/grid12/grid12_timeslice20160321–20180619_justEXPvariogram.png\" alt=\"justEXPvariogram\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we visualize the total experimental variogram corresponding to the entire sampled time period for grid-cell 89 in the array (*grid89_timeslice20160321–20180619_justEXPvariogram.eps*). \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex4/variograms/grid89/grid89_timeslice20160321–20180619_justEXPvariogram.png\" alt=\"justEXPvariogram\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded experimental variogram range labelled with grid-cells of select variograms from above (*ZTD_grid_range_labelled.png*).\n",
    "\n",
    "<img src=\"support_docs/ZTD_grid_range_labelled.png\" alt=\"ZTD_grid_range\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Example 5. Generate seasonal phase/amplitude plots, as listed under [section #5](#overview_5)<a id='example_5'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "#### Seasonal phase/amplitude analysis overview\n",
    "\n",
    "The following sinusoidal form is used in this program to derive seasonal metrics:\n",
    "\n",
    "**A * sin(ωt + p) + c**\n",
    "\n",
    "\n",
    "Where *A* = amplitude (m), *ω* = angular frequency (rad/s), *p* = phase shift (days, positive to the left of origin), *c* = offset (days)\n",
    "\n",
    "And from which we derive and pass frequency (1/years) = (*ω*/2π) * 31556952, and period = 1/frequency (scaled to years)\n",
    "\n",
    "Finally the conventional origin *t0* = January 1, 1970"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "<img src=\"support_docs/Sinusoidal_form.png\" alt=\"footprint\" width=\"700\">\n",
    "<center><b>Fig. 2</b> Sinusoidal form used to derive seasonal metrics.</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Produce plots illustrating seasonal phase/amplitude/period fits for each individual station (**`--station_seasonal_phase`**) and averaged across each grid-cell (**`--grid_seasonal_phase`**). The standard deviation is also plotted across each grid-cell.\n",
    "\n",
    "Control the prerequisite minimum time-series span (years) a given station TS must span, and the prerequisite minimum fractional observations in span (fraction across specified `--timeinterval TIMEINTERVAL`, by default the entire span of input dataset). Here, we will specify a minimum time-series span of 2 years, and minimum fraction observation of 0.1 (i.e. **`--min_span 2 0.1`**). **Note** the latter constraint is intentionally loose as we have only downloaded data for every 5 days of the years, such that 5/365 ~ 0.014 would be expected for 'full' coverage through the specified time interval. Thus, the denser the temporal sampling, the more conservative you can be with the fractional observation control.\n",
    "\n",
    "Save figures of curve-fits vs data per station with  **`--phaseamp_per_station`** for debugging purposes. Not recommended for large-scale runs in the interest of speed/efficiency.\n",
    "\n",
    "Also save gridded arrays as raster files with **`--grid_to_raster`** so as to more conveniently replot in successive script calls (recommended).\n",
    "\n",
    "Finally, use the maximum number of allowed CPUs to save some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file products/UNRcombinedGPS_ztd.csv --workdir maps_ex5 --bounding_box '32 42 -125 -114' --grid_seasonal_phase --station_seasonal_phase --min_span 2 0.1 --phaseamp_per_station --grid_to_raster --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated rasters (i.e. the TIF files in the specified output directory):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated debug figures illustrating the curve-fits vs data (**`--phaseamp_per_station`**):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex5/phaseamp_per_station"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we visualize the time-series and curve-fit corresponding for one of the stations.\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/phaseamp_per_station/stationP335.png\" alt=\"stationCACH\" width=\"700\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!ls maps_ex5/figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Now we can take a look at the generated plots:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the seasonal phase of tropospheric zenith delay by station (*ZTD_station_seasonal_phase.png*) with a **`hot`** colorbar.\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_station_seasonal_phase.png\" alt=\"ZTD_station_seasonal_phase\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded mean of the station-wise seasonal phase of tropospheric zenith delay (*ZTD_grid_seasonal_phase.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_grid_seasonal_phase.png\" alt=\"ZTD_grid_seasonal_phase\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the seasonal amplitude of tropospheric zenith delay by station (*ZTD_station_seasonal_amplitude.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_station_seasonal_amplitude.png\" alt=\"ZTD_station_seasonal_amplitude\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded mean of the station-wise seasonal amplitude of tropospheric zenith delay (*ZTD_grid_seasonal_amplitude.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_grid_seasonal_amplitude.png\" alt=\"ZTD_grid_seasonal_amplitude\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the seasonal period of tropospheric zenith delay by station (*ZTD_station_delay_period.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_station_delay_period.png\" alt=\"ZTD_station_delay_period\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded period of the station-wise seasonal period of tropospheric zenith delay (*ZTD_grid_seasonal_period.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_grid_seasonal_period.png\" alt=\"ZTD_grid_seasonal_period\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded period standard deviation of the station-wise seasonal phase of tropospheric zenith delay (*ZTD_grid_seasonal_period_stdev.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex5/figures/ZTD_grid_seasonal_period_stdev.png\" alt=\"ZTD_grid_seasonal_period_stdev\" width=\"700\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true,
    "hidden": true
   },
   "source": [
    "### Example 6. Generate weather model/GNSS residual plots <a id='example_6'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Produce plots illustrating the residual between tropospheric zenith delay at specified GNSS stations and collocated weather-model delay nodes.\n",
    "\n",
    "GNSS data will again be downloaded with **`raiderDownloadGNSS.py`**, and GMAO weather-model derived delay will be computed with **`raiderDelay.py`**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Virtually access GNSS station location and zenith delay information for the year '2019', for every 12 days, at a UTC time of day 'HH:MM:SS' of '00:00:00', and across a geographic bounding box '32 42 -125 -114' spanning over California.\n",
    "\n",
    "The footprint of the specified geographic bounding box is again depicted in **Fig. 1**.\n",
    "\n",
    "In addition to querying for multiple years, we will also experiment with using the maximum number of allowed CPUs to save some time! Recall again that the default number of CPUs used for parallelization is 8.\n",
    "\n",
    "Note these features and similar examples are outlined in more detail in the companion notebook **`raiderDownloadGNSS/raiderDownloadGNSS_tutorial.ipynb`**.\n",
    "\n",
    "**NOTE** the following two sections will take ~12 minutes to complete. It is *recommended* to take a coffee break!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderDownloadGNSS.py --out GNSS_2019 --date 20190101 20191231 12 --returntime '00:00:00' --bounding_box '32 42 -125 -114' --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Compute tropospheric zenith delay from the GMAO weather-model for the year '2019', for every 12 days, at a UTC time of day 'HH:MM:SS' of '00:00:00', at stations located across a geographic bounding box '32 42 -125 -114' spanning over California and captured in the `GNSS_2019/gnssStationList_overbbox.csv` list generated by the `raiderDownloadGNSS.py` above (i.e. `--station_file GNSS_2019/gnssStationList_overbbox.csv`), and with an integration height limit `--zref` of 30,000 m. \n",
    "\n",
    "The footprint of the specified geographic bounding box is again depicted in **Fig. 1**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!mkdir GMAO_2019\n",
    "!cd GMAO_2019\n",
    "!raiderDelay.py --model GMAO --date 20190101 20191231 12 --time 00:00 --station_file ../GNSS_2019/gnssStationList_overbbox.csv --zref 30000 -v\n",
    "!cd ../"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Combine delay files derived above from the GMAO weather-model (`--raider 'GMAO_2019/GMAO_Delay_*.csv' --raiderDir GMAO_2019 --raider_column totalDelay`) and GNSS stations (`GNSS_2019/UNRcombinedGPS_ztd.csv --column ZTD`) respectively, passing only data which are collocated in space and time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderCombine.py --gnss GNSS_2019/UNRcombinedGPS_ztd.csv --column ZTD --raider 'GMAO_2019/GMAO_Delay_*.csv' --raiderDir GMAO_2019 --raider_column totalDelay --out Combined_delays_GNSSandGMAO_2019.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Using the file **`Combined_delays_GNSSandGMAO_2019.csv`** generated by **`raiderCombine.py`** as input and passing the weather-model/GNSS residual values (`--column_name ZTD_minus_RAiDER`), produce plots illustrating mean tropospheric zenith delay by station + across each grid-cell, and standard deviation of tropospheric zenith delay by station + across each grid-cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "!raiderStats.py --file Combined_delays_GNSSandGMAO_2019.csv --column_name ZTD_minus_RAiDER --workdir maps_ex6 --bounding_box '32 42 -125 -114' --station_delay_mean --station_delay_stdev --grid_delay_mean --grid_delay_stdev --grid_to_raster --cpus all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the mean tropospheric zenith delay by station (*ZTD_station_delay_mean.png*) with a **`hot`** colorbar. \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex6/figures/ZTD_minus_RAiDER_station_delay_mean.png\" alt=\"ZTD_minus_RAiDER_station_delay_mean\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file Combined_delays_GNSSandGMAO_2019.csv --column_name ZTD_minus_RAiDER  --workdir maps_ex6 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --station_delay_mean\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the standard deviation of tropospheric zenith delay by station (*ZTD_minus_RAiDER_station_delay_stdev.png*).  \n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex6/figures/ZTD_minus_RAiDER_station_delay_stdev.png\" alt=\"ZTD_minus_RAiDER_station_delay_stdev\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file Combined_delays_GNSSandGMAO_2019.csv --column_name ZTD_minus_RAiDER  --workdir maps_ex6 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --station_delay_stdev\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded mean tropospheric zenith delay (*ZTD_minus_RAiDER_grid_delay_mean.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex6/figures/ZTD_minus_RAiDER_grid_delay_mean.png\" alt=\"ZTD_minus_RAiDER_grid_delay_mean\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file Combined_delays_GNSSandGMAO_2019.csv --column_name ZTD_minus_RAiDER  --workdir maps_ex6 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --grid_delay_mean\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Below is the gridded standard deviation of tropospheric zenith delay (*ZTD_minus_RAiDER_grid_delay_stdev.png*).\n",
    "\n",
    "<img src=\"support_docs/maps/maps_ex6/figures/ZTD_minus_RAiDER_grid_delay_stdev.png\" alt=\"ZTD_minus_RAiDER_grid_delay_stdev\" width=\"700\">\n",
    "\n",
    "To generate this figure alone, run:\n",
    "```\n",
    "!raiderStats.py --file Combined_delays_GNSSandGMAO_2019.csv --column_name ZTD_minus_RAiDER  --workdir maps_ex6 --bounding_box '32 42 -125 -114' --timeinterval '2016-01-01 2018-12-31' --grid_delay_stdev\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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": false,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}