This notebook serves as an introduction to the OPERA Land Surface Disturbance (DIST) product and visualizing relevant raster layers for wildfire applications.
Download the provisional data at https://www.jpl.nasa.gov/go/opera/products/dist-product-suite
# Notebook dependencies
import xarray as xr
import hvplot.xarray
import geoviews as gv
import pyproj
from pyproj import Proj
import holoviews as hv
from bokeh.models import FixedTicker
hv.extension('bokeh')
import warnings
warnings.filterwarnings('ignore')
The DIST product maps per pixel vegetation disturbance (specifically, vegetation cover loss) from the Harmonized Landsat-8 and Sentinel-2 A/B (HLS) scenes. Vegetation disturbance is mapped when there is an indicated decrease in vegetation cover within an HLS pixel. This notebook focuses on relevant layers within the DIST_ALERT product for wildfire applications, which is released at the cadence of HLS imagery.
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 kilometers squared divided into 3660 rows and 3660 columns at 30 meter pixel spacing. Each tile overlaps neighbors by 4900 meters in each direction.
The DIST_ALERT product is distributed as a set of Cloud-Optimized GeoTIFF (COG) files to enable download of only particular layers of interest to a given user. All L3 DIST layers are stored in files following GeoTIFF format specifications.
# User-specified input
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-ANOM-MAX', 'VEG-DIST-DATE', 'VEG-DIST-STATUS']
bandpath = f"{data_dir}%s.tif"
# 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
# Creates geocube of stacked bands
da, crs = stack_bands(bandpath, bandlist)
# Creates basemap
base = gv.tile_sources.EsriTerrain.opts(width=1000, height=1000, padding=0.1)
Data Type: UInt8
Description: Difference between historical and current year observed vegetation cover at the date of maximum decrease, measured on scale from 0-100%
base = gv.tile_sources.EsriNatGeo.opts(width=1000, height=1000, padding=0.1)
veg_anom_max = da.z.sel({'band':1})
veg_anom_max.hvplot.image(x='longitude',
y='latitude',
crs=crs,
rasterize=True,
dynamic=True,
aspect='equal',
frame_width=500,
frame_height=500,
cmap='hot_r',
clim=(0,100), alpha=0.8).opts(title=f"VEG_ANOM_MAX", xlabel='Longitude', ylabel='Latitude') * base
Data Type: Int16
Description: Day of first loss anomaly detection in the last year, denoted as the number of days since December 31st, 2020.
veg_dist_date = da.z.sel({'band':2})
veg_dist_date.hvplot.image(x='longitude',
y='latitude',
crs=crs,
rasterize=True,
dynamic=True,
aspect='equal',
frame_width=500,
frame_height=500,
cmap='hot_r',
alpha=0.8).opts(title=f"VEG_DIST_DATE", xlabel='Longitude', ylabel='Latitude',clim=(0, 592)) * base
Data Type: UInt8
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."
veg_dist_status = da.z.sel({'band':3})
veg_dist_status.hvplot.image(x='longitude',
y='latitude',
crs=crs,
rasterize=True,
dynamic=True,
aspect='equal',
frame_width=500,
frame_height=500,
cmap='hot_r',
alpha=0.8).opts(title=f"VEG_DIST_STATUS", clim=(0,4), colorbar_opts={'ticker': FixedTicker(ticks=[0, 1, 2, 3, 4])}, xlabel='Longitude', ylabel='Latitude') * base
Layer Values: