Raster Analytics - Calculate wildfire landslide risk¶
In October 2017, wildfires raged through Sonoma and Napa counties, devastating surrounding communities. In the wake of these fires, the burn scars could cause further risk to public safety from a different kind of disaster: landslides. Post-fire landslides are particularly hazardous because there is more erosion and weaker soil in burned areas without vegetation to anchor the topsoil.
Groups handling rehabilitation, emergency planning and mitigation after a wildfire need to assess the vulnerability of the landscape to landslides. In this notebook, we will provide local emergency management teams a summary of post-wildfire landslide risk, so officials can target mitigation efforts to the most vulnerable watershed basins.
We will use the imagery layers to assess landslide risk per watershed within the burn area. We will create a landslide risk map and then summarize the landslide risk based on watershed sub-basins. We will use raster function chains to derive a burn severity map, a topographic slope map, and a landcover index map. These individual processing chains will be combined into one processing chain for distributed processing on the Raster Analytics server and then be summarized by watershed sub-basins.
Import required libraries¶
import arcgis
from arcgis.gis import GIS
from arcgis.raster.functions import *
from ipywidgets import *
gis = GIS(url='https://pythonapi.playground.esri.com/portal', username='arcgis_python', password='amazing_arcgis_123')
arcgis.raster.analytics.is_supported(gis)
Get data¶
For this analysis we need the following datasets
- a Landsat 8 imagery for before (Before_L8)
- a Landast 8 imagery for after (After_L8) the wildfire
- a DEM (digital elevation model) showing the elevation of the terrain
- a NLCD (National Landcover Dataset) showing land use and predominant vegetation type
- a watershed basin dataset
In the cells below, we access these datasets from the GIS
before_l8 = gis.content.search('title:Before_L8 owner:api_data_owner',
item_type = "Image Service",
outside_org=True)[0].layers[0]
after_l8 = gis.content.search('title:After_L8 owner:api_data_owner',
item_type = "Image Service",
outside_org=True)[0].layers[0]
before_l8
dem = gis.content.search('title:Sonoma_DEM owner:api_data_owner',
item_type = "Image Service",
outside_org=True)[0].layers[0]
nlcd = gis.content.search('title:Sonoma_NLCD2011 owner:api_data_owner',
item_type = "Image Service",
outside_org=True)[0].layers[0]
basins = gis.content.search('title:Sonoma_Basins owner:api_data_owner',
item_type = "Image Service",
outside_org=True)[0].layers[0]
# A preview of National Landcover Dataset layer
nlcd
Create a burn severity map¶
To compare the burn scars on the before and after Landsat imagery, we’ll choose the multispectral bands 5
,3
,2
to be displayed. The [5,3,2] band combination improves visibility of fire and burn scars. Healthy vegetation is shown in bright red, while stressed vegetation is displayed as dull red. Nonvegetated features such as bare and urban areas are displayed in various shades of gray and blue.
Below, we apply the same bands combination to the before_l8 and after_l8 layers.
infrared_before = extract_band(before_l8,
band_names = ['sr_band5','sr_band3','sr_band2'])
infrared_after = extract_band(after_l8,
band_names = ['sr_band5','sr_band3','sr_band2'])
Visual Assessment¶
Below, in order to visually compare the burn effects, we create two maps and load the extracted bands of before and after imagery.
# Create two maps to compare before and after imageries side by side
map1 = gis.map(location='-122.58, 38.45', zoomlevel=10)
map2 = gis.map(location='-122.58, 38.45', zoomlevel=10)
map1.layout = Layout(flex='1 1', height='500px', padding='10px')
map2.layout = Layout(flex='1 1', height='500px', padding='10px')
map1.add_layer(infrared_before)
map2.add_layer(infrared_after)
box = HBox([map1, map2])
box
From the maps above, we are able to visually observe the burn scars. Next, let us repeat this process, but this time, we will try to quantify the extent of forest fire.
Quantitative Assessment¶
A Normalized Burn Ratio (NBR) can be used to delineate the burned areas and identify the severity of the fire. The formula for NBR is very similar to that of NDVI except that it uses near-infrared band 5 and the short-wave infrared band 7:
$$ NBR = \frac{B_5 - B_7}{B_5 + B_7} $$The NBR equation was designed to be calculated from reflectance, but it can be calculated from radiance and digital_number (dn) with changes to the burn severity (discussed in the table below). For a given area, an NBR is calculated from an image just prior to the burn and a second NBR is calculated for an image immediately following the burn. Burn extent and severity is evaluated by taking the difference between these two index layers:
$$ \Delta NBR = NBR_{prefire} - NBR_{postfire} $$The meaning of the ∆NBR values can vary by scene, and interpretation in specific instances should always be based on some field assessment. However, the following table from the USGS FireMon program can be useful as a first approximation for interpreting the NBR difference:
\begin{align}{\Delta \mathbf{NBR}} \end{align} | Burn Severity |
---|---|
-2.0 to 0.1 | Regrowth and Unburned |
0.1 to 0.27 | Low severity burn |
0.27 to 0.44 | Medium severity burn |
0.44 to 0.66 | Moderate severity burn |
> 0.66 | High severity burn |
Source: http://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:normalized_burn_ratio
# Calculate before/after NBR indices and their difference
nbr_prefire = band_arithmetic(before_l8,
band_indexes = "(b5 - b7) / (b5 + b7)")
nbr_postfire = band_arithmetic(after_l8,
band_indexes = "(b5 - b7) / (b5 + b7)")
nbr_diff = nbr_prefire - nbr_postfire
# Use Remap function to reclassify the NBR difference score to 1-5
nbr_diff_remap = remap(nbr_diff,
input_ranges=[-2.0, 0.1, # Regrowth and Unburned
0.1, 0.27, # Low Severity burn
0.27, 0.44, # Medium Severity burn
0.44, 0.66, # Moderate Severity
0.66, 2.00], # High Severity
output_values=[1, 2, 3, 4, 5],
astype='u8')
# Create a colormap to show reclassified NBR indices with different color
burn_severity = colormap(nbr_diff_remap,
colormap=[[1, 56, 168, 0], [2, 141, 212, 0],
[3, 255, 255, 0], [4, 255, 128, 0],
[5, 255, 0, 0]])
To view the raster function chain visually, we install graphviz
Python library.
! conda install graphviz -y
# Overview of what raster functions have been applied to
# create burn_serverity layer
burn_severity.draw_graph()
# Visualize burnt areas
burn_severity