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.
A NASA Earthdata Login account is required to download the data used in this tutorial. You can create an account at the link provided.
Note 1: This notebook uses provisional products, which may differ slightly from operational products. Please refer to DIST product specification for more information.
Note2: DIST products are distributed via NASA's Distributed Active Archive Centers (DAACs), specifically the LP DAAC. 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.
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:
| Name | Alarm Date | Containment Date | Centroid (lat) | Centroid (lon) | Official Area Impacted (acres) |
|---|---|---|---|---|---|
| Gulch | 3/12/2022 | 3/14/2022 | 41.882 | -121.849 | 113 |
| Mountain | 9/2/2022 | 10/2/2022 | 41.462 | -122.66 | 13440 |
| McKinney | 7/29/2022 | 11/11/2022 | 41.814 | -122.841 | 60102 |
| Alex | 7/31/2022 | 8/3/2022 | 41.939 | -122.952 | 150 |
| Mill | 9/2/2022 | 9/11/2022 | 41.476 | -122.393 | 3939 |
| Coyote | 9/7/2022 | 9/8/2022 | 41.834 | -121.793 | 296 |
Data from the California Data Portal
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.
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. Note: The DIST-ANN product used here is a provisional product with an observation date exceeding the operational 1-year period.
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.
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.
DIST product data are stored and distributed via NASA's Distributed Active Archive Centers (DAACs), specifically the LP DAAC. 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.
First we import the necessary Python libraries. These come pre-installed in the opera_app anaconda environement within the OPERA Applications Gitub repository. We also import a collection of custom DIST-specific functions from a source file called dist_utils.py.
# Notebook dependencies
import hvplot.xarray
import geoviews as gv
import matplotlib.pyplot as plt
import geopandas as gpd
from datetime import datetime
import ipyleaflet
import leafmap.leafmap as leafmap
from bokeh.models import FixedTicker
import holoviews as hv
hv.extension('bokeh')
import warnings
warnings.filterwarnings('ignore')
import sys
sys.path.append('../../')
from src.dist_utils import *
The next cell configures gdal and provideds necessary authentication to successfully access LP DAAC cloud-hosted assets. 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.
# Check for valid Earthdata credentials. If entering for the first time, be very careful to enter username/password correctly!
check_netrc()
# Set GDAL configs to successfully access LP DAAC Cloud Assets via vsicurl
gdal.SetConfigOption("GDAL_HTTP_UNSAFESSL", "YES")
gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/cookies.txt')
gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/cookies.txt')
gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','FALSE')
gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF')
The OPERA Disturbance product is derived from HLS. The time series of HLS is what provides the quantification of vegetation disturbance.
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.
# Path to data location on LP DAAC and provide list of desired bands
hls1 = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T10TEM.2022202T185058.v2.0/HLS.L30.T10TEM.2022202T185058.v2.0.'
hls2 = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T10TEM.2022266T185119.v2.0/HLS.L30.T10TEM.2022266T185119.v2.0.'
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'
bandlist = ['B05','B04','B03','B02']
# Make true color image
make_hls_true_color(hls1, bandlist, 'true_color_7_21_22_NorCal.tif')
make_hls_true_color(hls2, bandlist, 'true_color_9_23_22_NorCal.tif')
# Make false color image
make_hls_false_color(hls1, bandlist, 'false_color_7_21_22_NorCal.tif')
make_hls_false_color(hls2, bandlist, 'false_color_9_23_22_NorCal.tif')
# Make ndvi image
make_hls_ndvi(hls1, bandlist, 'ndvi_7_21_22_NorCal.tif')
make_hls_ndvi(hls2, bandlist, 'ndvi_9_23_22_NorCal.tif')
# Make VEG-DIST-STATUS image
make_veg_dist_status_visual(veg_dist_status, 'veg_dist_status_NorCal.tif')
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.
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.
m = leafmap.Map(center=[42.0, -122.3], zoom=10)
m.split_map(
left_layer='tifs/true_color_7_21_22_NorCal.tif',
right_layer='tifs/true_color_9_23_22_NorCal.tif',
left_label="7-21-2022",
right_label="9-23-2022"
)
m
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).
m = leafmap.Map(center=[42.0, -122.3], zoom=10)
m.split_map(
left_layer='tifs/false_color_7_21_22_NorCal.tif',
right_layer='tifs/false_color_9_23_22_NorCal.tif',
left_label='7-21-2022',
right_label="9-23-2022"
)
m
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. 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.
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.
m = leafmap.Map(center=[42.0, -122.3], zoom=10)
m.split_map(
left_layer='tifs/ndvi_7_21_22_NorCal.tif',
right_layer='tifs/ndvi_9_23_22_NorCal.tif',
left_label="7-21-2022",
right_label="9-23-2022"
)
m
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?
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.
m = leafmap.Map(center=[42.0, -122.3], zoom=10)
m.split_map(
left_layer='tifs/ndvi_9_23_22_NorCal.tif',
right_layer='tifs/veg_dist_status_NorCal.tif',
left_label="9-23-2022 NDVI",
right_label="2022 DIST-ANN VEG-DIST-STATUS"
)
m
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.
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.
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.
### Open file containing paths to DIST-ALERT for Northern California ###
# Specify location of data and open as a list
dist_alert_tiles_local = '../links_to_data/10TEM_DIST-ALERT_links.txt'
file = open(dist_alert_tiles_local, "r")
lines = file.readlines()
file.close()
# Get only the VEG-DIST-STATUS and VEG-IND paths
veg_dist_status = []
veg_ind = []
for line in lines:
fp = line[:-1]
# Check if the file path ends with 'VEG-DIST-STATUS.tif' or 'VEG-IND' and add to corresponding list
if fp.endswith('VEG-DIST-STATUS.tif'):
veg_dist_status.append(line.strip())
elif fp.endswith('VEG-IND.tif'):
veg_ind.append(line.strip())
# 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
veg_dist_status = veg_dist_status[0::10]
veg_ind = veg_ind[0::10]
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 and ipyleaflet docs.
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.
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.
### Plot VEG-DIST-STATUS on ipyleaflet map ###
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'
# Make leaflet map
m = leafmap.Map(basemap=ipyleaflet.basemaps.Esri.WorldImagery,
center=(42.0, -122.3),
zoom=9,
crs=ipyleaflet.projections.EPSG3857,
draw_control=False)
# Add raster and point data and draw functionality
m.add_cog_layer(dist_url, name="VEG-DIST-STATUS")
m.add_points_from_xy('../supplemental_data/ca_fires_2022.csv', x ='Centroid (Lon)', y = 'Centroid (Lat)', layer_name="CaliforniaFires")
dc = ipyleaflet.DrawControl(
polygon={},
rectangle={"shapeOptions": {"color": "blue"}},
circlemarker={},
circle = {},
polyline={}
)
# Draw an AOI on an interactive map
print('Select an Area of Interest using the tools on the left side of the map.')
dc.on_draw(handle_draw)
m.add_control(dc)
display(m)
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.
### Compute the area of each class within the bounding box through time ###
# User-Defined Parameters
from shapely import box
ll = dc.last_draw['geometry']['coordinates'][0][0]
ur = dc.last_draw['geometry']['coordinates'][0][2]
aoi = box(ll[0], ll[1], ur[0], ur[1])
# Make GeoDataFrame from AOI
gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[aoi])
# Compute the affected areas through time
affected_areas_through_time = []
veg_ind_through_time = []
dates = []
for i,tile in enumerate(veg_dist_status[1:]):
date = extract_date_from_string(tile)
date = datetime.strptime(date, '%Y%j')
area_stats = leafmap.zonal_stats(gdf, tile, stats=['count'], categorical=True)[0]
ind_stats = leafmap.zonal_stats(gdf, veg_ind[i], stats=['mean'], categorical=False)[0]
if ind_stats['mean'] is not None:
dates.append(date)
veg_ind_through_time.append(ind_stats['mean'])
del area_stats['count']
for status in range(5):
if status not in area_stats:
area_stats[status] = 0
area_stats = dict(sorted(area_stats.items()))
pixel_area = 30 * 30
affected_areas_through_time.append(compute_areas([area_stats], pixel_area, 'alert', date))
# Make Pandas dataframes for the computed statistics
combined_areas = pd.concat(affected_areas_through_time, axis=0)
veg_indicators = pd.DataFrame({'Date':dates, 'Vegetation Indicator':veg_ind_through_time})
veg_indicators.set_index('Date', inplace=True)
Now we can plot the results graphically using plot functionality within the open-source pandas library, whose docs are available here.
# Group the area data by "Date" and "Class"
grouped_data = combined_areas.groupby(['Date', 'VEG-DIST-STATUS Class'])['Area (km2)'].sum()
# Reset the index of the grouped data to turn it back into a DataFrame
grouped_df = grouped_data.reset_index()
# Pivot the data to have "Class" values as columns and "Date" as index
pivot_df = grouped_df.pivot(index='Date', columns='VEG-DIST-STATUS Class', values='Area (km2)')
# Plot the data
fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(15,5))
pivot_df.plot(ax=ax1,marker='.')
veg_indicators.plot(ax=ax2,marker='.')
# Add labels and title
ax1.set_xlabel('Date')
ax1.set_ylabel('Area (km2)')
ax2.set_xlabel('Date')
ax2.set_ylabel('Vegetation Indicator')
plt.show()
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.
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.
# Specify filepath location of OPERA DIST_ANN tile and desired bands (Note: bandlist is not comprehensive of all available layers)
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_'
bandlist = ['VEG-DIST-STATUS', 'VEG-HIST', 'VEG-IND-MAX','VEG-ANOM-MAX', 'VEG-DIST-CONF', 'VEG-DIST-DATE',
'VEG-DIST-COUNT', 'VEG-DIST-DUR', 'VEG-LAST-DATE', 'GEN-DIST-STATUS', 'GEN-ANOM-MAX', 'GEN-DIST-CONF',
'GEN-DIST-DATE', 'GEN-DIST-COUNT', 'GEN-DIST-DUR', 'GEN-LAST-DATE']
bandpath = f"{data_dir}%s.tif"
# Create geocube of stacked bands
da_ann, crs = stack_bands(bandpath, bandlist)
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.
base = gv.tile_sources.EsriNatGeo.opts(width=1000, height=1000, padding=0.1)
veg_dist_status = da_ann.z.where(da_ann['z']!=255).sel({'band':1})
color_key = {
"Confirmed, <50% ongoing": "#ffffb2",
"Confirmed, <50% completed": "#f45629",
"Confirmed, ≥50% ongoing": "#feb751",
"Confirmed, ≥50% completed": "#bd0026",
}
levels = 4
ticks = [2.5, 3.5, 4.5, 5.5]
ticker = FixedTicker(ticks=ticks)
labels = dict(zip(ticks, color_key))
veg_dist_status.where(veg_dist_status!=0).hvplot.image(x='longitude',
y='latitude',
crs=crs,
rasterize=True,
dynamic=True,
aspect='equal',
frame_width=500,
frame_height=500,
clim=(2,6), alpha=0.8).opts(title=f"VEG_DIST_STATUS", xlabel='Longitude',
ylabel='Latitude', color_levels = levels,
cmap=tuple(color_key.values()),
colorbar_opts={'ticker': ticker, 'major_label_overrides':labels}) * base
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.
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.
### Plot VEG-DIST-STATUS on ipyleaflet map ###
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'
# Make leaflet map
m = leafmap.Map(basemap=ipyleaflet.basemaps.Esri.WorldImagery,
center=(42.0, -122.3),
zoom=9,
crs=ipyleaflet.projections.EPSG3857,
draw_control=False)
# Add raster and point data and draw functionality
m.add_cog_layer(dist_url, name="VEG-DIST-STATUS")
m.add_points_from_xy('../supplemental_data/ca_fires_2022.csv', x ='Centroid (Lon)', y = 'Centroid (Lat)', layer_name="CaliforniaFires")
dc = ipyleaflet.DrawControl(
polygon={},
rectangle={"shapeOptions": {"color": "blue"}},
circlemarker={},
circle = {},
polyline={}
)
# Draw an AOI on an interactive map
print('Select an Area of Interest using the tools on the left side of the map.')
dc.on_draw(handle_draw)
m.add_control(dc)
display(m)
Now we can compute cumulative area statistics, below.
# Use user-defined AOI for statistics
from shapely import box
ll = dc.last_draw['geometry']['coordinates'][0][0]
ur = dc.last_draw['geometry']['coordinates'][0][2]
aoi = box(ll[0], ll[1], ur[0], ur[1])
# Create GeoDataFrame from user-specified bounding box
gdf = gpd.GeoDataFrame(index=[0], crs='epsg:4326', geometry=[aoi])
# Compute statistics within bounding box
stats = leafmap.zonal_stats(gdf, dist_url, stats=['count'], categorical=True)[0]
del stats['count']
for i in [0,5,2,4,6]:
if i not in stats:
stats[i] = 0
affected_areas = compute_areas([stats], pixel_area, product='ann', date=None)
affected_areas
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?
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.