#!/usr/bin/env python
# coding: utf-8

# # Visualizing the evolution of a wildfire using the OPERA DIST product
# ---
# 
# **This notebook serves as a visualization tool using the OPERA Land Surface Disturbance (DIST) product to illustrate the progression of an active wildfire in McKinney, Klamath National Forest, western Siskiyou County, California. **
# 
# <font color='red'>Note: This notebook uses provisional OPERA DIST products. Download the provisional data at https://www.jpl.nasa.gov/go/opera/products/dist-product-suite </font>

# In[ ]:


# Notebook dependencies
import math
import xarray as xr
import hvplot.xarray
import geoviews as gv
import pyproj
import numpy as np
import geopandas as gpd
import holoviews as hv
import panel.widgets as pnw
import datetime
from datetime import date, timedelta, datetime

from bokeh.models import FixedTicker
hv.extension('bokeh')

import warnings
warnings.filterwarnings('ignore')


# ---
# ## **Data Information Input**

# In the code cell below, the user should specify the:
# * Start and end days of interest <br>
# * Date iteration step<br>
# * Data directory<br>
# * Band list<br>
# * Anomaly threshold<br><br>
# 
# **<font color='red'>Note: The cell below is the only code in the notebook that should be modified. </font>**
# 

# In[ ]:


# Filter dates to focus only the relevant timeframe, e.g. McKinney Fire start date: Day 575 (July 29, 2022)
start_date = datetime(2022,8,1) # 29Jul2022 McKinney started
end_date = datetime(2022,8,15)  # 15Aug2022 Acquisition date of the input HLS
n = 2  # Interval day for the slider. e.g., every 2 days

data_dir ='https://opera-provisional-products.s3.us-west-2.amazonaws.com/DIST/DIST_HLS/WG/DIST-ALERT/McKinney_Wildfire/OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1/OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_'
bandlist = ['VEG-DIST-STATUS', 'VEG-DIST-DATE', 'VEG-ANOM-MAX']

# Maximum Anomaly threshold 
anom_threshold = 15 # filter pixels with maximum anomaly < 15%

# Focus the wildfire area extent calculation on a small region within the DIST product tile
bounds = [2000,3000,0,1000]   # upper-left y-pixel, lower-left y-pixel, upper-left x-pixel, lower-right x-pixel


# For the McKinney Wildfire event, we want to focus solely on the time period when the fire started (July 29, 2022). The wildfire's evolution is clear if we look over every 2 days between August 01 and 15, 2022. <br>
# *** 

# ### **<font color='red'> -- Do not modify any of the code below -- </font>**

# In[ ]:


# Global variables
total_pixel_count = 3660 * 3660
pixel_area = 30 * 30
ref_date = datetime(2020, 12, 31)
starting_day = (start_date - ref_date).days
ending_day = (end_date - ref_date).days + 1

bandpath = f"{data_dir}%s.tif"


# In[ ]:


# Functions
def stack_bands(bandpath:str, bandlist:list): 
    '''
    Returns geocube with three bands stacked into one multi-dimensional array.
            Parameters:
                    bandpath (str): Path to bands that should be stacked
                    bandlist (list): Three bands that should be stacked
            Returns:
                    bandStack (xarray Dataset): Geocube with stacked bands
                    crs (int): Coordinate Reference System corresponding to bands
    '''
    bandStack = []; bandS = []; bandStack_ = [];
    for i,band in enumerate(bandlist):
        if i==0:
            bandStack_ = xr.open_rasterio(bandpath%band)
            crs = pyproj.CRS.to_epsg(pyproj.CRS.from_proj4(bandStack_.crs))
            bandStack_ = bandStack_ * bandStack_.scales[0]
            bandStack = bandStack_.squeeze(drop=True)
            bandStack = bandStack.to_dataset(name='z')
            bandStack.coords['band'] = i+1
            bandStack = bandStack.rename({'x':'longitude', 'y':'latitude', 'band':'band'})
            bandStack = bandStack.expand_dims(dim='band')  
        else:
            bandS = xr.open_rasterio(bandpath%band)
            bandS = bandS * bandS.scales[0]
            bandS = bandS.squeeze(drop=True)
            bandS = bandS.to_dataset(name='z')
            bandS.coords['band'] = i+1
            bandS = bandS.rename({'x':'longitude', 'y':'latitude', 'band':'band'})
            bandS = bandS.expand_dims(dim='band')
            bandStack = xr.concat([bandStack, bandS], dim='band')
    return bandStack, crs

def time_and_area_cube(dist_status, dist_date, veg_anom_max, bounds, step=3):
    '''
    Returns geocube with time and area dimensions.
            Parameters:
                    dist_status (xarray DataArray): Disturbance Status band
                    anom_max (xarray DataArray): Maximum Anomaly band
                    dist_date (xarray DataArray): Disturbance Date band
                    step (int): Increment between each day in time series
            Returns:
                    wildfire_extent (xarray Dataset): Geocube with time and area dimensions
    '''
    lats = np.array(dist_status.latitude)
    lons = np.array(dist_status.longitude)
    expanded_array1 = []
    expanded_array2 = []
    respective_areas = {}
    
    for i in range(starting_day, ending_day, step):
        vg = dist_status.where((veg_anom_max > anom_threshold) & (dist_date > starting_day) & (dist_date <= i))
        extent_area = compute_area(vg.data,bounds)
        date = standard_date(str(i))
        coords =  {'lat': lats, 'lon': lons, 'time': date, 'area':extent_area}
        time_and_area = xr.DataArray(vg.data, coords=coords, dims=['lat', 'lon'])
        expanded_time_and_area = xr.concat([time_and_area], 'time')
        expanded_time_and_area = expanded_time_and_area.to_dataset(name='z')
        expanded_array2.append(expanded_time_and_area)
    area_extent = xr.concat(expanded_array2[:], dim='time')
    return area_extent

def compute_area(data_,bounds):
    '''
    Returns area of wildfire extent for single day.
            Parameters:
                    data (numpy array): Dist Status values (1.0-4.0)
            Returns:
                    fire_area (str): Wildfire extent area in kilometers squared
    '''
    data = data_[bounds[0]:bounds[1], bounds[2]:bounds[3]]
    fire_pixel_count = len(data[np.where(data>0)])
    fire_area = fire_pixel_count * pixel_area * pow(10, -6)
    fire_area = str(math.trunc(fire_area)) + " kilometers squared"
    return fire_area

def standard_date(day):
    '''
    Returns the inputted day number as a standard date.
            Parameters:
                    day (str): Day number that should be converted
            Returns:
                    res (str): Standard date corresponding to inputted day
    '''
    day.rjust(3 + len(day), '0')
    res_date = ref_date + timedelta(days=int(day))
    res = res_date.strftime("%m-%d-%Y")
    return res


# ### **<font color='red'> -- Do not modify any of the code above -- </font>**

# ---
# ## **1. Prepare the Geocube: Create and visualize the multidimensional dataset**

# In[ ]:


# Stack the bands to create a multidimensional "geocube"
da, crs = stack_bands(bandpath, bandlist)
da  # Inspect the geocube dataset created


# In[ ]:


# Visualize all the band stacked together
da.hvplot(x='longitude', 
          y='latitude', 
          crs=crs, 
          rasterize=True, 
          dynamic=True, 
          cmap='hot_r', 
          colorbar=False, 
          project=True, 
          frame_height=250, 
          shared_axes=False, 
          hover=True).layout()


# ---
# ### **Band 1: Vegetation Disturbance Status (VEG-DIST-STATUS)**
# 
# **Data Type:** UInt8<br>
# **Description:** Indication of vegetation cover loss (vegetation disturbance); "provisional" is used from the first detection until vegetation disturbance is detected for consecutive number of HLS scenes, when it is then labeled "confirmed."<br>
# 
# **Layer Values:**<br> 
# * **0:** No disturbance<br>
# * **1:** Provisional (**first detection**) Disturbance with vegetation cover change <50% <br>
# * **2:** Confirmed (**recurrent detection**) Disturbance with vegetation cover change < 50% <br> 
# * **3:** Provisional Disturbance with vegetation cover change ≥ 50% <br>
# * **4:** Confirmed Disturbance with vegetation cover change ≥ 50%  <br> 
# 
# ### **Band 2: Date of Initial Vegetation Disturbance (VEG-DIST-DATE)**
# 
# **Data Type:** Int16<br>
# **Description:** Day of first loss anomaly detection in the last year, defined as the 365 day period before the current date (366 for leap years). Day denoted as the number of days since December 31, 2020.<br>
# 
# ### **Band 3: Maximum Vegetation Anomaly (VEG-ANOM-MAX)** 
# 
# **Data Type:** UInt8<br>
# **Description:** Difference between historical and current year observed vegetation cover at the date of maximum decrease, measured on scale from 0-100%<br>
# 
# ---

# In[ ]:


# Visualize one band from the stack and overlay on a basemap
base = gv.tile_sources.EsriNatGeo.opts(width=1000, height=1000, padding=0.1)
da.z.sel({'band':1}).hvplot(x='longitude', 
                            y='latitude', 
                            crs=crs, 
                            rasterize=True, 
                            dynamic=True, 
                            cmap='hot_r', 
                            clim=(0,4), 
                            colorbar=True, 
                            project=True, 
                            frame_height=400, 
                            hover=True).opts(title='VEG_DIST_STATUS',
                                            colorbar_opts={'ticker': FixedTicker(ticks=[0, 1, 2, 3, 4])}) * base


# ---
# ## **2. Visualize the evolution of the McKinney Wildfire**

# In[ ]:


# Create a multidimensional array as a function of time and calculate the area of the burn as it evolves through time
area_coverage = time_and_area_cube(da.z.sel({'band':1}), da.z.sel({'band':2}), da.z.sel({'band':3}), bounds, step=n)
area_coverage


# In[ ]:


## Visualize the evolution of the fire overlaid on a basemap

# Create a basemap
base = gv.tile_sources.EsriNatGeo()

# Plot wildfire area extent time series
area_coverage_slider = area_coverage.z.interactive.sel(time=pnw.DiscreteSlider).hvplot(x='lon', 
                                                                                       y='lat', 
                                                                                       crs=crs, 
                                                                                       kind='image', 
                                                                                       rasterize=True, 
                                                                                       dynamic=True,
                                                                                       cmap='hot_r', 
                                                                                       aspect='equal', 
                                                                                       frame_width=600, 
                                                                                       frame_height=600).opts(active_tools=['wheel_zoom'],
                                                                                                              xlabel='Longitude', 
                                                                                                              ylabel='Latitude',  
                                                                                                              clim=(0,4), 
                                                                                                              colorbar_opts={'ticker': FixedTicker(ticks=[0, 1, 2, 3, 4])}, 
                                                                                                              alpha=0.9)
area_coverage_slider * base


# The plot above displays a time series **area extent** of the McKinney Wildfire. The slider iterates over every 2 days from August 01, 2022 through Aug 15, 2022.

# In[ ]:


## Plot the latest wildfire extent and overlay the latest perimter shapefile from National Interagency Fire Center (NIFC)
# https://data-nifc.opendata.arcgis.com/maps/wfigs-current-wildland-fire-perimeters

# Load shapefile of the McKinney perimeter from NIFC
shp_ = gpd.read_file("./aux_files/McKinney_NIFC.shp")
mckinney_nifc = gv.Polygons(shp_.geometry[0]).opts(line_color='red', color='None', line_width=2)

# Plot the shapefile the last time-step and the basemap together
base*mckinney_nifc*area_coverage.z.interactive.sel(time='08-15-2022').hvplot(x='lon', 
                                                                             y='lat', 
                                                                             crs=crs, 
                                                                             geo='True', 
                                                                             kind='image', 
                                                                             rasterize=True, 
                                                                             cmap='hot_r', 
                                                                             aspect='equal', 
                                                                             frame_height=600, 
                                                                             frame_width=600, 
                                                                             alpha=0.9).opts(active_tools=['wheel_zoom'], 
                                                                                             xlabel='Longitude', 
                                                                                             ylabel='Latitude',  
                                                                                             clim=(1,4), 
                                                                                             colorbar_opts={'ticker': FixedTicker(ticks=[1, 2, 3, 4])}, 
                                                                                             alpha=0.9)


# **Layer Values:**
# * **0:** No disturbance<br>
# * **1:** Provisional (**first detection**) Disturbance with vegetation cover change <50% <br>
# * **2:** Confirmed (**recurrent detection**) Disturbance with vegetation cover change < 50% <br> 
# * **3:** Provisional Disturbance with vegetation cover change ≥ 50% <br>
# * **4:** Confirmed Disturbance with vegetation cover change ≥ 50%  <br> 

# ---
# ## **Summary** 
# 
# The McKinney wildfire started on Friday July 29th, 2022 approx. 02:15 PM in the Klamath National Forest, Siskiyou County California. As of September 1, 2022, the fire is already 99% contained with estimated 60,138 acres (~243 km<sup>2</sup>) of area burned. <br><br>
# In this Jupyter notebook, we showcase how end-users can leverage the **OPERA Land Surface Disturbance (DIST)** product **to interactively visualize the evolution and calculate the extent of a wildfire**.
# <br>
# ***

# In[ ]: