GDV Tool 0.5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Disclaimer: This tool is still being developed and should be used as only an experimental exploratory aid at this point in time. Do not base policy decisions on these analysis results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "code_folding": []
   },
   "source": [
    "# Groundwater dependent vegetation detection, trend and change detection tool\n",
    "***\n",
    "### Background and summary\n",
    "Dewatering of mine pits contributes to groundwater depletion, potentially placing proximal groundwater-dependent vegetation at risk of water stress. While earth observation data has been trialled to monitor the effects of dewatering on vegetation health, a lack of ground truthing has limited its uptake. This project has developed a validated remote sensing approach for detecting deep-rooted vegetation species (e.g. <i>Eucalyptus camaldulensis, E. victrix, Melaleuca argentea</i>) and measuring/monitoring declining heath symptoms (e.g. reduced greenness/leaf moisture). This approach has been developed into a basic tool that leverages the Digital Earth Australia Platform (DEA), providing access to temporal Landsat and Sentinel 2 Imagery for analysis. \n",
    "\n",
    "A spatial multicriteria evaluation (SMCE) model is used to integrate a temporal series of remotely-sensed vegetation indices to estimate groundwater dependent vegetation likelihood on a scale of 0-1. Detected high-likelihood areas can be explored further for health increase/decrease via Mann-Kendall trend analyses (original and seasonal methods), change point detection (dynamic programming and pelt methods), and/or the change detection change vector analysis. Model accuracy has been calibrated and validated using over 500 ground-truthed vegetation health assessment sites captured in 2019. Results are presented below using interactive Jupyter notebook controls. \n",
    "\n",
    "### Problem statement\n",
    "Mine dewatering is the removal of unwanted groundwater from open pit mines to ensure continual operation and allow mineral extraction. In some cases, this process can result in groundwater depletion, affecting the health of vegetation that draws most of its water requirements from near the water-table, potentially resulting in mortality. This has elicited a requirement for industry to continually monitor and report on vegetation health. \n",
    "\n",
    "This monitoring has traditionally been undertaken through floristic surveys to visually record the status of vegetation, with studies often done annually or seasonally (3-6 months). However, it is poorly targeted and annual observations may miss early identification of changes that could otherwise allow proactive management. An improved strategy involves directing resources at locations where vegetation health has declined relative to historical observations. Remote sensing-based monitoring can compare current observations to a long historical bank of imagery to detect anomalous changes with <i>c.</i> monthly frequency.\n",
    "\n",
    "### Aims\n",
    "This project/tool seeks to determine if vegetation health monitoring can be achieved through remote sensing techniques leveraging the growing online and accessible archive of earth observational data. There are four aims:\n",
    "1. Develop GDV likelihood model for mapping potential GDV areas to assist survey targeting;\n",
    "2. Validate and refine GDV likelihood model via field visits and ground-truthing;\n",
    "3. Implement time-series trend and change analysis to detect GDV decline and change dynamics; and\n",
    "4. Develop into interactive Jupyter Notebook as a foundation for future tool development and deliver to project partners.\n",
    "\n",
    "### Methodology\n",
    "Please see associated technical report for a deeper explanation of the methodology underpinning this tool.\n",
    "\n",
    "### Datasets\n",
    "Satellite imagery from Landsat 5, 7, 8 and Sentinel-2 is freely available and have a reveist time over Australia of 15 days and 3-5 days, respectively. The expansive historical archive, dating back to 1986, makes Landsat ideal for monitoring long-term, gradual changes to vegetation. Sentinel-2 has a smaller archive (mid-2015 onwards only), but can provide a higher resolution image (10 metre) to improve vegetation detection. Vegetation health observation points captured during the validation phase of this project are also provided in the data folder, which can be used to validate the model for several study areas (optional).\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tool section overview\n",
    "This tool will take users through numerous processing steps that are required to detect vegetation and calculate health trends and changes. An overview of the key sections of this tool are outlied below:<br></br>\n",
    "<u>Section 00: Initialise code</u>:\n",
    "1. Initialise tool libraries, storage, files, etc.\n",
    "\n",
    "<u>Section 01: Generate GDV likelihood</u>:\n",
    "1. Load Landsat 5-7 and Sentinel-2 images for an area of interest\n",
    "2. Compute vegetation and moisture indices to highlight presence of vegetation (various indices offered)\n",
    "3. Interpolate missing data (optional) and standardise (via invariant site detection, orthogonal polynomial coefficients and fuzzy sigmoidals)\n",
    "4. Calculate, threshold and visualise groundwater dependent vegetation likelihood data via pairwise comparison weighting and ROC curve/standard deviation\n",
    "\n",
    "<u>Section 02: Perform Mann-Kendall trend analyses</u>: \n",
    "1. Load Landsat/Sentinel data, compute vegetation index, interpolate and standardise, apply a statistical test to remove outliers\n",
    "2. Set trend analysis parameters for mann-kendall analysis type (e.g. original or seasonal, quarterly or monthly, trend significance and direction)\n",
    "3. Visualise trend analysis results\n",
    "\n",
    "<u>Section 03: Perform spectral change vector analysis</u>: \n",
    "1. Load Landsat/Sentinel data, compute vegetation index, interpolate and standardise, apply a statistical test to remove outliers\n",
    "2. Set baseline and comparison year date ranges and perform change vector analysis over time period\n",
    "3. Visualise change vector analysis results\n",
    "***\n",
    "### Getting started\n",
    "Simply start from 'Section 00' below and run each cell in a downward direction. Some steps may take less than a second to process, others may take minutes. Missing a step (or 'cell') will eventually cause an error, so it is vital all steps are run and in order. If you make a mistake, the tool can be quickly reset by clicking the 'Restart the Kernal' icon in the top menu (i.e. the refresh icon). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 00: Initialise code\n",
    "This is a short (but vital) section of the tool, as it loads all external code into the Jupyter Notebook. This code allows you to perform analyses, and the GDV Tool cannot be used without it. All code loaded by the two code cells in this section exist in folder 'gdv_scripts' and 'dea_scripts' in the project folder.\n",
    "\n",
    "Please run cells 00, 00A and 00B below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 00: Initialise required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install pymannkendall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install ruptures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 00A: Initialise modules and scripts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# import python modules\n",
    "import warnings\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import HTML, Button, Output\n",
    "from IPython.display import clear_output\n",
    "\n",
    "# initialise user-defined modules and scripts\n",
    "from gdv_scripts import helpers         # general helper functions used throughout tool\n",
    "from gdv_scripts import parameters      # storage classes for holding user control parameter changes\n",
    "from gdv_scripts import controls        # methods for initialising and observing ui widget controls\n",
    "from gdv_scripts import maps            # methods for initialising and observing leaflet maps\n",
    "from gdv_scripts import data            # methods for fetching and preparing data cube data\n",
    "from gdv_scripts import indices         # methods for fetching and preparing vege/moist indices\n",
    "from gdv_scripts import interps         # methods for interpolation\n",
    "from gdv_scripts import stats           # methods for generating statistics and analyses\n",
    "from gdv_scripts import plotting        # methods for plotting various things!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 00B: Create instances of required parameter storage classes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# create instances of param storage classes\n",
    "prms_general = parameters.GeneralParams()     # hold general parameters (e.g. project name, output folder, etc.)\n",
    "prms_like = parameters.LikelihoodParams()     # hold parameters for gdv likelihood (e.g. veg index, stability, thresholds)\n",
    "prms_trend = parameters.TrendParams()         # hold parameters for trend analysis (e.g. veg index, trend type, time slice) \n",
    "prms_change = parameters.ChangeParams()       # hold parameters for change detection analysis (e.g. veg index, time slice, thresholding)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 01: Generate GDV likelihood layers\n",
    "The first key section of this tool involves the use of spatial multicriteria evaluation modelling techniques to differentiate GDV from non-GDV to an acceptable level of accuracy using freely available remote sensing satellite imagery. It is hypothesised that GDV can be differentiated from non-GDV on satellite imagery by combining metrics of greenness, leaf moisture and stability with the construction of a weighting scheme disproportionally in favour of the dry season (September-November) vs the wet season (January-March). \n",
    "\n",
    "The key processing steps applied in this section are as follows:\n",
    "1. Fetch a ‘cube’ of Landsat or Sentinel satellite data for a user-defined study area and timespan using the DEA ODC;\n",
    "2. Generate vegetation and moisture indices based on user-selection;\n",
    "3. Reduce data to wet (Jan-Mar) and dry (Sep-Nov) annual medians and interpolate missing data (optional);\n",
    "4. Standardise vegetation and moisture using an invariant target site and orthogonal polynomial coefficients approach;\n",
    "5. Re-scale values using fuzzy sigmoidal functions;\n",
    "6. Generate vegetation and moisture stability for each year using wet and dry seasons;\n",
    "7. Produce a GDV likelihood layers via Analytical Hierarchy Process (AHP) Pairwise Comparison weighting approach;\n",
    "8. Slice GDV likelihood layers into GDV/non-GDV areas using ground-truth point data (optional) or basic standard deviation;\n",
    "9. Export results as GDV likelihood layers (median GDV likelihood of all time or per year in data 'cube').\n",
    "\n",
    "Proceed to step 01A when ready."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01A: Set satellite imagery requirements and define study area boundary\n",
    "The following control panel is used to set the type of satellite data (e.g. Landsat or Sentinel), the range of dates you require for analysis, and how the tool will handle clouds in satellite imagery. This section also allows you to select (or draw) a study area for which to generate GDV likelihood layers.\n",
    "\n",
    "The control panel parameters are:\n",
    "* `platform`: Satellite platform to fetch data from (i.e. Landsat 5-8, Sentinel 2AB). Recommended: Landsat.\n",
    "Landsat offers 35+ years of satellite data at a pixel resolution of 30m. Sentinel 2 offers a superior pixel resolution of 10m, but currently offers < 5 years of satellite data. Simply speaking, Sentinel 2 is more likely to produce more accurate GDV likelihood layers due to its improved pixel resolution, but far less robust trend detection. Landsat is, thus, just the opposite. Despite the improved pixel resolution of Sentinel 2, Landsat is recommended as using 10+ years of data is typically ideal in all sections of this tool.\n",
    "* `product`: Satellite data DEA ODC analysis ready product to use (i.e. NBAR, NBART). Recommended: NBART.\n",
    "Satellite imagery requires significant atmospheric correction for time-series analysis. DEA ODC offers atmospherically-corrected products in the form of NBAR (Nadir-corrected BRDF Adjusted Reflectance) and NBART (Nadir-corrected BRDF Adjusted Reflectance with Terrain Illumination Reflectance Correction). For information, see: https://docs.dea.ga.gov.au/data/data.html.\n",
    "* `query dates`: Start and end dates for which satellite imagery will be fetched (e.g. 1999 to 2019). Recommended: range of >= 10 years.\n",
    "Satellite imagery (and thus analysis) will be limited to the start and end dates set here. It is recommended >= 10 years is used to reduce the influence of seasonal fluctuations during GDV likelihood layer generation. Note: end date is inclusive. A minimum of three years required.\n",
    "* `max cloud cover`: Maximum % of clouds allowed in satellite images. Recommended: between 2% - 15%.\n",
    "Clouds can result in inaccurate vegetation and moisture indices. If the percentage of clouds is greater than the value set here, the image will not be used in analysis. \n",
    "* `mask clouds`: Mask out any satellite image pixels that contain clouds if they are found to exist. Recommended: Yes.\n",
    "Cloud pixels are converted into 'nan' if this is set to Yes. Otherwise, the cloud pixel will be left as is in the satellite imagery.\n",
    "\n",
    "Under the control panel is an interactive map for selecting/drawing a study area boundary. Numerous pre-defined study areas (see red rectangles) have already been drawn for use during development. Feel free to use these by clicking to select. You may also draw a new study area using the 'draw new boundary' icon (the black square). Map controls (from top left to bottom right) include:\n",
    "* `full screen mode`: Make the interactive map full screen size.\n",
    "* `measurement tool`: Take linear and path-based distance and area measurements.\n",
    "* `zoom in/out`: Zoom map in (+) and out (-).\n",
    "* `draw new boundary`: Draw a new study area boundary. Click, drag and let go of mouse to set a new study area, then proceed to next step.\n",
    "* `edit new boundary`: Make edits to or delete newly drawn study area boundary.\n",
    "* `base map layers`: Choose between ESRI satellite imagery or Open Street Maps topographic basemap layers.\n",
    "\n",
    "Tips:\n",
    "* If drawing a new study area, an area between 15000 - 35000 ha is recommended. Small study areas may result in processing problems, while very large study areas will cause memory errors.\n",
    "* A study area that covers a range of land types is recommended. Try and capture disturbance, bare-ground, vegetation, drainage, hills.\n",
    "* Regarding query dates, 10+ years is recommended for analysis. This will reduce seasonal and bush fire variation in imagery. A minimum of 3 years required.\n",
    "\n",
    "Run the code cell below to enable the control panel and interactive map. Set your required parameters. Proceed to step 01B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build satellite data fetch dashboard and study area selection map\n",
    "def build_sat_fetch_dash_and_map(params, study_center):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params:\n",
    "        print('Likelihood parameters not set. Did you run the above cells?')\n",
    "        return None, None\n",
    "    elif not all(study_center):\n",
    "        print('No study center. Did you run the above cells?')\n",
    "        return None, None\n",
    "       \n",
    "    # build dash and map\n",
    "    sat_fetch_dash = controls.build_sat_fetch_dashboard(params=params)                  # dashboard\n",
    "    study_area_map = maps.load_sat_fetch_map(params=params, study_center=study_center)  # map\n",
    "    \n",
    "    # return dash and map\n",
    "    return sat_fetch_dash, study_area_map\n",
    "\n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_like = parameters.LikelihoodParams()\n",
    "\n",
    "# get sat fetch dash and map\n",
    "sat_fetch_dash, study_area_map = build_sat_fetch_dash_and_map(prms_like, prms_general.study_center)\n",
    "display(sat_fetch_dash, study_area_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01B: Fetch satellite imagery from Open Data Cube\n",
    "Fetch satellite data from the Digital Earth Australia (DEA) Open Data Cube (ODC) for the satellite options and study area set in the above control panel. For context, DEA is a digital platform that catalogues large amounts of Earth observation data covering Australia. It is underpinned by the ODC, which is an open source software package that provides a Python-based API for high performance querying and data access. For more information, visit (https://www.ga.gov.au/dea/odc).\n",
    "\n",
    "This step is easy. Simply run the code cell beneath these notes to begin fetching satellite data. Note: this can take a few minutes, depending on the study area size and date range size.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Fetch raw satellite data into an xarray dataset;\n",
    "2. Rename bands for conformity reasons;\n",
    "3. Drop unrequired bands (e.g. metadata);\n",
    "4. Reduce satellite data to wet and dry Pilbara seasons (Jan-Mar and Sep-Nov, respectively)\n",
    "\n",
    "Run the code cell below to fetch satellite data. This may take a few minutes. Proceed to step 01B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# fetch cube data for gdv likelihood analysis and mapping\n",
    "def fetch_gdv_cube_data(params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params.check_general_params():\n",
    "        print('Not all required parameters set above. Please choose study area and satellite parameters.')\n",
    "        return None\n",
    "    elif helpers.valid_year_range(params.query_dates) < 3:\n",
    "        print('Date range is too small, three or more years required. We suggest >= 10 years. Please adjust above.')\n",
    "        return None \n",
    "\n",
    "    # fetch data based on user parameters\n",
    "    ds = data.fetch_cube_data(platform=params.platform, product=params.product, q_dates=params.query_dates, \n",
    "                              cloud_cover=params.cloud_cover, cloud_mask=params.cloud_mask, \n",
    "                              coords_nw=params.alb_nw, coords_se=params.alb_se)    \n",
    "\n",
    "    if ds:\n",
    "        # reduce dataset down to wet (jfm) and dry (son) months\n",
    "        ds = data.reduce_ds_to_wet_and_dry(ds)\n",
    "\n",
    "        # remove above dashboard and map\n",
    "        print('\\nDisabling above dashboard and map. Please re-run the tool if you require a new study area.')\n",
    "        controls.disable_controls(sat_fetch_dash)\n",
    "        study_area_map.close()\n",
    " \n",
    "        return ds\n",
    "\n",
    "    else:\n",
    "        return None\n",
    "        \n",
    "ds = fetch_gdv_cube_data(prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01C: Set vegetation and moisture indices and options\n",
    "The following control panel is used to set the required vegetation and moisture indices (and associated options) that underpin the spatial multicriteria evaluation modelling techniques implemented in this tool. Ultimately, they form the foundation that is used to generate GDV likelihood layers. \n",
    "\n",
    "Remote sensing of vegetation in semi-arid regions like the Pilbara is challenged by significant soil background reflectance. Thus, choice of vegetation index is an important consideration, as model outputs are likely to be vary depending on chosen index. A suite of remote sensing indices were tested during the development of this tool, and those found to produce optimal results across multiple tests sites have been provided to users. Based on our preliminary findings, <i>SLAVI</i>, <i>MAVI</i>, <i>Tasselled Cap Greenness</i>, and <i>NDVI</i> produced consistently better GDV likelihood models. Regardless, the available vegetation and moisture indices include:\n",
    "\n",
    "Vegetation indices:\n",
    "* NDVI (Normalised Difference Vegetation Index) &emsp; &ensp; $\\frac{NIR - RED}{NIR + RED}$ &emsp; &emsp; &emsp;\n",
    "\n",
    "* SAVI (Soil-adjusted Vegetation Index): &emsp; &emsp; &emsp; &emsp; $\\frac{(1 + L)(NIR - RED)}{NIR + RED + L}$\n",
    "\n",
    "* STVI-3 (Stress-related Vegetation Index): &emsp; &emsp; &emsp; &emsp; $\\frac{(NIR)}{RED + MIR}$\n",
    "\n",
    "* SLAVI (Specific Leaf Area Vegetation Index): &emsp; &emsp; &emsp; $\\frac{(NIR)}{RED + SWIR}$\n",
    "\n",
    "* MSVI-1 (Mid-infrared Stress Vegetation Index 1): &emsp; &emsp; $\\frac{(NIR)}{MIR}$\n",
    "\n",
    "* MSVI-2 (Mid-infrared Stress Vegetation Index 2): &emsp; &emsp; $\\frac{(NIR)}{SWIR}$\n",
    "\n",
    "* MSVI-3 (Mid-infrared Stress Vegetation Index 3): &emsp; $\\frac{(NIR)}{MIR + SWIR}$\n",
    "\n",
    "* MAVI (Moisture-adjusted Vegetation Index): &emsp; &ensp; $\\frac{(NIR - RED)}{NIR + RED + SWIR}$\n",
    "\n",
    "* Greenness (Taselled Cap derived): &emsp; &emsp; &emsp; &emsp; &emsp; &emsp; &emsp; &ensp; $\\text{N/A}$\n",
    "\n",
    "Moisture indices:\n",
    "* NDMI (Normalised Difference Moisture Index): &emsp; &ensp; $\\frac{NIR - SWIR}{NIR + SWIR}$\n",
    "\n",
    "* GVMI (Global Vegetation Moisture Index): &emsp; &ensp; $\\frac{(NIR + 0.1) - (SWIRII + 0.02)}{(NIR + 0.1) + (SWIRII + 0.02)}$\n",
    "\n",
    "The control parameters on the below control panel are:\n",
    "* `vegetation index`: Vegetation index to be used in GDV likelihood modelling. Recommended: MAVI, SLAVI, TCAP Greenness, NDVI. \n",
    "One of nine vegetation indices can be used to represent greenness (Tucker 1979) in the spatial multicriteria evaluation modelling techniques used to generate GDV likelihood layers. See the bove formulae for more information. Further information these various indices can be found at the Index Database website (https://www.indexdatabase.de).\n",
    "* `moisture index`: Moisture index to be used in GDV likelihood modelling. Recommended: NDMI. \n",
    "one of two moisture indices can be used to represent leaf moisture in the spatial multicriteria evaluation modelling techniques used to generate GDV likelihood layers. Moisture has been found to be useful in vegetation monitoring within Australia's rangelands (e.g. O’Neill 1996; Andrew and Warrener 2017). See the bove formulae for more information. Further information these various indices can be found at the Index Database website (https://www.indexdatabase.de).\n",
    "* `soil adjustment (L) value`: Used to minimise soil brightness influences when using the SAVI index. Only applicable when using SAVI.\n",
    "An L value of 0.5 in reflectance space has been found to minimize soil brightness variations and eliminate the need for additional calibration for different soils when using SAVI. Useful to calibrate SAVI, if used.\n",
    "* `rescale values`: Rescales index value range from -1to1 to 0to2 where applicable. Recommended: Yes\n",
    "Used during index standardisation process and fuzzy sigmoidal rescaling. It is highly recommended users rescale indices.\n",
    "* `show figures`: Show a slide of all vegetation and moisture images once calculated. Recommended: Yes\n",
    "Setting this to Yes will print out a figure of all vegetation and moisture indices available in the ODC dataset. Useful to turn off when you are familiar with the study area.\n",
    "\n",
    "Tips:\n",
    "* Based on prior model validation tests, optimal indices for the Yandi, Roy Hill and Ophthalmia study areas include MAVI, SLAVI, NDVI or TCAP Greenness. One of these is recommended.\n",
    "* Hill-dense study areas trend to benefit from NDVI, SAVI or MAVI indices. Flat, stony plain-dense study areas benefit moe from SLAVI and TCAP.\n",
    "* It is highly recommended index values are rescaled (i.e. set Rescale values control to 'Yes'\n",
    "* Processing can take quite awhile depending on study area size. Please wait!\n",
    "\n",
    "Run the code cell below to enable the control panel. Set your required parameters. Proceed to step 01D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build vegetation and moisture index generation dashboard\n",
    "def build_indices_dash(params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params.check_raw_idx_params():\n",
    "        print('Vegetation/moisture parameters not set. Did you run the above cells?')\n",
    "        return None\n",
    "        \n",
    "    # build dash\n",
    "    indices_dash = controls.build_indices_dashboard(params=params)\n",
    "    \n",
    "    # return dash\n",
    "    return indices_dash\n",
    "    \n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_like.reset_raw_idx_params()\n",
    "\n",
    "# get indices dash\n",
    "indices_dash = build_indices_dash(prms_like)\n",
    "display(indices_dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01D: Generate vegetation and moisture indices from satellite imagery\n",
    "Generate vegetation and moisture indices for every date returned from the ODC using the parameters set in the above control panel. \n",
    "\n",
    "This step is easy. Simply run the code cell beneath these notes to begin generating vegetation and moisture indices. Note: this can take a few minutes, depending on the study area size and date range size.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Calculate vegetation and moisture indices from xarray dataset;\n",
    "2. Drop unnecessary bands;\n",
    "3. Resample data down to wet and dry medians for every year in dataset;\n",
    "4. Compute dataset into memory;\n",
    "5. Calculate and present summary statistics;\n",
    "6. Perform a preliminary z-test to highlight any potential data outliers;\n",
    "7. Produce graph of mean annual vegetation and moisture indices for wet and dry season; and\n",
    "8. Show figures of all indices (if requested)\n",
    "\n",
    "Run the code cell below to generate vegetation and moisture indices. This may take a few minutes. Proceed to step 01E."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     4,
     7,
     10,
     18,
     21,
     24
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# generate vegetation and moisture indices, resample down to wet and dry, provide summary stats and graph\n",
    "def generate_gdv_indices_wet_and_dry(ds, params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not ds:\n",
    "        print('No cube dataset provided. Please ensure all cells above have been run.')\n",
    "        return None\n",
    "    elif not indices.check_required_bands_exist(ds=ds):\n",
    "        print('Not all required bands available in dataset. Please re-run the tool again.')\n",
    "        return ds\n",
    "    elif not params.check_raw_idx_params():\n",
    "        print('No vegetation or moisture index set. Please ensure they are set above above.')\n",
    "        return ds\n",
    "    \n",
    "    # create vegetation and moisture indices and resample to wet and dry season\n",
    "    ds = indices.create_resample_idxs_wet_and_dry(ds=ds, params=params)\n",
    "    \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    elif len(ds.groupby('time.month')) != 2:\n",
    "        print('Num of months of data in dataset not equal to 2. Please include more years in your query above.')\n",
    "        return ds\n",
    "        \n",
    "    # remove above dashboard and map\n",
    "    print('Disabling above dashboard. Please re-run the tool if you require a different index.\\n')\n",
    "    controls.disable_controls(indices_dash)\n",
    "        \n",
    "    # calculate some summary stats and do quick z-score test to help user\n",
    "    stats.calc_summary_stats_wet_and_dry(ds=ds)\n",
    "    ds = stats.zscore_flag_show_reset_wet_dry(ds=ds, z_val=1.96, reset=False)  \n",
    "    \n",
    "    # show graphs\n",
    "    fig_raw_idx_stats = plotting.plot_raw_idx_stats(ds=ds, fignum=1)\n",
    "    \n",
    "    # plot figs if requested\n",
    "    if params.plot_raw_idx:\n",
    "        fig_raw_idx, fig_raw_mst = plotting.plot_raw_idx(ds=ds, stretch_min=0.01, stretch_max=0.99)\n",
    "\n",
    "    return ds\n",
    "\n",
    "ds = generate_gdv_indices_wet_and_dry(ds=ds, params=prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01E: Set vegetation and moisture index correction and interpolation options\n",
    "The following control panel is used to set optional vegetation and moisture index correction and interpolation options. \n",
    "\n",
    "Satellite imagery is prone to data gaps (usually due to the presence of considerable cloud cover), missing pixels and occasional outliers (e.g. post-bushfire greening, climate events). Some of these issues can be corrected through linear interpolation and statistical outlier testing (e.g. z-test). The methods provided here attempt to correct these issues. Please note, interpolation and outlier removal is completely optional (although recommended). \n",
    "\n",
    "The control parameters are:\n",
    "* `handle missing data`: Choose how to handle missing satellite images. Default: Interpolate.\n",
    "Interpolation is applied using linear interpolation and polynomials. Images are split into seasons, and given a missing image, the images at years before and after are used to perform linear interpolation. If any dates at the start and end year are missing, the whole year will be dropped. The alternative is to drop the image (i.e. the whole year) completely.\n",
    "* `fill missing pixels`: Fill in missing pixels for partially missing images. Default: No.\n",
    "Sometimes satellite images may contain missing pixels due to cloud masking processes. These can be 'filled in' via linear interpolation similar to above. It is recommended this is not applied.\n",
    "* `perform z-test`: Perform a z-test to detect and remove significant image outliers. Default: Yes.\n",
    "A z-test is used to validate the hypothesis that a sample (e.g. a vegetation or moisture image) belongs to the same population (i.e. the rest of the dataset). This is done by comparing each date (the sample) to the populaton mean and standard deviation, and is useful for determining significant outliers. See here for more information: https://en.wikipedia.org/wiki/Z-test.\n",
    "* `z-score value`: Set the z-score (or critical value) for the z-test, if applied. Default: 1.95 (equivalent to p-value < 0.05).\n",
    "Z-scores are measures of standard deviation. For example, 1.96 is interpreted as 1.96 standard deviations away from the mean. We can relate standard deviation with probabilities (p-values), allowing for the ability to attach z-scores to p-values. A z-score of 1.96, when using a 95% confidence level (p-value 0.05), are -1.96 and +1.96 standard deviations. In other words, any image z-score outside of this range will be dropped. To remove anything outside p-value 0.01 (i.e. very significant outliers), use a z-score of 2.58. Conversely, to use a p-value 0.1 (i.e. less significant outliers), use a z-score of 1.65.\n",
    "* `show figures`: Show a slide of all corrected and interpolated vegetation and moisture images once calculated.\n",
    "Setting this to Yes will print out a figure of all corrected/interpolated vegetation and moisture indices available in the ODC dataset. Useful to turn off when you are familiar with the study area.\n",
    "\n",
    "Tips:\n",
    "* Interpolation is optional.\n",
    "* A higher z-score value will remove more image outliers, a lower z-score value will remove less image outliers. A z-score of 1.96 is recommended.\n",
    "\n",
    "Run the code cell below to enable the control panel. Set your required parameters. Proceed to step 01F."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build clean and interpolate vegetation and moisture index dashboard\n",
    "def build_interpolation_dash(params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params:\n",
    "        print('Likelihood parameters not set. Did you run the above cells?')\n",
    "        return None\n",
    "        \n",
    "    # build dash\n",
    "    interpolation_dash = controls.build_interpolation_dashboard(params=params)\n",
    "    \n",
    "    # return dash\n",
    "    return interpolation_dash\n",
    "    \n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_like.reset_interpolate_idx_params()\n",
    "\n",
    "# get indices dash\n",
    "interpolation_dash = build_interpolation_dash(prms_like)\n",
    "display(interpolation_dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01F: Correct and interpolate vegetation and moisture indices\n",
    "Correct and interpolate vegetation and moisture image gaps using the parameters set in the above control panel.\n",
    "\n",
    "Simply run the code cell beneath these notes to begin correcting and interpolating vegetation and moisture indices, if requested.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Perform a z-test to highlight any potential data outliers;\n",
    "2. Check if first and last years in dataset are missing any seasons (must be dropped if so);\n",
    "3. Check for any missing images and perform linear interpolation (if requested), else drop the year;\n",
    "4. Calculate and present summary statistics;\n",
    "7. Produce graph of mean annual vegetation and moisture indices for wet and dry season; and\n",
    "8. Show figures of all indices (if requested)\n",
    "\n",
    "Run the code cell below to generate vegetation and moisture indices. This may take a few minutes. Proceed to step 01G."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     4,
     7,
     10,
     13,
     24,
     45,
     66,
     69,
     102
    ],
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# do zscore to flag outliers, remove first and last if missing data, \n",
    "def interpolate_gdv_indices_wet_and_dry(ds, params):\n",
    "        \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_interpolate_idx_params():\n",
    "        print('No cleaning/interpolation parameters set. Please ensure they are set above.')\n",
    "        return ds\n",
    "    elif len(ds['time.year']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    elif len(ds.groupby('time.month')) != 2:\n",
    "        print('Num of months of data in dataset not equal to 2. Please include more years in your query above.')\n",
    "        return ds\n",
    "        \n",
    "    # do zscore and set any years to nan if it contains a flagged outlier\n",
    "    if params.zscore:\n",
    "        ds = stats.zscore_flag_show_reset_wet_dry(ds=ds, z_val=params.zscore_value, reset=True)\n",
    "    else:\n",
    "        print('Skipping Z-score test, not requested by user.')\n",
    "        \n",
    "    # drop first year if any months within are (excluding all nan dates), else skip. same for final year, except drop if any all nan too.\n",
    "    ds = interps.drop_incomplete_first_year_wet_dry(ds=ds)\n",
    "    ds = interps.drop_incomplete_final_year_wet_dry(ds=ds)\n",
    "    \n",
    "    # basic checks\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # flag whole years if any all nan months exist\n",
    "    ds = interps.flag_years_with_missing_wet_or_dry(ds=ds)\n",
    "\n",
    "    # interpolate and extrapolate, or drop missing data\n",
    "    if params.miss_data_handle == 'interpolate':\n",
    "        ds = interps.interp_wet_dry(ds=ds)\n",
    "    elif params.miss_data_handle == 'drop':\n",
    "        print('Warning, you have chosen to drop missing data. This may result in significant data loss. Deleting data.')\n",
    "        ds = ds.dropna('time', how='all')\n",
    "        \n",
    "    # if user requested, fill missing pixels\n",
    "    if params.fill_pixels:\n",
    "        ds = interps.fill_pixels_wet_dry(ds=ds)\n",
    "\n",
    "    # remove above dashboard and map\n",
    "    print('\\nDisabling above dashboard. Please re-run the tool if you require a different index.\\n')\n",
    "    controls.disable_controls(interpolation_dash)\n",
    "        \n",
    "    # calculate some summary stats and do quick z-score test to help user\n",
    "    stats.calc_summary_stats_wet_and_dry(ds=ds)\n",
    "    \n",
    "    # show graphs\n",
    "    fig_interp_idx_stats = plotting.plot_raw_idx_stats(ds=ds, fignum=2)\n",
    "    \n",
    "    # plot figs if requested\n",
    "    if params.plot_interp_idx:\n",
    "        fig_interp_idx, fig_interp_mst = plotting.plot_raw_idx(ds=ds, stretch_min=0.01, stretch_max=0.99)\n",
    "\n",
    "    return ds\n",
    "\n",
    "ds = interpolate_gdv_indices_wet_and_dry(ds=ds, params=prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01G: Set invariant target sites parameters prior to standardisation\n",
    "The following control panel is used to set vegetation and moisture standardisation parameters. It is intended to iterative with step 01H below.\n",
    "\n",
    "Broad scale time-series monitoring of vegetation and moisture cover depends on relationships between field observations of cover and the spectral information measured by satellite imagery sensors. Variation in measured spectral information caused by factors other than variation in vegetation and moisture cover on the ground can reduce the accuracy at which vegetation cover estimates can be made from the imagery. Reduction of these variations is carried out via standardisation. \n",
    "\n",
    "Here, vegetation and moisture indices are standardised using an 'invariant target site' approach. This process involves determining the most green, moist and steady (i.e. consistently un-changing) image pixels across all dates within the dry season. The most green and moist pixels are obtained from vegetation/moisture median baseline images with values above the percentile set by the user below. Steadiness is determined from a vegetation/moisture trend slope images calculated via orthogonal polynomial regression (Snedecor 1956) for dry season images. Steadiest pixels can be obtained where slope equals (or is close to) 0. Finally, these invariant target sites, along with fuzzy sigmoidal functions (Robertson et al. 2004), are used to standardise and re-scale image values to a range from 0-1 to make all layers commensurate.\n",
    "\n",
    "Users have a lot of flexability when generating invariant target sites. The control parameters on the below control panel are:\n",
    "* `greennest/moistest percentile`: the percentile value used to threshold vegetation and moisture images into 'greennest' and 'moistest' sites. Recommended: 95-99%.\n",
    "The most green and moist pixels across all dry season images are required for standardisation. These pixels are derived from vegetation/moisture median baseline images that have a vegetation/moisture value greater than the percentile value set here. For example, a percentile value of 99.5% will designate vegetation pixels as 'greennest' where they are higher than 99.5% of the rest of the image. \n",
    "* `steadiest slope percentile`: the percentile value used to threshold the orthogonal polynomial regression slope into 'steadiest' pixels. Recommended: 4-6%.\n",
    "The steadiest pixels derived from the orthogonal polynomial regression slope (i.e. slope values equal or close to 0) are required for standardisation. For example, a percentile of 5% will select slope pixels that within 5% either side (i.e. +/- values) of 0 (i.e. the closest values to 0). These are the steadiest pixels across all time. \n",
    "* `num max sites`: the maximum number of invariant sites to be utlimately used in standardisation. Recommended: 25-75.\n",
    "Occasionally too many invariant target sites are generated from the above values. Too many invariant target sites may reduce the seperation of high and very high vegetation/moisture values. As a general rule of thumb, approximately 50 invariant targets produces realistic results for moderate study area sizes (i.e. mine site scale). Reduce to around 25 sites for smaller study areas (sub-mine site size).\n",
    "* `show figures`: Show a slide of various images used in the invariant target site and standardisation process.\n",
    "Setting this to Yes will print out a figure of various images used in the invariant target site and standardisation. Useful to turn off when you are familiar with the study area.\n",
    "\n",
    "Tips:\n",
    "* For greennest/moistest/steadiest slope percentile, increase by a % if too many pixels are returned, reduce by a % if too few are returned.\n",
    "* For steadiest slope percentile, decrease by a % if too many pixels are returned, increase by a few % if too few are returned.\n",
    "* Attempt to achieve approximately 25 sites for smaller study areas (sub-mine site scale) and around 50 sites for moderate sized study areas (mine site and surrounds scale).\n",
    "* This step and step 01H below are iterative. If you do not obtain desirable results in step 01H, change the control panel parameters here and re-run the code cell at step 01H.\n",
    "\n",
    "Run the code cell below to enable the control panel. Set your required parameters. Proceed to step 01H. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build clean and interpolate vegetation and moisture index dashboard\n",
    "def build_prepare_standardisation_dash(params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params:\n",
    "        print('Likelihood parameters not set. Did you run the above cells?')\n",
    "        return None\n",
    "        \n",
    "    # build dash\n",
    "    prepare_standardisation_dash = controls.build_prepare_standardisation_dashboard(params=params)\n",
    "    \n",
    "    # return dash\n",
    "    return prepare_standardisation_dash\n",
    "    \n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_like.reset_prepare_standardisation_params()\n",
    "\n",
    "# get indices dash\n",
    "prepare_standardisation_dash = build_prepare_standardisation_dash(prms_like)\n",
    "display(prepare_standardisation_dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01H: Generate and visualise invariant target sites parameters prior to standardisation\n",
    "Generate and visualise the invariant target sites generated using the parameters set in the above control panel. These invariant target sites will be used to apply standardisation in the next step (step 01I). The resulting figures will assist you in determining if the number and distribution of invariant target sites is optimal. See the tips in step 01G above for advice on what to aim for (i.e. typically 25-75 sites).\n",
    "\n",
    "Simply run the code cell beneath these notes to generate and visualise invariant target sites.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. All time median vegetation and moisture index images for the dry season generated;\n",
    "2. Obtain greennest/moistest pixels where vegetation/moisture values greater than user-defined 'greennest/moistest percentile';\n",
    "3. Perform orthogonal polynomial regression to obtain trend slope across all dry season images;\n",
    "4. Obtain steadiest slope pixels where slope values less than/greater than user-defined 'steadiest slope percentile';\n",
    "7. Create invariant target site mask from above outputs;\n",
    "8. Reduce number of invariant target sites to maximum number allowed (or skip if less than already); and\n",
    "8. Show various figures of main stages of invariant target site generation process (if requested)\n",
    "\n",
    "Run the code cell below to generate and visualise invariant target sites indices. This shouldn't take long. If too few invariant target sites were returned, modify the parameters in the control above and re-run. When adequate, Proceed to step 01I."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# prepare invariant green/moist target standardisation parameters and sites\n",
    "def prepare_standardisation_gdv_indices_dry(ds, params):\n",
    "        \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_prepare_standardisation_params():\n",
    "        print('No standardisation parameters set. Please ensure they are set above.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return None\n",
    "    elif len(ds.groupby('time.month')) != 2:\n",
    "        print('Num of months of data in dataset not equal to 2. Please include more years in your query above.')\n",
    "        return None\n",
    "            \n",
    "    # tell user\n",
    "    print('Generating invariant target sites for vegetation/moisture for dry season. Please wait.\\n') \n",
    "    \n",
    "    # get median all time for veg, mst, wet, dry and display\n",
    "    med_alltime_dry = ds.where(ds['time.month'] > 6, drop=True).median('time')   \n",
    "    \n",
    "    # get n% percentile of veg, mst, create igt masks, give user basic info and display\n",
    "    print('Getting highest valued vegetation and moisture pixels. Please wait.') \n",
    "    mask_green = stats.mask_via_percentile(da=med_alltime_dry['veg_idx'], percentile_val=params.green_moist_percentile)\n",
    "    mask_moist = stats.mask_via_percentile(da=med_alltime_dry['mst_idx'], percentile_val=params.green_moist_percentile)\n",
    "    \n",
    "    # get opc slope for veg, wet and display\n",
    "    print('\\nGenerating orthogonal polynomial coefficient and slope for vegetation and moisture data. Please wait.') \n",
    "    slope_veg = stats.build_opc_slope_dry(ds, var='veg_idx')\n",
    "    slope_mst = stats.build_opc_slope_dry(ds, var='mst_idx')\n",
    "    \n",
    "    # check if data valid\n",
    "    if not slope_veg.any() or not slope_mst.any():\n",
    "        print('No slope dataset returned. Aborting.')\n",
    "        return None\n",
    "    \n",
    "    # get steadiest veg, mst from slope pixels and display\n",
    "    print('> Getting steadiest parts of slope for vegetation and moisture data. Please wait.') \n",
    "    mask_steady_veg = stats.get_steadiest_opc_slope(da=slope_veg, slope_steadiness_perecentile=params.slope_steady_percentile)\n",
    "    mask_steady_mst = stats.get_steadiest_opc_slope(da=slope_mst, slope_steadiness_perecentile=params.slope_steady_percentile)\n",
    "    \n",
    "    # check if data valid\n",
    "    if not mask_steady_veg.any() or not mask_steady_mst.any():\n",
    "        print('No slope mask dataset returned. Aborting.')\n",
    "        return None\n",
    "    \n",
    "    # get veg, mst invariant sites from masks and plot\n",
    "    print('\\nBuilding invariant green/moist sites. Please wait.') \n",
    "    mask_invariant_greenest = stats.make_invariant_sites(mask_highest=mask_green, mask_steadiest=mask_steady_veg)\n",
    "    mask_invariant_moistest = stats.make_invariant_sites(mask_highest=mask_moist, mask_steadiest=mask_steady_mst)\n",
    "\n",
    "    # reduce invariant veg, mst sites down to user defined if < user defined, plot\n",
    "    print('\\nBuilding reduced invariant green/moist sites. Please wait.') \n",
    "    mask_max_invariant_greenest = stats.reduce_invariant_sites_mask(da=med_alltime_dry['veg_idx'], mask=mask_invariant_greenest, max_num_sites=params.max_num_sites)\n",
    "    mask_max_invariant_moistest = stats.reduce_invariant_sites_mask(da=med_alltime_dry['mst_idx'], mask=mask_invariant_moistest, max_num_sites=params.max_num_sites)\n",
    "    \n",
    "    # get mask where nan are - sigmoidal converts to 1 so must mask back out\n",
    "    mask_nan = xr.where(np.isnan(ds), True, False).max('time')\n",
    "    slope_veg = slope_veg.where(~mask_nan['veg_idx'])\n",
    "    slope_mst = slope_mst.where(~mask_nan['mst_idx'])\n",
    "    \n",
    "    if params.plot_standard_idx:\n",
    "        print('Preparing various standardisation figures. Please wait.') \n",
    "        fig_med_alltime_idx, axs_med_alltime_idx = plotting.plot_median_alltime_idx_dry(med_alltime_dry, stretch_min=0.001, stretch_max=0.999, fignum=6)\n",
    "        fig_med_mask_idx, axs_med_mask_idx = plotting.plot_mask_alltime_idx_dry(mask_green, mask_moist, fignum=7)\n",
    "        fig_opc_slope_idx, axs_opc_slope_idx  = plotting.plot_opc_slope_idx_dry(slope_veg=slope_veg, slope_mst=slope_mst, fignum=8)\n",
    "        fig_steady_slope_idx, axs_steady_slope_idx = plotting.plot_mask_steady_slope_dry(mask_steady_veg, mask_steady_mst, fignum=9)\n",
    "        fig_invariant_targets_idx, axs_invariant_targets_idx = plotting.plot_invariant_targets_dry(mask_invariant_greenest, mask_invariant_moistest, fignum=10)\n",
    "        fig_reduced_invariants_idx, axs_reduced_invariants_idx = plotting.plot_reduced_invariant_targets_dry(mask_max_invariant_greenest, mask_max_invariant_moistest, fignum=11)\n",
    "    \n",
    "    # tlel user and return\n",
    "    print('\\nSuccessfully generated invariant green/moisture sites for standardisation.\\n') \n",
    "    return mask_invariant_greenest, mask_invariant_moistest\n",
    "    \n",
    "mask_inv_veg, mask_inv_mst = prepare_standardisation_gdv_indices_dry(ds=ds, params=prms_like)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01I: Perform vegetation and moisture standardisation using invariant target sites\n",
    "Perform vegetation and moisture index standardisation using the invariant target sites generated in step 01G and step 01H above. This process makes use of increasing fuzzy sigmoidal functions (Robertson et al. 2004) to re-scale images to a new range of 0-1 to ensure images commensurate across seasons and years.\n",
    "\n",
    "Simply run the code cell beneath these notes to perform standardisation. Note: if you selected 'show figures' in the control panel above (step 01G), a slide of figures will show in the results below too.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. For each season of each year, obtain low and high inflection points for vegetation and moisture within invariant target sites; and\n",
    "2. Re-scale values using increasing fuzzy sigmoidal function;\n",
    "3. Show various standardised vegetation and moisture resulting from the process (if requested)\n",
    "\n",
    "Run the code cell below to perform standardisation. This shouldn't take long. Proceed to step 01J."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# standardise veg, mst da to dry season using invariant targets, low and high inflection. used in apply func\n",
    "def do_gdv_standard_indices(ds, mask_inv_veg, mask_inv_mst, low_percentile, high_percentile, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not mask_inv_veg.any() or not mask_inv_mst.any():\n",
    "        print('No invariant target masks exist or all empty. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return None\n",
    "    elif len(ds.groupby('time.month')) != 2:\n",
    "        print('Num of months of data in dataset not equal to 2. Please include more years in your query above.')\n",
    "        return None\n",
    "            \n",
    "    # tell user\n",
    "    print('Standardising raw vegetation and moisture indices using invariant targets and increasing sigmoidal function. Please wait.') \n",
    "    \n",
    "    # get mask where nan are - sigmoidal converts to 1 so must mask back out\n",
    "    mask_nan = xr.where(np.isfinite(ds), 1, 0)\n",
    "    \n",
    "    # get lowest veg, mst inflection from median all time dry\n",
    "    low_inflect_veg = np.nanpercentile(ds['veg_idx'].median('time'), q=low_percentile)\n",
    "    low_inflect_mst = np.nanpercentile(ds['mst_idx'].median('time'), q=low_percentile)\n",
    "    \n",
    "    # perform increase sigmoidal on every image in ds\n",
    "    ds = ds.groupby('time').apply(stats.perform_increase_sigmoid, mask_veg=mask_inv_veg,  mask_mst=mask_inv_mst, \n",
    "                                  high_percentile=high_percentile, low_inflect_veg=low_inflect_veg, low_inflect_mst=low_inflect_mst)\n",
    "    \n",
    "    # bring back those nan areas based on above mask. if no nan were wanted these wouldnt exist by now\n",
    "    ds = ds.where(mask_nan)\n",
    "    \n",
    "    if params.plot_standard_idx:\n",
    "        print('Preparing standardised indices figures. Please wait.') \n",
    "        fig_standard_idx, fig_standard_mst = plotting.plot_standard_idx(ds=ds, stretch_min=0.0001, stretch_max=0.9999)\n",
    "\n",
    "    # tell user and return\n",
    "    print('Successfully standardised raw vegetation and moisture indices.') \n",
    "    print('\\nDisabling above dashboard. Please re-run the tool if you require different invariant target sites.')\n",
    "    controls.disable_controls(prepare_standardisation_dash)\n",
    "    return ds\n",
    "    \n",
    "ds = do_gdv_standard_indices(ds=ds, mask_inv_veg=mask_inv_veg, mask_inv_mst=mask_inv_mst, low_percentile=5, high_percentile=95, params=prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01J: Generate and standardise vegetation and moisture stability\n",
    "Generate additional vegetation and stability model input layers that underpin the spatial multicriteria evaluation modelling techniques used to derive GDV likelihood layers using standardised vegetation and moisture images developed above. \n",
    "\n",
    "GDV has been found to maintain certain characteristics that permit its detection from remotely sensed imagery, including perennial high greenness and perennial moist foliage (O’Grady et al. 2011; Barron et al. 2014). As we have already generated standardised vegetation and moisture, persistent green and moisture can simply be obtained by subtracting wet season vegetation/moisture from dry season standardised images. Stability, which improves the likelihood of a GDV-site when coupled with high vegetation and moisture layers, can be inferred from small differences (ideally 0) between these two extremes. A increasing/decreasing fuzzy sigmoidal function (Robertson et al. 2004) is also applied to re-scale values from -1 to 1 to 0-1 to make stability images commensurate across years.\n",
    "\n",
    "This step does not require user input and therefore presents no control panel or control parameters. \n",
    "\n",
    "Simply run the code cell beneath these notes to generate standardised vegetation and moisture stability. Note: if you selected 'show figures' in the control panel above (step 01G), a slide of figures will show in the results below too.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Subtract standardised vegetation and moisture images in the wet season from the dry season each season of each year;\n",
    "2. Re-scale values from a range of -1 to 1 to 0 to 1 using increasing/decreasing fuzzy sigmoidal function; and\n",
    "3. Show standardised vegetation and moisture stability from the process (if requested)\n",
    "\n",
    "Run the code cell below to generate standardised vegetation and moisture stability layers. This shouldn't take long. Proceed to step 01K."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def do_gdv_standard_stability(ds, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return None\n",
    "    elif len(ds.groupby('time.month')) != 2:\n",
    "        print('Num of months of data in dataset not equal to 2. Please include more years in your query above.')\n",
    "        return None\n",
    "        \n",
    "    # tell user\n",
    "    print('Generating vegetation and moisture stability. Please wait.') \n",
    "            \n",
    "    # generate stability (wet - dry) across each year, populate stability list with das\n",
    "    stability_list = []\n",
    "    ds.groupby('time.year').apply(stats.calc_seasonal_stability, stability_list=stability_list)\n",
    "    \n",
    "    if len(stability_list) == 0:\n",
    "        print('No stability data generated. Aborting.')\n",
    "        return None\n",
    "        \n",
    "    # concat all arrays into dataset and rename vars\n",
    "    stability = xr.concat(stability_list, dim='time').sortby('time')\n",
    "    stability = stability.rename({'veg_idx': 'veg_stb', 'mst_idx': 'mst_stb'})\n",
    "    \n",
    "    # plot\n",
    "    #if params.plot_standard_idx:\n",
    "        #fig_stability_veg, fig_stability_mst = plotting.plot_raw_stability_idx(ds=stability)\n",
    "    \n",
    "    # standardise stability using increase/decreasing sigmoidal, rescale -1to1 to 0to2\n",
    "    print('Standardising vegetation and moisture stability. Please wait.') \n",
    "    stability = stability + 1\n",
    "    stability = stability.groupby('time').apply(stats.perform_inc_dec_sigmoid)\n",
    "\n",
    "    # mask out low veg, mst areas (these come back as very stable)\n",
    "    stability = stats.mask_seasonal_stability(stability=stability, ds=ds, percentile=90)\n",
    "    \n",
    "    if params.plot_standard_idx:\n",
    "        fig_stand_stable_veg, fig_stand_stable_mst = plotting.plot_stand_stable_idx(ds=stability)\n",
    "        \n",
    "    # tell user and return\n",
    "    print('Successfully standardised raw vegetation and moisture indices.')      \n",
    "    return stability\n",
    "    \n",
    "stability = do_gdv_standard_stability(ds=ds, params=prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01K: Weight vegetation, moisture and stability images and generate GDV likelihood\n",
    "Apply pre-developed weights to vegetation, moisture and stability images and generate annual GDV likelihood model layers.\n",
    "\n",
    "Weights for each vegetation, moisture and stability image (i.e. evidence layers) have been developed using pairwise comparison in the context of the Analytical Hierarchy Process (Saatay 1980). These weights are provided in the table here:\n",
    "\n",
    "| Variable | Season | Weight |\n",
    "| --- | --- | --- |\n",
    "| Vegetation | Dry | 0.362481346 |\n",
    "| Vegetation | Wet | 0.162613224 |\n",
    "| Vegetation Stability | N/A | 0.057825868 |\n",
    "| Moisture | Dry | 0.284883443 |\n",
    "| Moisture | Wet | 0.104621642 |\n",
    "| Moisture Stability | N/A | 0.027574478 |\n",
    "\n",
    "These weights are applied to each standardised vegetation, moisture and stability layer and summed for each year. Thus, annual GDV likelihood models are produced, each on a continuous scale of 0-1, with 0 representing locations highly unlikely to be GDV and 1 representing high likelihood. \n",
    "\n",
    "This step does not require user input and therefore presents no control panel or control parameters. \n",
    "\n",
    "Simply run the code cell beneath these notes to generate annual GDV likelihood layers. Note: if you selected 'show figures' in the control panel above (step 01G), a slide of figures will show in the results below too.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Apply pairwise comparison-derived weighting values to each vegetation, moisture and stability image for each year;\n",
    "2. Sum the weight evidence layers for each year; and\n",
    "3. Show GDV likelihood layers from the process (if requested)\n",
    "\n",
    "Run the code cell below to generate annual GDV likelihood layers. This shouldn't take long. Proceed to step 01L."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def do_gdv_likelihood(ds, stability, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds or not stability:\n",
    "        print('Error occured, no datasets exists. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) <= 2 or len(stability['time']) <= 2:\n",
    "        print('Less than three years of data in datasets. Please include more years in your query above.')\n",
    "        return ds\n",
    "    elif len(ds.groupby('time.month')) != 2:\n",
    "        print('Num of months of data in dataset not equal to 2. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # tell user\n",
    "    print('Generating groundwater dependent vegetation likelihood maps. Please wait.') \n",
    "    \n",
    "    # generate stability (wet - dry) across each year, populate stability list with das\n",
    "    likelihood_list = []\n",
    "    ds.groupby('time.year').apply(stats.generate_likelihood_map, stability=stability, likelihood_list=likelihood_list, params=params)\n",
    "    \n",
    "    if len(likelihood_list) == 0:\n",
    "        print('No likelihood data generated. Aborting.')\n",
    "        return None\n",
    "    \n",
    "    # concat all arrays into dataset and rename vars\n",
    "    likelihood = xr.concat(likelihood_list, dim='time').sortby('time')\n",
    "    likelihood = likelihood.rename('likelihood').to_dataset()\n",
    "    \n",
    "    # tell user\n",
    "    print('Successfully generated groundwater dependent vegetation likelihood maps.') \n",
    "    \n",
    "    if params.plot_standard_idx:\n",
    "        fig_likelihood = plotting.plot_likelihood(ds=likelihood, stretch_min=0, stretch_max=1)\n",
    "    \n",
    "    return likelihood\n",
    "\n",
    "likelihood = do_gdv_likelihood(ds=ds, stability=stability, params=prms_like) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01L: Generate potential bushfire areas (optional)\n",
    "Generate potential bushfire (and infrastructure) areas across years to assist in explaining exceptionally low/high areas or years of GDV likelihood (if they occur above). This step is completely optional and can be skipped completely if burn areas are of no importance to your analysis.\n",
    "\n",
    "The process uses the nor normalised burn ratio (NBR), which is applied to annual mean satellite images. Following that, the delta NBR (∆NBR) is calculated to estimate burn severity using pre- and post-fire NBR images. These pre- and post-fire images are simply a particular year of NBR subtract the previous year NBR. An overview of this process is offered by Keeley (2009).\n",
    "\n",
    "This step does not require user input and therefore presents no control panel or control parameters. \n",
    "\n",
    "Simply run the code cell beneath these notes to generate potential burn areas for each GDV likelihood layer. \n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Fetch all available satellite imagery for the date range set way back in step 01A and 01B;\n",
    "1. Calculate normalised burn ratio (NBR) index for each image for each year;\n",
    "3. Resample data down to annual NBR image means;\n",
    "2. For each year in NBR dataset, minus the previous year to calculate delta NBR (i.e. bush fire severity; and\n",
    "3. Plot these results over GDV likelihood images.\n",
    "\n",
    "Run the code cell below to generate potential bush fire area layers. This shouldn't take long. Proceed to step 01M."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# fetch cube data for normalised burn ratio bushfire and infrastructure mapping\n",
    "def fetch_nbr_cube_data(likelihood, params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params.check_general_params():\n",
    "        print('Not all required parameters set above. Please choose study area and satellite parameters.')\n",
    "        return None\n",
    "    elif helpers.valid_year_range(params.query_dates) < 3:\n",
    "        print('Date range is too small, three or more years required. We suggest >= 10 years. Please adjust above.')\n",
    "        return None \n",
    "\n",
    "    # fetch data based on user parameters\n",
    "    nbr = data.fetch_cube_data(platform=params.platform, product=params.product, q_dates=params.query_dates, cloud_cover=params.cloud_cover, \n",
    "                               cloud_mask=params.cloud_mask, coords_nw=params.alb_nw, coords_se=params.alb_se)\n",
    "    \n",
    "    if nbr:\n",
    "        with warnings.catch_warnings():\n",
    "            warnings.simplefilter(\"ignore\", category=RuntimeWarning)\n",
    "        \n",
    "            # calculate nbr and resample to annual means\n",
    "            nbr = indices.calc_nbr_index(ds=nbr, index='nbr')\n",
    "            nbr = indices.create_resample_nbr_annual(ds=nbr)\n",
    "\n",
    "            # masking out non high intensity fire\n",
    "            nbr = xr.where(nbr > 0.1, 1, np.nan)\n",
    "\n",
    "            # apply majority filter to remove salt n pepper\n",
    "            nbr = stats.apply_major_filter(ds=nbr)\n",
    "\n",
    "            # convert to mask of 1 else nan\n",
    "            nbr = nbr.where(np.isfinite(nbr))\n",
    "\n",
    "            # drop any dates where all nan\n",
    "            nbr = nbr.dropna(dim='time', how='all')\n",
    "            \n",
    "            # transform year string to datetime\n",
    "            like = likelihood.copy(deep=True)\n",
    "            like['time'] = like['time'].data.astype(np.datetime64)\n",
    "            \n",
    "            # get datasets for plotting (annual like, nbr means)\n",
    "            like_burn_years = like['likelihood'].where(like['time'] == nbr['time']).groupby('time.year').mean()\n",
    "            nbr_burn_years = nbr['nbr_idx'].groupby('time.year').mean()\n",
    "\n",
    "            # generat bush fire figure\n",
    "            print('Preparing bush fire and infrastructure figure. Please wait.')\n",
    "            fig_nbr_idx = plotting.plot_nbr_idx(like_years=like_burn_years, nbr_years=nbr_burn_years)\n",
    "        \n",
    "            # clear memory\n",
    "            nbr = None\n",
    "                        \n",
    "fetch_nbr_cube_data(likelihood=likelihood, params=prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01M: Set threshold method for GDV likelihood layers\n",
    "The following control panel is used to set parameters for GDV likelihood presence/absence thresholding. It is intended to be iterative with step 01N below.\n",
    "\n",
    "Although the GDV likelihood layers generated above can indeed be quite insightful in their current format, clarity can be improved by seperating GDV and non-GDV areas via thresholding techniques. Here, thresholding can be applied to dichotomise GDV locations from non-GDV locations so as to formulate contiguous areas. This process is accomplished in one of two ways including:\n",
    "1. ROC-slicing (e.g. Robinson et al. 2010) using known GDV locations captured via field surveys; or\n",
    "2. application of standard deviation calculation to select high GDV likelihood above a certain threshold.\n",
    "\n",
    "The ROC-slicing approach is limited in that it requires a shapefile of GDV tree species (e.g. <i>Eucalyptus victrix</i>) locations. This data must be uploaded to JupyterLab in a specific format (see below) and then selected in the control below. Regardless, this approach will provide an output that contains a ROC graph showing GDV likelihood model true positive/negative and area under the curve (AUC) metrics, which is obviously beneficial. If a shapefile is unavailable, thresholding via standard deviation has worked well in our tests.\n",
    "\n",
    "Users have a lot of flexability when thresholding GDV likelihood. The control parameters on the below control panel are:\n",
    "* `select groundtruth points`: Select a shapefile of ground-truthed GDV tree species locations.\n",
    "The locations in this shapefile will be used to perform ROC-slicing to dichotamise GDV likelihood into presence/absence and provide an indicator of model accuracy. Click 'Select', navigate to the file (must end in .shp), then click 'Change'. Please ensure your shapefile adheres to the rules listed below. Note: ignore this control if you would rather use standard deviation to threshold likelihood layers.\n",
    "* `threshold type`: Set the thresholding method to apply to the GDV likelihood layers.\n",
    "If 'groundtruth shapefile' is selected, a shapefile must be provided and ROC-slicing will be performed. If standard deviation is selected, ensure threshold date and standard deviation  slider value is correct.\n",
    "* `threshold date`: Select a GDV likelihood layer year to threshold and display on interactive map below. \n",
    "If thresholding is done via standard deviation, setting this to 'All' is recommended to reduce seasonal fluctuations. If a groundtruth shapefile is provided, it is recommended that the year the data was collected from is selected instead. \n",
    "* `standard deviation`: Set the standard deviation value which will be used to dichotomise GDV likelihood in presence/absence. Recommended: 1-2.\n",
    "Higher standard deviation values will return less thresholded GDV likelihood pixels. Conversely, lower standard deviation values will return a larger amount of GDV likelihood. This value can be used to fine-tune the dichotomisation of GDV and non-GDV based on experience in your study area.\n",
    "* `show map`: Show an interactive map once thresholding is complete. Recommended: Yes.\n",
    "If set to 'Yes', an interactive map will be displayed once thresholding is finished. This will contain both a GDV likelihood layer pre-thresholding and post-thresholding from the year set in the 'threshold date' control.\n",
    "* `export folder`: Set the folder where you would like to save thresholded GDV likelihood layers. \n",
    "If a folder location is set, a netcdf (.nc) file will be exported to this location. This .nc file can be imported in Section 02 and 03 of this tool to reduce the need to always perform this section of the tool.\n",
    "\n",
    "Shapefile requirements:\n",
    "* Shapefile must be of geometry type POINT and be of file type .shp;\n",
    "* Coordinate system must be in GDA94 Australian Albers (EPSG: 3577);\n",
    "* A column containing data on GDV tree presence and absence called GDV_ACT must be included in the shapefile;\n",
    "* The GDV_ACT column must only contain presence (value = 1) and absence (value = 0) information; and\n",
    "* It is highly recommended that the shapefile represents observations from a single year (i.e. not spanning a decade).\n",
    "\n",
    "Run the code cell below to enable the control panel. Set your required parameters. Proceed to step 01N. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Disclaimer: Shapefile is no longer available in workshop material. Please provide your own or use the Standard Deviation thresholding method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build gdv likelihood threshold interface and interactive map\n",
    "def build_likelihood_threshold_dash(likelihood, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not likelihood:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params and len(likelihood['time']) <= 2:\n",
    "        print('Likelihood parameters not set and no likelihood dataset available. Did you run the above cells?')\n",
    "        return None\n",
    "    \n",
    "    # get all years from likelihood dataset as list\n",
    "    years_list = [(year, year) for year in list(likelihood['time'].values)]\n",
    "    years_list.append(('All', 'all'))\n",
    "\n",
    "    if len(years_list) == 0:\n",
    "        print('Number of likelihood dates equals 0. Did you run the above cells?')\n",
    "        return None\n",
    "    \n",
    "    # build dash\n",
    "    likelihood_threshold_dash = controls.build_likelihood_threshold_dashboard(params=params, years_list=years_list)\n",
    "    \n",
    "    # return dash and map\n",
    "    return likelihood_threshold_dash\n",
    "    \n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_like.reset_likelihood_threshold_params()\n",
    "\n",
    "# get indices dash\n",
    "likelihood_threshold_dash = build_likelihood_threshold_dash(likelihood=likelihood, params=prms_like)\n",
    "display(likelihood_threshold_dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01N: Perform GDV likelihood layer thresholding and map the result\n",
    "Perform GDV likelihood thresholding and plot the results onto an interactive map using the method selected in step 01M above. This process makes use of ROC-slicing (e.g. Robinson et al. 2010) or standard deviation-based thresholding, depending on user-choice above. If ROC-slicing is required, a shapefile of GDV tree locations must be provided.\n",
    "\n",
    "Simply run the cell below. Note: step 01M and 01N (this step) are intended to be interative. If the threshold produced below does not match expectations, feel free to alter the parameters above (especially the standard deviation or threshold date) and re-run this code cell to modify results. Note: if you selected 'show map' in the control panel above (step 01M), an interactive map will show in the results below too.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. If threshold type is set to groundtruth data, provided shapefile will be checked and prepared;\n",
    "2. ROC-slicing will be performed if shapefile provided and ROC curve will be plotted and layers thresholded using optimal cutoff;\n",
    "3. If standard deviation selected, threshold calculated and threshold performed;\n",
    "4. Map the original and thresholded GDV likelihood layer results onto interactive map (if requested)\n",
    "\n",
    "Run the code cell below to perform standardisation. This shouldn't take long. Proceed to step 01O."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def do_gdv_likelihood_threshold(likelihood, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not likelihood:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif len(likelihood['time']) <= 2:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return None\n",
    "    \n",
    "    # tell user\n",
    "    print('Attempting to threshold groundwater dependent vegetation likelihood map. Please wait.')\n",
    "    \n",
    "    # get all data or specific year\n",
    "    if params.like_thresh_date == 'all':\n",
    "        scene = likelihood.median('time')\n",
    "    elif params.like_thresh_date != 'all':\n",
    "        scene = likelihood.sel(time=params.like_thresh_date)\n",
    "        \n",
    "    # threshold likelihood dependong on user selection\n",
    "    if params.like_thresh_type == 'stdv':\n",
    "        print('\\nThresholding groundwater dependent vegetation likelihood using standard deviation. Please wait.')\n",
    "        scene_thresh = stats.thresh_via_stdv(da=scene['likelihood'], num_of_std_devs=params.like_thresh_stdv)\n",
    "        scene_thresh = scene_thresh.to_dataset()\n",
    "        \n",
    "    elif params.like_thresh_type == 'shape':\n",
    "        if params.ground_sites_path:\n",
    "            scene_thresh, shp, fpr, tpr, auc, cutoff = stats.thresh_via_shape(da=scene['likelihood'], shp_path=params.ground_sites_path)\n",
    "            \n",
    "            if scene_thresh.any():\n",
    "                scene_thresh = scene_thresh.to_dataset()\n",
    "                fig_shape_roc, ax_shape_roc = plotting.plot_shape_and_roc_stats(da=scene['likelihood'], shp=shp, tpr=tpr, fpr=fpr, auc=auc, fignum=12)\n",
    "                plt.show()\n",
    "            else:\n",
    "                print('> No data was returned following likelihood thresholding. Aborting.')\n",
    "                return None, None\n",
    "        else:\n",
    "            print('> Groundtruth Shapefile threshold method selected but no shapefile provided. Please select a shapefile and try again.')\n",
    "            return None, None\n",
    "        \n",
    "    # basic check\n",
    "    if not scene.any() or not scene_thresh.any():\n",
    "        print('No likelihood data to display. Aborting.')\n",
    "        return None, None\n",
    "        \n",
    "    # write data if selected\n",
    "    if params.write_data_path:\n",
    "        data.write_like_data_path(da_raw=scene, da_thresh=scene_thresh, out_path=params.write_data_path)\n",
    "        \n",
    "    # plot onto map if selected\n",
    "    if params.map_likelihood:\n",
    "        print('Preparing interactive map. Please wait.')\n",
    "        likelihood_map = maps.load_likelihood_map(da_raw=scene['likelihood'], da_thresh=scene_thresh['likelihood'], params=prms_like)\n",
    "        display(likelihood_map)\n",
    "    \n",
    "    return scene, scene_thresh\n",
    "\n",
    "scene_like, scene_like_threshed = do_gdv_likelihood_threshold(likelihood=likelihood, params=prms_like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 01O: Finalise GDV likelihood layer\n",
    "Finalise GDV likelihood layers and processes. This performs basic data cleaning and information transferal to new variables, etc. This step is required before Section 02 and Section 03 can be used.\n",
    "\n",
    "This step does not require user input and therefore presents no control panel or control parameters. Simply run the code cell beneath these notes to finalise GDV likelihood layers.\n",
    "\n",
    "The following basic steps are undertaken:\n",
    "1. Clear original satellite imagery and un-needed variables;\n",
    "2. Clear various mask variables; and\n",
    "3. Transfer information to new metadata variable classes.\n",
    "\n",
    "Run the code cell below to finalise GDV likelihood layers. This shouldn't take long. You may now proceed to Section 02 or Section 03."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# clean up existing datasets, variables\n",
    "def clean_up_data(prms_like, prms_trend, prms_change):\n",
    "    \n",
    "    # tell user and delete\n",
    "    print('Clearing original dataset variable (ds).')\n",
    "    ds = None\n",
    "\n",
    "    print('Clearing original stability variable (stability).')\n",
    "    stability = None\n",
    "\n",
    "    print('Clearing original mask for vegetation invariant targets (mask_inv_veg).')\n",
    "    mask_inv_veg = None\n",
    "\n",
    "    print('Clearing original mask for moisture invariant targets (mask_inv_mst).')\n",
    "    mask_inv_mst = None\n",
    "    \n",
    "    print('Transferring metadata over to trend/change parameters.')\n",
    "    prms_change.crs = prms_trend.crs = prms_like.crs\n",
    "    prms_change.affine = prms_trend.affine = prms_like.affine\n",
    "    prms_change.gcs_nw = prms_trend.gcs_nw = prms_like.gcs_nw\n",
    "    prms_change.gcs_se = prms_trend.gcs_se = prms_like.gcs_se\n",
    "    prms_change.alb_nw = prms_trend.alb_nw = prms_like.alb_nw\n",
    "    prms_change.alb_se = prms_trend.alb_se = prms_like.alb_se\n",
    "    \n",
    "    # disable likelihood dashboard\n",
    "    print('Disabling likelihood dashboard. Rerun the tool if you require a different threshold.')\n",
    "    try:\n",
    "        controls.disable_controls(likelihood_threshold_dash)\n",
    "    except:\n",
    "        print('No need to disable likelihood dashboard. Has not been initialised. Proceeding.')\n",
    "            \n",
    "    # tell user and return\n",
    "    print('Cleaning process complete.')\n",
    "    return ds, stability, mask_inv_veg, mask_inv_mst, prms_change, prms_trend\n",
    "\n",
    "ds, stability, mask_inv_veg, mask_inv_mst, prms_change, prms_trend = clean_up_data(prms_like=prms_like, prms_trend=prms_trend, prms_change=prms_change)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 02: Monitor GDV health trends via Mann-Kendall trend analysis\n",
    "The second section of this tool involves the application of several trend analyses to detect anomalous, continuously decreasing (or increasing) trends in plant health. The intent is to provide an assistive, user-driven, interactive approach to exploring GDV health trends over a specific time period and study area boundary. The resulting trend layers are coupled with graphing tools that facilitate exploration of vegetation health history anywhere within the study area (i.e. both inside and outside of GDV trend areas). Finally, change point detection (e.g. Andersen 2009) is also incorporated to allow users to detect 'break points' in vegetation health history (i.e. dates when vegetation health signatures significantly shifted away from norm).\n",
    "\n",
    "Mann-Kendall (MK) tests (Mann 1945; Hipel and McLeod 2005) are used to perform trend analysis. A MK test determines if there is a monotonic increase, monotonic decrease or no significant change over time. The MK non-parametric trend test is a function of the ranks of the observations rather than their actual values and so does not require the assumption of normally distributed data (Hamed 2008) like, for example, linear regression. Furthermore, the test allows the mapping of p-values and thus provides a confidence level for interpretation. Mann-Kendall is applied using the pymannkendall python module (Hussain et al. 2019). Two variations of the MK test are included in this tool:\n",
    "1. `original mk`: is a nonparametric test which does not consider serial correlation or seasonal effects. \n",
    "This method is most appropriate when exploring health trends for a specific season or month across years. For example, exploring trends at a monitoring site for every Janurary from 2009 to 2019. \n",
    "2. `seasonal mk`: is a nonparametric test that considers seasonal effects and is suited for time-series data with seasonal trends (Hirsch et al. 1982). This method is most appropriate when exploring vegetation trends across every season or month across years. For example, exploring trends at a monitoring site across every season beginning 2009 and ending 2020. This approach is recommended as it allows more data to be incorporated in the MK test.\n",
    "\n",
    "Change point detection is used in time-series analysis to detect abupt shifts in time series trends (i.e. shifts in a time series' instantanous velocity). This can be applied to a time-series of vegetation health observations to detect anomalous sequences/states in the vegetation health history, and to detect dates where sudden changes occurred (e.g. Andersen 2009). This tool uses two offline, exact segmentation change point algorithms (Truong et al. 2019, Truong et al. 2020) to perform change point detection:\n",
    "1. `pruned exact linear time (pelt) search method`: is an exact search method that detects break points via the minimization of costs and is generally considered to produce highly optimal and consistent results (Wambui et al. 2014). Users- determine the number of break points returned using a 'penelty' value.\n",
    "2. `dynamic programming (dynpro) search method`: is also an exact method used for automated change point detection. It has a considerable computational cost of O(Qn^2), where Q is the max number of change points and n is the number of data points (Truong et al. 2020). Generally considered less strict than PELT and, thus, returns more break points.\n",
    "\n",
    "The key processing steps applied in this section are as follows:\n",
    "1. Fetch a ‘cube’ of Landsat or Sentinel satellite data for a user-defined study area and timespan using the DEA ODC;\n",
    "2. Generate vegetation index based on user-selection;\n",
    "3. Reduce data to selected season(s)/month(s) and interpolate missing data (optional);\n",
    "4. Standardise vegetation using a basic invariant target site approach;\n",
    "5. Perform MK test;\n",
    "6. Map trend results onto interactive map; and\n",
    "7. Graph vegetation history and perform change point detection on any pixel clicked by user.\n",
    "\n",
    "It is required that Section 00 and Section 01 has been run prior to this section. Proceed to step 02A when ready."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 02A: Set trend analysis parameters\n",
    "The following control panel is used to set the various trend analysis parameters including satellite data, timespan of trend analysis, Mann-Kendall analysistype, and interpolation and standardisation approaches. Most of these parameters have been explained in Section 01 and will not be repeated here. New control parameters not yet seen are:\n",
    "* `select existing likelihood data`: Select an existing thresholded GDV likelihood dataset and import. Recommended: No.\n",
    "Users can import previously generated thresholded GDV likelihood netcdf (.nc) datasets. This can be used to eliminate need to run Section 01 every time monitoring is udnertaken.\n",
    "* `mann-kendall type`: Choose whether to do an original or seasonal Mann-Kendall test. Recommended: Seasonal.\n",
    "The original MK test is most appropriate to find trends for a specific month or season for a user-defined range of years (e.g. health trends every March across 10 years). The seasonal MK test is most appropriate when every season or month of every year is required, allowing for seasonal fluctuation to be considered in the trend. Seasonal is recommended as more data can be incorporated into the MK test, potentially improving test accuracy.\n",
    "* `time slice`: Choose the slice of time used in the MK test. Recommended: Quarters if seasonal MK, anything for original MK.\n",
    "For seasonal MK, season to season (i.e. quarter to quarter) is recommended, as less data will need to be interpolated. If month to month is selected, missing months will be interpolated, but more information will be used to detect a more accurate trend. For original MK, choose a specific season (or quarter) or month to perform the MK test across.\n",
    "* `trend dates`: Start and End dates for analysis (e.g. 2009 to 2019). End date is inclusive.\n",
    "It is recommended 10+ years is selected to provide the MK and change point detection tests sufficient data during processing. A minimum of 5 years is required, which can limit the use of Sentinel 2 in this analysis.\n",
    "* `sig trend only`: Display only significant trends (i.e. MK test pixels with p-value < 0.05) on the resulting interactive map. Recommended: No.\n",
    "If set to 'No', all trends (increasing, decrease, no trend) will be shown on the interactive map. If 'Yes', only those trends that are statistically significant will be returned.\n",
    "* `show trend`: Choose which type of trend (e.g. increasing, decreasing, all) to show on the output interactive map.\n",
    "If set to 'All', all increasing, decreasing and no trend pixels will be shown on resulting map. Otherwise, specific trend types will be returned instead.\n",
    "\n",
    "Tips:\n",
    "* Landsat satellite must be used if GDV likelihood used Landsat imagery. Likewise for Sentinel 2.\n",
    "* A minimum of 5 years is required for trend analysis. It is recommended 10+ years are used.\n",
    "* Landsat satellite is recommended here, as Sentinel 2 has yet to reach >= 5 full years of data (will be available end of 2020).\n",
    "\n",
    "Run the code cell below to enable the control panel and interactive map. Set your required parameters. Proceed to step 01B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build mannkendall trend control panel\n",
    "def build_trend_dash(params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params:\n",
    "        print('Trend parameters not set. Did you run the above cells?')\n",
    "        return None, None\n",
    "       \n",
    "    # build dash and map\n",
    "    trend_dash = controls.build_trend_dashboard(params=params)\n",
    "    \n",
    "    # return dash and map\n",
    "    return trend_dash\n",
    "\n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_trend.reset_trend_params()\n",
    "\n",
    "# get sat fetch dash and map\n",
    "trend_dash = build_trend_dash(params=prms_trend)\n",
    "display(trend_dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 02B: Perform trend analysis\n",
    "Perform MK trend analysis on thresholded GDV likelihood areas and plot the results onto an interactive map using the parameters selected in step 02A above. \n",
    "This process makes use of one of two Mann-Kendall trend analysis techniques, as set by the user.\n",
    "\n",
    "Simply run the cell below. Note: step 02A and 02B (this step) are intended to be interative. If the MK test result produced below does not match expectations, feel free to alter the parameters above (especially the MK test type,  time slice and trend dates) and re-run this code cell to modify results.\n",
    "\n",
    "The basic steps are undertaken during this process are outlined in the Section 02 notes above.\n",
    "\n",
    "Run the code cell below to perform MK trend analysis. This can take quite awhile to complete, so please wait. Proceed to step 02C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# prepare for mannkendall trend analysis!\n",
    "def prepare_for_trend_analysis(scene_like_threshed, params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params.check_general_params():\n",
    "        print('Not all required parameters set above. Please choose study area and satellite parameters.')\n",
    "        return None\n",
    "    elif helpers.valid_year_range(params.query_dates) < 5:\n",
    "        print('Date range is too small, five or more years required. We suggest >= 10 years. Please adjust above.')\n",
    "        return None\n",
    "        \n",
    "    # enable xr apply progress bar (not implemented in xr... yet)\n",
    "    helpers.enable_xr_progressbar()\n",
    "\n",
    "    # fetch data based on user parameters\n",
    "    ds = data.fetch_cube_data(platform=params.platform, product=params.product, q_dates=params.query_dates, cloud_cover=params.cloud_cover, \n",
    "                              cloud_mask=params.cloud_mask, coords_nw=params.alb_nw, coords_se=params.alb_se)  \n",
    "        \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not indices.check_required_bands_exist(ds=ds):\n",
    "        print('Not all required bands available in dataset. Please re-run the tool again.')\n",
    "        return ds\n",
    "    elif not params.check_trend_idx_params():\n",
    "        print('No vegetation index set. Please ensure they are set above above.')\n",
    "        return ds\n",
    "    \n",
    "    # add an empty first date array if doesn't exist. this is important for resampling\n",
    "    ds = data.add_empty_first_date_array(ds=ds, time_slice=params.time_slice, first_query_date=params.query_dates[0])\n",
    "        \n",
    "    # reduce dataset down to specific quarter or month (original mk only, seasonal uses everything)\n",
    "    if params.mk_type == 'original':\n",
    "        print('\\nOrginal mannkendall selected, restricting to chosen quarter/month, hold tight...')\n",
    "        ds = data.reduce_to_specific_month_or_season(ds=ds, time_slice=params.time_slice)\n",
    "        \n",
    "    # create vegetation indices and resample to specified time slice\n",
    "    ds = indices.create_resample_idxs_trend_time_slice(ds=ds, params=params)\n",
    "            \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) < 5:\n",
    "        print('Less than five years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # calculate some summary stats and graph\n",
    "    stats.calc_summary_stats_trend(ds=ds)\n",
    "    fig_raw_idx_trend_stats = plotting.plot_raw_trend_stats(ds=ds, fignum=13)\n",
    "    plt.show()\n",
    "\n",
    "    # do z-score if time slice is original mannkendall (seasonal likes outliers, so ignore if seasonal)\n",
    "    if params.mk_type == 'original' and params.zscore:\n",
    "        ds = stats.zscore_flag_show_reset_trend(ds=ds, z_val=params.zscore_value, reset=True)\n",
    "    else:\n",
    "        print('Skipping z-score outlier test for seasonal mann-kendall. Seasonal mann-kendall will consider outliers.')\n",
    "        \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_trend_interpolate_idx_params():\n",
    "        print('No cleaning/interpolation parameters set. Please ensure they are set above.')\n",
    "        return ds\n",
    "    elif len(ds['time.year']) < 5:\n",
    "        print('Less than three years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "\n",
    "    # interpolate and extrapolate, or drop missing data\n",
    "    if params.miss_data_handle == 'interpolate':\n",
    "        if params.mk_type == 'original':\n",
    "            ds = interps.interp_original_trend(ds=ds)\n",
    "        elif params.mk_type == 'seasonal':\n",
    "            ds = interps.interp_seasonal_trend(ds=ds)\n",
    "    elif params.miss_data_handle == 'drop':\n",
    "        print('Warning, you have chosen to drop missing data. This may result in significant data loss. Deleting data.')\n",
    "        ds = ds.dropna('time', how='all')\n",
    "        \n",
    "    # if user requested, fill missing pixels\n",
    "    if params.fill_pixels:\n",
    "        ds = interps.fill_pixels_trend(ds=ds)\n",
    "    \n",
    "    # calculate some summary stats and show graphs\n",
    "    stats.calc_summary_stats_trend(ds=ds)   \n",
    "    fig_raw_idx_trend_clean_stats = plotting.plot_raw_trend_stats(ds=ds, fignum=14)\n",
    "    plt.show()\n",
    "    \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_trend_standardisation_params():\n",
    "        print('No standardisation parameters set. Please ensure they are set above.')\n",
    "        return ds\n",
    "    elif len(ds['time.year']) < 5:\n",
    "        print('Less than five years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # standardise using basic greennest site and / max approach\n",
    "    if params.mk_type == 'original':\n",
    "        ds = stats.do_standardisation_indices_trend(ds=ds, params=params)\n",
    "    elif params.mk_type == 'seasonal':\n",
    "        print('Standardisation not required for seasonal mannkendall. Skipping and proceeding.')\n",
    "        \n",
    "    return ds\n",
    "\n",
    "\n",
    "# perform dataset mask (to gdv), stack data into vectors, and do mannkendall\n",
    "def do_mannkendall(ds, scene_like_threshed, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds or not scene_like_threshed:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_mk_trend_params():\n",
    "\n",
    "        print('Not all required parameters set above. Please choose study area and satellite parameters.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) < 5:\n",
    "        print('Less than five years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "        \n",
    "    # prepare empty mk dataset, mask out non-gdv from ds using likelihood thresh, stack into pixel vectors\n",
    "    mk_result = stats.prepare_empty_mk_dataset(scene_like_threshed=scene_like_threshed)\n",
    "    masked = stats.mask_non_gdv_data(ds=ds, scene_like_threshed=scene_like_threshed)\n",
    "    stacked = stats.stack_and_prepare_pixels(ds=masked)\n",
    "    \n",
    "    # get seasonal period if seasonal mk being undertaken\n",
    "    period = None\n",
    "    if params.mk_type == 'seasonal':\n",
    "        if params.time_slice == 'qt':\n",
    "            period = 4\n",
    "        elif params.time_slice == 'mt':\n",
    "            period = 12\n",
    "  \n",
    "    # perform mannkendall analysis on each pxiel vector, unstack once done, transfer back to original result dataset\n",
    "    unstacked = stats.perform_mk(stacked=stacked, mk_type=params.mk_type, sig_only=params.mk_sig_only, trend_type=params.trend_type, period=period)\n",
    "    mk_result = stats.transfer_unstack_to_mk_result(mk_result=mk_result, unstacked=unstacked) \n",
    "    \n",
    "    # add crs back on incase its gone\n",
    "    mk_result.attrs['crs'] = params.crs\n",
    "    \n",
    "    # tell user, return\n",
    "    print('Mann-kendall analysis successfully completed.')\n",
    "    return mk_result  \n",
    "\n",
    "# check if change map exists, disable otherwise\n",
    "if 'change_map' in locals():\n",
    "    change_map.close()\n",
    "\n",
    "# check if prior work done to create likelihood data, if not set none\n",
    "if 'scene_like_threshed' not in locals() or not scene_like_threshed:\n",
    "    scene_like_threshed, prms_trend = data.load_and_prepare_like_netcdf(params=prms_trend)\n",
    "    \n",
    "# fetch data and do mannkendall\n",
    "if scene_like_threshed and prms_like and prms_trend:\n",
    "    plat = helpers.identify_platform(da=scene_like_threshed)\n",
    "    \n",
    "    if prms_trend.platform == plat:\n",
    "        prms_trend.correct_sat_launch_date()\n",
    "        ds = prepare_for_trend_analysis(scene_like_threshed=scene_like_threshed, params=prms_trend)\n",
    "        mk_result = do_mannkendall(ds=ds, scene_like_threshed=scene_like_threshed, params=prms_trend)\n",
    "    else:\n",
    "        print('Cannot perform trend analysis using platform: {0} when likelihood used: {1}. Platforms must match.'.format(prms_trend.platform, plat))\n",
    "        print('It is recommended running Section 01 again using the matching platform or changing your above selection.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 02C: Visualise trend analysis results, GDV health graphing tools and change point detection controls\n",
    "The following step is used to plot MK test result layers onto an interactive map to facilitate trend exploration. Users can also click on any pixel within the study area boundary to graph vegetation health signatures and perform change point detection to determine dates where vegetation health significantly shifted from the norm. This step is intended to be iterative with steps 02A and 02B above.\n",
    "\n",
    "The results of MK tests provide three layers that users can explore on the interactive map:\n",
    "* `tau value layer`: The rank correlation coefficient. Indicative of trend direction (increasing, decreasing, no trend). \n",
    "A value of 1 (very green pixels) indicates a perfect, continuously increasing trend, or in our case, continously increasing vegetation health. In contrast, a Tau value of -1 (very red pixels) is a perfect, continuously decreasing trend, or in our case, continuously declining vegetation health. A Tau of 0 (yellow-ish pixels) indicates no continuous trend in any particular direction. Values can occur anywhere between -1 to 1.\n",
    "* `p-value layer`: The probability of obtaining the observed trend, as in most statistical tests. \n",
    "Typical significance levels for given hypothesis test are 0.1, 0.05, and 0.01. Whiter pixels are very significant, black/red pixels are not. As a general rule of thumb, a p-value < 0.05 is a standard for determining statistical significance. As an example here, a pixel with a declining trend of -0.8 with a p-value of 0.05 can be considered as a significant decline trend.\n",
    "* `sen's slope layer`: A more robust trend layer that ignores trends that don’t occur for at least 30% of the time series. \n",
    "Shows robust, relatively longer-term trends and will ignore a short drought, La Nina, El-Nino cycle.\n",
    "\n",
    "The Tau, P-value and Sen's slope layers act as a guide to direct users to potential GDV decline or regrowth areas. More information can be obtained about these trends by simply clicking on the pixel of interest and running a change point detection analysis on it. This is achieved using the control panel underneat the interactive map. The control panel parameters include:\n",
    "* `model`: Cost function model. Users can select between Least absolute deviation (L1), Least squared deviation (L2), or Kernalised mean change (rbf). Recommended: rbf.\n",
    "These models dictate how the function detects changes within the vegetation health signal. L1 detects changes in the median of the signal. Overall it is a robust restimator of a shift in the central point of vegetation samples. L2 is the same, but detects mean-shifts in the vegetation health signal. Rbf detects changes in the mean of the embedded signal, based on Arlot et al. (2012). The Rbf model is recommended.\n",
    "* `min break distance`: Minimum distance between break points. One unit is equal to one date. Default is 2. Recommended: 1-3.\n",
    "Smaller value will allow break points to be detcted closer together, larger value will ensure gaps exists between breaks. Helpful for differentiating between short or long-term breaks.\n",
    "* `jump value`: Sets how large a 'jump' between vegetation health regimes need to be before they are detected. Recommended: 1.\n",
    "Smaller value will allow smaller shifts in vegetation health regimes to be detected. A larger value will only return very significant shifts in health regimes.\n",
    "* `num breaks (dynpro)`: Forces the Dynamic Programming change point detection method to return a set number of breaks points. Recommended: 2-4.\n",
    "Given a time-series of vegetation health values, tell the Dynamic Programming method how many break points to return. Even numbers only (start and end of breaks = 1 regime).\n",
    "* `break penelty (pelt)`: Set the penalty value for the PELT change point detection function. Recommend: 2-4.\n",
    "The penalty value is used in the PELT function to guard against model overfitting. Can be used to detect optimal breaks through minimisation of costs. Smaller values may return more breaks, but may suffer overfitting.\n",
    "* `show button`: Used to perform a change point detection with set parameters on the last clicked pixel. \n",
    "CLick this to generate a graph of vegetation health and perform a PELT and DynPro change point detection test for current pixel.\n",
    "\n",
    "Tip:\n",
    "* For the Tau layer, red pixels indicate some magnitude of negative trend. The more red a pixel is, the more negative that trend was. This indicates vegetation continously declined date after date; \n",
    "* For the Tau layer, green pixels indicate some magnitude of positive trend. The more green a pixel is, the more positive that trend was. This indicates vegetation continously increased date after date;\n",
    "* The Tau map can be somewhat decieving if not assessed alongside the p-value. Not all trends detected will be statistically significant. Jump between the Tau and P-value map to explore where trends are actually significant (or opt to return only significant trends using the control panel above);\n",
    "* Click on any pixel to get a reading of the Tau, P-Value and Sen's slope at that particular location;\n",
    "* Click on any pixel and click the 'Show' button in the Change Point Detection control panel to perform a change point detection and see vegetation graphs;\n",
    "* The change point detection default values are recommended; and\n",
    "* The change point detection processors are not limited to GDV likelihood and trend pixels - click anywhere in your study area and click 'Show' to see potential breaks in vegetation health.\n",
    "\n",
    "Run the code cell below to generate the interactive map and change point detection control panel. This shouldn't take long. This is the final step involved with trend analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# prepare interactive map \n",
    "def build_trend_map_and_cpd_panel(ds, mk_result, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds or not mk_result:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    \n",
    "    # globals\n",
    "    trend_x = None\n",
    "    trend_y = None\n",
    "    trend_msg = None\n",
    "    trend_info = HTML(value='<h3></h3>')\n",
    "\n",
    "    # reset cpd params\n",
    "    params.reset_cpd_params()\n",
    "    \n",
    "    # capture click\n",
    "    def handle_click(**kwargs):\n",
    "        global trend_x, trend_y\n",
    "        global trend_msg\n",
    "\n",
    "        if kwargs['type'] == 'click':\n",
    "            try:\n",
    "                lon, lat = kwargs['coordinates'][1], kwargs['coordinates'][0]\n",
    "                trend_x, trend_y = helpers.transform_point(lon=lon, lat=lat, in_epsg='epsg:4326', out_epsg='epsg:3577')\n",
    "                trend_msg = stats.get_mk_pixel_info(mk_result=mk_result, x=trend_x, y=trend_y)\n",
    "                trend_info.value = '<h3>{0}</h3>'.format(trend_msg)\n",
    "                \n",
    "            except:\n",
    "                print('Problem during trend_map_and_cpd_panel click event. Skipping.')\n",
    "            \n",
    "    # on button click\n",
    "    def handle_button(event):\n",
    "        global trend_x\n",
    "        global trend_y\n",
    "\n",
    "        with output:\n",
    "            clear_output()\n",
    "\n",
    "            # do change point detection and store results\n",
    "            vec, dates, result_pelt, result_dynp = stats.do_change_pd(ds=ds, x=trend_x, y=trend_y, model=params.model, min_size=params.min_size, jump=params.jump, \n",
    "                                                                      dynp_n_bkps=params.dynp_n_bkps, pelt_penelty=params.pelt_penelty)\n",
    "\n",
    "            if (str(vec) == 'None' or str(dates) == 'None') or (len(vec) == 0 or len(dates) == 0):\n",
    "                print('No vector returned. Did you click outside the study area? Try again.')\n",
    "                return\n",
    "            \n",
    "            # plot pelt cbd\n",
    "            if result_pelt:\n",
    "                fig_mk_pelt, fig_axs_pelt = plotting.plot_cpd(signal=vec, dates=dates, result=result_pelt, cbd_type='PELT', fignum=21)               \n",
    "                display(fig_mk_pelt)\n",
    "            else:\n",
    "                display('Change point (PELT) did not detect a break. Skipping.')\n",
    "\n",
    "            # plot dynp cbd\n",
    "            if result_dynp:\n",
    "                fig_mk_dynp, fig_axs_dynp = plotting.plot_cpd(signal=vec, dates=dates, result=result_dynp, cbd_type='DynPro', fignum=22)\n",
    "                display(fig_mk_dynp)\n",
    "            else:\n",
    "                display('Change point (DynPro) did not detect a break. Skipping.')\n",
    "\n",
    "    # load interactive map and add interactivity\n",
    "    trend_map = maps.load_trend_map(ds=ds, mk_result=mk_result, params=params)\n",
    "    trend_map.on_interaction(handle_click)\n",
    "\n",
    "    # load cpd dashboard\n",
    "    trend_cpd_dash = controls.build_trend_cpd_dashboard(params=params, handle_button=handle_button)\n",
    "                        \n",
    "    # load button and output and add interactivity\n",
    "    output = Output()\n",
    "\n",
    "    # return all controls\n",
    "    return trend_map, trend_info, trend_cpd_dash, output\n",
    "\n",
    "# get interactive map and dash\n",
    "trend_map, trend_info, trend_cpd_dash, trend_output = build_trend_map_and_cpd_panel(ds=ds, mk_result=mk_result, params=prms_trend)\n",
    "\n",
    "# show controls\n",
    "display(trend_map, trend_info, trend_cpd_dash, trend_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 03: Monitor GDV health changes via spectral change vector analysis\n",
    "The third and final section of this tool uses time-series spectral change vector analysis to detect different types of change and change magnitude of plant health. The intent is to provide an assistive, user-driven, interactive approach to exploring GDV health change over a specific time period and study area boundary. The resulting change layers are coupled with graphing tools that facilitate exploration of vegetation health history anywhere within the study area (i.e. both inside and outside of GDV change areas). Finally, as with Section 02, change point detection (e.g. Andersen 2009) is incorporated to allow users to detect 'break points' in vegetation health history (i.e. dates when vegetation health signatures significantly shifted away from norm).\n",
    "\n",
    "Spectral change vector analysis (CVA) has been successfully used to derive areas of vegetation and moisture changes, such as vegetation health increase or decrease, from satellite imagery (e.g. Malila 1980). CVA calculates angle and magnitude of change vectors derived from a two-band image, usually derived from the process of Tasselled Cap Transformation (Crist 1985). It is typically used to identify spectral changes between two identical scenes which were acquired at different times. For each pixel between the two images, CVA calculates the change vector in the two-dimensional spectral space. This results in 'angle' and 'magnitude' layer outputs, the former of which has been found to represent different types of vegetation or moisture change depending on angle range. A figure demonstrating these angles ranges is presented below (image from Zanchetta and Bitelli 2017): \n",
    "\n",
    "![alt text](imgs/cva_angles.png \"CVA Angles\")\n",
    "\n",
    "This image demonstrates the angle of change (bearing of arrow) and magnitude (distance of line) calculated from the difference between two images captured on different dates (e.g. a vegetation median baseline vs 2018 vegetation). In this example, the angle has a bearing of approx. 180 degrees. Based on the image, this falls into a quadrant that is indicative of vegetation decline (i.e. bare-soil increase). In contrast, an angle between 270-359 degrees would be indicative of vegetation increase. A threshold median factor can also be set in this CVA method. This calculates a threshold magnitude for which pixels are considered stable, i.e. no change, and can be used to control the sensitivity of what magnitudes are returned.\n",
    "\n",
    "Change point detection is used in time-series analysis to detect abupt shifts in time series trends (i.e. shifts in a time series' instantanous velocity). This can be applied to a time-series of vegetation health observations to detect anomalous sequences/states in the vegetation health history, and to detect dates where sudden changes occurred (e.g. Andersen 2009). This tool uses two offline, exact segmentation change point algorithms (Truong et al. 2019, Truong et al. 2020) to perform change point detection:\n",
    "1. `pruned exact linear time (pelt) search method`: is an exact search method that detects break points via the minimization of costs and is generally considered to produce highly optimal and consistent results (Wambui et al. 2014). Users- determine the number of break points returned using a 'penelty' value.\n",
    "2. `dynamic programming (dynpro) search method`: is also an exact method used for automated change point detection. It has a considerable computational cost of O(Qn^2), where Q is the max number of change points and n is the number of data points (Truong et al. 2020). Generally considered less strict than PELT and, thus, returns more break points.\n",
    "\n",
    "The key processing steps applied in this section are as follows:\n",
    "1. Fetch a ‘cube’ of Landsat or Sentinel satellite data for a user-defined study area and timespan using the DEA ODC;\n",
    "2. Generate Tasselled Cap index;\n",
    "3. Reduce data to baseline and comparison datasets for specified season(s) or month(s);\n",
    "4. Standardise vegetation using a basic invariant target site approach;\n",
    "5. Perform spectral CVA for every comparison year against baseline data;\n",
    "6. Map trend results onto interactive map; and\n",
    "7. Graph vegetation history and perform change point detection on any pixel clicked by user.\n",
    "\n",
    "It is required that Section 00 and Section 01 has been run prior to this section. Proceed to step 03A when ready."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 03A: Set change analysis parameters\n",
    "The following control panel is used to set the various spectral change vector analysis parameters including satellite data, timespan of baseline and comparison date ranges, change vector angle and magnitude values, and interpolation and standardisation approaches. Most of these parameters have been explained in Section 01 and will not be repeated here. New control parameters not yet seen are:\n",
    "\n",
    "* `select existing likelihood data`: Select an existing thresholded GDV likelihood dataset and import. Recommended: No.\n",
    "Users can import previously generated thresholded GDV likelihood netcdf (.nc) datasets. This can be used to eliminate need to run Section 01 every time monitoring is udnertaken.\n",
    "* `baseline dates`: Set the range of dates used to form the vegetation/soil baseline in which all other vegetation/soil images will be compared to. Recommended: 20+ years.\n",
    "CVA requires a baseline image to which all comparison images are compared against. This baseline is generated from median Tasselled Cap greeness and brightness bands of all data within the date range set here. Typically, 10-20 years forms a good baseline in which current vegetation/soil can be compared to.\n",
    "* `comparison dates`: Set the range of dates that will be compared to the baseline determine change. Recommended 5+ years.\n",
    "Each date in the comparison date range will be compared seperately against the baseline image. For example, a range of dry season vegetation images from 2010-2019 will produce a change map for dry season 2010, dry season 2011, dry season 2012, etc. These will be visualised later.\n",
    "* `time slice`: Choose to perform CVA across a particular month or season.\n",
    "Choose a specific season (quarter) or month to perform the CVA on for every year in your comparison date range. For example, setting Q4 (SON, September to November) will pull satellite imagery only within those months for every year.\n",
    "* `angle range`: Set the range of CVA angles to return. Recommended: 90-180 degrees (i.e. vegetation decline).\n",
    "Different angle ranges of angles represent different change types. Based on Malila (1980), angle ranges include: 0-90deg for moisture decline, 90-180deg for vegetation decline, 180-270deg for moisture increase, and 270-360deg for vegetation increase. Set the angle range to 0-360 for all angles.\n",
    "* `magnitude threshold`: Set the median threshold factor. Recommended: 1-2.\n",
    "The magnitude threshold factor value calculates a threshold magnitude for which pixels are considered stable, i.e. no change. This allows users to control the amount of magnitude (or change) that will be returned. A lower threshold factor (e.g. 1) returns more low-magnitude (as well as high) change pixels. Conversley, a higher magnitude factor (e.g. 3) returns less low-magnitude (as well as high) change pixels.\n",
    "\n",
    "Tips:\n",
    "* Landsat satellite must be used if GDV likelihood used Landsat imagery. Likewise for Sentinel 2.\n",
    "* The earliest comparison date must not be earlier than the earliest baseline date.\n",
    "* Both baseline and comparison date range is recommended to >= 5 years, although not a requirement.\n",
    "* A good approach to designating a date range for your baseline is a recent range of dates prior to mine operation. \n",
    "* Landsat satellite is recommended here, as Sentinel 2 has yet to reach >= 5 full years of data (will be available end of 2020).\n",
    "\n",
    "Run the code cell below to enable the control panel and interactive map. Set your required parameters. Proceed to step 03B."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# build change vector analysis control panel\n",
    "def build_change_dash(params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params:\n",
    "        print('Change parameters not set. Did you run the above cells?')\n",
    "        return None, None\n",
    "       \n",
    "    # build dash and map\n",
    "    change_dash = controls.build_change_dashboard(params=params)\n",
    "    \n",
    "    # return dash and map\n",
    "    return change_dash\n",
    "\n",
    "# reset likelihood params (if run again after setting)\n",
    "prms_change.reset_change_params()\n",
    "\n",
    "# get sat fetch dash and map\n",
    "change_dash = build_change_dash(params=prms_change)\n",
    "display(change_dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 03B: Perform change analysis\n",
    "Perform spectral change vector analysis on thresholded GDV likelihood areas and plot the results onto an interactive map using the parameters selected in step 03B above. \n",
    "\n",
    "Simply run the cell below. Note: step 03A and 03B (this step) are intended to be interative. If the change vector result produced below does not match expectations, feel free to alter the parameters above (especially the angle range, magnitude threshold, time slice and baseline date range) and re-run this code cell to modify results.\n",
    "\n",
    "The basic steps are undertaken during this process are outlined in the Section 03 notes above.\n",
    "\n",
    "Run the code cell below to perform spectral change vector analysis. This can take quite awhile to complete, so please wait. Proceed to step 03C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# prepare for change vector analysis!\n",
    "def prepare_for_change_analysis(scene_like_threshed, params):\n",
    "    \n",
    "    # basic checks\n",
    "    if not params.check_general_params():\n",
    "        print('Not all required parameters set above. Please choose study area and satellite parameters.')\n",
    "        return None\n",
    "        \n",
    "    # enable xr apply progress bar (not implemented in xr... yet)\n",
    "    helpers.enable_xr_progressbar()\n",
    "    \n",
    "    # get earliest an latest dates\n",
    "    earliest_date, latest_date = helpers.get_earliest_and_latest_date_range(params.base_dates, params.comp_dates)\n",
    "    \n",
    "    if not earliest_date and not latest_date:\n",
    "        print('Problem with change base and comparison dates. Please check. Aborting.')\n",
    "        return None\n",
    "    \n",
    "    # create query date tuple\n",
    "    query_dates = (earliest_date, latest_date)\n",
    "    \n",
    "    # fetch data based on user parameters\n",
    "    ds = data.fetch_cube_data(platform=params.platform, product=params.product, q_dates=query_dates, cloud_cover=params.cloud_cover, \n",
    "                              cloud_mask=params.cloud_mask, coords_nw=params.alb_nw, coords_se=params.alb_se)  \n",
    "    \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not indices.check_required_bands_exist(ds=ds):\n",
    "        print('Not all required bands available in dataset. Please re-run the tool again.')\n",
    "        return ds\n",
    "    elif not params.check_change_idx_params():\n",
    "        print('No vegetation index set. Please ensure they are set above above.')\n",
    "        return ds\n",
    "    \n",
    "    # add an empty first date array if doesn't exist. this is important for resampling\n",
    "    ds = data.add_empty_first_date_array(ds=ds, time_slice=params.time_slice, first_query_date=query_dates[0])\n",
    "        \n",
    "    # reduce dataset down to specific quarter or month\n",
    "    print('\\nRestricting dataset down to chosen quarter/month, hold tight.')\n",
    "    ds = data.reduce_to_specific_month_or_season(ds=ds, time_slice=params.time_slice)\n",
    "    \n",
    "    # create vegetation indices and resample to specific season\n",
    "    ds = indices.create_resample_idxs_change_time_slice(ds=ds, params=params)\n",
    "            \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) == 0:\n",
    "        print('No years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # calculate some summary stats and graph\n",
    "    stats.calc_summary_stats_change(ds=ds)\n",
    "    fig_raw_idx_change_stats = plotting.plot_raw_change_stats(ds=ds, fignum=30)\n",
    "    plt.show()\n",
    "    \n",
    "    # do z-score if requested\n",
    "    if params.zscore:\n",
    "        ds = stats.zscore_flag_show_reset_change(ds=ds, z_val=params.zscore_value, reset=True)\n",
    "    else:\n",
    "        print('Skipping z-score outlier test for seasonal mann-kendall. Seasonal mann-kendall will consider outliers.')\n",
    "        \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_change_interpolate_idx_params():\n",
    "        print('No cleaning/interpolation parameters set. Please ensure they are set above.')\n",
    "        return ds\n",
    "    elif len(ds['time.year']) == 0:\n",
    "        print('No years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # interpolate and extrapolate, or drop missing data\n",
    "    if params.miss_data_handle == 'interpolate':\n",
    "        ds = interps.interp_change_vector(ds=ds)\n",
    "    elif params.miss_data_handle == 'drop':\n",
    "        print('Warning, you have chosen to drop missing data. This may result in significant data loss. Deleting data.')\n",
    "        ds = ds.dropna('time', how='all')\n",
    "        \n",
    "    # if user requested, fill missing pixels\n",
    "    if params.fill_pixels:\n",
    "        ds = interps.fill_pixels_change(ds=ds)\n",
    "    \n",
    "    # calculate some summary stats and show graphs\n",
    "    stats.calc_summary_stats_change(ds=ds)   \n",
    "    fig_raw_idx_change_clean_stats = plotting.plot_raw_change_stats(ds=ds, fignum=31)\n",
    "    plt.show()\n",
    "    \n",
    "    # basic check\n",
    "    if not ds:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_change_standardisation_params():\n",
    "        print('No standardisation parameters set. Please ensure they are set above.')\n",
    "        return ds\n",
    "    elif len(ds['time.year']) == 0:\n",
    "        print('No years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # standardise using basic greennest/ brightest sites and / max approach\n",
    "    ds = stats.do_standardisation_indices_change(ds=ds, params=params)\n",
    "    \n",
    "    # get median all time for veg, brt and plot\n",
    "    med_alltime = ds.median('time')\n",
    "    \n",
    "    # tell user and plot baselines\n",
    "    print('\\nPlotting greeness and brightness baseline images. Hold tight.')\n",
    "    fig_change_baselines = plotting.plot_median_alltime_idx_change(ds=med_alltime, stretch_min=0.001, stretch_max=0.999, fignum=32)\n",
    "    plt.show()\n",
    "    \n",
    "    return ds\n",
    "    \n",
    "# perform dataset split (baseline to comparison) and do change vector analysis\n",
    "def do_cva(ds, scene_like_threshed, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds or not scene_like_threshed:\n",
    "        print('Error occured after compute, no dataset exists. Aborting.')\n",
    "        return None\n",
    "    elif not params.check_cva_params():\n",
    "        print('Not all required parameters set above. Please choose change vector analysis parameters.')\n",
    "        return None\n",
    "    elif len(ds['time.year']) == 0:\n",
    "        print('Less than five years of data in dataset. Please include more years in your query above.')\n",
    "        return ds\n",
    "    \n",
    "    # split data into baseline and comparison datasets\n",
    "    print('Splitting dataset into baseline and comparison sets. Please wait.')\n",
    "    base_ds, comp_ds = stats.seperate_baseline_and_compare_ds(ds=ds, base_dates=params.base_dates, comp_dates=params.comp_dates)\n",
    "    \n",
    "    # prepare comparison dataset for holding cva results\n",
    "    print('> Preparing comparison dataset to hold cva variables. Please wait.')\n",
    "    comp_ds = stats.prepare_cva_comparison_dataset(comp_ds=comp_ds)\n",
    "        \n",
    "    # do change vector as apply function\n",
    "    print('\\nPerforming change vector analysis! This should not take too long, please wait.')\n",
    "    cva_ds = comp_ds.groupby('time').apply(stats.apply_cva, base_da=base_ds, tmf=params.change_threshold)\n",
    "    \n",
    "    # drop unnecessary vars left over\n",
    "    print('> Dropping unnecessary variables. Please wait.')\n",
    "    cva_ds = cva_ds.drop(['magnitude', 'brt_diff', 'veg_diff', 'sel'], errors='ignore')\n",
    "    \n",
    "    # restrict dcva angles to user requested range\n",
    "    print('\\nReducing change vector analysis result to user specified angles, cleaning up magnitudes. Please wait.')\n",
    "    cva_ds = stats.reduce_cva_to_specific_angles(cva_ds=cva_ds, min_ang=params.change_angles_range[0], max_ang=params.change_angles_range[1])\n",
    "    \n",
    "    # convert magnitude units into percentages\n",
    "    print('> Converting magnitude units into percentages. Please wait.')\n",
    "    cva_ds['magout'] = cva_ds['magout'] * 100    \n",
    "    \n",
    "    # mask out non-gdv areas\n",
    "    print('\\nReducing change vector analysis angles and magnitudes to groundwater dependent vegetation likelihood areas. Please wait.')\n",
    "    cva_ds = stats.mask_cva_to_gdv_likelihood(cva_ds=cva_ds, likelihood=scene_like_threshed['likelihood'])\n",
    "    \n",
    "    # tell user and return\n",
    "    print('Change vector analysis successfully completed.')\n",
    "    return cva_ds\n",
    "\n",
    "\n",
    "# check if trend map exists, disable otherwise\n",
    "if 'trend_map' in locals():\n",
    "    trend_map.close()\n",
    "\n",
    "# check if prior work done to create likelihood data, if not set none\n",
    "if 'scene_like_threshed' not in locals() or not scene_like_threshed:\n",
    "    scene_like_threshed, prms_change = data.load_and_prepare_like_netcdf(params=prms_change)\n",
    "    \n",
    "# fetch data and do change vector analysis\n",
    "if scene_like_threshed and prms_like and prms_change:\n",
    "    plat = helpers.identify_platform(da=scene_like_threshed)\n",
    "    \n",
    "    if prms_change.platform == plat:\n",
    "        prms_change.correct_sat_launch_date()\n",
    "        ds = prepare_for_change_analysis(scene_like_threshed=scene_like_threshed, params=prms_change)\n",
    "        cva_ds = do_cva(ds=ds, scene_like_threshed=scene_like_threshed, params=prms_change)\n",
    "    else:\n",
    "        print('Cannot perform change analysis using platform: {0} when likelihood used: {1}. Platforms must match.'.format(prms_change.platform, plat))\n",
    "        print('It is recommended running Section 01 again using the matching platform or changing your above selection.')   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 03C: Visualise change analysis results, GDV health graphing tools and change point detection controls\n",
    "The following step is used to plot change vector analysis result layers onto an interactive map to facilitate change exploration. Users can 'scrub' across dates using the slider in the top right corner of the interactive map. Users can also click on any pixel within the study area boundary to graph vegetation health signatures and perform change point detection to determine dates where vegetation health significantly shifted from the norm. This step is intended to be iterative with steps 03A and 03B above.\n",
    "\n",
    "The results of spectral change vector analysis provide two layers that users can explore on the interactive map:\n",
    "* `angles layer`: angle of change between baseline and comparison date. Different angles represent different changes (see notes at start of Section 02).\n",
    "The default value is an angle range of 90-180 degrees, which represents vegetation decline. The colour of these pixels will vary due to the large range of potential values (i.e. 0-360).\n",
    "* `magnitudes layer`: magnitude of distance between baseline and comparison date. The higher the magnitude, the larger the distance (or difference) was between dates. Darker red pixels are higher magnitude amounts, orange colours are lower magnitude values. Magnitude ranges from 0-100.\n",
    "\n",
    "The angle and magnitude layers act as a guide to direct users to potential GDV decline areas. More information can be obtained about these layers by simply clicking on the pixel of interest and running a change point detection analysis on it. This is achieved using the control panel underneat the interactive map. The control panel parameters include:\n",
    "* `model`: Cost function model. Users can select between Least absolute deviation (L1), Least squared deviation (L2), or Kernalised mean change (rbf). Recommended: rbf.\n",
    "These models dictate how the function detects changes within the vegetation health signal. L1 detects changes in the median of the signal. Overall it is a robust restimator of a shift in the central point of vegetation samples. L2 is the same, but detects mean-shifts in the vegetation health signal. Rbf detects changes in the mean of the embedded signal, based on Arlot et al. (2012). The Rbf model is recommended.\n",
    "* `min break distance`: Minimum distance between break points. One unit is equal to one date. Default is 2. Recommended: 1-3.\n",
    "Smaller value will allow break points to be detcted closer together, larger value will ensure gaps exists between breaks. Helpful for differentiating between short or long-term breaks.\n",
    "* `jump value`: Sets how large a 'jump' between vegetation health regimes need to be before they are detected. Recommended: 1.\n",
    "Smaller value will allow smaller shifts in vegetation health regimes to be detected. A larger value will only return very significant shifts in health regimes.\n",
    "* `num breaks (dynpro)`: Forces the Dynamic Programming change point detection method to return a set number of breaks points. Recommended: 2-4.\n",
    "Given a time-series of vegetation health values, tell the Dynamic Programming method how many break points to return. Even numbers only (start and end of breaks = 1 regime).\n",
    "* `break penelty (pelt)`: Set the penalty value for the PELT change point detection function. Recommend: 2-4.\n",
    "The penalty value is used in the PELT function to guard against model overfitting. Can be used to detect optimal breaks through minimisation of costs. Smaller values may return more breaks, but may suffer overfitting.\n",
    "* `show button`: Used to perform a change point detection with set parameters on the last clicked pixel. \n",
    "CLick this to generate a graph of vegetation health and perform a PELT and DynPro change point detection test for current pixel.\n",
    "\n",
    "Tip:\n",
    "* For the angles layer, colour map is a rainbow. You must click on the pixel to determine the precise angle this colour represents;\n",
    "* For the magnitude layer, orange pixels indicate lower magnitude of change. The more red a pixel is, the more magnitude of change that was detected at that pixel;\n",
    "* Click on any pixel to get a reading of the angle or magnitude at that particular location;\n",
    "* Click on any pixel and click the 'Show' button in the Change Point Detection control panel to perform a change point detection and see vegetation graphs;\n",
    "* The change point detection default values are recommended; and\n",
    "* The change point detection processors are not limited to GDV likelihood and change vector pixels - click anywhere in your study area and click 'Show' to see potential breaks in vegetation health.\n",
    "\n",
    "Run the code cell below to generate the interactive map and change point detection control panel. This shouldn't take long. This is the final step involved with change vector analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# prepare interactive cva map \n",
    "def build_cva_map_and_cpd_panel(ds, cva_ds, params):\n",
    "    \n",
    "    # basic check\n",
    "    if not ds or not cva_ds:\n",
    "        print('No dataset exists. Aborting.')\n",
    "        return None\n",
    "    \n",
    "    # globals\n",
    "    change_x, change_y = None, None\n",
    "    change_msg = None\n",
    "    change_info = HTML(value='<h3></h3>')\n",
    "\n",
    "    # reset cpd params\n",
    "    params.reset_cpd_params()\n",
    "    \n",
    "    # capture click\n",
    "    def handle_click(**kwargs):\n",
    "        global change_x, change_y\n",
    "        global change_msg\n",
    "\n",
    "        if kwargs['type'] == 'click':\n",
    "            try:\n",
    "                lon, lat = kwargs['coordinates'][1], kwargs['coordinates'][0]\n",
    "                change_x, change_y = helpers.transform_point(lon=lon, lat=lat, in_epsg='epsg:4326', out_epsg='epsg:3577')\n",
    "                \n",
    "                date = [control.widget.value for control in change_map.controls if 'SelectionSlider' in str(control)]\n",
    "                date = date[0]\n",
    "                \n",
    "                change_msg = stats.get_cva_pixel_info(cva_da=cva_ds.sel(time=date), x=change_x, y=change_y)\n",
    "                change_info.value = '<h3>{0}</h3>'.format(change_msg)\n",
    "                \n",
    "            except:\n",
    "                print('Problem during build_cva_map_and_cpd_panel click event. Skipping.')\n",
    "            \n",
    "    # on button click\n",
    "    def handle_button(event):\n",
    "        global change_x, change_y\n",
    "\n",
    "        with output:\n",
    "            clear_output()\n",
    "            \n",
    "            # do change point detection and store results for veg\n",
    "            vec_veg, dates_veg, result_pelt_veg, result_dynp_veg = stats.do_cva_change_pd(ds=ds['veg_idx'].to_dataset(), index='veg_idx', x=change_x, y=change_y, \n",
    "                                                                                          model=params.model, min_size=params.min_size, jump=params.jump, \n",
    "                                                                                          dynp_n_bkps=params.dynp_n_bkps, pelt_penelty=params.pelt_penelty)\n",
    "            \n",
    "            # do change point detection and store results for brt\n",
    "            vec_brt, dates_brt, result_pelt_brt, result_dynp_brt = stats.do_cva_change_pd(ds=ds['brt_idx'].to_dataset(), index='brt_idx', x=change_x, y=change_y, \n",
    "                                                                                          model=params.model, min_size=params.min_size, jump=params.jump, \n",
    "                                                                                          dynp_n_bkps=params.dynp_n_bkps, pelt_penelty=params.pelt_penelty)           \n",
    "            \n",
    "            # check if anything empty\n",
    "            if (str(vec_veg) == 'None' or str(dates_veg) == 'None') or (len(vec_veg) == 0 or len(dates_veg) == 0):\n",
    "                print('No veg vector returned. Did you click outside the study area? Try again.')\n",
    "                return\n",
    "            if (str(vec_brt) == 'None' or str(dates_brt) == 'None') or (len(vec_brt) == 0 or len(dates_brt) == 0):\n",
    "                print('No brt vector returned. Did you click outside the study area? Try again.')\n",
    "                return\n",
    "            \n",
    "            # plot pelt cbd\n",
    "            if result_pelt_veg and result_pelt_brt:\n",
    "                fig_cva_pelt, fig_cva_axs_pelt = plotting.plot_cva_cpd(signal_veg=vec_veg, signal_bright=vec_brt, dates=dates_veg, \n",
    "                                                                       result_veg=result_pelt_veg, result_brt=result_pelt_brt,\n",
    "                                                                       cbd_type='PELT', fignum=41)\n",
    "                display(fig_cva_pelt)\n",
    "            else:\n",
    "                display('Change point (PELT) did not detect a break. Skipping.')\n",
    "\n",
    "            # plot dynp cbd\n",
    "            if result_dynp_veg and result_dynp_brt:\n",
    "                fig_cva_dynp, fig_cva_axs_dynp =  plotting.plot_cva_cpd(signal_veg=vec_veg, signal_bright=vec_brt, dates=dates_veg, \n",
    "                                                                        result_veg=result_dynp_veg, result_brt=result_dynp_brt,\n",
    "                                                                        cbd_type='DynPro', fignum=42)\n",
    "                display(fig_cva_dynp)\n",
    "            else:\n",
    "                display('Change point (DynPro) did not detect a break. Skipping.')\n",
    "\n",
    "    # load interactive map and add interactivity\n",
    "    change_map = maps.load_change_map(ds=ds, cva_ds=cva_ds, params=params)\n",
    "    change_map.on_interaction(handle_click)\n",
    "\n",
    "    # load cpd dashboard\n",
    "    change_cpd_dash = controls.build_change_cpd_dashboard(params=params, handle_button=handle_button)\n",
    "                    \n",
    "    # load button and output and add interactivity\n",
    "    output = Output()\n",
    "\n",
    "    # return all controls\n",
    "    return change_map, change_info, change_cpd_dash, output\n",
    "\n",
    "# get interactive map and dash\n",
    "change_map, change_info, change_cpd_dash, change_output = build_cva_map_and_cpd_panel(ds=ds, cva_ds=cva_ds, params=prms_change)\n",
    "\n",
    "# show controls\n",
    "display(change_map, change_info, change_cpd_dash, change_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 04: References\n",
    "\n",
    "* Robertson, M., Villet, M. and Palmer, A. (2004). A fuzzy classification technique for predicting species’ distributions: applications using invasive alien plants and indigenous insects. Diversity and Distributions, 10(5-6), pp.461-474.\n",
    "\n",
    "* Robinson, T., van Klinken, R. and Metternicht, G. (2010). Comparison of alternative strategies for invasive species distribution modeling. Ecological Modelling, 221(19), pp.2261-2269.\n",
    "\n",
    "* Saaty, T.L., 1980. “The Analytic Hierarchy Process.” McGraw-Hill, New York.\n",
    "\n",
    "* Snedecor, G.W. (1956). Statistical methods: Applied to experiments in agriculture and biology. The Iowa State College Press.\n",
    "\n",
    "* Mann, H.B. (1945), Nonparametric tests against trend, Econometrica, 13, 245-259.\n",
    "\n",
    "* Hipel, K.W. and McLeod, A.I., (2005). Time Series Modelling of Water Resources and Environmental Systems.\n",
    "\n",
    "* Hirsch, R. M., Slack, J. R., & Smith, R. A. (1982). Techniques of trend analysis for monthly water quality data. Water resources research, 18(1), 107:121.\n",
    "\n",
    "* C. Truong, L. Oudre, N. Vayatis. (2020) Selective review of offline change point detection methods. Signal Processing, 167:107299, 2020.\n",
    "\n",
    "* Hamed, K. (2008). Trend detection in hydrologic data: The Mann–Kendall trend test under the scaling hypothesis. Journal of Hydrology, 349(3-4), pp.350-363.\n",
    "\n",
    "* Wambui, G.D., Waititu, G.A. and Wanjoya, A. (2015). The power of the pruned exact linear time (PELT) test in multiple changepoint detection. American Journal of Theoretical and Applied Statistics, 4(6), pp.581-586.\n",
    "\n",
    "* Truong, C., Oudre, L. and Vayatis, N., 2019. Selective review of offline change point detection methods. Signal Processing, p.107299.\n",
    "\n",
    "* Hussain et al., (2019). pyMannKendall: a python package for non parametric Mann Kendall family of trend tests.. Journal of Open Source Software, 4(39), 1556,\n",
    "\n",
    "* Arlot, S., Celisse, A., and Z. Harchaoui. Kernel change-point detection. arXiv preprint arXiv:1202.3878, 1(0000):1–26, 2012.\n",
    "\n",
    "* Malila, W.A., 1980. Change vector analysis: an approach for detecting forest changes with Landsat. In LARS symposia (p. 385).\n",
    "\n",
    "* Zanchetta, A. and Bitelli, G., 2017. A combined change detection procedure to study desertification using opensource tools. Open Geospatial Data, Software and Standards, 2(1), p.10.\n",
    "\n",
    "* Crist, E.P., 1985. A TM tasseled cap equivalent transformation for reflectance factor data. Remote Sensing of Environment, 17(3), pp.301-306.\n",
    "\n",
    "* O’Neill, AL, 1996, Satellite-derived vegetation indices applied to semi-arid shrublands in Australia, Australian Geographer, 50, pp. 185-199.\n",
    "\n",
    "* Andrew, M. and Warrener, H. (2017). Detecting microrefugia in semi-arid landscapes from remotely sensed vegetation dynamics. Remote Sensing of Environment, 200, pp.114-124.\n",
    "\n",
    "* Tucker, CJ 1979, ‘Red and photographic infrared linear combinations for monitoring vegetation’, Remote Sensing of Environment, vol. 8, pp. 127-150.\n",
    "\n",
    "* O’Grady, A. P., Carter, J.L. and Bruce, J. (2011) Can we predict groundwater discharge from terrestrial ecosystems using existing eco-hydrological concepts? Hydrology and Earth System Science Discussions, 8: 8231-8253.\n",
    "\n",
    "* Barron, O. V., Emelyanova, I., Van Niel, T.G., Pollock, D., and Hodgson, G. (2014) Mapping groundwater-dependent ecosystems using remote sensing measures of vegetation and moisture dynamics, Hydrological Processes, 28: 272-385.\n",
    "\n",
    "* Andersen, T., Carstensen, J., Hernandez-Garcia, E. and Duarte, C.M., 2009. Ecological thresholds and regime shifts: approaches to identification. Trends in Ecology & Evolution, 24(1), pp.49-57.\n",
    "\n",
    "* Mann, H.B. (1945), Nonparametric tests against trend, Econometrica, 13, 245-259.\n",
    "\n",
    "* Hipel, K.W. and McLeod, A.I., (2005). Time Series Modelling of Water Resources and Environmental Systems.\n",
    "\n",
    "* Hamed, K. (2008). Trend detection in hydrologic data: The Mann–Kendall trend test under the scaling hypothesis. Journal of Hydrology, 349(3-4), pp.350-363.\n",
    "\n",
    "* Keeley, J. E. (2009). Fire intensity, fire severity and burn severity: A brief review and suggested usage. International Journal of Wildland Fire, 18(1), 116–126."
   ]
  }
 ],
 "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.6.9"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}