{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "20925cd3",
   "metadata": {},
   "source": [
    "# Visualizing and Analyzing OPERA DIST-ALERT and DIST-ANN Products to Visualize to Wildfire Impact in Northern California\n",
    "This notebook is designed to showcase the use of the OPERA DIST-ALERT and DIST-ANN products to visualize vegetation disturbance associated with a series of wildfires that affected northern California and southern Oregon during the 2022 calendar year. OPERA DIST-ALERT data enables timelapse examination of the extent and severity of vegetation loss, whereas the DIST-ANN data provide a yearly summary of confirmed vegetation loss in a single composite tile. A Jupyter Notebook specifically focusing on the DIST-ALERT product to investigate a wildfire is available [here](https://github.com/OPERA-Cal-Val/OPERA_Applications/blob/main/DIST/Wildfire/McKinney.ipynb).\n",
    "\n",
    "A [NASA Earthdata Login](https://urs.earthdata.nasa.gov/) account is required to download the data used in this tutorial. You can create an account at the link provided.\n",
    "\n",
    "*<font color='red'>Note 1: This notebook uses provisional products, which may differ slightly from operational products. Please refer to [DIST product specification](https://d2pn8kiwq2w21t.cloudfront.net/documents/ProductSpec_DIST_HLS.pdf) for more information. </font>*<br><br>\n",
    "*<font color='red'>Note2: DIST products are distributed via NASA's Distributed Active Archive Centers (DAACs), specifically the [LP DAAC](https://lpdaac.usgs.gov). However, the DIST-ANN data accessed and visualized in this particular notebook were produced specifically for instructional purposes and specific for this notebook, and are instead, stored and distributed by the Global Land Analysis and Discovery (GLAD) group at the University of Maryland.</font>*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acbd93f1",
   "metadata": {},
   "source": [
    "## Background\n",
    "A series of fires affected counties of northern California during the 2021 and 2022 wildfire seasons. We will be looking at DIST-ALERT and DIST-ANN tiles that captures the  vegetation disturbance associated with a subset of these fires located in the region of Klamath National Forest in western Siskiyou County:\n",
    "\n",
    "| Name | Alarm Date | Containment Date | Centroid (lat)| Centroid (lon) | Official Area Impacted (acres) |\n",
    "| ---  | --- | --- | --- | --- | --- |\n",
    "| Gulch | 3/12/2022 | 3/14/2022 | 41.882 | -121.849 | 113 |\n",
    "| Mountain | 9/2/2022 | 10/2/2022 | 41.462 | -122.66 | 13440 |\n",
    "| McKinney | 7/29/2022 | 11/11/2022 | 41.814 | -122.841 | 60102 |\n",
    "| Alex | 7/31/2022 | 8/3/2022 | 41.939 | -122.952 | 150 |\n",
    "| Mill | 9/2/2022 | 9/11/2022 | 41.476 | -122.393 | 3939 |\n",
    "| Coyote | 9/7/2022 | 9/8/2022 | 41.834 | -121.793 | 296 | \n",
    "\n",
    "*Data from the [California Data Portal](https://data.ca.gov/dataset/california-fire-perimeters-all)*\n",
    "\n",
    "We will be visualizing the spatial and temporal impacts of the fires on vegetation changes caused by these fires using the **OPERA DIST-ALERT and DIST-ANN products**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09046434",
   "metadata": {},
   "source": [
    "### DIST Product Suite Background\n",
    "---\n",
    "The land Disturbance product suite (**DIST**) maps vegetation disturbance from Harmonized Landsat-8 and Sentinel-2 A/B (HLS) scenes. Disturbance is detected when vegetation cover decreases or spectral variation is outside a historical norm within an HLS pixel. Two DIST products compose the DIST product suite: 1) the **DIST-ALERT** product, capturing vegetation disturbance at the cadence of HLS sampling (2-3 days); and 2) the **DIST-ANN** product, summarizing the confirmed changes of the DIST-ALERT products from previous calendar year.  *<font color='red'>Note: The DIST-ANN product used here is a provisional product with an observation date exceeding the operational 1-year period. </font>*\n",
    "\n",
    "This notebook provides a step-by-step workflow visualizing **DIST-ANN** raster layers for the 2022 calendar year. An analogous notebook for the **DIST-ALERT** product may be accessed in the [OPERA Applications Github repository](https://github.com/OPERA-Cal-Val/OPERA_Applications/blob/main/DIST/Wildfire/Intro_To_DIST.ipynb).\n",
    "\n",
    "### Metadata\n",
    "---\n",
    "\n",
    "HLS products provide surface reflectance (SR) data from the Operational Land Imager (OLI) aboard the Landsat-8 remote sensing satellite and the Multi-Spectral Instrument (MSI) aboard the Sentinel-2 A/B remote sensing satellite. HLS products are distributed over projected map coordinates aligned with the Military Grid Reference System (MGRS). Each tile covers 109.8 square kilometers divided into 3660 rows and 3660 columns at 30 meter pixel spacing. Each tile overlaps neighbors by 4900 meters in each direction. The **DIST-ANN** product is stored distributed as a set of 16 Cloud-Optimized GeoTIFF (COG) files.  Details specific to the available raster layers and their properties are available in the [OPERA DIST Product Specifications Document](https://d2pn8kiwq2w21t.cloudfront.net/documents/ProductSpec_DIST_HLS.pdf).\n",
    "\n",
    "### Data Distribution\n",
    "DIST product data are stored and distributed via NASA's Distributed Active Archive Centers (DAACs), specifically the [LP DAAC](https://lpdaac.usgs.gov).  However, the DIST-ANN data accessed and visualized in this particular notebook were produced specifically for instructional purposes and specific for this notebook, and are instead, stored and distributed by the Global Land Analysis and Discovery (GLAD) group at the University of Maryland.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a65b1d51",
   "metadata": {},
   "source": [
    "## Library Imports and Configurations\n",
    "First we import the necessary Python libraries. These come pre-installed in the `opera_app` anaconda environement within the [OPERA Applications Gitub repository](https://github.com/OPERA-Cal-Val/OPERA_Applications). We also import a collection of custom DIST-specific functions from a source file called `dist_utils.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ecde9f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notebook dependencies\n",
    "import hvplot.xarray\n",
    "import geoviews as gv\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "from datetime import datetime\n",
    "import ipyleaflet\n",
    "import leafmap.leafmap as leafmap\n",
    "from bokeh.models import FixedTicker\n",
    "\n",
    "import holoviews as hv\n",
    "hv.extension('bokeh')\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import sys\n",
    "sys.path.append('../../')\n",
    "from src.dist_utils import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6040454f",
   "metadata": {},
   "source": [
    "The next cell configures gdal and provideds necessary authentication to successfully access LP DAAC cloud-hosted assets.\n",
    "We first determine that valid Earthdata credentials are present in a .netrc file, used for accessing LP DAAC Cloud Assets. We also set the configuration options for GDAL to access the data successfully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec10efd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check for valid Earthdata credentials. If entering for the first time, be very careful to enter username/password correctly!\n",
    "check_netrc()\n",
    "\n",
    "# Set GDAL configs to successfully access LP DAAC Cloud Assets via vsicurl\n",
    "gdal.SetConfigOption(\"GDAL_HTTP_UNSAFESSL\", \"YES\")\n",
    "gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')\n",
    "gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')\n",
    "gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','FALSE')\n",
    "gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be589981",
   "metadata": {},
   "source": [
    "## Create Products from Harmonized Landsat and Sentinel (HLS) Data\n",
    "The OPERA Disturbance product is derived from HLS. The time series of HLS is what provides the quantification of vegetation disturbance. \n",
    "\n",
    "In the next cell, we specify and access HLS data from Earthdata to provide some background and insight into how the DIST products are derived. We have chosen an HLS tile that covers the fire-affected region. We retreive this tile for two different time periods, namely July 21, 2022 and September 23, 2022. The tiles have the exact same spatial extent. The next cell will use custom DIST functions contained in `dist_utils.py` to produce three HLS products: (1) True color; (2) False color; and (3) Normalized Difference Vegetation Index (NDVI). These files will be available in a subdirectory on the user's file system called `tifs`. The false color and NDVI products help to visualize changes in vegetation that are further quantified with the HLS DIST-ALERT product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e6ba4c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to data location on LP DAAC and provide list of desired bands\n",
    "hls1 = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T10TEM.2022202T185058.v2.0/HLS.L30.T10TEM.2022202T185058.v2.0.'\n",
    "hls2 = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T10TEM.2022266T185119.v2.0/HLS.L30.T10TEM.2022266T185119.v2.0.'\n",
    "veg_dist_status = 'https://glad.umd.edu/projects/opera/SEP/DIST-ANN/10/T/E/M/2022/OPERA_L3_DIST-ANN-HLS_T10TEM_2022_2023136T210508Z_30_v0_VEG-DIST-STATUS.tif'\n",
    "bandlist = ['B05','B04','B03','B02']\n",
    "\n",
    "# Make true color image\n",
    "make_hls_true_color(hls1, bandlist, 'true_color_7_21_22_NorCal.tif')\n",
    "make_hls_true_color(hls2, bandlist, 'true_color_9_23_22_NorCal.tif')\n",
    "\n",
    "# Make false color image\n",
    "make_hls_false_color(hls1, bandlist, 'false_color_7_21_22_NorCal.tif')\n",
    "make_hls_false_color(hls2, bandlist, 'false_color_9_23_22_NorCal.tif')\n",
    "\n",
    "# Make ndvi image\n",
    "make_hls_ndvi(hls1, bandlist, 'ndvi_7_21_22_NorCal.tif')\n",
    "make_hls_ndvi(hls2, bandlist, 'ndvi_9_23_22_NorCal.tif')\n",
    "\n",
    "# Make VEG-DIST-STATUS image\n",
    "make_veg_dist_status_visual(veg_dist_status, 'veg_dist_status_NorCal.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c03375c",
   "metadata": {},
   "source": [
    "## Visualize Pre- and Post-Fire HLS Products\n",
    "We now visualize the true color, false color, and NDVI images for the area of interest. First let's visualize the pre- and post-fire true color HLS data. This visualization workflow uses the open-source `leafmap` Python library to create informative timelapse visualizations of the area of interst. `leafmap` provides numerous powerful capabilities for geospatial data visualization and analysis in the Jupyer environment. For additional details on the `leafmap` library, see the [leafmap docs](https://leafmap.org).\n",
    "\n",
    "Below we create a series of 'split maps' which show the pre-fire image on the left hand side and the post-fire image on the right hand side. The map that appears after executing the cell is interactive. Move the slider to the right to view the July 21, 2022 (pre-fire) image and to the left to view the September 23, 2022 (post-fire) image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc2c4254",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[42.0, -122.3], zoom=10)\n",
    "m.split_map(\n",
    "    left_layer='tifs/true_color_7_21_22_NorCal.tif',\n",
    "    right_layer='tifs/true_color_9_23_22_NorCal.tif',\n",
    "    left_label=\"7-21-2022\",\n",
    "    right_label=\"9-23-2022\"\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb841f13",
   "metadata": {},
   "source": [
    "Close examination reveals a difference in the visual appearance of the landscape in the left side of the post-fire true color image. This is the burn scar associated with the McKinney Fire. You may notice other areas that appear visually different between the two scenes. However, visually the fire-affected area is not very easy to delineate in the true color image. We can leverage both the false color and NDVI images to more clearly see the affected region. Let's first look at a false color image of the same area (explained further below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1ae3776",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[42.0, -122.3], zoom=10)\n",
    "m.split_map(\n",
    "    left_layer='tifs/false_color_7_21_22_NorCal.tif',\n",
    "    right_layer='tifs/false_color_9_23_22_NorCal.tif',\n",
    "    left_label='7-21-2022',\n",
    "    right_label=\"9-23-2022\"\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "915bddfd",
   "metadata": {},
   "source": [
    "Above is what is called a 'false color' image. In this false color image, the near infrared (NIR) band replaces the red band. The NIR band is sensitive to the chlorophyll in leafy vegetation, making the false color image a useful tool for investigating vegatation change between two images.\n",
    "Here, red color is a strong indicator of vegetation cover, whereas blue and green colors are more indicative of non-vegetative cover. Notice the regions which appear more red in the left (pre-fire) image and more blue-green in the second (post-fire) image. These regions have undergone vegetation loss between the two HLS scenes. Zoom in on the image to view the changes in more detail. Also note the changes in the closely clustered circular regions in the eastern and central regions of the tile  – these are agricultural changes.\n",
    "\n",
    "Yet another way to highlight vegetation pixels is through the use of the normalized difference vegetation index (NDVI). This is a band ratio between the NIR and red bands, which produces an image with values ranging from 0-1. Values nearer 0 are unlikely to be vegetation, whereas values near 1 are very likely to represent vegetation. Let's have a look at the NDVI images we have produced. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e18fe35",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[42.0, -122.3], zoom=10)\n",
    "m.split_map(\n",
    "    left_layer='tifs/ndvi_7_21_22_NorCal.tif',\n",
    "    right_layer='tifs/ndvi_9_23_22_NorCal.tif',\n",
    "    left_label=\"7-21-2022\",\n",
    "    right_label=\"9-23-2022\"\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f75d148a",
   "metadata": {},
   "source": [
    "In the above interactive map, darker colors represent NDVI values nearer to 0 (non-vegetation) while lighter colors represent NDVI values nearer to 1 (most likely vegetation). Pan between the pre- and post-fire NDVI tiles. Is the burn scar easier to delineate now?\n",
    "\n",
    "Let's now turn to look at the DIST-HLS product, and specifically the DIST-ANN yearly summary of vegetation change for the region of interest. The next cell pans between the NDVI image we have produced and the co-located DIST-HLS data. While the DIST-HLS product is not derived from the NDVI data directly, there is a clear relationship between the vegetation decrease indicated in the NDVI and the pixels classified as 'disturbance' in the DIST-ANN data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a650589",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[42.0, -122.3], zoom=10)\n",
    "m.split_map(\n",
    "    left_layer='tifs/ndvi_9_23_22_NorCal.tif',\n",
    "    right_layer='tifs/veg_dist_status_NorCal.tif',\n",
    "    left_label=\"9-23-2022 NDVI\",\n",
    "    right_label=\"2022 DIST-ANN VEG-DIST-STATUS\"\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40dfbc5e",
   "metadata": {},
   "source": [
    "In the sections below, we will explore the OPERA DIST-ALERT and DIST-ANN product suite to see how these qualitative changes visible in the HLS may be more further quantified."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecaf2f4d",
   "metadata": {},
   "source": [
    "## Load DIST-ALERT data\n",
    "First we will explore the OPERA DIST-ALERT product, which provides a measure of vegetation disturbance at the cadence of the HLS data (2-3 days). The next two cells path to and create a list containing a subset of DIST-ALERT tiles spanning the 2022-2023 calendar years. The data explored here are currently stored and distributed by the Global Land Analysis and Discovery (GLAD) group at the University of Maryland, and are accessed through a series of url links contained in the provided `10TEM_DIST-ALERT_links.txt`. This file is available in the `DIST` repository by default. \n",
    "\n",
    "We explore two layers of the DIST-ALERT product, namely the (1) `VEG-DIST-STATUS` and (2) `VEG-IND` layers. `VEG-DIST-STATUS` tracks pixels exhibiting provisional and confirmed vegetation disturbance, whereas `VEG-IND` tracks the pixel-wise vegetation indicator value. These layers are helpful for depicting the spatial extent and severity of vegetation loss.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4e56e3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Open file containing paths to DIST-ALERT for Northern California ###\n",
    "\n",
    "# Specify location of data and open as a list\n",
    "dist_alert_tiles_local = '../links_to_data/10TEM_DIST-ALERT_links.txt'\n",
    "file = open(dist_alert_tiles_local, \"r\")\n",
    "lines = file.readlines()\n",
    "file.close()\n",
    "\n",
    "# Get only the VEG-DIST-STATUS and VEG-IND paths\n",
    "veg_dist_status = []\n",
    "veg_ind = []\n",
    "for line in lines:\n",
    "    fp = line[:-1]\n",
    "    # Check if the file path ends with 'VEG-DIST-STATUS.tif' or 'VEG-IND' and add to corresponding list\n",
    "    if fp.endswith('VEG-DIST-STATUS.tif'):\n",
    "        veg_dist_status.append(line.strip())\n",
    "    elif fp.endswith('VEG-IND.tif'):\n",
    "        veg_ind.append(line.strip())\n",
    "    \n",
    "# There are ~250 VEG-DIST-STATUS and VEG-HIST tiles, lets keep every 10th tile. This results in a subset of ~25 tiles for exploration\n",
    "veg_dist_status = veg_dist_status[0::10]\n",
    "veg_ind = veg_ind[0::10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0259636c",
   "metadata": {},
   "source": [
    "## Visualize spatial extent and severity of vegetation change through time with the OPERA DIST-ALERT product\n",
    "Let's first examine the spatial extent of the fire. This notebook uses the open-source `leafmap` and `ipyleaflet` libraries for visualization and for enabling custom user-defined areas of interest on interactive maps. For more information on these libraries, see the [leafmap](https://leafmap.org) and [ipyleaflet](https://ipyleaflet.readthedocs.io/en/latest/) docs.\n",
    "\n",
    "Throughout this notebook, the user will encounter several interactive maps displaying a raster layer a series of tools on the left-hand side for zooming and drawing bounding boxes with `ipyleaflet`. The next cell is the first example of this functionality. The interactive map will include a raster file and several point locations with information related to the wildfires. This data is pulled from the `ca_fires_2022.csv` contained in the `/supplemental_data/` directory.\n",
    "\n",
    "The McKinney Fire is the largest white pseudo-circular area in the western portion of the raster. Use the `rectangle` tool to draw a bounding box over this region to retreive information from the underlying raster data within the bounding box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af04e900",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot VEG-DIST-STATUS on ipyleaflet map ###\n",
    "dist_url ='https://glad.umd.edu/projects/opera/SEP/DIST-ANN/10/T/E/M/2022/OPERA_L3_DIST-ANN-HLS_T10TEM_2022_2023136T210508Z_30_v0_VEG-DIST-STATUS.tif'\n",
    "\n",
    "# Make leaflet map\n",
    "m = leafmap.Map(basemap=ipyleaflet.basemaps.Esri.WorldImagery,\n",
    "    center=(42.0, -122.3),\n",
    "    zoom=9,\n",
    "    crs=ipyleaflet.projections.EPSG3857,\n",
    "    draw_control=False)\n",
    "\n",
    "# Add raster and point data and draw functionality\n",
    "m.add_cog_layer(dist_url, name=\"VEG-DIST-STATUS\")\n",
    "m.add_points_from_xy('../supplemental_data/ca_fires_2022.csv', x ='Centroid (Lon)', y = 'Centroid (Lat)', layer_name=\"CaliforniaFires\")\n",
    "dc = ipyleaflet.DrawControl(\n",
    "    polygon={},\n",
    "    rectangle={\"shapeOptions\": {\"color\": \"blue\"}},\n",
    "    circlemarker={},\n",
    "    circle = {},\n",
    "    polyline={}\n",
    ")\n",
    "\n",
    "# Draw an AOI on an interactive map\n",
    "print('Select an Area of Interest using the tools on the left side of the map.')\n",
    "dc.on_draw(handle_draw)\n",
    "m.add_control(dc)\n",
    "display(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1c53b18",
   "metadata": {},
   "source": [
    "Below, we use the built-in `leafmap` zonal statistics tool to compute the area for each class in the `VEG-DIST-STATUS` tile and the mean vegetation indicator value in the `VEG-IND` tile within this rectangular region. We then plot the area and mean vegetation indicator as a time-series to visualize how vegetation disturbance extent and severity has evolved through time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f87a797",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Compute the area of each class within the bounding box through time ###\n",
    "# User-Defined Parameters\n",
    "from shapely import box\n",
    "ll = dc.last_draw['geometry']['coordinates'][0][0]\n",
    "ur = dc.last_draw['geometry']['coordinates'][0][2]\n",
    "aoi = box(ll[0], ll[1], ur[0], ur[1])\n",
    "\n",
    "# Make GeoDataFrame from AOI\n",
    "gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[aoi])\n",
    "\n",
    "# Compute the affected areas through time\n",
    "affected_areas_through_time = []\n",
    "veg_ind_through_time = []\n",
    "dates = []\n",
    "for i,tile in enumerate(veg_dist_status[1:]):\n",
    "    \n",
    "    date = extract_date_from_string(tile)\n",
    "    date = datetime.strptime(date, '%Y%j')\n",
    "    area_stats = leafmap.zonal_stats(gdf, tile, stats=['count'], categorical=True)[0]\n",
    "    ind_stats = leafmap.zonal_stats(gdf, veg_ind[i], stats=['mean'], categorical=False)[0]\n",
    "    if ind_stats['mean'] is not None:\n",
    "        dates.append(date)\n",
    "        veg_ind_through_time.append(ind_stats['mean'])\n",
    "    del area_stats['count']\n",
    "    for status in range(5):\n",
    "        if status not in area_stats:\n",
    "            area_stats[status] = 0\n",
    "    area_stats = dict(sorted(area_stats.items()))\n",
    "    pixel_area = 30 * 30\n",
    "    affected_areas_through_time.append(compute_areas([area_stats], pixel_area, 'alert', date))\n",
    "\n",
    "# Make Pandas dataframes for the computed statistics\n",
    "combined_areas = pd.concat(affected_areas_through_time, axis=0)\n",
    "veg_indicators = pd.DataFrame({'Date':dates, 'Vegetation Indicator':veg_ind_through_time})\n",
    "veg_indicators.set_index('Date', inplace=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "420a3e03",
   "metadata": {},
   "source": [
    "Now we can plot the results graphically using plot functionality within the open-source `pandas` library, whose docs are available [here](https://pandas.pydata.org/docs/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e89d33d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Group the area data by \"Date\" and \"Class\"\n",
    "grouped_data = combined_areas.groupby(['Date', 'VEG-DIST-STATUS Class'])['Area (km2)'].sum()\n",
    "\n",
    "# Reset the index of the grouped data to turn it back into a DataFrame\n",
    "grouped_df = grouped_data.reset_index()\n",
    "\n",
    "# Pivot the data to have \"Class\" values as columns and \"Date\" as index\n",
    "pivot_df = grouped_df.pivot(index='Date', columns='VEG-DIST-STATUS Class', values='Area (km2)')\n",
    "\n",
    "# Plot the data\n",
    "fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(15,5))\n",
    "pivot_df.plot(ax=ax1,marker='.')\n",
    "veg_indicators.plot(ax=ax2,marker='.')\n",
    "\n",
    "# Add labels and title\n",
    "ax1.set_xlabel('Date')\n",
    "ax1.set_ylabel('Area (km2)')\n",
    "ax2.set_xlabel('Date')\n",
    "ax2.set_ylabel('Vegetation Indicator')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4be44ecd",
   "metadata": {},
   "source": [
    "If the McKinney Fire region was selected with the bounding box, the above plots should indicate a notable decrease in area corresponding to `VEG-DIST-STATUS` class 0 (No disturbance) and a corresponding increase in `VEG-DIST-STATUS` class 4 (Confirmed; ≥50% disturbance) that occurs around August 2022. This corresponds to the time to the occurance of the McKinney Fire in early August, when disturbance increased significantly. Likewise, the plot of vegetation indicator through time indicates a noteable decrease at the same time, inicating the severity of vegetation loss."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1535e10e",
   "metadata": {},
   "source": [
    "## Visualize and Analyze the DIST-ANN Data\n",
    "We will now explore the annual summary of disturbance associated with the McKinney Fire. Below we load in the layers of the `DIST-ANN` product to visualize a summary of the spatial extent of disturbance for the 2022 calendar year. The second cell below stacks the layers to create a 'geocube' of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c5645ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify filepath location of OPERA DIST_ANN tile and desired bands (Note: bandlist is not comprehensive of all available layers)\n",
    "data_dir = 'https://glad.umd.edu/projects/opera/SEP/DIST-ANN/10/T/E/M/2022/OPERA_L3_DIST-ANN-HLS_T10TEM_2022_2023136T210508Z_30_v0_'\n",
    "bandlist = ['VEG-DIST-STATUS', 'VEG-HIST', 'VEG-IND-MAX','VEG-ANOM-MAX', 'VEG-DIST-CONF', 'VEG-DIST-DATE', \n",
    "            'VEG-DIST-COUNT', 'VEG-DIST-DUR', 'VEG-LAST-DATE', 'GEN-DIST-STATUS', 'GEN-ANOM-MAX', 'GEN-DIST-CONF',\n",
    "            'GEN-DIST-DATE', 'GEN-DIST-COUNT', 'GEN-DIST-DUR', 'GEN-LAST-DATE']\n",
    "bandpath = f\"{data_dir}%s.tif\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "337f9168",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create geocube of stacked bands\n",
    "da_ann, crs = stack_bands(bandpath, bandlist)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9111026d",
   "metadata": {},
   "source": [
    "### Visualize the DIST-ANN VEG_DIST_STATUS Layer\n",
    " Any pixel which registers as confirmed disturbance in the DIST-ALERT `VEG-DIST-STATUS` data throughout the calendar year will be added to the yearly DIST-ANN `VEG-DIST-STATUS` layer for that year. The DIST-ANN `VEG-DIST-STATUS` layer tracks only confirmed changes, and whether the disturbance was greater or less than 50% when compared to historical vegetation cover. Below we visualize the DIST-ANN `VEG-DIST-STATUS` layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f59379a",
   "metadata": {},
   "outputs": [],
   "source": [
    "base = gv.tile_sources.EsriNatGeo.opts(width=1000, height=1000, padding=0.1)\n",
    "veg_dist_status = da_ann.z.where(da_ann['z']!=255).sel({'band':1})\n",
    "\n",
    "color_key = {\n",
    "    \"Confirmed, <50% ongoing\": \"#ffffb2\",\n",
    "    \"Confirmed, <50% completed\": \"#f45629\",\n",
    "    \"Confirmed, ≥50% ongoing\": \"#feb751\",\n",
    "    \"Confirmed, ≥50% completed\": \"#bd0026\",\n",
    "}\n",
    "\n",
    "levels = 4\n",
    "ticks = [2.5, 3.5, 4.5, 5.5]\n",
    "ticker = FixedTicker(ticks=ticks)\n",
    "labels = dict(zip(ticks, color_key))\n",
    "\n",
    "veg_dist_status.where(veg_dist_status!=0).hvplot.image(x='longitude', \n",
    "                             y='latitude', \n",
    "                             crs=crs, \n",
    "                             rasterize=True,\n",
    "                             dynamic=True, \n",
    "                             aspect='equal', \n",
    "                             frame_width=500, \n",
    "                             frame_height=500,\n",
    "                             clim=(2,6), alpha=0.8).opts(title=f\"VEG_DIST_STATUS\", xlabel='Longitude', \n",
    "                                                         ylabel='Latitude', color_levels = levels, \n",
    "                                                         cmap=tuple(color_key.values()),\n",
    "                                                         colorbar_opts={'ticker': ticker, 'major_label_overrides':labels}) * base"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37d1895b",
   "metadata": {},
   "source": [
    "### Compute Cumulative Areas of Disturbance \n",
    "We can use the DIST-ANN product to compute the total area of disturbance for a given region. To do so, we can use an interactive map and draw functionality provided by the open-source `ipyleaflet` library, just as we did for the DIST-ALERT analysis above.\n",
    "\n",
    "Below, use the draw tools to draw a rectangular region over an area that appears to show disturbance (we recommend the McKinney Fire area). Using custom functions in `dist_utils` we compute cumulative disturbance for each class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c0ef010",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Plot VEG-DIST-STATUS on ipyleaflet map ###\n",
    "dist_url ='https://glad.umd.edu/projects/opera/SEP/DIST-ANN/10/T/E/M/2022/OPERA_L3_DIST-ANN-HLS_T10TEM_2022_2023136T210508Z_30_v0_VEG-DIST-STATUS.tif'\n",
    "\n",
    "# Make leaflet map\n",
    "m = leafmap.Map(basemap=ipyleaflet.basemaps.Esri.WorldImagery,\n",
    "    center=(42.0, -122.3),\n",
    "    zoom=9,\n",
    "    crs=ipyleaflet.projections.EPSG3857,\n",
    "    draw_control=False)\n",
    "\n",
    "# Add raster and point data and draw functionality\n",
    "m.add_cog_layer(dist_url, name=\"VEG-DIST-STATUS\")\n",
    "m.add_points_from_xy('../supplemental_data/ca_fires_2022.csv', x ='Centroid (Lon)', y = 'Centroid (Lat)', layer_name=\"CaliforniaFires\")\n",
    "dc = ipyleaflet.DrawControl(\n",
    "    polygon={},\n",
    "    rectangle={\"shapeOptions\": {\"color\": \"blue\"}},\n",
    "    circlemarker={},\n",
    "    circle = {},\n",
    "    polyline={}\n",
    ")\n",
    "\n",
    "# Draw an AOI on an interactive map\n",
    "print('Select an Area of Interest using the tools on the left side of the map.')\n",
    "dc.on_draw(handle_draw)\n",
    "m.add_control(dc)\n",
    "display(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "355045fc",
   "metadata": {},
   "source": [
    "Now we can compute cumulative area statistics, below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c20b3e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use user-defined AOI for statistics\n",
    "from shapely import box\n",
    "ll = dc.last_draw['geometry']['coordinates'][0][0]\n",
    "ur = dc.last_draw['geometry']['coordinates'][0][2]\n",
    "aoi = box(ll[0], ll[1], ur[0], ur[1])\n",
    "\n",
    "# Create GeoDataFrame from user-specified bounding box\n",
    "gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[aoi])\n",
    "\n",
    "# Compute statistics within bounding box\n",
    "stats = leafmap.zonal_stats(gdf, dist_url, stats=['count'], categorical=True)[0]\n",
    "del stats['count']\n",
    "for i in [0,5,2,4,6]:\n",
    "    if i not in stats:\n",
    "        stats[i] = 0\n",
    "\n",
    "affected_areas = compute_areas([stats], pixel_area, product='ann', date=None)\n",
    "affected_areas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "618d3408",
   "metadata": {},
   "source": [
    "Above, you should see a dataframe containing rows for each `VEG-DIST-STATUS` class, their description, and the cumulative areas in square kilometers and hectares. How does this area compare to the official area reported in the table at the top of this notebook?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b428ce08",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "This notebook provides a tool for accessing and visulizing data of the OPERA DIST-ALERT and DIST-ANN products and quantifying the spatiotemporal impact of a wildfire on vegetation change. The DIST-ALERT data track the temporal record of vegetation disturbance, wheras the DIST-ANN data provide a summary of the disturbance that occured over the course of the calendar year. This notebook demonstrates how the `leaflet` and `ipyleaflet` libraries may be used to create enhanced visualizations over user-specified areas. For additional workflows, see the [OPERA Applications Github repository](https://github.com/OPERA-Cal-Val/OPERA_Applications)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}