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.

Table of Contents¶

Import required libraries
Get data
Create-a-burn-severity-map
- Visual Assessment
- Quantitative Assessment
Create a slope index map
Create the landslide risk map
Summarize landslide risk by subbasin
Conclusion

Import required libraries¶

In [1]:

import arcgis
from arcgis.gis import GIS
from arcgis.raster.functions import *
from ipywidgets import *

In [2]:

gis = GIS('https://pythonapi.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')
arcgis.raster.analytics.is_supported(gis)

Out[2]:

True

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

In [3]:

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

Out[3]:

In [4]:

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

Out[4]:

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.

In [5]:

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.

In [6]:

# 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

Out[6]:

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

In [7]:

# 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

In [8]:

# 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.

In [9]:

! conda install graphviz -y

Solving environment: done

## Package Plan ##

  environment location: /opt/conda

  added / updated specs: 
    - graphviz


The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    graphviz-2.40.1            |       h21bd128_2         6.9 MB  defaults
    ca-certificates-2018.12.5  |                0         123 KB  defaults
    ------------------------------------------------------------
                                           Total:         7.1 MB

The following packages will be UPDATED:

    ca-certificates: 2018.11.29-ha4d7672_0    conda-forge --> 2018.12.5-0             defaults
    cryptography:    2.3.1-py36hb7f436b_1000  conda-forge --> 2.4.2-py36h1ba5d50_0    defaults
    curl:            7.63.0-h646f8bb_1000     conda-forge --> 7.63.0-hbc83047_1000    defaults
    graphviz:        2.40.1-h21bd128_2        anaconda    --> 2.40.1-h21bd128_2       defaults
    grpcio:          1.16.0-py36h4f00d22_1000 conda-forge --> 1.16.1-py36hf8bcb03_1   defaults
    libcurl:         7.63.0-h01ee5af_1000     conda-forge --> 7.63.0-h20c2e04_1000    defaults
    libpq:           10.6-h13b8bad_1000       conda-forge --> 11.1-h20c2e04_0         defaults
    openssl:         1.0.2p-h14c3975_1002     conda-forge --> 1.1.1a-h7b6447c_0       defaults
    pycurl:          7.43.0.2-py36hb7f436b_0  defaults    --> 7.43.0.2-py36h1ba5d50_0 defaults
    python:          3.6.6-hd21baee_1003      conda-forge --> 3.6.8-h0371630_0        defaults
    qt:              5.9.6-h8703b6f_2         defaults    --> 5.9.7-h5867ecd_1        defaults

The following packages will be DOWNGRADED:

    certifi:         2018.11.29-py36_1000     conda-forge --> 2018.11.29-py36_0       defaults
    krb5:            1.16.2-hc83ff2d_1000     conda-forge --> 1.16.1-h173b8e3_7       defaults
    libssh2:         1.8.0-h1ad7b7a_1003      conda-forge --> 1.8.0-h1ba5d50_4        defaults


Downloading and Extracting Packages
graphviz-2.40.1      | 6.9 MB    | ##################################### | 100% 
ca-certificates-2018 | 123 KB    | ##################################### | 100% 
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

In [10]:

# Overview of what raster functions have been applied to 
# create burn_serverity layer
burn_severity.draw_graph()

Out[10]:

In [11]:

# Visualize burnt areas
burn_severity

Out[11]: