{ "cells": [ { "cell_type": "markdown", "id": "4686ee13", "metadata": {}, "source": [ "# Visualizing the evolution of a wildfire using the OPERA DIST product\n", "---\n", "\n", "**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. **\n", "\n", "Note: This notebook uses provisional OPERA DIST products. Download the provisional data at https://www.jpl.nasa.gov/go/opera/products/dist-product-suite " ] }, { "cell_type": "code", "execution_count": null, "id": "a3dfcb53", "metadata": {}, "outputs": [], "source": [ "# Notebook dependencies\n", "import math\n", "import xarray as xr\n", "import hvplot.xarray\n", "import geoviews as gv\n", "import pyproj\n", "import numpy as np\n", "import geopandas as gpd\n", "import holoviews as hv\n", "import panel.widgets as pnw\n", "import datetime\n", "from datetime import date, timedelta, datetime\n", "\n", "from bokeh.models import FixedTicker\n", "hv.extension('bokeh')\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "markdown", "id": "be05d52a", "metadata": {}, "source": [ "---\n", "## **Data Information Input**" ] }, { "cell_type": "markdown", "id": "2adcbb74", "metadata": {}, "source": [ "In the code cell below, the user should specify the:\n", "* Start and end days of interest
\n", "* Date iteration step
\n", "* Data directory
\n", "* Band list
\n", "* Anomaly threshold

\n", "\n", "**Note: The cell below is the only code in the notebook that should be modified. **\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1b7e6115", "metadata": {}, "outputs": [], "source": [ "# Filter dates to focus only the relevant timeframe, e.g. McKinney Fire start date: Day 575 (July 29, 2022)\n", "start_date = datetime(2022,8,1) # 29Jul2022 McKinney started\n", "end_date = datetime(2022,8,15) # 15Aug2022 Acquisition date of the input HLS\n", "n = 2 # Interval day for the slider. e.g., every 2 days\n", "\n", "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_'\n", "bandlist = ['VEG-DIST-STATUS', 'VEG-DIST-DATE', 'VEG-ANOM-MAX']\n", "\n", "# Maximum Anomaly threshold \n", "anom_threshold = 15 # filter pixels with maximum anomaly < 15%\n", "\n", "# Focus the wildfire area extent calculation on a small region within the DIST product tile\n", "bounds = [2000,3000,0,1000] # upper-left y-pixel, lower-left y-pixel, upper-left x-pixel, lower-right x-pixel" ] }, { "cell_type": "markdown", "id": "e5dd4b17", "metadata": {}, "source": [ "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.
\n", "*** " ] }, { "cell_type": "markdown", "id": "355272c8", "metadata": {}, "source": [ "### ** -- Do not modify any of the code below -- **" ] }, { "cell_type": "code", "execution_count": null, "id": "5d075d54", "metadata": {}, "outputs": [], "source": [ "# Global variables\n", "total_pixel_count = 3660 * 3660\n", "pixel_area = 30 * 30\n", "ref_date = datetime(2020, 12, 31)\n", "starting_day = (start_date - ref_date).days\n", "ending_day = (end_date - ref_date).days + 1\n", "\n", "bandpath = f\"{data_dir}%s.tif\"" ] }, { "cell_type": "code", "execution_count": null, "id": "ae2f7ce4", "metadata": {}, "outputs": [], "source": [ "# Functions\n", "def stack_bands(bandpath:str, bandlist:list): \n", " '''\n", " Returns geocube with three bands stacked into one multi-dimensional array.\n", " Parameters:\n", " bandpath (str): Path to bands that should be stacked\n", " bandlist (list): Three bands that should be stacked\n", " Returns:\n", " bandStack (xarray Dataset): Geocube with stacked bands\n", " crs (int): Coordinate Reference System corresponding to bands\n", " '''\n", " bandStack = []; bandS = []; bandStack_ = [];\n", " for i,band in enumerate(bandlist):\n", " if i==0:\n", " bandStack_ = xr.open_rasterio(bandpath%band)\n", " crs = pyproj.CRS.to_epsg(pyproj.CRS.from_proj4(bandStack_.crs))\n", " bandStack_ = bandStack_ * bandStack_.scales[0]\n", " bandStack = bandStack_.squeeze(drop=True)\n", " bandStack = bandStack.to_dataset(name='z')\n", " bandStack.coords['band'] = i+1\n", " bandStack = bandStack.rename({'x':'longitude', 'y':'latitude', 'band':'band'})\n", " bandStack = bandStack.expand_dims(dim='band') \n", " else:\n", " bandS = xr.open_rasterio(bandpath%band)\n", " bandS = bandS * bandS.scales[0]\n", " bandS = bandS.squeeze(drop=True)\n", " bandS = bandS.to_dataset(name='z')\n", " bandS.coords['band'] = i+1\n", " bandS = bandS.rename({'x':'longitude', 'y':'latitude', 'band':'band'})\n", " bandS = bandS.expand_dims(dim='band')\n", " bandStack = xr.concat([bandStack, bandS], dim='band')\n", " return bandStack, crs\n", "\n", "def time_and_area_cube(dist_status, dist_date, veg_anom_max, bounds, step=3):\n", " '''\n", " Returns geocube with time and area dimensions.\n", " Parameters:\n", " dist_status (xarray DataArray): Disturbance Status band\n", " anom_max (xarray DataArray): Maximum Anomaly band\n", " dist_date (xarray DataArray): Disturbance Date band\n", " step (int): Increment between each day in time series\n", " Returns:\n", " wildfire_extent (xarray Dataset): Geocube with time and area dimensions\n", " '''\n", " lats = np.array(dist_status.latitude)\n", " lons = np.array(dist_status.longitude)\n", " expanded_array1 = []\n", " expanded_array2 = []\n", " respective_areas = {}\n", " \n", " for i in range(starting_day, ending_day, step):\n", " vg = dist_status.where((veg_anom_max > anom_threshold) & (dist_date > starting_day) & (dist_date <= i))\n", " extent_area = compute_area(vg.data,bounds)\n", " date = standard_date(str(i))\n", " coords = {'lat': lats, 'lon': lons, 'time': date, 'area':extent_area}\n", " time_and_area = xr.DataArray(vg.data, coords=coords, dims=['lat', 'lon'])\n", " expanded_time_and_area = xr.concat([time_and_area], 'time')\n", " expanded_time_and_area = expanded_time_and_area.to_dataset(name='z')\n", " expanded_array2.append(expanded_time_and_area)\n", " area_extent = xr.concat(expanded_array2[:], dim='time')\n", " return area_extent\n", "\n", "def compute_area(data_,bounds):\n", " '''\n", " Returns area of wildfire extent for single day.\n", " Parameters:\n", " data (numpy array): Dist Status values (1.0-4.0)\n", " Returns:\n", " fire_area (str): Wildfire extent area in kilometers squared\n", " '''\n", " data = data_[bounds[0]:bounds[1], bounds[2]:bounds[3]]\n", " fire_pixel_count = len(data[np.where(data>0)])\n", " fire_area = fire_pixel_count * pixel_area * pow(10, -6)\n", " fire_area = str(math.trunc(fire_area)) + \" kilometers squared\"\n", " return fire_area\n", "\n", "def standard_date(day):\n", " '''\n", " Returns the inputted day number as a standard date.\n", " Parameters:\n", " day (str): Day number that should be converted\n", " Returns:\n", " res (str): Standard date corresponding to inputted day\n", " '''\n", " day.rjust(3 + len(day), '0')\n", " res_date = ref_date + timedelta(days=int(day))\n", " res = res_date.strftime(\"%m-%d-%Y\")\n", " return res" ] }, { "cell_type": "markdown", "id": "61c1d66c", "metadata": {}, "source": [ "### ** -- Do not modify any of the code above -- **" ] }, { "cell_type": "markdown", "id": "dae1e499", "metadata": {}, "source": [ "---\n", "## **1. Prepare the Geocube: Create and visualize the multidimensional dataset**" ] }, { "cell_type": "code", "execution_count": null, "id": "45a194a5", "metadata": {}, "outputs": [], "source": [ "# Stack the bands to create a multidimensional \"geocube\"\n", "da, crs = stack_bands(bandpath, bandlist)\n", "da # Inspect the geocube dataset created" ] }, { "cell_type": "code", "execution_count": null, "id": "bba53fbb", "metadata": {}, "outputs": [], "source": [ "# Visualize all the band stacked together\n", "da.hvplot(x='longitude', \n", " y='latitude', \n", " crs=crs, \n", " rasterize=True, \n", " dynamic=True, \n", " cmap='hot_r', \n", " colorbar=False, \n", " project=True, \n", " frame_height=250, \n", " shared_axes=False, \n", " hover=True).layout()\n" ] }, { "cell_type": "markdown", "id": "38814edc", "metadata": {}, "source": [ "---\n", "### **Band 1: Vegetation Disturbance Status (VEG-DIST-STATUS)**\n", "\n", "**Data Type:** UInt8
\n", "**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.\"
\n", "\n", "**Layer Values:**
\n", "* **0:** No disturbance
\n", "* **1:** Provisional (**first detection**) Disturbance with vegetation cover change <50%
\n", "* **2:** Confirmed (**recurrent detection**) Disturbance with vegetation cover change < 50%
\n", "* **3:** Provisional Disturbance with vegetation cover change ≥ 50%
\n", "* **4:** Confirmed Disturbance with vegetation cover change ≥ 50%
\n", "\n", "### **Band 2: Date of Initial Vegetation Disturbance (VEG-DIST-DATE)**\n", "\n", "**Data Type:** Int16
\n", "**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.
\n", "\n", "### **Band 3: Maximum Vegetation Anomaly (VEG-ANOM-MAX)** \n", "\n", "**Data Type:** UInt8
\n", "**Description:** Difference between historical and current year observed vegetation cover at the date of maximum decrease, measured on scale from 0-100%
\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": null, "id": "77307aea", "metadata": {}, "outputs": [], "source": [ "# Visualize one band from the stack and overlay on a basemap\n", "base = gv.tile_sources.EsriNatGeo.opts(width=1000, height=1000, padding=0.1)\n", "da.z.sel({'band':1}).hvplot(x='longitude', \n", " y='latitude', \n", " crs=crs, \n", " rasterize=True, \n", " dynamic=True, \n", " cmap='hot_r', \n", " clim=(0,4), \n", " colorbar=True, \n", " project=True, \n", " frame_height=400, \n", " hover=True).opts(title='VEG_DIST_STATUS',\n", " colorbar_opts={'ticker': FixedTicker(ticks=[0, 1, 2, 3, 4])}) * base\n" ] }, { "cell_type": "markdown", "id": "519c7207", "metadata": {}, "source": [ "---\n", "## **2. Visualize the evolution of the McKinney Wildfire**" ] }, { "cell_type": "code", "execution_count": null, "id": "466c265b", "metadata": {}, "outputs": [], "source": [ "# Create a multidimensional array as a function of time and calculate the area of the burn as it evolves through time\n", "area_coverage = time_and_area_cube(da.z.sel({'band':1}), da.z.sel({'band':2}), da.z.sel({'band':3}), bounds, step=n)\n", "area_coverage" ] }, { "cell_type": "code", "execution_count": null, "id": "5df98bc2", "metadata": { "scrolled": false }, "outputs": [], "source": [ "## Visualize the evolution of the fire overlaid on a basemap\n", "\n", "# Create a basemap\n", "base = gv.tile_sources.EsriNatGeo()\n", "\n", "# Plot wildfire area extent time series\n", "area_coverage_slider = area_coverage.z.interactive.sel(time=pnw.DiscreteSlider).hvplot(x='lon', \n", " y='lat', \n", " crs=crs, \n", " kind='image', \n", " rasterize=True, \n", " dynamic=True,\n", " cmap='hot_r', \n", " aspect='equal', \n", " frame_width=600, \n", " frame_height=600).opts(active_tools=['wheel_zoom'],\n", " xlabel='Longitude', \n", " ylabel='Latitude', \n", " clim=(0,4), \n", " colorbar_opts={'ticker': FixedTicker(ticks=[0, 1, 2, 3, 4])}, \n", " alpha=0.9)\n", "area_coverage_slider * base" ] }, { "cell_type": "markdown", "id": "b630585a", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "id": "25f79dc4", "metadata": {}, "outputs": [], "source": [ "## Plot the latest wildfire extent and overlay the latest perimter shapefile from National Interagency Fire Center (NIFC)\n", "# https://data-nifc.opendata.arcgis.com/maps/wfigs-current-wildland-fire-perimeters\n", "\n", "# Load shapefile of the McKinney perimeter from NIFC\n", "shp_ = gpd.read_file(\"./aux_files/McKinney_NIFC.shp\")\n", "mckinney_nifc = gv.Polygons(shp_.geometry[0]).opts(line_color='red', color='None', line_width=2)\n", "\n", "# Plot the shapefile the last time-step and the basemap together\n", "base*mckinney_nifc*area_coverage.z.interactive.sel(time='08-15-2022').hvplot(x='lon', \n", " y='lat', \n", " crs=crs, \n", " geo='True', \n", " kind='image', \n", " rasterize=True, \n", " cmap='hot_r', \n", " aspect='equal', \n", " frame_height=600, \n", " frame_width=600, \n", " alpha=0.9).opts(active_tools=['wheel_zoom'], \n", " xlabel='Longitude', \n", " ylabel='Latitude', \n", " clim=(1,4), \n", " colorbar_opts={'ticker': FixedTicker(ticks=[1, 2, 3, 4])}, \n", " alpha=0.9)\n" ] }, { "cell_type": "markdown", "id": "3dec64a4", "metadata": {}, "source": [ "**Layer Values:**\n", "* **0:** No disturbance
\n", "* **1:** Provisional (**first detection**) Disturbance with vegetation cover change <50%
\n", "* **2:** Confirmed (**recurrent detection**) Disturbance with vegetation cover change < 50%
\n", "* **3:** Provisional Disturbance with vegetation cover change ≥ 50%
\n", "* **4:** Confirmed Disturbance with vegetation cover change ≥ 50%
" ] }, { "cell_type": "markdown", "id": "afe9cbfb", "metadata": {}, "source": [ "---\n", "## **Summary** \n", "\n", "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 km2) of area burned.

\n", "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**.\n", "
\n", "***" ] }, { "cell_type": "code", "execution_count": null, "id": "d033d7d1", "metadata": {}, "outputs": [], "source": [] } ], "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.7.12" } }, "nbformat": 4, "nbformat_minor": 5 }